© 2010 W.Ehrhardt Last update Feb. 24, 2010

Basic information about the AMath package

Contents

Introduction

The AMath unit implements accurate mathematical functions, it makes many routines of Delphi's math unit available to other supported Pascal versions and fixes bugs and inaccuracies of Delphi.

The elementary mathematical functions include: exponential, logarithmic, trigonometric, hyperbolic, inverse circular / hyperbolic functions. Then there are polynomial, vector, statistic operations as well as floating point and FPU control functions.

All one argument elementary transcendental functions have peak relative errors less than 1.9e-19, the value for power(x,y) is 2.8e-18 (for |x|,|y| < 1000) and 4.8e-18 (for |x|,|y| < 2000).

Unit QuadSolv accurately solves quadratic equations with double coefficients. The functions implement ideas of G.E. Forsythe, W. Kahan, P.H. Sterbenz (high precision calculation of discriminant, scaling by powers of two etc).

The units SpecFun, SpecFunX, SFCommon implement special functions for extended and double precision; the following function groups are available: gamma and related, AGM and Jacobian elliptic functions, elliptic integrals, error function with standard normal distribution, exponential integrals and logarithmic integral.

Currently all functions have double and extended versions (with name suffix x), eg. erfc vs. erfcx. Generally the extended versions have larger relative errors (measured in corresponding machine epsilons eps_x or eps_d) than their double counterparts, especially gammax and betax for large arguments. Note that some functions are sensitive to small changes in the argument; therefore in high precision comparisons argument values should be used, that are representable in both calculations.


AMath functions

