How do you interoperate with other applications using double-precision floating-point?

General FreeBASIC programming questions.
cbruce
Posts: 163
Joined: Sep 12, 2007 19:13
Location: Dallas, Texas

How do you interoperate with other applications using double-precision floating-point?

Post by cbruce »

Since double-precision floating-point values/digits/rounding can (do) vary between programming languages and platforms... how do you write code that uses doubles that needs to interoperate with other people's applications?

Thanks!
Bruce
MrSwiss
Posts: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: How do you interoperate with other applications using double-precision floating-point?

Post by MrSwiss »

Seems that you are somewhat confused, since a Double (FB or C, ...) is a defined
quantity, in terms of being 64 bit / 53 bit mantissa (52 + 1 implied), as defined by:
IEEE 754, binary64. (in FB/C called: Double)
Therefore, Double in FB = Double in any other language, the Name can be different
(the format remains the same, if standard compliant).
cbruce
Posts: 163
Joined: Sep 12, 2007 19:13
Location: Dallas, Texas

Re: How do you interoperate with other applications using double-precision floating-point?

Post by cbruce »

.
I've already had problems with the number of significant digits MrSwiss.

Here's why... (I think)... because of the number of significant digits available to Doubles in different compilers, runtimes, etc.:

Code: Select all

 FreeBASIC (15 digits):
 4.940656458412465   e-324 to  1.797693134862316   e+308
-4.940656458412465   e-324 to -1.797693134862316   e+308

 Visual Basic (17 digits):
 4.94065645841246544 e-324 to  1.79769313486231570 e+308
-4.94065645841246544 e-324 to -1.79769313486231570 e+308

 VBA (13 digits):
 4.94065645841247    e-324 to  1.79769313486232    e+308
-4.94065645841247    e-324 to -1.79769313486231    e+308

 .NET (13 digits)
-1.79769313486232    e+308 to  1.79769313486232    e+308
(don't ask me why this is weird - it's from the Microsoft .NET 4.5 doc.)

 IEEE 754 (16 digits):
 4.9406564584124654  e-324 to  1.7976931348623157  e+308
-4.9406564584124654  e-324 to -1.7976931348623157  e+308 
I have exactly the same function coded in FreeBASIC that is in a paper I'm working from where the function was coded in MatLab. Here are the results of the two:

Code: Select all

FreeBASIC result:
z[n+1] =  0.007268783820012809     + 0.05340708046876657 * i

MatLab result:
z[n+1] =  0.0072687838200128062515 + 0.05340708046876657438 * i
"if FB = MatLab" returns FALSE... even though it's the same function and everyone should be using IEEE 754 Doubles.
.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: How do you interoperate with other applications using double-precision floating-point?

Post by srvaldez »

MathWorks wrote:MATLAB constructs the double-precision (or double) data type according to IEEE® Standard 754 for double precision.
so does FB, please show some math calculations (not jus the results) in MATLAB and the same in FB
Last edited by srvaldez on May 21, 2019 2:25, edited 1 time in total.
cbruce
Posts: 163
Joined: Sep 12, 2007 19:13
Location: Dallas, Texas

Re: How do you interoperate with other applications using double-precision floating-point?

Post by cbruce »

IEEE 754 64-bit
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: How do you interoperate with other applications using double-precision floating-point?

Post by srvaldez »

my apologies cbruce, I edited my post after searching the web.
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: How do you interoperate with other applications using double-precision floating-point?

Post by jj2007 »

srvaldez wrote:please show some math calculations (not jus the results) in MATLAB and the same in FB
Right, results only don't make much sense. For your inspiration:

Code: Select all

Testl2e	REAL10 1.442695040888963407
Testl2t	REAL10 3.321928094887362348
TestPI 	REAL10 3.141592653589793238

  Init
  fld Testl2e
  fld Testl2t
  fld TestPI
  fmul
  fmul
  deb 4, "Log2(e)*Log2(10)*PI", ST(0)
EndOfCode
Output: Log2(e)*Log2(10)*PI ST(0) 15.05617449128342534

This is what the hardware allows. In case of doubt: Wolfram Alpha has a much better precision to offer.
cbruce
Posts: 163
Joined: Sep 12, 2007 19:13
Location: Dallas, Texas

Re: How do you interoperate with other applications using double-precision floating-point?

Post by cbruce »

.
First - in case you are having problems sleeping... here is a little light reading:

The pitfalls of verifying floating-point computations - Monniaux - May 23, 2008 - 0701192.pdf
(This one talks specifically about the CPU register and memory issues)
https://arxiv.org/pdf/cs/0701192.pdf

What Every Computer Scientist Should Know About Floating-Point Arithmetic
https://docs.oracle.com/cd/E19957-01/80 ... dberg.html

What you never wanted to know about floating point but will be forced to find out
https://www.volkerschatz.com/science/float.html

The main takeaway in all of this turns out to be:
  1. Rounding errors due to compiler optimizations and re-sequencing of instructions.
  2. Floating-point values can change when they are stored, as occurs when they are simply spilled to the stack. This occurs because some processors and C++ implementations may store floating-point values with extended precision in registers but less precision in memory.
