' ******************************************************************
'                          FBMath include file
' ******************************************************************

#INCLIB "fbmath"

#INCLUDE "string.bi"

' ******************************************************************
'                         Mathematical constants
' ******************************************************************

#DEFINE Pi         3.141592653589793#   ' Pi
#DEFINE Ln2        0.6931471805599453#  ' Ln(2)
#DEFINE Ln10       2.302585092994046#   ' Ln(10)
#DEFINE LnPi       1.1447298858494#     ' Ln(Pi)
#DEFINE InvLn2     1.442695040888963#   ' 1/Ln(2)
#DEFINE InvLn10    0.4342944819032518#  ' 1/Ln(10)
#DEFINE TwoPi      6.283185307179586#   ' 2*Pi
#DEFINE PiDiv2     1.570796326794897#   ' Pi/2
#DEFINE SqrtPi     1.772453850905516#   ' Sqrt(Pi)
#DEFINE Sqrt2Pi    2.506628274631001#   ' Sqrt(2*Pi)
#DEFINE InvSqrt2Pi 0.3989422804014327#  ' 1/Sqrt(2*Pi)
#DEFINE LnSqrt2Pi  0.9189385332046728#  ' Ln(Sqrt(2*Pi))
#DEFINE Ln2PiDiv2  0.9189385332046728#  ' Ln(2*Pi)/2
#DEFINE Sqrt2      1.414213562373095#   ' Sqrt(2)
#DEFINE Sqrt2Div2  0.7071067811865476#  ' Sqrt(2)/2
#DEFINE Gold       1.618033988749895#   ' Golden Mean (1 + Sqrt(5))/2
#DEFINE CGold      0.3819660112501052#  ' 2 - Gold

' ------------------------------------------------------------------
' Machine-dependent constants
' ------------------------------------------------------------------

#DEFINE MachEp    2.220446049250313D-16  ' Floating point precision: 2^(-52)
#DEFINE MaxNum    1.797693134862315D+308 ' Max. floating point number: 2^1024
#DEFINE MinNum    2.225073858507202D-308 ' Min. floating point number: 2^(-1022)
#DEFINE MaxLog  709.782712893384#        ' Max. argument for Exp = Ln(MaxNum)
#DEFINE MinLog -708.3964185322641#       ' Min. argument for Exp = Ln(MinNum)
#DEFINE MaxFac  170                      ' Max. argument for Factorial
#DEFINE MaxGam  171.624376956302#        ' Max. argument for Gamma
#DEFINE MaxLgm    2.556348D+305          ' Max. argument for LnGamma

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

#IFNDEF True
  #DEFINE True -1
#ENDIF

#IFNDEF False
  #DEFINE False 0
#ENDIF

' ------------------------------------------------------------------
' 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
#DEFINE FPLoss     -5  ' Partial loss of precision

' ------------------------------------------------------------------
' Error codes for matrix computations
' ------------------------------------------------------------------

#DEFINE MatOk       0  ' No error
#DEFINE MatNonConv -1  ' Non-convergence
#DEFINE MatSing    -2  ' Quasi-singular matrix
#DEFINE MatErrDim  -3  ' Non-compatible dimensions
#DEFINE MatNotPD   -4  ' Matrix not positive definite

' ------------------------------------------------------------------
' Error codes for optimization and nonlinear equations
' ------------------------------------------------------------------

#DEFINE OptOk         0  ' No error
#DEFINE OptNonConv   -1  ' Non-convergence
#DEFINE OptSing      -2  ' Quasi-singular hessian matrix
#DEFINE OptBigLambda -3  ' Too high Marquardt parameter

' ******************************************************************
'                         Complex numbers
' ******************************************************************

TYPE Complex
  X AS DOUBLE
  Y AS DOUBLE
END TYPE

' ******************************************************************
'                     Error-handling function
' ******************************************************************

DECLARE FUNCTION MathErr () AS INTEGER
' Returns the error code from the last procedure call

' ******************************************************************
'                 Elementary mathematical functions
' ******************************************************************

DECLARE FUNCTION Min (X AS DOUBLE, Y AS DOUBLE) AS DOUBLE
' Minimum of 2 reals

DECLARE FUNCTION Max (X AS DOUBLE, Y AS DOUBLE) AS DOUBLE
' Maximum of 2 reals

DECLARE FUNCTION IMin (X AS INTEGER, Y AS INTEGER) AS INTEGER
' Minimum of 2 integers

DECLARE FUNCTION IMax (X AS INTEGER, Y AS INTEGER) AS INTEGER
' Maximum of 2 integers

DECLARE FUNCTION Expo (X AS DOUBLE) AS DOUBLE
' Exponential

DECLARE FUNCTION Exp2 (X AS DOUBLE) AS DOUBLE
' 2^X

DECLARE FUNCTION Exp10 (X AS DOUBLE) AS DOUBLE
' 10^X

DECLARE FUNCTION Ln (X AS DOUBLE) AS DOUBLE
' Natural log

DECLARE FUNCTION Log2 (X AS DOUBLE) AS DOUBLE
' Log, base 2

DECLARE FUNCTION Log10 (X AS DOUBLE) AS DOUBLE
' Decimal log

DECLARE FUNCTION LogA (X AS DOUBLE, A AS DOUBLE) AS DOUBLE
' Log, base A

DECLARE FUNCTION FixAngle (Theta AS DOUBLE) AS DOUBLE
' Set Theta in [-Pi..Pi]

DECLARE FUNCTION Pythag (X AS DOUBLE, Y AS DOUBLE) AS DOUBLE
' Sqrt(X^2 + Y^2)

' ******************************************************************
'                        Rounding functions
' ******************************************************************

DECLARE FUNCTION Round(X           AS DOUBLE,_
                       Digit_Count AS INTEGER = 0) AS DOUBLE
' Rounds X to Digit_Count decimal places

DECLARE FUNCTION Ceil(X AS DOUBLE) AS INTEGER
' Ceiling function