Elementary numerical functions
cbrt               Return the cube root of x
ceil               Return the smallest integer >= x; |x| <= MaxLongint
ceilx              Return the smallest integer >= x
floor              Return the largest integer <= x; |x| <= MaxLongint
floorx             Return the largest integer <= x
fmod               Return x mod y, y<>0, sign(result) = sign(x)
hypot              Return sqrt(x*x + y*y)
intpower           Return x^n; via binary exponentiation
modf               Return frac(x) and trunc(x) in ip, |x| <= MaxLongint
nroot              Return the nth root of x; n<>0, x >= 0 if n is even
sqrt1pm1           Return sqrt(1+x)-1 with increased accuracy for x near 0, x>=-1
Elementary transcendental functions
arccos             Return the inverse circular cosine of x, |x| <= 1
arccosh            Return the inverse hyperbolic cosine, x >= 1. For x near 1 use arccosh1p(x-1) to reduce cancellation errors!
arccosh1p          Calculate arccosh(1+x), x>=0, with increased accuracy for x near 0
arccot             Return the sign symmetric inverse circular cotangent; arccot(x) = arctan(1/x), x <> 0
arccotc            Return the continuous inverse circular cotangent; arccotc(x) = Pi/2 - arctan(x)
arccoth            Return the inverse hyperbolic cotangent of x, |x| > 1
arccsc             Return the inverse cosecant of x, |x| >= 1
arccsch            Return the inverse hyperbolic cosecant of x, x <> 0
arcsec             Return the inverse secant of x, |x| >= 1
arcsech            Return the inverse hyperbolic secant of x, 0 < x <= 1
arcsin             Return the inverse circular sine of x, |x| <= 1
arcsinh            Return the inverse hyperbolic sine of x
arctan2            Return arctan(y/x); result in [-Pi..Pi] with correct quadrant
arctanh            Return the inverse hyperbolic tangent of x, |x| < 1
cos                Accurate version of circular cosine, uses system.cos for |x| <= Pi/4
cosh               Return the hyperbolic cosine of x
cosPi              Return cos(Pi*x), result will be 1 for abs(x) >= 2^64
cot                Return the circular cotangent of x, x mod Pi <> 0
coth               Return the hyperbolic cotangent of x, x<>0
csc                Return the circular cosecant of x, x mod Pi <> 0
csch               Return the hyperbolic cosecant of x, x<>0
exp                Accurate exp, result good to extended precision
exp10              Return 10^x
exp2               Return 2^x
expm1              Return exp(x)-1 with increased accuracy for x near 0
expx2              Return exp(x*|x|) with damped error amplification in computing exp of the product
ln1p               Return ln(1+x) with increased accuracy for x near 0
lnxp1              Delphi alias for ln1p
log10              Return base 10 logarithm of x
log2               Return base 2 logarithm of x
logbase            Return base b logarithm of x
logN               Delphi alias for logbase
power              Return x^y; if frac(y)<>0 then x must be > 0
sec                Return the circular secant of x, x mod Pi <> Pi/2
sech               Return the hyperbolic secant of x
sin                Accurate version of circular sine, uses system.sin for |x| <= Pi/4
sincos             Return accurate values s=sin(x), c=cos(x)
sincosPi           Return s=sin(Pi*x), c=cos(Pi*x); (s,x)=(0,1) for abs(x) >= 2^64
sinh               Return the hyperbolic sine of x with increased accuracy for x near 0
sinPi              Return sin(Pi*x), result will be 0 for abs(x) >= 2^64
tan                Return the circular tangent of x, x mod Pi <> Pi/2
tanh               Return the hyperbolic tangent of x with increased accuracy for x near 0
Floating point functions
copysign           Return abs(x)*sign(y)
copysignd          Return abs(x)*sign(y)
frexp              Return the mantissa m and exponent e of x with x = m*2^e, 0.5 <= abs(m) < 1; if x is 0, +-INF, NaN or denormal, return m=x, e=0
frexpd             Return the mantissa m and exponent e of d with d = m*2^e, 0.5 <= abs(m) < 1; if d is 0, +-INF, NaN or denormal, return m=d, e=0
ilogb              Return base 2 exponent of x, -MaxLongint for x=0/denormal, MaxLongint if INF or Nan
IsInf              Return true if x is +INF or -INF
IsInfD             Return true if d is +INF or -INF
IsNaN              Return true if x is a NaN
IsNaND             Return true if d is a NaN
IsNaNorInf         Return true if x is a NaN or infinite
IsNaNorInfD        Return true if d is a NaN or infinite
ldexp              Return x*2^e
ldexpd             Return d*2^e
predd              Return next representable double after d in the direction -Inf
predx              Return next representable extended after x in the direction -Inf
scalbn             Return x*2^e
succd              Return next representable double after d in the direction +Inf
succx              Return next representable extended after x in the direction +Inf
FPU control functions
Get8087CW          Return the FPC control word
GetExceptionMask   Return the  current exception mask
GetPrecisionMode   Return the current precision control mode
GetRoundMode       Return the current rounding mode
Set8087CW          Set new FPU control word
SetExceptionMask   Set new exception mask
SetPrecisionMode   Set new precision control mode and return the old precision
SetRoundMode       Set new rounding mode and return the old mode
Polynomial, Vector, Statistic Operations
dot2               Accurate dot product sum(x[i]*y[i], i=0..n-1) of two double vectors
dot2x              Accurate dot product sum(x[i]*y[i], i=0..n-1) of two extended vectors
mean               Compute accurate mean = sum(a[i], i=0..n-1)/n of double vector
MeanAndStdDev      Accurate mean and sample standard deviation of double vector
norm2              Calculate the 2-norm = sqrt(sum(a[i]^2, i=0..n-1)) of an double vector
norm2x             Calculate the 2-norm = sqrt(sum(a[i]^2, i=0..n-1)) of an extended vector
PolEval            Evaluate polynomial; return a[0] + a[1]*x + ... + a[n-1]*x^(n-1)
PolEvalX           Evaluate polynomial; return a[0] + a[1]*x + ... + a[n-1]*x^(n-1)
ssqd               Calculate sum(a[i]^2, i=0..n-1) = scale^2*sumsq, scale>=0, sumsq>0
ssqx               Calculate sum(a[i]^2, i=0..n-1) = scale^2*sumsq, scale>=0, sumsq>0
sum2               Compute accurate sum(a[i], i=0..n-1) of double vector
sum2x              Compute accurate sum(a[i], i=0..n-1) of extended vector
sumsqr             Calculate sum(a[i]^2, i=0..n-1) of an double vector
sumsqrx            Calculate sum(a[i]^2, i=0..n-1) of an extended vector
Argument reduction for trigonometric functions
rem_pio2_cw        Cody/Waite  reduction of x:  z = x - n*Pi/2, |z| <= Pi/4, result = n mod 8
rem_pio2_ph        Payne/Hanek reduction of x:  z = x - n*Pi/2, |z| <= Pi/4, result = n mod 8
rem_pio2           Argument    reduction of x:  z = x - n*Pi/2, |z| <= Pi/4, result = n mod 8.
                   Uses Payne/Hanek if |x| > ph_cutoff, Cody/Waite otherwise