A friend of mine sent me these links and told me that for their software they require third parties to certify interoperability of their packages, both in software and on platform, in the areas where they utilize floating point calculations - specifically because of these problems.

(I'll give you a code stub next).
.
Last edited by cbruce on May 21, 2019 3:36, edited 1 time in total.
cbruce
Posts: 163
Joined: Sep 12, 2007 19:13
Location: Dallas, Texas

Re: How do you interoperate with other applications using double-precision floating-point?

Post by cbruce »

You'll need the fbcomplex library ...

Code: Select all

' ******************************************************************
' fbcomplex.bas : Complex number library for FreeBASIC
' ------------------------------------------------------------------
' Based on ComplexMath Delphi library by E. F. Glynn
' http://www.efg2.com/Lab/Mathematics/Complex/index.html
' ------------------------------------------------------------------
' Compile with: fbc -c fbcomplex.bas 
'         then: ar r libfbcomplex.a fbcomplex.o
' Place libfbcomplex.a in the FB library subdirectory
' ******************************************************************

' ------------------------------------------------------------------
' Mathematical constants
' ------------------------------------------------------------------

#DEFINE MaxNum    1.797693134862315E+308  ' Max. floating point number: 2^1024
#DEFINE MaxLog  709.782712893384          ' Max. argument for Exp
#DEFINE MinLog -708.3964185322641         ' Min. argument for Exp
#DEFINE TwoPi     6.283185307179586       ' 2*Pi
#DEFINE PiDiv2    1.570796326794897       ' Pi/2
#DEFINE Ln2PiDiv2 0.9189385332046728      ' Ln(2*Pi)/2

' ------------------------------------------------------------------
' Error codes for function evaluations
' ------------------------------------------------------------------

#DEFINE FOk         0  ' No error
#DEFINE FDomain    -1  ' Argument domain error
#DEFINE FSing      -2  ' Singularity
#DEFINE FOverflow  -3  ' Overflow range error
#DEFINE FUnderflow -4  ' Underflow range error

' ------------------------------------------------------------------
' Boolean constants
' ------------------------------------------------------------------

'''#DEFINE True -1
'''#DEFINE False 0

' ------------------------------------------------------------------
' Type definition
' ------------------------------------------------------------------

TYPE Complex
  X AS DOUBLE
  Y AS DOUBLE
END TYPE

' ------------------------------------------------------------------
' Global variables
' ------------------------------------------------------------------

DIM SHARED AS INTEGER ErrCode 
' Error code from the latest function evaluation

DIM SHARED AS Complex C_zero   = (0, 0)
DIM SHARED AS Complex C_one    = (1, 0)
DIM SHARED AS Complex C_i      = (0, 1)     
DIM SHARED AS Complex C_MaxNum = (MaxNum, MaxNum)

' ------------------------------------------------------------------
' Function returning error code
' ------------------------------------------------------------------

FUNCTION CMathErr () AS INTEGER
  CMathErr = ErrCode
END FUNCTION

' ------------------------------------------------------------------
' General functions
' ------------------------------------------------------------------

FUNCTION Cmplx(BYVAL X AS DOUBLE, BYVAL Y AS DOUBLE) AS Complex
' Returns the Complex number X + iY

  ErrCode = FOk
  RETURN TYPE<Complex> (X, Y)
END FUNCTION

FUNCTION Polar(BYVAL R AS DOUBLE, BYVAL Theta AS DOUBLE) AS Complex
' Returns the Complex number R * (cos(Theta) + i * sin(Theta))

  ErrCode = FOk
  
  IF R < 0 THEN
    ErrCode = FDomain
    RETURN C_zero
  ELSEIF R = 0 THEN
    RETURN C_zero
  END IF
 
  RETURN TYPE<Complex> (R * COS(Theta), R * SIN(Theta))
END FUNCTION

FUNCTION CReal(BYREF Z AS Complex) AS DOUBLE
' Real part of Z

  ErrCode = FOk
  RETURN Z.X
END FUNCTION

FUNCTION CImag(BYREF Z AS Complex) AS DOUBLE
' Imaginary part of Z

  ErrCode = FOk
  RETURN Z.Y
END FUNCTION

FUNCTION CSgn(BYREF Z AS Complex) AS INTEGER
' Complex sign

  ErrCode = FOk
  IF Z.X > 0 THEN
    RETURN 1
  ELSEIF Z.X < 0 THEN
    RETURN -1
  ELSEIF Z.Y > 0 THEN
    RETURN 1
  ELSEIF Z.Y < 0 THEN
    RETURN -1
  ELSE
    RETURN 0
  END IF
END FUNCTION

SUB CSwap(BYREF X AS Complex, BYREF Y AS Complex)
' Exchanges two Complexes

  DIM AS Complex Temp
  
  ErrCode = FOk

  Temp = X
  X = Y
  Y = Temp
END SUB

FUNCTION CAbs(BYREF Z AS Complex) AS DOUBLE
' Modulus of Z

  ErrCode = FOk
  
  IF Z.X = 0 THEN RETURN ABS(Z.Y)
  IF Z.Y = 0 THEN RETURN ABS(Z.X)

  DIM AS DOUBLE AbsX, AbsY, R

  AbsX = ABS(Z.X)
  AbsY = ABS(Z.Y)

  IF AbsX > AbsY THEN
    R = AbsY / AbsX
    RETURN AbsX * SQR(1 + R * R)
  ELSEIF AbsY = 0 THEN
    RETURN 0
  ELSE
    R = AbsX / AbsY
    RETURN AbsY * SQR(1 + R * R)
  END IF
END FUNCTION

FUNCTION CArg(BYREF Z AS Complex) AS DOUBLE
' Argument of Z

  ErrCode = FOk
  RETURN ATAN2(Z.Y, Z.X)
END FUNCTION

FUNCTION CConj(BYREF Z AS Complex) AS Complex
' Complex conjugate
  
  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X, - Z.Y)    