DECLARE FUNCTION Floor(X AS DOUBLE) AS INTEGER
' Floor function

' ******************************************************************
'                       Hyperbolic functions
' ******************************************************************

DECLARE FUNCTION Sinh (X AS DOUBLE) AS DOUBLE
' Hyperbolic sine

DECLARE FUNCTION Cosh (X AS DOUBLE) AS DOUBLE
' Hyperbolic cosine

DECLARE FUNCTION Tanh (X AS DOUBLE) AS DOUBLE
' Hyperbolic tangent

DECLARE FUNCTION ASinh (X AS DOUBLE) AS DOUBLE
' Inverse hyperbolic sine

DECLARE FUNCTION ACosh (X AS DOUBLE) AS DOUBLE
' Inverse hyperbolic cosine

DECLARE FUNCTION ATanh (X AS DOUBLE) AS DOUBLE
' Inverse hyperbolic tangent

DECLARE SUB SinhCosh (      X     AS DOUBLE,_
                      BYREF SinhX AS DOUBLE,_
                      BYREF CoshX AS DOUBLE)
' Hyperbolic sine and cosine

' ******************************************************************
'               Gamma function and related functions
' ******************************************************************

DECLARE FUNCTION Fact (N AS INTEGER) AS DOUBLE
' Factorial

DECLARE FUNCTION Gamma (X AS DOUBLE) AS DOUBLE
' Gamma function

DECLARE FUNCTION SgnGamma (X AS DOUBLE) AS INTEGER
' Sign of Gamma function

DECLARE FUNCTION LnGamma (X AS DOUBLE) AS DOUBLE
' Ln of Gamma function

DECLARE FUNCTION DiGamma (X AS DOUBLE) AS DOUBLE
' Digamma or Psi function: d Ln Gamma(X) / dX

DECLARE FUNCTION TriGamma (X AS DOUBLE) AS DOUBLE
' Trigamma or Psi-prime function: d^2 Ln Gamma(X) / dX^2

DECLARE FUNCTION IGamma (A AS DOUBLE, X AS DOUBLE) AS DOUBLE
' Incomplete Gamma function

DECLARE FUNCTION JGamma (A AS DOUBLE, X AS DOUBLE) AS DOUBLE
' Complement of incomplete Gamma function

DECLARE FUNCTION InvGamma(A AS DOUBLE, Y AS DOUBLE) AS DOUBLE
' Returns X such that IGamma(A, X) = Y

DECLARE FUNCTION Erf (X AS DOUBLE) AS DOUBLE
' Error function

DECLARE FUNCTION Erfc (X AS DOUBLE) AS DOUBLE
' Complement of error function

DECLARE FUNCTION Stirling (X AS DOUBLE) AS DOUBLE
' Approximates Gamma function with Stirling's formula, for X >= 30

DECLARE FUNCTION StirLog (X AS DOUBLE) AS DOUBLE
' Approximates Ln(Gamma) with Stirling's formula, for X >= 13

' ******************************************************************
'                          Beta function
' ******************************************************************

DECLARE FUNCTION Beta (X AS DOUBLE,_
                       Y AS DOUBLE) AS DOUBLE
' Beta function

DECLARE FUNCTION IBeta (A AS DOUBLE,_
                        B AS DOUBLE,_
                        X AS DOUBLE) AS DOUBLE
' Incomplete Beta function

DECLARE FUNCTION InvBeta(A AS DOUBLE,_
                         B AS DOUBLE,_
                         Y AS DOUBLE) AS DOUBLE
' Returns X such that IBeta(A, B, X) = Y

' ******************************************************************
'                      Lambert's function
' ******************************************************************

DECLARE FUNCTION LambertW(X       AS DOUBLE,_
                          UBranch AS INTEGER = True,_
                          Offset  AS INTEGER = False) AS DOUBLE

' ******************************************************************
'                       Binomial distribution
' ******************************************************************

DECLARE FUNCTION Binomial(N AS INTEGER, K AS INTEGER) AS DOUBLE
' Binomial coefficient C(N,K)

DECLARE FUNCTION PBinom(N AS INTEGER, P AS DOUBLE,_
                        K AS INTEGER) AS DOUBLE
' Prob(X = K) for binomial distribution B(N,P)

DECLARE FUNCTION FBinom(N AS INTEGER, P AS DOUBLE,_
                        K AS INTEGER) AS DOUBLE
' Prob(X <= K) for binomial distribution B(N,P)

' ******************************************************************
'                        Poisson distribution
' ******************************************************************

DECLARE FUNCTION PPoisson(Mu AS DOUBLE, K AS INTEGER) AS DOUBLE
' Prob(X = K) for Poisson distribution P(Mu)

DECLARE FUNCTION FPoisson(Mu AS DOUBLE, K AS INTEGER) AS DOUBLE
' Prob(X <= K) for Poisson distribution P(Mu)

' ******************************************************************
'                        Normal distribution
' ******************************************************************

DECLARE FUNCTION DNorm(X AS DOUBLE) AS DOUBLE
' Density of standard normal distribution

DECLARE FUNCTION FNorm(X AS DOUBLE) AS DOUBLE
' Prob(U <= X) for the standard normal distribution

DECLARE FUNCTION PNorm(X AS DOUBLE) AS DOUBLE
' Prob(|U| > X) for the standard normal distribution

DECLARE FUNCTION InvNorm(P AS DOUBLE) AS DOUBLE
' Returns X such that FNorm(X) = P

' ******************************************************************
'                        Student distribution
' ******************************************************************

DECLARE FUNCTION DStudent(Nu AS INTEGER, X AS DOUBLE) AS DOUBLE
' Density of Student distribution with Nu d.o.f.

DECLARE FUNCTION FStudent(Nu AS INTEGER, X AS DOUBLE) AS DOUBLE
' Prob(t <= X) for Student dist. with Nu d.o.f.

DECLARE FUNCTION PStudent(Nu AS INTEGER, X AS DOUBLE) AS DOUBLE
' Prob(|t| > X) for Student dist. with Nu d.o.f.