rem_int2           Argument reduction of x: z*Pi = x*Pi - n*Pi/2, |z|<=1/4, result = n mod 8.
                   Used for argument reduction in sin(Pi*x) and cos(Pi*x).

Other functions
Dbl2Hex            return d as a big-endian hex string
Ext2Hex            return x as a big-endian hex string
fisEQd             return true if x and y are bit-identical
fisEQx             return true if x and y are bit-identical
fisNEd             return true if x and y are not bit-identical
fisNEx             return true if x and y are not bit-identical

QuadSolv functions

squad(a,b,c: double; var x1,y1,x2,y2: double): Solve the quadratic equation a*x^2 + b*x + c = 0. Result is the number of different solutions: 0 (if a=b=0), 1 (x1), or 2 (x1,x2). If the result is = -2, x1+i*y1 and x2+i*y2 are the two complex solutions. No precautions against over/underflow, NAN/INF coefficients return 0.

squadx(a,b,c: double; var x1,y1,x2,y2: double): Solve the quadratic equation a*x^2 + b*x + c = 0. Result is the number of different solutions: 0 (if a=b=0 or INF/NAN), 1 (x1), or 2 (x1,x2). If the result is = -2, then x1 + i*y1 and x2 + i*y2 are the two complex solutions. Uses scaling by powers of two to minimize over/underflows.

squad_selftest: Selftest for squad core routine, result=0 if OK, else number of the failing test case.


Special functions

Note that the extended versions have x suffixes, eg. erfcx vs. erfc.

AGM and Jacobian elliptic functions
agm(x,y)           Return the arithmetic-geometric mean of |x| and |y|
sncndn(x,mc...)    Calculate the Jacobian elliptic functions sn,cn,dn for argument x and complementary parameter mc.
Elliptic integrals
ell_rc(x,y)        Return Carlson's degenerate elliptic integral RC; x>=0, y<>0
ell_rf(x,y,z)      Return Carlson's elliptic integral of the first kind; x,y,z >=0, at most one =0
ell_rd(x,y,z)      Return Carlson's elliptic integral of the 2nd kind; z>0; x,y >=0, at most one =0
ell_rj(x,y,z,r)    Return Carlson's elliptic integral of the 3rd kind; r<>0; x,y,z >=0, at most one =0

cel1(kc)           Return the complementary complete elliptic integral of the 1st kind, kc<>0
cel2(kc,a,b)       Return the general complete elliptic integral of the 2nd kind, kc<>0
cel(kc,p,a,b)      Return the general complete elliptic integral, kc<>0, Cauchy principle value if p<0
el1(x,kc)          Return the complementary complete elliptic integral of the 1st kind
el2(x,kc,a,b)      Return the general complete elliptic integral of the 2nd kind, kc<>0
el3(x,kc,p)        Return the incomplete elliptic integral of the 3rd kind, 1+p*x^2<>0

EllipticF(z,k)     Return the incomplete elliptic integral of the 1st kind; |z|<=1, |k*z|<1
EllipticK(k)       Return the complete elliptic integral of the 1st kind, |k| < 1
EllipticCK(kc)     Return the complementary complete elliptic integral of the 1st kind, kc<>0
EllipticE(z,k)     Return the incomplete elliptic integrals of the 2nd kind, |z|<=1, |k|<=1
EllipticEC(k)      Return the complete elliptic integral of the 2nd kind, |k| < 1
EllipticCE(kc)     Return the complementary complete elliptic integral of the 2nd kind
EllipticPi(z,nu,k) Return the incomplete elliptic integral of the 3rd kind, |z|<=1, |k*z|<1
EllipticPiC(nu,k)  Return the complete elliptic integral of the 3rd kind, |k|<1
EllipticCPi(nu,k)  Return the complementary complete elliptic integral of the 3rd kind, k<>0
Error function and standard normal distribution
erf(x)             Return the error function erf(x)
erfc(x)            Return complementary error function erfc(x) = 1-erf(x)
normstd_pdf(x)     Return the std. normal probability density function exp(-x^2/2)/sqrt(2*Pi)
normstd_cdf(x)     Return the standard normal distribution function
normstd_inv(y)     Return the inverse standard normal distribution function, 0 < y < 1.
                   For x=normstd_inv(y) and y from (0,1), normstd_cdf(x) ~ y
