Basic information about the AMath package

• Introduction
• AMath functions
  • Elementary numerical functions
  • Elementary transcendental functions
  • Floating point functions
  • FPU control functions
  • Polynomial, Vector, Statistic Operations
  • Argument reduction for trigonometric functions
  • Other functions
• QuadSolv functions
• Special functions
  • AGM and Jacobian elliptic functions
  • Elliptic integrals
  • Error function and standard normal distribution
  • Exponential integrals and logarithmic integral
  • Gamma function and related
• References

Introduction

AMath 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

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

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

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

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

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).

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

Special 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.

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

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

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(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/