DECLARE FUNCTION InvStudent(Nu AS INTEGER, P AS DOUBLE) AS DOUBLE
' Returns X such that FStudent(Nu, X) = P

' ******************************************************************
'                        Snedecor distribution
' ******************************************************************

DECLARE FUNCTION DSnedecor(Nu1 AS INTEGER, Nu2 AS INTEGER,_
                           X AS DOUBLE) AS DOUBLE
' Density of Fisher-Snedecor distribution with Nu1 and Nu2 d.o.f.

DECLARE FUNCTION FSnedecor(Nu1 AS INTEGER, Nu2 AS INTEGER,_
                           X AS DOUBLE) AS DOUBLE
' Prob(F <= X) for Fisher-Snedecor dist. with Nu1 and Nu2 d.o.f.

DECLARE FUNCTION PSnedecor(Nu1 AS INTEGER, Nu2 AS INTEGER,_
                           X AS DOUBLE) AS DOUBLE
' Prob(F > X) for Fisher-Snedecor dist. with Nu1 and Nu2 d.o.f.

DECLARE FUNCTION InvSnedecor(Nu1 AS INTEGER, Nu2 AS INTEGER,_
                             P AS DOUBLE) AS DOUBLE
' Returns X such that FSnedecor(Nu1, Nu2, X) = P

' ******************************************************************
'                        Khi-2 distribution
' ******************************************************************

DECLARE FUNCTION DKhi2(Nu AS INTEGER, X AS DOUBLE) AS DOUBLE
' Density of Khi-2 distribution with Nu d.o.f.

DECLARE FUNCTION FKhi2(Nu AS INTEGER, X AS DOUBLE) AS DOUBLE
' Prob(Khi2 <= X) for the Khi-2 distribution with Nu d.o.f.

DECLARE FUNCTION PKhi2(Nu AS INTEGER, X AS DOUBLE) AS DOUBLE
' Prob(Khi2 > X) for the Khi-2 distribution with Nu d.o.f.

DECLARE FUNCTION InvKhi2(Nu AS INTEGER, P AS DOUBLE) AS DOUBLE
' Returns X such that FKhi2(Nu, X) = P

' ******************************************************************
'                     Exponential distribution
' ******************************************************************

DECLARE FUNCTION DExpo(A AS DOUBLE, X AS DOUBLE) AS DOUBLE
' Density of exponential distribution with parameter A

DECLARE FUNCTION FExpo(A AS DOUBLE, X AS DOUBLE) AS DOUBLE
' Prob(Z <= X) for exponential dist. with parameter A

' ******************************************************************
'                        Beta distribution
' ******************************************************************

DECLARE FUNCTION DBeta(A AS DOUBLE, B AS DOUBLE,_
                       X AS DOUBLE) AS DOUBLE
' Density of Beta distribution with parameters A and B

DECLARE FUNCTION FBeta(A AS DOUBLE, B AS DOUBLE,_
                       X AS DOUBLE) AS DOUBLE
' Prob(Z <= X) for Beta dist. with parameters A and B

' ******************************************************************
'                        Gamma distribution
' ******************************************************************

DECLARE FUNCTION DGamma(A AS DOUBLE, B AS DOUBLE,_
                        X AS DOUBLE) AS DOUBLE
' Density of Gamma dist. with parameters A and B

DECLARE FUNCTION FGamma(A AS DOUBLE, B AS DOUBLE,_
                        X AS DOUBLE) AS DOUBLE
' Prob(Z <= X) for Gamma dist. with parameters A and B

' ******************************************************************
'                  Matrices and linear equations
' ******************************************************************

DECLARE SUB GaussJordan(A() AS DOUBLE,_
                        BYREF Det AS DOUBLE)
' Transforms matrix A by the Gauss-Jordan method

DECLARE SUB LinEq(A() AS DOUBLE,_
                  B() AS DOUBLE,_
                  BYREF Det AS DOUBLE)
' Solves linear system A * X = B by Gauss-Jordan method

DECLARE SUB LU_Decomp(A() AS DOUBLE)
' LU decomposition

DECLARE SUB LU_Solve(A() AS DOUBLE, B() AS DOUBLE, X() AS DOUBLE)
' Solution of linear system from LU decomposition

DECLARE SUB QR_Decomp(A() AS DOUBLE, R() AS DOUBLE)
' QR decomposition

DECLARE SUB QR_Solve(Q() AS DOUBLE, R() AS DOUBLE, _
                     B() AS DOUBLE, X() AS DOUBLE)
' Solution of linear system from QR decomposition

DECLARE SUB SV_Decomp(A() AS DOUBLE, S() AS DOUBLE, V() AS DOUBLE)
' Singular value decomposition

DECLARE SUB SV_Solve(U() AS DOUBLE, S() AS DOUBLE, V() AS DOUBLE, _
                     B() AS DOUBLE, X() AS DOUBLE)
' Solution of linear system from SVD

DECLARE SUB SV_SetZero(S() AS DOUBLE, Tol AS DOUBLE)

' Set singular values to zero with relative tolerance Tol

DECLARE SUB SV_Approx(U() AS DOUBLE, S() AS DOUBLE, _
                      V() AS DOUBLE, A() AS DOUBLE)
' Matrix approximation from SVD

DECLARE SUB Cholesky(A() AS DOUBLE, L() AS DOUBLE)
' Cholesky factorization of a positive definite symmetric matrix

DECLARE SUB Jacobi(A()      AS DOUBLE,_
                   MaxIter  AS INTEGER,_
                   Tol      AS DOUBLE,_
                   Lambda() AS DOUBLE,_
                   V()      AS DOUBLE)
' Eigenvalues and eigenvectors of a symmetric matrix

DECLARE SUB EigenVals(A()      AS DOUBLE,_
                      Lambda() AS Complex)
' Eigenvalues of a general square matrix

DECLARE SUB EigenVect(A()      AS DOUBLE,_
                      Lambda() AS complex,_
                      V()      AS DOUBLE)