END FUNCTION

' ------------------------------------------------------------------
' Logarithm, Exponential, Roots
' ------------------------------------------------------------------

FUNCTION CLog(BYREF Z AS Complex) AS Complex
' Principal part of Complex logarithm

  DIM AS DOUBLE R, LnR

  ErrCode = FOk
  
  R = CAbs(Z)
  
  IF R < 0 THEN
    ErrCode = FDomain
    LnR = - MaxNum
  ELSEIF R = 0 THEN
    ErrCode = FSing
    LnR = - MaxNum
  ELSE
    LnR = LOG(R)
  END IF
  
  RETURN TYPE<Complex> (LnR, CArg(Z))
END FUNCTION  

FUNCTION CExp(BYREF Z AS Complex) AS Complex
' Complex exponential

  DIM AS DOUBLE ExpX

  ErrCode = FOk

  IF Z.X < MinLog THEN
    ErrCode = FUnderflow
    ExpX = 0
  ELSEIF Z.X > MaxLog THEN
    ErrCode = FOverflow
    ExpX = MaxNum
  ELSE
    ExpX = EXP(Z.X)
  END IF
  
  RETURN TYPE<Complex> (ExpX * COS(Z.Y), ExpX * SIN(Z.Y))
END FUNCTION  

FUNCTION CRoot(BYREF Z AS Complex, BYVAL K AS INTEGER, BYVAL N AS INTEGER) AS Complex
' CRoot can calculate all 'N' roots of 'Z' by varying 'K' from 0 to N-1
' This is an application of DeMoivre's theorem.

  IF N <= 0 OR K < 0 OR K >= N THEN
    ErrCode = FDomain
    RETURN C_zero  
  END IF

  ErrCode = FOk
  RETURN Polar(CAbs(Z) ^ (1 / N), (CArg(Z) + K * TwoPi) / N)
END FUNCTION

FUNCTION CSqrt(BYREF Z AS Complex) AS Complex
  ErrCode = FOk
  RETURN Polar(SQR(CAbs(Z)), 0.5 * CArg(Z))
END FUNCTION

' ------------------------------------------------------------------
' Operators
' ------------------------------------------------------------------

OPERATOR = (BYREF Z1 AS Complex, BYREF Z2 AS Complex) AS INTEGER
  ErrCode = FOk
  RETURN (Z1.X = Z2.X AND Z1.Y = Z2.Y) 
END OPERATOR

OPERATOR <> (BYREF Z1 AS Complex, BYREF Z2 AS Complex) AS INTEGER
  ErrCode = FOk
  RETURN (Z1.X <> Z2.X OR Z1.Y <> Z2.Y) 
END OPERATOR

OPERATOR + (BYREF Z1 AS Complex, BYREF Z2 AS Complex) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (Z1.X + Z2.X, Z1.Y + Z2.Y) 
END OPERATOR

OPERATOR + (BYVAL A AS DOUBLE, BYREF Z AS Complex) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X + A, Z.Y) 
END OPERATOR

OPERATOR + (BYREF Z AS Complex, BYVAL A AS DOUBLE) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X + A, Z.Y) 
END OPERATOR

OPERATOR - (BYREF Z1 AS Complex, BYREF Z2 AS Complex) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (Z1.X - Z2.X, Z1.Y - Z2.Y) 
END OPERATOR

OPERATOR - (BYVAL A AS DOUBLE, BYREF Z AS Complex) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (A - Z.X, - Z.Y) 
END OPERATOR

OPERATOR - (BYREF Z AS Complex, BYVAL A AS DOUBLE) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X - A, Z.Y) 
END OPERATOR

OPERATOR - (BYREF Z AS Complex) AS Complex
' Unary minus

  ErrCode = FOk
  RETURN TYPE<Complex> (- Z.X, - Z.Y) 
END OPERATOR

OPERATOR * (BYREF Z1 AS Complex, BYREF Z2 AS Complex) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (Z1.X * Z2.X - Z1.Y * Z2.Y, _
                        Z1.X * Z2.Y + Z1.Y * Z2.X) 