Exponential integrals and logarithmic integral
e1(x)              Return the exponential integral E_1(x) = integral(exp(-x*t)/t, t=1..Inf), x > 0
ei(x)              Return the exponential integral Ei(x)  = PV-integral(exp(t)/t, t=-Inf..x)
en(n,x)            Return the exponential integral E_n(x) = integral(exp(-x*t)/t^n, t=1..Inf), x > 0
li(x)              Return the logarithmic integral li(x)  = PV-integral(1/ln(t), t=0..x), x > 1
Gamma function and related
gamma(x)           Return gamma(x), |x| <= MAXGAM; invalid if x is non-positive integer
invgamma(x)        Return 1/gamma(x)
signgamma(x)       Return sign(gamma(x)), useless for 0 or negative integer
lngamma(x)         Return ln(|gamma(x)|), |x| <= MAXLGM, invalid if x is non-positive integer.
psi(x)             Return the psi (digamma) function of x, INF if x is non-positive integer
beta(x,y)          Return beta(x,y)=gamma(x)*gamma(y)/gamma(x+y)
lnbeta(x,y)        Return the logarithm of |beta(x,y)|=|gamma(x)*gamma(y)/gamma(x+y)|
zeta(s)            Return the Riemann zeta function at s, s<>1
zeta1p(x)          Return the Riemann zeta function at 1+x, x<>0


References

  1. [HMF]: M. Abramowitz, I.A. Stegun. Handbook of Mathematical Functions. New York, 1970, http://www.math.sfu.ca/~cbm/aands/
  2. Intel, IA-32 Architecture Software Developer's Manual Volume 2A: Instruction Set Reference, A-M
  3. D. Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, 1991, http://docs.sun.com/source/806-3568/ncg_goldberg.html; extended and edited reprint via http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.22.6768
  4. ISO/IEC 10967-2, Information technology: Language independent arithmetic, Part 2: Elementary numerical functions, http://standards.iso.org/ittf/PubliclyAvailableStandards/c024427_ISO_IEC_10967-2_2001(E).zip
  5. FDLIBM 5.3 (Freely Distributable LIBM), developed at Sun Microsystems, see http://www.netlib.org/fdlibm/ or http://www.validlab.com/software/fdlibm53.tar.gz
  6. K.C. Ng, Argument Reduction for Huge Arguments: Good to the Last Bit, Technical report, SunPro, 1992. Available from http://www.validlab.com/arg.pdf
  7. Cephes Mathematical Library, Version 2.8, http://www.moshier.net/#Cephes or http://www.netlib.org/cephes/
  8. T. Ogita, S.M. Rump, and S. Oishi, Accurate sum and dot product, SIAM J. Sci. Comput., 26 (2005), pp. 1955-1988. Available as http://www.ti3.tu-harburg.de/paper/rump/OgRuOi05.pdf
  9. N.J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., Philadelphia, 2002, http://www.maths.manchester.ac.uk/~higham/asna/
  10. R. Bulirsch, Numerical Calculation of Elliptic Integrals and Elliptic Functions, Numerische Mathematik 7, 78-90, 1965
  11. R. Bulirsch, Numerical Calculation of Elliptic Integrals and Elliptic Functions, part III. Numerische Mathematik 13, 305-315, 1969
  12. B.C. Carlson, Computing Elliptic Integrals by Duplication, Numerische Mathematik 33, 1-16, 1979
  13. W.H. Press et al, Numerical Recipes in C, 2nd ed., Cambridge, 1992
  14. SLATEC Common Mathematical Library, Version 4.1, July 1993 (general purpose mathematical and statistical routines written in Fortran 77), http://www.netlib.org/slatec
  15. W.Kahan, On the Cost of Floating-Point Computation Without Extra-Precise Arithmetic, 2004 http://www.eecs.berkeley.edu/~wkahan/Qdrtcs.pdf
  16. G.E. Forsythe, How do you solve a quadratic equation? Stanford University Technical Report no. CS40, 1966 ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/66/40/CS-TR-66-40.pdf
  17. P.H. Sterbenz, Floating-Point Computation. Englewood Cliffs, 1974. Chap.9.3: Carefully written programs/quadratic equation, p.246ff
  18. W.Y. Sit, Quadratic Programming? 1997. http://www.mmrc.iss.ac.cn/ascm/ascm03/sample.pdf
  19. Boost C++ Libraries, Release 1.42.0, 2010. http://www.boost.org/