' Eigenvalues and eigenvectors of a general square matrix

' ******************************************************************
'                         Optimization
' ******************************************************************

DECLARE SUB MinBrack(      Func AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                     BYREF A    AS DOUBLE,_
                     BYREF B    AS DOUBLE,_
                     BYREF C    AS DOUBLE,_
                     BYREF Fa   AS DOUBLE,_
                     BYREF Fb   AS DOUBLE,_
                     BYREF Fc   AS DOUBLE)
' Brackets a minimum of function Func between 3 points A, B, C

DECLARE SUB GoldSearch(      Func    AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                             A       AS DOUBLE,_
                             B       AS DOUBLE,_
                             MaxIter AS INTEGER,_
                             Tol     AS DOUBLE,_
                       BYREF Xmin    AS DOUBLE,_
                       BYREF Ymin    AS DOUBLE)
' Minimizes function Func by the golden search method

DECLARE SUB LinMin(      Func     AS FUNCTION(X() AS DOUBLE) AS DOUBLE,_
                         X()      AS DOUBLE,_
                         DeltaX() AS DOUBLE,_
                   BYREF R        AS DOUBLE,_
                         MaxIter  AS INTEGER,_
                         Tol      AS DOUBLE,_
                   BYREF Fmin     AS DOUBLE)
' Minimizes function Func from X() in the direction specified by DeltaX()

DECLARE SUB NumGrad(Func AS FUNCTION(X() AS DOUBLE) AS DOUBLE,_
                    X()  AS DOUBLE,_
                    Eta  AS DOUBLE,_
                    G()  AS DOUBLE)
' Numerical Gradient

DECLARE SUB NumHessGrad(Func AS FUNCTION(X() AS DOUBLE) AS DOUBLE,_
                        X()  AS DOUBLE,_
                        Eta  AS DOUBLE,_
                        G()  AS DOUBLE,_
                        H()  AS DOUBLE)
' Numerical Gradient and Hessian

DECLARE SUB Newton(      Func     AS FUNCTION(X() AS DOUBLE) AS DOUBLE,_
                         HessGrad AS SUB(X() AS DOUBLE, G() AS DOUBLE, H() AS DOUBLE),_
                         X()      AS DOUBLE,_
                         MaxIter  AS INTEGER,_
                         Tol      AS DOUBLE,_
                   BYREF Fmin     AS DOUBLE,_
                         G()      AS DOUBLE,_
                         Hinv()   AS DOUBLE,_
                   BYREF Det      AS DOUBLE)
' Minimizes a function of several variables by Newton's method.

DECLARE SUB Marquardt(      Func     AS FUNCTION(X() AS DOUBLE) AS DOUBLE,_
                            HessGrad AS SUB(X() AS DOUBLE, G() AS DOUBLE, H() AS DOUBLE),_
                            X()      AS DOUBLE,_
                            MaxIter  AS INTEGER,_
                            Tol      AS DOUBLE,_
                      BYREF Fmin     AS DOUBLE,_
                            G()      AS DOUBLE,_
                            Hinv()   AS DOUBLE,_
                      BYREF Det      AS DOUBLE)
' Minimizes a function of several variables by Marquardt's method.

DECLARE SUB BFGS(      Func     AS FUNCTION(X() AS DOUBLE) AS DOUBLE,_
                       Gradient AS SUB(X() AS DOUBLE, G() AS DOUBLE),_
                       X()      AS DOUBLE,_
                       MaxIter  AS INTEGER,_
                       Tol      AS DOUBLE,_
                 BYREF Fmin     AS DOUBLE,_
                       G()      AS DOUBLE,_
                       Hinv()   AS DOUBLE)
' Minimizes a function of several variables by the BFGS method.

DECLARE SUB Simplex(      Func    AS FUNCTION(X() AS DOUBLE) AS DOUBLE,_
                          X()     AS DOUBLE,_
                          MaxIter AS INTEGER,_
                          Tol     AS DOUBLE,_
                    BYREF Fmin    AS DOUBLE)
' Minimizes a function of several variables by the simplex method.

DECLARE FUNCTION CompVec(X()    AS DOUBLE,_
                         Xref() AS DOUBLE,_
                         Tol    AS DOUBLE) AS INTEGER
' Compares two vectors

' ******************************************************************
'                      Nonlinear equations
' ******************************************************************

DECLARE SUB RootBrack(      Func AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                      BYREF X    AS DOUBLE,_
                      BYREF Y    AS DOUBLE,_
                      BYREF FX   AS DOUBLE,_
                      BYREF FY   AS DOUBLE)
' Brackets a root of function Func between X and Y

DECLARE SUB Bisect(      Func    AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                   BYREF X       AS DOUBLE,_
                   BYREF Y       AS DOUBLE,_
                         MaxIter AS INTEGER,_
                         Tol     AS DOUBLE,_
                   BYREF F       AS DOUBLE)
' Bisection method for nonlinear equation

DECLARE SUB Secant(      Func    AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                   BYREF X       AS DOUBLE,_
                   BYREF Y       AS DOUBLE,_
                         MaxIter AS INTEGER,_
                         Tol     AS DOUBLE,_
                   BYREF F       AS DOUBLE)
' Secant method for nonlinear equation

DECLARE SUB NewtEq(      Func    AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                         Deriv   AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                   BYREF X       AS DOUBLE,_
                         MaxIter AS INTEGER,_
                         Tol     AS DOUBLE,_
                   BYREF F       AS DOUBLE)
' Newton-Raphson solver for nonlinear equation

DECLARE SUB NumJacobian(Equations AS SUB(X() AS DOUBLE, F() AS DOUBLE),_
                        X()       AS DOUBLE,_
                        Eta       AS DOUBLE,_
                        D()       AS DOUBLE)
' Numerical Jacobian