END OPERATOR

OPERATOR * (BYVAL A AS DOUBLE, BYREF Z AS Complex) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (A * Z.X, A * Z.Y) 
END OPERATOR

OPERATOR * (BYREF Z AS Complex, BYVAL A AS DOUBLE) AS Complex
  ErrCode = FOk
  RETURN TYPE<Complex> (A * Z.X, A * Z.Y)
END OPERATOR

OPERATOR / (BYREF Z1 AS Complex, BYREF Z2 AS Complex) AS Complex

  IF Z2.X = 0 AND Z2.Y = 0 THEN
    ErrCode = FOverflow
    RETURN C_MaxNum
  END IF

  DIM AS DOUBLE Temp = 1 / (Z2.X * Z2.X + Z2.Y * Z2.Y)
  
  ErrCode = FOk
  RETURN TYPE<Complex> ((Z1.X * Z2.X + Z1.Y * Z2.Y) * Temp, _
                        (Z1.Y * Z2.X - Z1.X * Z2.Y) * Temp)
END OPERATOR

OPERATOR / (BYVAL A AS DOUBLE, BYREF Z AS Complex) AS Complex

  IF Z.X = 0 AND Z.Y = 0 THEN
    ErrCode = FOverflow
    RETURN C_MaxNum
  END IF

  DIM AS DOUBLE Temp = A / (Z.X * Z.X + Z.Y * Z.Y)
  
  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X * Temp, - Z.Y * Temp)
END OPERATOR

OPERATOR / (BYREF Z AS Complex, BYVAL A AS DOUBLE) AS Complex

  IF A = 0 THEN
    ErrCode = FOverflow
    RETURN C_MaxNum
  END IF

  ErrCode = FOk
  RETURN TYPE<Complex> (Z.X / A, Z.Y / A)
END OPERATOR

OPERATOR ^ (BYREF A AS Complex, BYVAL N AS INTEGER) AS Complex
' Computes A^N by repeated multiplications

  DIM AS INTEGER M
  DIM AS Complex B, Res
  
  ErrCode = FOk

  IF N < 0 THEN 
    M = - N  
    B = 1 / A
    IF ErrCode = FOverflow THEN RETURN C_MaxNum
  ELSE
    M = N
    B = A
  END IF  
  
  ' Legendre's algorithm for minimizing the number of multiplications
  
  Res = C_one
  DO WHILE M > 0
    IF (M AND 1) THEN Res = Res * B
    B = B * B
    M = M SHR 1
  LOOP
  
  RETURN Res
END OPERATOR

OPERATOR ^ (BYREF A AS Complex, BYREF B AS Complex) AS Complex

  ErrCode = FOk
  
  IF A.X = 0 AND A.Y = 0 THEN
    IF B.X = 0 AND B.Y = 0 THEN
      RETURN C_one   ' lim A^A = 1 when A -> 0
    ELSE
      RETURN C_zero  ' 0^B = 0
    END IF  
  END IF
  
  RETURN CExp(B * CLog(A))
END OPERATOR

OPERATOR ^ (BYREF A AS Complex, BYVAL B AS DOUBLE) AS Complex

  IF FRAC(B) = 0 THEN RETURN A ^ CINT(B)

  ErrCode = FOk
  
  IF A.X = 0 AND A.Y = 0 THEN
    IF B = 0 THEN
      RETURN C_one
    ELSEIF B > 0 THEN
      RETURN C_zero
    ELSE
      ErrCode = FSing
      RETURN C_MaxNum
    END IF
  END IF

  RETURN Polar(CAbs(A) ^ B, CArg(A) * B)
END OPERATOR

OPERATOR ^ (BYVAL A AS DOUBLE, BYREF B AS Complex) AS Complex
  
  IF A < 0 THEN
    ErrCode = FSing
    RETURN C_zero
  END IF
  
  ErrCode = FOk
  
  IF A = 0 THEN
    IF B.X = 0 AND B.Y = 0 THEN
      RETURN C_one
    ELSE
      RETURN C_zero
    END IF
  END IF
    
  RETURN CExp(B * LOG(A))
END OPERATOR

' ------------------------------------------------------------------
' Polynomials and rational fractions
' ------------------------------------------------------------------

DECLARE FUNCTION CPoly OVERLOAD (BYREF Z   AS Complex, _ 
                                 Coef()    AS Complex, _ 
                                 BYVAL Deg AS INTEGER) AS Complex

DECLARE FUNCTION CPoly (BYREF Z   AS Complex, _ 
                        Coef()    AS DOUBLE,  _ 
                        BYVAL Deg AS INTEGER) AS Complex
                       
DECLARE FUNCTION CFrac OVERLOAD (BYREF Z    AS Complex, _
                                 Coef()     AS Complex, _ 
                                 BYVAL Deg1 AS INTEGER, _ 
                                 BYVAL Deg2 AS INTEGER) AS Complex
                                 
