How do you interoperate with other applications using double-precision floating-point?
How do you interoperate with other applications using double-precision floating-point?
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
Thanks!
Bruce
Re: How do you interoperate with other applications using double-precision floating-point?
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).
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).
Re: How do you interoperate with other applications using double-precision floating-point?
.
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.:
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:
"if FB = MatLab" returns FALSE... even though it's the same function and everyone should be using IEEE 754 Doubles.
.
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
Code: Select all
FreeBASIC result:
z[n+1] = 0.007268783820012809 + 0.05340708046876657 * i
MatLab result:
z[n+1] = 0.0072687838200128062515 + 0.05340708046876657438 * i
.
Re: How do you interoperate with other applications using double-precision floating-point?
so does FB, please show some math calculations (not jus the results) in MATLAB and the same in FBMathWorks wrote:MATLAB constructs the double-precision (or double) data type according to IEEE® Standard 754 for double precision.
Last edited by srvaldez on May 21, 2019 2:25, edited 1 time in total.
Re: How do you interoperate with other applications using double-precision floating-point?
my apologies cbruce, I edited my post after searching the web.
Re: How do you interoperate with other applications using double-precision floating-point?
Right, results only don't make much sense. For your inspiration:srvaldez wrote:please show some math calculations (not jus the results) in MATLAB and the same in FB
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
This is what the hardware allows. In case of doubt: Wolfram Alpha has a much better precision to offer.
Re: How do you interoperate with other applications using double-precision floating-point?
.
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:
(I'll give you a code stub next).
.
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:
- Rounding errors due to compiler optimizations and re-sequencing of instructions.
- 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.
(I'll give you a code stub next).
.
Last edited by cbruce on May 21, 2019 3:36, edited 1 time in total.
Re: How do you interoperate with other applications using double-precision floating-point?
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
Re: How do you interoperate with other applications using double-precision floating-point?
.
And here is a test code stub that will compile and run...
And this should be the result:
.
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
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
Re: How do you interoperate with other applications using double-precision floating-point?
.
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.
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.
Re: How do you interoperate with other applications using double-precision floating-point?
.
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.
.
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.
.
Re: How do you interoperate with other applications using double-precision floating-point?
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.
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.
Re: How do you interoperate with other applications using double-precision floating-point?
.
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).
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).
Re: How do you interoperate with other applications using double-precision floating-point?
hello cbruce
I performed the calculations using mpfr, and it looks like Matlab is a bit more accuratecbruce 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:
"if FB = MatLab" returns FALSE... even though it's the same function and everyone should be using IEEE 754 Doubles.Code: Select all
FreeBASIC result: z[n+1] = 0.007268783820012809 + 0.05340708046876657 * i MatLab result: z[n+1] = 0.0072687838200128062515 + 0.05340708046876657438 * i
.
Code: Select all
z[n+1] = 0.007268783820012805476843273672420147277562822265626 + 0.05340708046876657175587196067215500961693672363281*I