DECLARE SUB NewtEqs(Equations AS SUB(X() AS DOUBLE, F() AS DOUBLE),_
                    Jacobian  AS SUB(X() AS DOUBLE, D() AS DOUBLE),_
                    X()       AS DOUBLE,_
                    MaxIter   AS INTEGER,_
                    Tol       AS DOUBLE,_
                    F()       AS DOUBLE)
' Newton-Raphson solver for system of equations

DECLARE SUB Broyden(Equations AS SUB(X() AS DOUBLE, F() AS DOUBLE),_
                    X()       AS DOUBLE,_
                    MaxIter   AS INTEGER,_
                    Tol       AS DOUBLE,_
                    F()       AS DOUBLE)
' Broyden's solver for system of equations

' ******************************************************************
'                        Polynomials
' ******************************************************************

DECLARE FUNCTION Poly(X      AS DOUBLE,_
                      Coef() AS DOUBLE,_
                      Deg    AS INTEGER) AS DOUBLE
' Evaluates a polynomial

DECLARE FUNCTION RFrac(X      AS DOUBLE,_
                       Coef() AS DOUBLE,_
                       Deg1   AS INTEGER,_
                       Deg2   AS INTEGER) AS DOUBLE
' Evaluates a rational fraction

DECLARE FUNCTION RootPol1(      A AS DOUBLE,_
                                B AS DOUBLE,_
                          BYREF X AS DOUBLE) AS INTEGER
' Linear equation

DECLARE FUNCTION RootPol2(Coef() AS DOUBLE,_
                          Z()    AS Complex) AS INTEGER
' Quadratic equation

DECLARE FUNCTION RootPol3(Coef() AS DOUBLE,_
                          Z()    AS Complex) AS INTEGER
' Cubic equation

DECLARE FUNCTION RootPol4(Coef() AS DOUBLE,_
                          Z()    AS Complex) AS INTEGER
' Quartic equation

DECLARE FUNCTION RootPol(Coef() AS DOUBLE,_
                         Deg    AS INTEGER,_
                         Z()    AS Complex) AS INTEGER
' Roots of polynomial from the companion matrix

DECLARE FUNCTION SetRealRoots(Deg AS INTEGER,_
                              Z() AS Complex,_
                              Tol AS DOUBLE) AS INTEGER
' Set imaginary parts to 0

DECLARE SUB SortRoots(Deg AS INTEGER, Z() AS Complex)
' Sort the roots of a polynomial

' ******************************************************************
'                        Numerical integration
' ******************************************************************

DECLARE FUNCTION GausLeg(Func AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                         A    AS DOUBLE,_
                         B    AS DOUBLE) AS DOUBLE
' Integral of function Func from A to B by Gauss-Legendre method

DECLARE FUNCTION GausLeg0(Func AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                          B    AS DOUBLE) AS DOUBLE
' Integral of function Func from 0 to B by Gauss-Legendre method

DECLARE FUNCTION Convol(Func1 AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                        Func2 AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                        T     AS DOUBLE) AS DOUBLE
' Convolution product of two functions Func1 and Func2 at time T

DECLARE FUNCTION TrapInt(X() AS DOUBLE, Y() AS DOUBLE) AS DOUBLE
' Integral of a curve (X, Y) by trapezoidal rule

DECLARE SUB ConvTrap(Func1 AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                     Func2 AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
                     T()   AS DOUBLE,_
                     Y()   AS DOUBLE)
' Convolution product of 2 functions by the trapezoidal rule

DECLARE SUB RKF45(      F      AS SUB (T AS DOUBLE, Y() AS DOUBLE, Yp() AS DOUBLE), _
                        Y()    AS DOUBLE, _
                        Yp()   AS DOUBLE, _
                  BYREF T      AS DOUBLE, _
                        Tout   AS DOUBLE, _
                        RelErr AS DOUBLE, _
                        AbsErr AS DOUBLE, _
                  BYREF Flag   AS INTEGER)
' Integration of a system of differential equations

' ******************************************************************
'                     Fast Fourier Transform
' ******************************************************************

DECLARE SUB FFT(InArray() AS Complex, OutArray() AS Complex)
' Computes the FFT of data in InArray - Result in OutArray

DECLARE SUB IFFT(InArray() AS Complex, OutArray() AS Complex)
' Computes the inverse FFT of data in InArray - Result in OutArray

DECLARE FUNCTION CalcFrequency(FrequencyIndex AS INTEGER, _
                               InArray() AS Complex) AS Complex
' Returns the FFT for the frequency specified by FrequencyIndex

' ******************************************************************
'                         Random Numbers
' ******************************************************************

DECLARE FUNCTION IRanGen() AS UINTEGER
' Returns a 32-bit random integer in [0, 2^32-1)

DECLARE FUNCTION IRanGen31() AS LONG
' Returns a 31-bit random integer

DECLARE FUNCTION RanGen1() AS DOUBLE
' Returns a random number on [0,1]-real-interval

DECLARE FUNCTION RanGen2() AS DOUBLE
' Returns a random number on [0,1)-real-interval

DECLARE FUNCTION RanGen3() AS DOUBLE
' Returns a random number on (0,1)-real-interval

DECLARE FUNCTION RanGen53() AS DOUBLE
' Returns a random number on [0,1) with 53-bit resolution

' ******************************************************************
'               Gaussian Random Number Generators
' ******************************************************************

DECLARE FUNCTION RanGaussStd AS DOUBLE
' Random number from the standard normal distribution

DECLARE FUNCTION RanGauss(Mu AS DOUBLE, Sigma AS DOUBLE) AS DOUBLE
' Random number from normal dist (mean Mu, SD Sigma)

DECLARE SUB RanMult(M() AS DOUBLE, L() AS DOUBLE, X() AS DOUBLE)
' Random vector X from multinormal dist. (mean M, Cholesky factor L)

DECLARE SUB RanMultIndep(M() AS DOUBLE, S() AS DOUBLE, X() AS DOUBLE)
' Random vector X from multinormal dist. (uncorrelated)

' ******************************************************************
'                     Markov Chain Monte-Carlo
' ******************************************************************