DECLARE FUNCTION CFrac (BYREF Z    AS Complex, _
                        Coef()     AS DOUBLE, _ 
                        BYVAL Deg1 AS INTEGER, _ 
                        BYVAL Deg2 AS INTEGER) AS Complex
                        
' ------------------------------------------------------------------                        
                                 
FUNCTION CPoly(BYREF Z   AS Complex, _ 
               Coef()    AS Complex, _ 
               BYVAL Deg AS INTEGER) AS Complex

' Polynomial with complex coefficients
' P(Z) = Coef(0) + Coef(1) * Z + Coef(2) * Z^2 +... + Coef(Deg) * Z^Deg

  DIM AS INTEGER I
  DIM AS Complex P

  P = Coef(Deg)
  FOR I = Deg - 1 TO 0 STEP -1
    P = P * Z + Coef(I)
  NEXT I
  
  RETURN P
END FUNCTION

FUNCTION CPoly(BYREF Z   AS Complex, _ 
               Coef()    AS DOUBLE,  _ 
               BYVAL Deg AS INTEGER) AS Complex

' Polynomial with real coefficients

  DIM AS INTEGER I
  DIM AS Complex P = Cmplx(Coef(Deg), 0)
  
  FOR I = Deg - 1 TO 0 STEP -1
    P = P * Z + Coef(I)
  NEXT I
  
  RETURN P
END FUNCTION

FUNCTION CFrac(BYREF Z    AS Complex, _
               Coef()     AS Complex, _ 
               BYVAL Deg1 AS INTEGER, _ 
               BYVAL Deg2 AS INTEGER) AS Complex

' Rational fraction with complex coefficients
'
'               Coef(0) + Coef(1) * Z + ... + Coef(Deg1) * Z^Deg1
' F(Z) = --------------------------------------------------------------------
'        Coef(Deg1 + 1) + Coef(Deg1+2) * Z + ... + Coef(Deg1+Deg2+1) * Z^Deg2


  DIM AS INTEGER I
  DIM AS Complex P = Coef(Deg1)
  DIM AS Complex Q = Coef(Deg1 + Deg2 + 1)

  FOR I = Deg1 - 1 TO 0 STEP -1
    P = P * Z + Coef(I)
  NEXT I

  FOR I = Deg1 + Deg2 TO Deg1 + 1 STEP -1
    Q = Q * Z + Coef(I)
  NEXT I
  
  RETURN P / Q
END FUNCTION 

FUNCTION CFrac(BYREF Z    AS Complex, _
               Coef()     AS DOUBLE,  _ 
               BYVAL Deg1 AS INTEGER, _ 
               BYVAL Deg2 AS INTEGER) AS Complex

' Rational fraction with real coefficients

  DIM AS INTEGER I
  DIM AS Complex P = Cmplx(Coef(Deg1), 0)
  DIM AS Complex Q = Cmplx(Coef(Deg1 + Deg2 + 1), 0)

  FOR I = Deg1 - 1 TO 0 STEP -1
    P = P * Z + Coef(I)
  NEXT I

  FOR I = Deg1 + Deg2 TO Deg1 + 1 STEP -1
    Q = Q * Z + Coef(I)
  NEXT I
  
  RETURN P / Q
END FUNCTION 

' ------------------------------------------------------------------
' Trigonometric functions
' ------------------------------------------------------------------

SUB SinhCosh(BYVAL X AS DOUBLE, BYREF SinhX AS DOUBLE, BYREF CoshX AS DOUBLE)
' Compute Sinh(X) and Cosh(X) with only one exponential     
  
  IF X < MinLog OR X > MaxLog THEN
    ErrCode = FOverflow
    CoshX = MaxNum
    SinhX = SGN(X) * CoshX
    EXIT SUB
  END IF
  
  DIM AS DOUBLE ExpX, ExpMinusX
  
  ErrCode = FOk
  ExpX = EXP(X)
  ExpMinusX = 1 / ExpX
  SinhX = 0.5 * (ExpX - ExpMinusX)
  CoshX = 0.5 * (ExpX + ExpMinusX)
END SUB  

FUNCTION CSin(BYREF Z AS Complex) AS Complex
  
  DIM AS DOUBLE SinhY, CoshY
  
  SinhCosh Z.Y, SinhY, CoshY
  RETURN TYPE<Complex> (SIN(Z.X) * CoshY, COS(Z.X) * SinhY)
END FUNCTION

FUNCTION CCos(BYREF Z AS Complex) AS Complex
  
  DIM AS DOUBLE SinhY, CoshY
  
  SinhCosh Z.Y, SinhY, CoshY
  RETURN TYPE<Complex> (COS(Z.X) * CoshY, - SIN(Z.X) * SinhY)
END FUNCTION

SUB CSinCos(BYREF Z AS Complex, BYREF SinZ AS Complex, BYREF CosZ AS Complex)
    
  DIM AS DOUBLE SinX, CosX, SinhY, CoshY
  
  SinX = SIN(Z.X)
  CosX = COS(Z.X)
  
  SinhCosh Z.Y, SinhY, CoshY
  
  SinZ = TYPE<Complex> (SinX * CoshY, CosX * SinhY)
  CosZ = TYPE<Complex> (CosX * CoshY, - SinX * SinhY)