DECLARE SUB InitMHParams(NCycles  AS INTEGER, _
                         MaxSim   AS INTEGER, _
                         SavedSim AS INTEGER)
' Initialize Metropolis-Hastings parameters

DECLARE SUB Hastings(      Func   AS FUNCTION(X() AS DOUBLE) AS DOUBLE, _
                           T      AS DOUBLE, _
                           X()    AS DOUBLE, _
                           V()    AS DOUBLE, _
                           Xmat() AS DOUBLE, _
                           Xmin() AS DOUBLE, _
                     BYREF Fmin   AS DOUBLE)
' Simulation of probability density by Metropolis-Hastings algorithm

' ******************************************************************
'                       Simulated Annealing
' ******************************************************************

DECLARE SUB InitSAParams(NT      AS INTEGER = 5,   _
                         NS      AS INTEGER = 15,  _
                         RT      AS DOUBLE  = 0.9, _
                         NCycles AS INTEGER = 1)
' Initialize Simulated Annealing parameters

DECLARE SUB SA_CreateLogFile(LogFileName AS STRING = "simann.txt")
' Initialize log file

DECLARE SUB SimAnn(      Func   AS FUNCTION(X() AS DOUBLE) AS DOUBLE, _
                         X()    AS DOUBLE, _
                         Xmin() AS DOUBLE, _
                         Xmax() AS DOUBLE, _
                   BYREF Fmin   AS DOUBLE)
' Minimization of a function of several variables by simulated annealing

' ******************************************************************
'                         Genetic Algorithm
' ******************************************************************

DECLARE SUB InitGAParams(NP AS INTEGER = 200, _
                         NG AS INTEGER = 40,  _
                         SR AS DOUBLE  = 0.5, _
                         MR AS DOUBLE  = 0.1, _
                         HR AS DOUBLE  = 0.5)
' Initialize Genetic Algorithm parameters

DECLARE SUB GA_CreateLogFile(LogFileName AS STRING = "genalg.txt")
' Initialize log file

DECLARE SUB GenAlg(      Func   AS FUNCTION(X() AS DOUBLE) AS DOUBLE, _
                         X()    AS DOUBLE, _
                         Xmin() AS DOUBLE, _
                         Xmax() AS DOUBLE, _
                   BYREF Fmin   AS DOUBLE)
' Minimization of a function of several variables by genetic algorithm

' ******************************************************************
'                            Statistics
' ******************************************************************

TYPE StatClass    ' Statistical class
  Inf AS DOUBLE   ' Lower bound
  Sup AS DOUBLE   ' Upper bound
  N   AS INTEGER  ' Number of values
  F   AS DOUBLE   ' Frequency
  D   AS DOUBLE   ' Density
END TYPE

DECLARE SUB QuickSort(SortArray() AS DOUBLE, _
                      Low         AS INTEGER, _
                      High        AS INTEGER)
' Sort array

DECLARE FUNCTION Median(X()    AS DOUBLE, _
                        Sorted AS INTEGER = True) AS DOUBLE
' Median of sample X

DECLARE FUNCTION Mean(X() AS DOUBLE) AS DOUBLE
' Mean of sample X

DECLARE FUNCTION StDev(X() AS DOUBLE, M AS DOUBLE) AS DOUBLE
' Standard deviation estimated from sample X

DECLARE FUNCTION StDevP(X() AS DOUBLE, M AS DOUBLE) AS DOUBLE
' Standard deviation of population

DECLARE FUNCTION Skewness(X() AS DOUBLE, _
                          M AS DOUBLE, Sigma AS DOUBLE) AS DOUBLE
' Skewness of sample X

DECLARE FUNCTION Kurtosis(X() AS DOUBLE, _
                          M AS DOUBLE, Sigma AS DOUBLE) AS DOUBLE
' Kurtosis of sample X

DECLARE SUB Interval(      X1     AS DOUBLE,  _
                           X2     AS DOUBLE,  _
                           MinSub AS INTEGER, _
                           MaxSub AS INTEGER, _
                     BYREF XMin   AS DOUBLE,  _
                     BYREF XMax   AS DOUBLE,  _
                     BYREF XStep  AS DOUBLE)
' Compute an appropriate interval for a set of values

DECLARE SUB StudIndep(      N1  AS INTEGER, N2 AS INTEGER, _
                            M1  AS DOUBLE,  M2 AS DOUBLE,  _
                            S1  AS DOUBLE,  S2 AS DOUBLE,  _
                      BYREF T   AS DOUBLE,                 _
                      BYREF DoF AS INTEGER)
' Student's test for independent samples

DECLARE SUB StudPaired(      X() AS DOUBLE, _
                             Y() AS DOUBLE, _
                       BYREF T   AS DOUBLE, _
                       BYREF DoF AS INTEGER)
' Student's test for paired samples

DECLARE SUB AnOVa1(      N()   AS INTEGER, _
                         M()   AS DOUBLE,  _
                         S()   AS DOUBLE,  _
                   BYREF V_f   AS DOUBLE,  _
                   BYREF V_r   AS DOUBLE,  _
                   BYREF F     AS DOUBLE,  _
                   BYREF DoF_f AS INTEGER, _
                   BYREF DoF_r AS INTEGER)
' One-way analysis of variance

DECLARE SUB AnOVa2(N     AS INTEGER, _
                   M()   AS DOUBLE,  _
                   S()   AS DOUBLE,  _
                   V()   AS DOUBLE,  _
                   F()   AS DOUBLE,  _
                   DoF() AS INTEGER)
' Two-way analysis of variance

DECLARE SUB Snedecor(      N1   AS INTEGER, N2 AS INTEGER, _
                           S1   AS DOUBLE,  S2 AS DOUBLE,  _
                     BYREF F    AS DOUBLE,                 _
                     BYREF DoF1 AS INTEGER,                _
                     BYREF DoF2 AS INTEGER)
' Snedecor's test (comparison of 2 variances)

DECLARE SUB Bartlett(      N()  AS INTEGER, _
                           S()  AS DOUBLE,  _
                     BYREF Khi2 AS DOUBLE,  _
                     BYREF Dof  AS INTEGER)
' Bartlett's test (comparison of several variances)

DECLARE SUB Mann_Whitney(      X1() AS DOUBLE, _
                               X2() AS DOUBLE, _
                         BYREF U    AS DOUBLE, _
                         BYREF Eps  AS DOUBLE)
' Mann-Whitney test

DECLARE SUB Wilcoxon(      X()   AS DOUBLE,  _
                           Y()   AS DOUBLE,  _
                           Ndiff AS INTEGER, _
                     BYREF T     AS DOUBLE,  _
                     BYREF Eps   AS DOUBLE)
' Wilcoxon test

DECLARE SUB Kruskal_Wallis(      N() AS INTEGER, _
                                 X() AS DOUBLE,  _
                           BYREF H   AS DOUBLE,  _
                           BYREF DoF AS INTEGER)
' Kruskal-Wallis test

DECLARE SUB Khi2_Conform(      Obs()   AS INTEGER, _
                               Calc()  AS DOUBLE,  _
                               N_estim AS INTEGER, _
                         BYREF Khi2    AS DOUBLE,  _
                         BYREF DoF     AS INTEGER)
' Khi-2 test for conformity

DECLARE SUB Khi2_Indep(      Obs() AS INTEGER, _
                       BYREF Khi2  AS DOUBLE,  _
                       BYREF DoF   AS INTEGER)
' Khi-2 test for independence

DECLARE SUB Woolf_Conform(      Obs()   AS INTEGER, _
                                Calc()  AS DOUBLE,  _
                                N_estim AS INTEGER, _
                          BYREF G       AS DOUBLE,  _
                          BYREF DoF     AS INTEGER)
' Woolf test for conformity

DECLARE SUB Woolf_Indep(      Obs() AS INTEGER, _
                        BYREF G     AS DOUBLE,  _
                        BYREF DoF   AS INTEGER)
' Woolf test for independence

DECLARE SUB Distrib(      X()  AS DOUBLE,  _
                          A    AS DOUBLE,  _
                          H    AS DOUBLE,  _
                    BYREF Ntot AS INTEGER, _
                          C()  AS StatClass)
' Generates a statistical distribution

' ******************************************************************
'                            Curve fitting
' ******************************************************************

TYPE TRegTest
  Vr  AS DOUBLE   ' Residual variance
  R2  AS DOUBLE   ' Coefficient of determination
  R2a AS DOUBLE   ' Adjusted coeff. of determination
  F   AS DOUBLE   ' Variance ratio (explained/residual)
  Nu1 AS INTEGER  ' Number of d.o.f. (numerator)
  Nu2 AS INTEGER  ' Number of d.o.f. (denominator)
END TYPE

DECLARE SUB LinFit(X() AS DOUBLE, _
                   Y() AS DOUBLE, _
                   B() AS DOUBLE, _
                   V() AS DOUBLE)
' Linear regression: Y = B(0) + B(1) * X

DECLARE SUB WLinFit(X() AS DOUBLE, _
                    Y() AS DOUBLE, _
                    S() AS DOUBLE, _
                    B() AS DOUBLE, _
                    V() AS DOUBLE)
' Weighted linear regression: Y = B(0) + B(1) * X

DECLARE SUB MulFit(X() AS DOUBLE, _
                   Y() AS DOUBLE, _
                   B() AS DOUBLE, _
                   V() AS DOUBLE)
' Multiple linear regression

DECLARE SUB WMulFit(X() AS DOUBLE, _
                    Y() AS DOUBLE, _
                    S() AS DOUBLE, _
                    B() AS DOUBLE, _
                    V() AS DOUBLE)
' Weighted multiple linear regression

DECLARE SUB PolFit(X() AS DOUBLE,  _
                   Y() AS DOUBLE,  _
                   B() AS DOUBLE,  _
                   V() AS DOUBLE)
' Polynomial regression

DECLARE SUB WPolFit(X() AS DOUBLE, _
                    Y() AS DOUBLE, _
                    S() AS DOUBLE, _
                    B() AS DOUBLE, _
                    V() AS DOUBLE)
' Weighted polynomial regression

DECLARE SUB RegTest(      Y()     AS DOUBLE, _
                          Ycalc() AS DOUBLE, _
                          V()     AS DOUBLE, _
                    BYREF Test    AS TRegTest)
' Test of unweighted regression

DECLARE SUB WRegTest(      Y()     AS DOUBLE, _
                           Ycalc() AS DOUBLE, _
                           S()     AS DOUBLE, _
                           V()     AS DOUBLE, _
                     BYREF Test    AS TRegTest)
' Test of weighted regression

DECLARE SUB VecMean(X() AS DOUBLE, M() AS DOUBLE)
' Mean vector M() from matrix X()

DECLARE SUB VecSD(X() AS DOUBLE, M() AS DOUBLE, S() AS DOUBLE)
' Vector of standard deviations S() from matrix X() and mean vector M()

DECLARE SUB MatVarCov(X() AS DOUBLE, M() AS DOUBLE, V() AS DOUBLE)
' Variance-covariance matrix V() from matrix X() and mean vector M()

DECLARE SUB MatCorrel(V() AS DOUBLE, R() AS DOUBLE)
' Correlation matrix R() from variance-covariance matrix V()

DECLARE SUB PCA(R() AS DOUBLE, Lambda() AS DOUBLE, C() AS DOUBLE, Rc() AS DOUBLE)
' Principal component analysis of correlation matrix

DECLARE SUB ScaleVar(X() AS DOUBLE, M() AS DOUBLE, S() AS DOUBLE, Z() AS DOUBLE)
' Scales a set of variables by subtracting means and dividing by SD

DECLARE SUB PrinFac(Z() AS DOUBLE, C() AS DOUBLE, F() AS DOUBLE)
' Computes principal factors