END SUB    
    
FUNCTION CTan(BYREF Z AS Complex) AS Complex

  DIM AS DOUBLE X2, Y2, SinX2, CosX2, SinhY2, CoshY2, Temp
  
  X2 = 2 * Z.X
  Y2 = 2 * Z.Y
  
  SinX2 = SIN(X2)
  CosX2 = COS(X2)
  
  SinhCosh Y2, SinhY2, CoshY2
  
  IF ErrCode = FOk THEN Temp = CosX2 + CoshY2 ELSE Temp = CoshY2
  
  IF Temp <> 0 THEN
    RETURN TYPE<Complex> (SinX2 / Temp, SinhY2 / Temp)
  ELSE
    ' Z = Pi/2 + k*Pi 
    ErrCode = FSing  
    RETURN TYPE<Complex> (MaxNum, 0)
  END IF
END FUNCTION

FUNCTION CASin(BYREF Z AS Complex) AS Complex
    
  DIM AS DOUBLE  Rp, Rm, S, T, X2, XX, YY
  DIM AS Complex Z1
  
  Z1 = TYPE<Complex> (Z.Y, - Z.X)
  
  X2 = 2 * Z.X
  XX = Z.X * Z.X
  YY = Z.Y * Z.Y
  S  = XX + YY + 1
  Rp = 0.5 * SQR(S + X2)
  Rm = 0.5 * SQR(S - X2)
  T  = Rp + Rm
   
  ErrCode = FOk
  
  RETURN TYPE<Complex> (ASIN(Rp - Rm), _
                        CSgn(Z1) * LOG(T + SQR(T * T - 1))) 
END FUNCTION

FUNCTION CACos(BYREF Z AS Complex) AS Complex
  ErrCode = FOk
  RETURN PiDiv2 - CASin(Z)
END FUNCTION

FUNCTION CATan(BYREF Z AS Complex) AS Complex
    
  IF (Z.X = 0) AND (ABS(Z.Y) = 1) THEN  ' Z = +/- i 
    ErrCode = FSing
    RETURN TYPE<Complex> (0, SGN(Z.Y) * MaxNum)
  END IF

  DIM AS DOUBLE XX, YY, Yp1, Ym1
  
  XX  = Z.X * Z.X
  YY  = Z.Y * Z.Y
  Yp1 = Z.Y + 1
  Ym1 = Z.Y - 1

  ErrCode = FOk
  
  RETURN TYPE<Complex> (0.5 * (ATAN2(Z.X, - Ym1) - ATAN2(- Z.X, Yp1)), _
                        0.25 * LOG((XX + Yp1 * Yp1) / (XX + Ym1 * Ym1)))  
END FUNCTION

' ------------------------------------------------------------------
' Hyperbolic functions
' ------------------------------------------------------------------

FUNCTION CSinh(BYREF Z AS Complex) AS Complex
  
  DIM AS DOUBLE SinhX, CoshX, SinY, CosY
  
  SinY = SIN(Z.Y)
  CosY = COS(Z.Y)
  
  SinhCosh Z.X, SinhX, CoshX
  
  ErrCode = FOk
  RETURN TYPE<Complex> (SinhX * CosY, CoshX * SinY) 
END FUNCTION

FUNCTION CCosh(BYREF Z AS Complex) AS Complex
    
  DIM AS DOUBLE SinhX, CoshX, SinY, CosY
  
  SinY = SIN(Z.Y)
  CosY = COS(Z.Y)
  
  SinhCosh Z.X, SinhX, CoshX
  
  ErrCode = FOk
  RETURN TYPE<Complex> (CoshX * CosY, SinhX * SinY)  
END FUNCTION

SUB CSinhCosh(BYREF Z AS Complex, BYREF SinhZ AS Complex, BYREF CoshZ AS Complex)
    
  DIM AS DOUBLE SinY, CosY, SinhX, CoshX
  
  SinY = SIN(Z.Y)
  CosY = COS(Z.Y)
  
  SinhCosh Z.X, SinhX, CoshX
  
  SinhZ = TYPE<Complex> (SinhX * CosY, CoshX * SinY)
  CoshZ = TYPE<Complex> (CoshX * CosY, SinhX * SinY)
END SUB    

FUNCTION CTanh(BYREF Z AS Complex) AS Complex

  DIM AS DOUBLE X2, Y2, SinY2, CosY2, SinhX2, CoshX2, Temp
  
  X2 = 2.0 * Z.X
  Y2 = 2.0 * Z.Y
  
  SinY2 = SIN(Y2)
  CosY2 = COS(Y2)
  
  SinhCosh X2, SinhX2, CoshX2
  
  IF ErrCode = FOk THEN Temp = CoshX2 + CosY2 ELSE Temp = CoshX2

  IF Temp = 0 THEN  ' Z = i * (Pi/2 + k*Pi) 
    ErrCode = FSing
    RETURN TYPE<Complex> (0, MaxNum)
  END IF  

  ErrCode = FOk   
  RETURN TYPE<Complex> (SinhX2 / Temp, SinY2 / Temp)