' ******************************************************************
'                     Graphics constants and type
' ******************************************************************

#DEFINE LinScale  0     ' Linear scale
#DEFINE LogScale  1     ' Logarithmic scale
#DEFINE MaxDec   10     ' Max. number of decimal places for labels

#DEFINE HorizGrid 1     ' Horizontal grid
#DEFINE VertiGrid 2     ' Vertical grid
#DEFINE BothGrid  3     ' Horizontal + vertical

TYPE TCurvParam         ' Curve plot parameters
  Col     AS INTEGER    ' Color
  Symbol  AS INTEGER    ' Point type
  Size    AS INTEGER    ' Point size
  Connect AS INTEGER    ' Points connected ?
  Text    AS STRING*20  ' Curve legend
END TYPE

' ******************************************************************
'                         Scale of graphics
' ******************************************************************

DECLARE SUB AutoScale(      X()   AS DOUBLE,  _
                            Scale AS INTEGER, _
                      BYREF XMin  AS DOUBLE,  _
                      BYREF XMax  AS DOUBLE,  _
                      BYREF XStep AS DOUBLE)
' Automatic scale determination

' ******************************************************************
'                           Screen Graphics
' ******************************************************************

DECLARE FUNCTION InitScreen(Mode    AS INTEGER = 12, _
                            FgColor AS INTEGER = 0, _
                            BgColor AS INTEGER = 15) AS INTEGER
' Initializes the graphic mode

DECLARE SUB InitGraph(X1     AS INTEGER = 15, _
                      X2     AS INTEGER = 85, _
                      Y1     AS INTEGER = 15, _
                      Y2     AS INTEGER = 85, _
                      Border AS INTEGER = -1)
' Initializes the graphic window

DECLARE SUB SetOxScale(Scale  AS INTEGER = LinScale, _
                       OxMin  AS DOUBLE  = 0, _
                       OxMax  AS DOUBLE  = 1, _
                       OxStep AS DOUBLE  = 0.2)
' Sets the scale on the Ox axis

DECLARE SUB SetOyScale(Scale  AS INTEGER = LinScale, _
                       OyMin  AS DOUBLE  = 0, _
                       OyMax  AS DOUBLE  = 1, _
                       OyStep AS DOUBLE  = 0.2)
' Sets the scale on the Oy axis

DECLARE FUNCTION Xpixel(X AS DOUBLE) AS INTEGER
' Converts user abscissa X to screen coordinate

DECLARE FUNCTION Ypixel(Y AS DOUBLE) AS INTEGER
' Converts user ordinate Y to screen coordinate

DECLARE SUB PlotOxAxis(Y AS DOUBLE = 0, Dec AS INTEGER = 2)
' Plots the Ox axis

DECLARE SUB PlotOyAxis(X AS DOUBLE = 0, Dec AS INTEGER = 2)
' Plots the Oy axis

DECLARE SUB PlotGrid(Grid AS INTEGER = BothGrid)
' Plots a grid

DECLARE SUB WriteTitles(GraphTitle AS STRING, _
                        XTitle     AS STRING, _
                        YTitle     AS STRING)
' Writes graph and axes titles

DECLARE SUB PlotPoint(X      AS INTEGER, _
                      Y      AS INTEGER, _
                      Symbol AS INTEGER = 1, _
                      Size   AS INTEGER = 2)
' Plots a point

DECLARE SUB PlotCurve(X()     AS DOUBLE, _
                      Y()     AS DOUBLE, _
                      Symbol  AS INTEGER = 1, _
                      Size    AS INTEGER = 2, _
                      Connect AS INTEGER = 0)
' Plots a set of points

DECLARE SUB PlotCurveWithErrorBars(X()     AS DOUBLE, _
                                   Y()     AS DOUBLE, _
                                   S()     AS DOUBLE, _
                                   Ns      AS INTEGER = 1, _
                                   Symbol  AS INTEGER = 1, _
                                   Size    AS INTEGER = 2, _
                                   Connect AS INTEGER = 0)
' Plots a set of points with error bars

DECLARE SUB ConRec(X()   AS INTEGER, _
                   Y()   AS INTEGER, _
                   Z()   AS DOUBLE, _
                   D()   AS DOUBLE, _
                   Col() AS INTEGER)
' Contour plot

DECLARE SUB SetClipping
' Limits the plots at the current viewport boundaries

DECLARE SUB PlotFunc(Func AS FUNCTION(X AS DOUBLE) AS DOUBLE)
' Plots a function

DECLARE SUB WriteLegend(CurvParam() AS TCurvParam)
' Writes curve legends

DECLARE SUB CloseGraph
' Quits the graphic mode

' ******************************************************************
'                         String functions
' ******************************************************************

DECLARE FUNCTION RFill(S AS STRING, L AS INTEGER) AS STRING
' Completes string S with trailing blanks for a total length L

DECLARE FUNCTION LFill(S AS STRING, L AS INTEGER) AS STRING
' Completes string S with leading blanks for a total length L

DECLARE FUNCTION CFill(S AS STRING, L AS INTEGER = 80) AS STRING
' Completes string S with leading blanks to center the string

DECLARE SUB Replace(S AS STRING, C1 AS STRING, C2 AS STRING)
' Replaces in string S all the occurences of character C1 by character C2

DECLARE SUB Parse(S AS STRING, Delim AS STRING, Elem() AS STRING, BYREF N AS INTEGER)
' Parses a string into its constitutive fields (separated by Delim)

DECLARE FUNCTION FloatToStr(X    AS DOUBLE,       _
                            Ntot AS INTEGER = 10, _
                            Ndec AS INTEGER = 4,  _
                            E    AS INTEGER = 0) AS STRING
' Converts a floating point number to a string

DECLARE FUNCTION RemZero(S AS STRING) AS STRING
' Remove non-significant zeroes

DECLARE FUNCTION IntToStr(N AS INTEGER, Ntot AS INTEGER = 10) AS STRING
' Converts an integer to a string