END FUNCTION

FUNCTION CASinh(BYREF Z AS Complex) AS Complex
    
  ErrCode = FOk
  RETURN - C_i * CASin(C_i * Z)
END FUNCTION

FUNCTION CACosh(BYREF Z AS Complex) AS Complex
  
  ErrCode = FOk
  RETURN CSgn(Cmplx(Z.Y, 1 - Z.X)) * C_i * CACos(Z)
END FUNCTION

FUNCTION CATanh(BYREF Z AS Complex) AS Complex
  
  IF ABS(Z.X) = 1 AND Z.Y = 0 THEN  ' Z = +/- 1
    ErrCode = FSing
    RETURN TYPE<complex> (SGN(Z.X) * MaxNum, 0)
  END IF

  ErrCode = FOk
  RETURN - C_i * CAtan(C_i * Z) 
END FUNCTION

' ------------------------------------------------------------------
' Gamma function
' ------------------------------------------------------------------

FUNCTION CApproxLnGamma(BYREF Z AS Complex) AS Complex
'  This is the approximation used in the National Bureau of
'  Standards "Table of the Gamma Function for Complex Arguments,"
'  Applied Mathematics Series 34, 1954. The NBS table was created
'  using this approximation over the area 9 < Re(z) < 10 and
'  0 < Im(z) < 10. Other table values were computed using the
'  relationship:
'        _                   _
'    ln | (z+1) = ln z + ln | (z) 

  DIM AS DOUBLE C(1 TO 8) = _
    {8.33333333333333E-02, -2.77777777777778E-03, _
     7.93650793650794E-04, -5.95238095238095E-04, _
     8.41750841750842E-04, -1.91752691752692E-03, _
     6.41025641025641E-03, -2.95506535947712E-02}
  
  DIM AS Complex Powers(1 TO 8), Temp, Sum
  
  DIM AS INTEGER I
  
  Powers(1) = 1 / Z
  
  Temp = Powers(1) * Powers(1)
  
  FOR I = 2 TO 8
    Powers(I) = Powers(I - 1) * Temp
  NEXT I  
  
  Sum = (Z - 0.5) * CLog(Z) - Z + Ln2PiDiv2 
  
  FOR I = 8 TO 1 STEP -1
    Sum = Sum + C(i) * Powers(i)
  NEXT I
  
  RETURN Sum
END FUNCTION

FUNCTION CLnGamma(BYREF Z AS Complex) AS Complex
  
  ErrCode = FOk
  
  ' Negative integer
  IF Z.X <= 0 AND FRAC(Z.X) = 0 AND Z.Y = 0 THEN  
    ErrCode = FSing
    RETURN C_MaxNum
  END IF
  
  ' 3rd or 4th quadrant
  IF Z.Y < 0 THEN                     
    RETURN CConj(CLnGamma(CConj(Z))) 
  END IF  

  ' "left" of NBS table range
  IF Z.X < 9 THEN  
    RETURN CLnGamma(Z + 1) - CLog(Z)
  END IF

  ' NBS table range: 9 < Re(z) < 10 
  RETURN CApproxLnGamma(Z)  
END FUNCTION    

' ------------------------------------------------------------------
' Printing functions
' ------------------------------------------------------------------

SUB CPrint(BYREF Z AS Complex, BYREF Mask AS STRING = "####.####")
' Prints a complex in rectangular form

  PRINT USING Mask; Z.X;
  IF Z.Y >= 0 THEN PRINT " + "; ELSE PRINT " - ";
  PRINT USING Mask; ABS(Z.Y); 
  PRINT " * i";
END SUB

SUB CPrintPolar(BYREF Z AS Complex, BYREF Mask AS STRING = "####.####")
' Prints a complex in polar form: Z = |Z| * exp(Arg(Z) * i)

  PRINT USING Mask; CAbs(Z); 
  PRINT " * exp(";
  PRINT USING Mask; CArg(Z); 
  PRINT " * i)"
END SUB

cbruce
Posts: 163
Joined: Sep 12, 2007 19:13
Location: Dallas, Texas

Re: How do you interoperate with other applications using double-precision floating-point?

Post by cbruce »

.
And here is a test code stub that will compile and run...

Code: Select all

' =======================================================================
' This stub program evaluates the Mann iteration of a Mandelbrot function.
' =======================================================================

#include "fbcomplex.bas"

function CStr(BYREF Z AS Complex) as string
    ' Convert a complex to a string in rectangular form.
    dim tmp as string
    tmp = str(Z.X)
    IF Z.Y >= 0 THEN tmp += " + " ELSE tmp += " - "
    tmp += str(ABS(Z.Y))
    tmp += " * i"
    return tmp
END function

SUB Mandelbrot (xt AS DOUBLE, yt AS DOUBLE)
    DIM AS Complex c, z, zn
    c = Cmplx(xt, yt)
    z = c
    dim as string tmp
    dim as integer n = 4      'A static iteration value for tesing.
    dim as double  s = 0.6    'A static s value for testing.
    dim as complex e = type <complex> (0.015923, -0.03179523) 'A static e value for testing.
    print
    print "      c = "; cstr(c)
    print "      s = "; str(s)
    print
    print "      n = "; str(n)
    print "      e = "; cstr(e)
    print
    ' Mandelbrot equation - Mann iteration method.
    ' f(z[n]) = z[n]*c*e              'Where c,z,e are Complex and z[0] = c
    ' z[n+1] = s*f(z[n])+(1-s)*z[n]   'Where 0<s<1
    for iter as integer = 0 to (n - 1)
        zn = z * c * e
        print "   iter = "; str(iter)
        print "      z = "; cstr(z)
        print "f(z[n]) = "; cstr(zn)
        '
        z  = (s * zn) + ((1-s) * z)
        print " z[n+1] = "; cstr(z)
        print
    next
    print
END SUB


Mandelbrot(0.325, 1.5125)  'pass a static c value for testing.
END

And this should be the result:

Code: Select all

      c = 0.325 + 1.5125 * i
      s = 0.6

      n = 4
      e = 0.015923 - 0.03179523 * i

   iter = 0
      z = 0.325 + 1.5125 * i
f(z[n]) = -0.0034857981 + 0.08503248483593751 * i
 z[n+1] = 0.12790852114 + 0.6560194909015625 * i

   iter = 1
      z = 0.12790852114 + 0.6560194909015625 * i
f(z[n]) = -0.002207244882903518 + 0.03670180238359441 * i
 z[n+1] = 0.0498390615262579 + 0.2844288777907816 * i

   iter = 2
      z = 0.0498390615262579 + 0.2844288777907816 * i
f(z[n]) = -0.001256231439215502 + 0.01583546970305873 * i
 z[n+1] = 0.01918188574697386 + 0.1232728329381479 * i

   iter = 3
      z = 0.01918188574697386 + 0.1232728329381479 * i
f(z[n]) = -0.0006732841312945596 + 0.006829912155845688 * i
 z[n+1] = 0.007268783820012809 + 0.05340708046876657 * i
.
cbruce
Posts: 163
Joined: Sep 12, 2007 19:13
Location: Dallas, Texas

Re: How do you interoperate with other applications using double-precision floating-point?

Post by cbruce »

.
Oh... and the MatLab code is just the couple of equations that I've commented into the code stub.

I can't show proper math equations here, so that's the best I could do.

Note:
The results are technically correct - it's just that there is a difference in the precision of the results in the FreeBASIC vs. MatLab instances.
cbruce
Posts: 163
Joined: Sep 12, 2007 19:13
Location: Dallas, Texas

Re: How do you interoperate with other applications using double-precision floating-point?

Post by cbruce »

.
Oh again...

My friend recommended that, if possible, when trying to implement interactive applications (with double precision needs) with multiple code bases and on multiple platforms... "IF" it is just the current running instances that need to pass accurate data back and forth for immediate needs (no storage and later retrieval of data by other applications)...

Have the apps on each end spin up something like the stub I just implemented - programmatically check their results and set the number of significant digits that they will utilize in the current instances based on several of these checks. If you only need to determine equality at a couple of points in your code - this can be a good solution that will allow you to keep the maximum number of significant digits in your answers.

Otherwise, just set a "stupid" (low) threshold rounding limit on everything and try that.
.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: How do you interoperate with other applications using double-precision floating-point?

Post by srvaldez »

hello cbruce
I suspect that it could be a precision difference in the complex math calculations between that of FB and Matlab, that is, one or the other does the complex math more accurately, unfortunately I can't verify this.
cbruce
Posts: 163
Joined: Sep 12, 2007 19:13
Location: Dallas, Texas

Re: How do you interoperate with other applications using double-precision floating-point?

Post by cbruce »

.
Thanks srvaldez, but it looks like no one can verify it between different compilers/platforms without a LOT of testing. So now I'm just trying to figure out which way I want to work around the issue in my application (as described above).
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: How do you interoperate with other applications using double-precision floating-point?

Post by srvaldez »

hello cbruce
cbruce wrote: I have exactly the same function coded in FreeBASIC that is in a paper I'm working from where the function was coded in MatLab. Here are the results of the two:

Code: Select all

FreeBASIC result:
z[n+1] =  0.007268783820012809     + 0.05340708046876657 * i

MatLab result:
z[n+1] =  0.0072687838200128062515 + 0.05340708046876657438 * i
"if FB = MatLab" returns FALSE... even though it's the same function and everyone should be using IEEE 754 Doubles.
.
I performed the calculations using mpfr, and it looks like Matlab is a bit more accurate

Code: Select all

z[n+1] =   0.007268783820012805476843273672420147277562822265626 +  0.05340708046876657175587196067215500961693672363281*I
Post Reply