© 2014 W.Ehrhardt Last update Dec. 07, 2014

Basic information about the AMath/DAMath packages

Contents

Introduction

The AMath package contains Pascal/Delphi source for accurate mathematical methods without using multi precision arithmetic; it is designed for the 80-bit extended data type. DAMath makes most of the AMath functions available for 64-bit systems without extended precision or 387-FPU.

Please note that the high accuracy can only be achieved with the rmNearest rounding mode; it decreases if other modes are used. The test suites run without 'failure warnings about relative errors' on Intel CPUs on Win98, Win2000, WinXP, and Win7. There may be some sporadic warnings with other processors or operating systems; normally, these are not AMath bugs but features of the CPU (and can be avoided by using slightly increased error levels). Note that this does not mean, that the Intel CPUs are 'better', it only reflects the development environment where the error bounds are defined; if developed on other systems, some warnings would occur for Intel CPUs.

The main parts of the packages are the AMath/DAMath units, the Tools units, the Special Functions units, and the complex units:

AMath and DAMath functions

Note: In DAMath the default functions are double precision and explicit extended versions are missing (e.g. succx).

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)
hypot3             Return sqrt(x*x + y*y + z*z)
intpower           Return x^n; via binary exponentiation (no overflow detection)
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
remainder          Return the IEEE754 remainder x REM y = x - rmNearest(x/y)*y
sqrt1pm1           Return sqrt(1+x)-1, accurate even for x near 0, x>=-1
Elementary transcendental functions
arccos             Return the inverse circular cosine of x, |x| <= 1
arccosd            Return the inverse circular cosine of x, |x| <= 1, result in degrees
arccosh            Return the inverse hyperbolic cosine, x >= 1. For x near 1 use arccosh1p(x-1) to reduce cancellation errors!
arccosh1p          Return arccosh(1+x), x>=0, accurate even 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)
arccotcd           Return the continuous inverse circular cotangent; arccotcd(x) = 90 - arctand(x), result in degrees
arccotd            Return the sign symmetric inverse circular cotangent, arccotd(x) = arctand(1/x), x <> 0, result in degrees
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
arcgd              Return the inverse Gudermannian function arcgd(x), |x| < Pi/2
archav             Return the inverse haversine archav(x), 0 <= x <= 1
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
arcsind            Return the inverse circular sine of x, |x| <= 1, result in degrees
arcsinh            Return the inverse hyperbolic sine of x
arctan2            Return arctan(y/x); result in [-Pi..Pi] with correct quadrant
arctand            Return the inverse circular tangent of x, result in degrees
arctanh            Return the inverse hyperbolic tangent of x, |x| < 1
compound           Return (1+x)^n; accurate version of Delphi/VP internal function
cos                Accurate version of circular cosine, uses system.cos for |x| <= Pi/4
cosd               Return cos(x), x in degrees
cosh               Return the hyperbolic cosine of x
coshm1             Return cosh(x)-1, accurate even for x near 0
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
cotd               Return cot(x), x in degrees
coth               Return the hyperbolic cotangent of x, x<>0
covers             Return the coversine covers(x) = 1 - sin(x)
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
exp10m1            Return 10^x - 1; special code for small x
exp2               Return 2^x
exp2m1             Return 2^x-1, accurate even for x near 0
exp3               Return 3^x
exp5               Return 5^x
exp7               Return 7^x
expm1              Return exp(x)-1, accurate even for x near 0
expmx2h            Return exp(-0.5*x^2) with damped error amplification
exprel             Return exprel(x) = (exp(x) - 1)/x,  1 for x=0
expx2              Return exp(x*|x|) with damped error amplification
gd                 Return the Gudermannian function gd(x)
hav                Return the haversine hav(x) = (1 - cos(x))/2
langevin           Return the Langevin function L(x) = coth(x) - 1/x, L(0) = 0
ln1mexp            Return ln(1-exp(x)), x<0
ln1p               Return ln(1+x), accurate even for x near 0
ln1pmx             Return ln(1+x)-x, accurate even for -0.5 <= x <= 0.5
lnxp1              Delphi alias for ln1p
log10              Return base 10 logarithm of x
log10p1            Return log10(1+x), accurate even for x near 0
log2               Return base 2 logarithm of x
log2p1             Return log2(1+x), accurate even for x near 0
logbase            Return base b logarithm of x
logistic           Return logistic(x) = 1/(1+exp(-x))
logit              Return logit(x) = ln(x/(1.0-x)), accurate near x=0.5
logN               Delphi alias for logbase
pow1p              Return (1+x)^y, x > -1
pow1pm1            Return (1+x)^y - 1; special code for small x,y
power              Return x^y; if frac(y)<>0 then x must be > 0
powm1              Return x^y - 1; special code for small x,y
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
sinc               Return the cardinal sine sinc(x) = sin(x)/x
sincos             Return accurate values s=sin(x), c=cos(x)
sincosd            Return sin(x) and cos(x), x in degrees
sincosPi           Return s=sin(Pi*x), c=cos(Pi*x); (s,c)=(0,1) for abs(x) >= 2^64
sincPi             Return the normalised cardinal sine sincPi(x) = sin(Pi*x)/(Pi*x)
sind               Return sin(x), x in degrees
sinh               Return the hyperbolic sine of x, accurate even for x near 0
sinhcosh           Return s=sinh(x) and c=cosh(x)
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
tand               Return tan(x), x in degrees
tanh               Return the hyperbolic tangent of x, accurate even for x near 0
vers               Return the versine vers(x) = 1 - cos(x)
Floating point functions
copysign           Return abs(x)*sign(y)
copysignd          Return abs(x)*sign(y)
copysigns          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, 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, return m=d, e=0
frexps             Return the mantissa m and exponent e of s with s = m*2^e, 0.5 <= abs(m) < 1; if s is 0, +-INF, NaN, return m=s, e=0
ilogb              Return base 2 exponent of x. For finite x ilogb = floor(log2(|x|)), otherwise -MaxLongint for x = 0 or MaxLongint if x = +-INF or Nan.
IsInf              Return true if x is +INF or -INF
IsInfD             Return true if d is +INF or -INF
IsInfS             Return true if s 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
IsNaNorInfS        Return true if s is a NaN or infinite
IsNaNS             Return true if s is a NaN
ldexp              Return x*2^e
ldexpd             Return d*2^e
ldexps             Return s*2^e
predd              Return next representable double after d in the direction -Inf
preds              Return next representable single after s in the direction -Inf
predx              Return next representable extended after x in the direction -Inf
rint               Return the integral value nearest x for the current rounding mode
scalbn             Return x*2^e
succd              Return next representable double after d in the direction +Inf
succs              Return next representable single after s in the direction +Inf
succx              Return next representable extended after x in the direction +Inf
ulpd               Return the 'unit in the last place': ulpd(d)=|d|-predd(|d|) for finite d
ulps               Return the 'unit in the last place': ulps(s)=|s|-preds(|s|) for finite s
ulpx               Return the 'unit in the last place': ulpx(x)=|x|-predx(|x|) for finite x
FPU control functions (Not in DAMath)
Get8087CW          Return the FPU 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
CSEvalX            Evaluate Chebyshev sum a[0]/2 + a[1]*T_1(x) +..+ a[n-1]*T_(n-1)(x) using Clenshaw algorithm
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 a double vector
MeanAndStdDev      Accurate mean and sample standard deviation of a double vector
MeanAndStdDevX     Accurate mean and sample standard deviation of an extended vector
meanx              Compute accurate mean = sum(a[i], i=0..n-1)/n of an extended vector
moment             Return the first 4 moments, skewness, and kurtosis of a double vector
momentx            Return the first 4 moments, skewness, and kurtosis of an extended vector
mssqd              Calculate mean mval and ssqd sum((a[i]-mval)^2) of a double vector
mssqx              Calculate mean mval and ssqx sum((a[i]-mval)^2) of an extended vector
norm2              Calculate the 2-norm = sqrt(sum(a[i]^2, i=0..n-1)) of a 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)
PolEvalS           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)
PolEvalEE          Evaluate polynomial with dynamic absolute error estimate (double version)
PolEvalEEX         Evaluate polynomial with dynamic absolute error estimate (extended version)
PolEvalCHE         Evaluate polynomial with dynamic absolute error estimate using compensated Horner scheme
PolEvalCHEX        Evaluate polynomial with dynamic absolute error estimate using compensated Horner scheme
psdev              Return the population standard deviation of a double vector
psdevx             Return the population standard deviation of an extended vector
pvar               Return the population variance of a double vector
pvarx              Return the population variance of an extended vector
rms                Calculate the RMS value sqrt(sum(a[i]^2, i=0..n-1)/n) of a double vector
rmsx               Calculate the RMS value sqrt(sum(a[i]^2, i=0..n-1)/n) of an extended vector
ssdev              Return the sample standard deviation of a double vector
ssdevx             Return the sample standard deviation of an extended vector
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 a 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 a double vector
sumsqrx            Calculate sum(a[i]^2, i=0..n-1) of an extended vector
svar               Return the sample variance of a double vector
svarx              Return the sample variance of an extended vector
Argument reduction for trigonometric functions
rem_2pi            Return x mod 2*Pi
rem_2pi_sym        Return x mod 2*Pi, -Pi <= result <= Pi
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).
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_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.
Basic double-extended (double-double) functions
Note that the double-double versions have d in their names, e.g. mul21d vs. mul21x
add21x             Return x = a+b
add2x              Return x = a+b
div21x             Return x = a/b,  b<>0
div2x              Return x = a/b,  b<>0
inv2x              Return x = 1/b,  b<>0
mul21x             Return x = a*b
mul2x              Return x = a*b
pow2xi             Return y = a^n,  frac(n)=0, a<>0 if n<0
sqr2x              Return x = a^2
sqrt2x             Return x = sqrt(a),  a >= 0
sub2x              Return x = a-b
xto2x              Return x = a
xxto2x             Return x = a+b
Other functions
angle2             Return the accurate angle between the vectors (x1,x2) and (y1,y2)
area_triangle      Return the area of the triangle defined by the points (xi,yi)
Dbl2Hex            Return d as a big-endian hex string
DegToRad           Convert angle x from degrees to radians
Ext2Dbl            Return x as double, or +-Inf if too large
Ext2Hex            Return x as a big-endian hex string
fisEQd             Return true if x and y are bit-identical
fisEQs             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
fisNEs             Return true if x and y are not bit-identical
fisNEx             Return true if x and y are not bit-identical
Hex2Dbl            Convert big-endian hex string to double,   OK if code=0, inverse of Dbl2Hex
Hex2Ext            Convert big-endian hex string to extended, OK if code=0, inverse of Ext2Hex
Hex2Sgl            Convert big-endian hex string to single,   OK if code=0, inverse of Sgl2Hex
in_triangle        Return true if the point (x,y) lies strictly inside the triangle defined  by the three points (xi,yi), false if it lies on a side or outside.
in_triangle_ex     Return +1 if the point (x,y) lies strictly inside the triangle defined by the three points (xi,yi), -1 if it lies strictly outside, 0 otherwise.
isign              Return the sign of x, 0 if x=0 or NAN
maxx               Return the maximum of two extendeds; x,y <> NAN
maxd               Return the maximum of two doubles; x,y <> NAN
maxs               Return the maximum of two singles; x,y <> NAN
mind               Return the minimum of two doubles; x,y <> NAN
mins               Return the minimum of two singles; x,y <> NAN
minx               Return the minimum of two extendeds; x,y <> NAN
orient2d           Return the mathematical orientation of the three points (xi,yi): >0 if the order is counterclockwise, <0 if clockwise, =0 if they are collinear.
RadToDeg           Convert angle x from radians to degrees
RandG              Random number from Gaussian (normal) distribution with given mean and standard deviation |StdDev|
RandG01            Random number from standard normal distribution (Mean=0, StdDev=1)
Sgl2Hex            Return s as a big-endian hex string

Special functions

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

Bessel functions and related
bessel_i0(x)         Return I0(x), the modified Bessel function of the 1st kind, order zero
bessel_i1(x)         Return I1(x), the modified Bessel function of the 1st kind, order one
bessel_j0(x)         Return J0(x), the Bessel function of the 1st kind, order zero
bessel_j1(x)         Return J1(x), the Bessel function of the 1st kind, order one
bessel_k0(x)         Return K0(x), the modified Bessel function of the 2nd kind, order zero, x>0
bessel_k1(x)         Return K1(x), the modified Bessel function of the 2nd kind, order one, x>0
bessel_y0(x)         Return Y0(x), the Bessel function of the 2nd kind, order zero; x>0
bessel_y1(x)         Return Y1(x), the Bessel function of the 2nd kind, order one; x>0

bessel_i0e(x)        Return I0(x)*exp(-|x|), the exponentially scaled modified Bessel function of the 1st kind, order zero
bessel_i1e(x)        Return I1(x)*exp(-|x|), the exponentially scaled modified Bessel function of the 1st kind, order one
bessel_k0e(x)        Return K0(x)*exp(x),    the exponentially scaled modified Bessel function of the 2nd kind, order zero, x>0
bessel_k1e(x)        Return K1(x)*exp(x),    the exponentially scaled modified Bessel function of the 2nd kind, order one, x>0

bessel_in(n,x)       Return I_n(x), the modified Bessel function of the 1st kind, order n.
bessel_jn(n,x)       Return J_n(x), the Bessel function of the 1st kind, order n.
bessel_kn(n,x)       Return K_n(x), the modified Bessel function of the 2nd kind, order n, x > 0.
bessel_yn(n,x)       Return Y_n(x), the Bessel function of the 2nd kind, order n, x > 0.

bessel_iv(v,x)       Return I_v(x), the modified Bessel function of the 1st kind, order v.
bessel_jv(v,x)       Return J_v(x), the Bessel function of the 1st kind, order v.
bessel_kv(v,x)       Return K_v(x), the modified Bessel function of the 2nd kind, order v, x > 0.
bessel_yv(v,x)       Return Y_v(x), the Bessel function of the 2nd kind, order v; x > 0.

bessel_ive(v,x)      Return I_v(x)*exp(-|x|), the exponentially scaled modified Bessel function of the 1st kind, order v.
bessel_kve(v,x)      Return K_v(x)*exp(x), the exponentially scaled modified Bessel function of the 2nd kind, order v, x>0

sph_bessel_in(n,x)   Return i_n(x), the modified spherical Bessel function of the 1st/2nd kind, order n
sph_bessel_jn(n,x)   Return j_n(x), the spherical Bessel function of the 1st kind, order n.
sph_bessel_kn(n,x)   Return k_n(x), the modified spherical Bessel function of the 3rd kind, order n, x>0
sph_bessel_yn(n,x)   Return y_n(x), the spherical Bessel function of the 2nd kind, order n >=0 , x<>0

sph_bessel_ine(n,x)  Return i_n(x)*exp(-|x|), the exponentially scaled modified spherical Bessel function of the 1st/2nd kind, order n
sph_bessel_kne(n,x)  Return k_n(x)*exp(x), the exponentially scaled modified spherical Bessel function of the 3rd kind, order n, x>0

airy_ai(x)           Return the Airy function Ai(x)
airy_aip(x)          Return the Airy function Ai'(x)
airy_bi(x)           Return the Airy function Bi(x)
airy_bip(x)          Return the Airy function Bi'(x)
airy_gi(x)           Return the Airy/Scorer function Gi(x) = 1/Pi*integral(sin(x*t+t^3/3), t=0..INF)
airy_hi(x)           Return the Airy/Scorer function Hi(x) = 1/Pi*integral(exp(x*t-t^3/3), t=0..INF)

kelvin_ber(x)        Return the Kelvin function ber(x)
kelvin_bei(x)        Return the Kelvin function bei(x)
kelvin_ker(x)        Return the Kelvin function ker(x), x > 0
kelvin_kei(x)        Return the Kelvin function kei(x), x >= 0
kelvin_kerkei(x,..)  Return the Kelvin functions kr=ker(x), ki=kei(x), x > 0
kelvin_berbei(x,..)  Return the Kelvin functions br=ber(x), bi=bei(x)

struve_h0(x)         Return H0(x), the Struve function of order 0
struve_h1(x)         Return H1(x), the Struve function of order 1
struve_l0(x)         Return L0(x), the modified Struve function of order 0
struve_l1(x)         Return L1(x), the modified Struve function of order 1
Elliptic integrals, elliptic / theta functions
comp_ellint_1(k)     Return the complete elliptic integral of the 1st kind, |k| < 1, same as EllipticK
comp_ellint_2(k)     Return the complete elliptic integral of the 2nd kind, |k| <=1, same as EllipticEC
comp_ellint_3(nu,k)  Return the complete elliptic integral of the 3rd kind, |k| < 1, same as EllipticPiC, nu<>1
ellint_1(phi,k)      Return the Legendre elliptic integral F(phi,k) of the 1st kind = integral(1/sqrt(1-k^2*sin(x)^2),x=0..phi), |k*sin(phi)| <= 1
ellint_2(phi,k)      Return the Legendre elliptic integral E(phi,k) of the 2nd kind = integral(sqrt(1-k^2*sin(x)^2),x=0..phi), |k*sin(phi)| <= 1
ellint_3(phi,nu,k)   Return the Legendre elliptic integral PI(phi,nu,k) of the 3rd kind = integral(1/sqrt(1-k^2*sin(x)^2)/(1-nu*sin(x)^2),x=0..phi), |k*sin(phi)| <= 1
heuman_lambda(phi,k) Return Heuman's function Lambda_0(phi,k) = F(phi,k')/K(k') + 2/Pi*K(k)*Z(phi,k'), |k|<=1
jacobi_zeta(phi,k)   Return the Jacobi Zeta function Z(phi,k) = E(phi,k) - E(k)/K(k)*F(phi,k), |k|<=1

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 1st 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 Bulirsch's complete elliptic integral of the 1st kind, k<>0
cel2(kc,a,b)         Return Bulirsch's complete elliptic integral of the 2nd kind, kc<>0
cel(kc,p,a,b)        Return Bulirsch's general complete elliptic integral, kc<>0, Cauchy principle value if p<0
el1(x,kc)            Return Bulirsch's incomplete elliptic integral of the 1st kind
el2(x,kc,a,b)        Return Bulirsch's incomplete elliptic integral of the 2nd kind, kc<>0
el3(x,kc,p)          Return Bulirsch's 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
EllipticKim(k)       Return K(i*k), the complete elliptic integral of the 1st kind with imaginary modulus = integral(1/sqrt(1-x^2)/sqrt(1+k^2*x^2),x=0..1)
EllipticCK(k)        Return the complementary complete elliptic integral of the  1st kind, k<>0
EllipticE(z,k)       Return the incomplete elliptic integrals of the 2nd kind, |z|<=1, |k*z|<=1
EllipticEC(k)        Return the complete elliptic integral of the 2nd kind, |k| < 1
EllipticECim(k)      Return E(i*k), the complete elliptic integral of the 2nd kind with imaginary modulus = integral(sqrt(1+k^2*x^2)/sqrt(1-x^2),x=0..1)
EllipticCE(k)        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, nu<>1
EllipticCPi(nu,k)    Return the complementary complete elliptic integral of the 3rd kind, k<>0, nu<>1

EllipticModulus(q)   Return the elliptic modulus k(q) = theta_2(q)^2/theta_3(q)^2, 0 <= q <= 1
EllipticNome(k)      Return the elliptic nome q(k) = exp(-Pi*EllipticCK(k)/EllipticK(k)), |k| < 1

jacobi_am(x,k)       Return the Jacobi amplitude am(x,k)
jacobi_arccn(x,k)    Return the inverse Jacobi elliptic function arccn(x,k), |x| <= 1, x >= sqrt(1 - 1/k^2) if k >= 1
jacobi_arccd(x,k)    Return the inverse Jacobi elliptic function arccd(x,k); |x| <= 1 if |k| < 1; |x| >= 1 if |k| > 1
jacobi_arccs(x,k)    Return the inverse Jacobi elliptic function arccs(x,k), |x| >= sqrt(k^2-1) for |k|>1}
jacobi_arcdc(x,k)    Return the inverse Jacobi elliptic function arcdc(x,k); |x| >= 1 if |k| < 1; |x| <= 1 if |k| > 1
jacobi_arcdn(x,k)    Return the inverse Jacobi elliptic function arcdn(x,k), 0 <= x <= 1, x^2 + k^2 > 1 if |k| < 1;  |x| <= 1 if |k| > 1
jacobi_arcds(x,k)    Return the inverse Jacobi elliptic function arcds(x,k), x^2 + k^2 >= 1}
jacobi_arcnc(x,k)    Return the inverse Jacobi elliptic function arcnc(x,k), x >= 1, x^2 <= k^2/(k^2-1) for |k|>1
jacobi_arcnd(x,k)    Return the inverse Jacobi elliptic function arcnd(x,k), x >= 1, x^2 <= k^2/(1-k^2) if k < 1
jacobi_arcns(x,k)    Return the inverse Jacobi elliptic function arcns(x,k), |x| >= 1, |x| >= k if k >= 1
jacobi_arcsc(x,k)    Return the inverse Jacobi elliptic function arcsc(x,k), |x| <= 1/sqrt(k^2-1) for |k| > 1
jacobi_arcsd(x,k)    Return the inverse Jacobi elliptic function arcsd(x,k), x^2*(1-k^2) <= 1
jacobi_arcsn(x,k)    Return the inverse Jacobi elliptic function arcsn(x,k), |x| <= 1 and |x*k| <= 1
jacobi_cd(x,k)       Return the Jacobi elliptic function cd(x,k)
jacobi_cn(x,k)       Return the Jacobi elliptic function cn(x,k)
jacobi_cs(x,k)       Return the Jacobi elliptic function cs(x,k)
jacobi_dc(x,k)       Return the Jacobi elliptic function dc(x,k)
jacobi_dn(x,k)       Return the Jacobi elliptic function dn(x,k)
jacobi_ds(x,k)       Return the Jacobi elliptic function ds(x,k)
jacobi_nc(x,k)       Return the Jacobi elliptic function nc(x,k)
jacobi_nd(x,k)       Return the Jacobi elliptic function nd(x,k)
jacobi_ns(x,k)       Return the Jacobi elliptic function ns(x,k)
jacobi_sc(x,k)       Return the Jacobi elliptic function sc(x,k)
jacobi_sd(x,k)       Return the Jacobi elliptic function sd(x,k)
jacobi_sn(x,k)       Return the Jacobi elliptic function sn(x,k)
jacobi_theta(n,x,q)  Return the Jacobi theta function theta_n(x,q), n=1..4, 0 <= q < 1
sncndn(x,mc,...)     Return the Jacobi elliptic functions sn,cn,dn for argument x and complementary parameter mc.
theta1p(q)           Return the derivative d/dx(theta_1(x,q)) at x=0, theta1p(q) = 2*q^(1/4)*sum((-1)^n*(2n+1)*q^(n*(n+1)),n=0..Inf), 0 <= q < 1
theta2(q)            Return Jacobi theta_2(q) = 2*q^(1/4)*sum(q^(n*(n+1)),n=0..Inf) 0 <= q < 1
theta3(q)            Return Jacobi theta_3(q) = 1 + 2*sum(q^(n*n)),n=1..Inf); |q| < 1
theta4(q)            Return Jacobi theta_4(q) = 1 + 2*sum((-1)^n*q^(n*(n+1)),n=1..Inf); |q| < 1
Error function and related
dawson(x)            Return Dawson's integral: dawson(x) = exp(-x^2)*integral(exp(t^2), t=0..x)
dawson2(p,x)         Return the generalized Dawson integral F(p,x) = exp(-x^p)*integral(exp(t^p), t=0..x); x,p >= 0
erf(x)               Return the error function erf(x) = 2/sqrt(Pi)*integral((exp(-t^2), t=0..x)
erfc(x)              Return the complementary error function erfc(x) = 1-erf(x)
erfce(x)             Return the exponentially scaled complementary error function erfce(x) = exp(x^2)*erfc(x)
erfc_inv(x)          Return the inverse function of erfc, erfc(erfc_inv(x)) = x, 0 < x < 2
erfg(p,x)            Return the generalized error function integral(exp(-t^p), t=0..x); x,p >= 0
erfi(x)              Return the imaginary error function erfi(x) = erf(ix)/i
erf_inv(x)           Return the inverse function of erf, erf(erf_inv(x)) = x, -1 < x < 1
erf_p(x)             Return the probability function erf_p = integral(exp(-t^2/2)/sqrt(2*Pi), t=-Inf..x)
erf_q(x)             Return the probability function erf_q = integral(exp(-t^2/2)/sqrt(2*Pi), t=x..Inf)
erf_z(x)             Return the probability function erf_z = exp(-x^2/2)/sqrt(2*Pi)
expint3(x)           Return the integral(exp(-t^3), t=0..x), x >= 0
Fresnel(x,s,c)       Return the Fresnel integrals S(x)=integral(sin(Pi/2*t^2),t=0..x) and C(x)=integral(cos(Pi/2*t^2),t=0..x)
FresnelC(x)          Return the Fresnel integral  C(x)=integral(cos(Pi/2*t^2),t=0..x)
FresnelS(x)          Return the Fresnel integral  S(x)=integral(sin(Pi/2*t^2),t=0..x)
gsi(x)               Return the Goodwin-Staton integral gsi(x) = integral(exp(-t*t)/(t+x), t=0..Inf), x > 0
inerfc(n,x)          Return the repeated integral of erfc, n >= -1; scaled with exp(x^2) for x>0
Exponential integrals and related
chi(x)               Return the hyperbolic cosine integral, chi(x) = EulerGamma + ln(|x|) + integral((cosh(t)-1)/t, t=0..|x|)
ci(x)                Return the cosine integral, ci(x) = EulerGamma + ln(|x|) + integral((cos(t)-1)/t, t=0..|x|)
cin(x)               Return the entire cosine integral, cin(x) = integral((1-cos(t))/t, t=0..x)
cinh(x)              Return the entire hyperbolic cosine integral, cinh(x) = integral((cosh(t)-1)/t, t=0..x)
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)
ein(x)               Return the entire exponential integral ein(x) = integral((1-exp(-t))/t, t=0..x)
en(n,x)              Return the exponential integral E_n(x) = integral(exp(-x*t)/t^n, t=1..Inf), x >= 0
gei(p,x)             Return the generalized exponential integral E_p(x) = integral(exp(-x*t)/t^p, t=1..Inf), x >= 0
li(x)                Return the logarithmic integral li(x)  = PV-integral(1/ln(t), t=0..x), x >= 0, x <> 1
shi(x)               Return the hyperbolic sine integral, shi(x) = integral(sinh(t)/t, t=0..x)
si(x)                Return the sine integral, si(x) = integral(sin(t)/t, t=0..x)
ssi(x)               Return the shifted sine integral, ssi(x) = si(x) - Pi/2
Gamma function and related
BatemanG(x)          Return the Bateman function G(x) = psi((x+1)/2) - psi(x/2); x <> 0,-1,-2,...
beta(x,y)            Return beta(x,y)=gamma(x)*gamma(y)/gamma(x+y)
beta3(a,b,x)         Return the non-normalised incomplete beta function B_x(a,b) = integral(t^(a-1)*(1-t)^(b-1), t=0..x);  0 <= x <= 1
binomial(n,k)        Return the binomial coefficient 'n choose k'
dfac(n)              Return the double factorial n!!, n<=MAXDFAC; INF for even n<0
fac(n)               Return the factorial n!, n < MAXGAM-1; INF if n<0
gamma(x)             Return gamma(x), x <= MAXGAM; invalid if x is a non-positive integer
gamma1pm1(x)         Return gamma(1+x)-1, accurate even for x near 0
gammastar(x)         Return Temme's gammastar(x) = gamma(x)/(sqrt(2*Pi)*x^(x-0.5)*exp(-x)), x>0.
gamma_delta_ratio(.) Return gamma(x)/gamma(x+d), accurate even for |d| << |x|
gamma_ratio(x,y)     Return gamma(x)/gamma(y)
ibeta(a,b,x)         Return the normalised incomplete beta function I_x(a,b), a>0, b>0, 0 <= x <= 1
ibeta_inv(a,b,y)     Return the functional inverse of the normalised incomplete beta function with a > 0, b > 0, and 0 <= y <= 1.
igamma(a,x)          Return the non-normalised upper incomplete gamma function GAMMA(a,x) = integral(exp(-t)*t^(a-1), t=x..Inf), x>=0
igammal(a,x)         Return the non-normalised lower incomplete gamma function gamma(a,x) = integral(exp(-t)*t^(a-1), t=0..x); x>=0, a<>0,-1,-2,..
igammap(a,x)         Return the normalised lower incomplete gamma function P(a,x), a >= 0, x >= 0
igammap_inv(a,p)     Inverse incomplete gamma: return x with P(a,x)=p, a>=0, 0 <= p < 1
igammaq(a,x)         Return the normalised upper incomplete gamma function Q(a,x), a >= 0, x >= 0
igammaq_inv(a,q)     Inverse complemented incomplete gamma: return x with Q(a,x)=q, a >= 0, 0 < q <= 1
igammat(a,x)         Return Tricomi's entire incomplete gamma function igammat(a,x) = igammal(a,x)/gamma(a)/x^a = P(a,x)/x^a
igamma_inv(a,p,q)    Return the inverse normalised incomplete gamma function; a > 0, p >= 0, q > 0 and p+q=1.
incgamma(a,x,p,q)    Return the normalised incomplete gamma functions P and Q, a >= 0, x >= 0
incgamma_inv(...)    General procedure for the inverse normalised incomplete gamma function
lnbeta(x,y)          Return the logarithm of |beta(x,y)|=|gamma(x)*gamma(y)/gamma(x+y)|
lnfac(n)             Return ln(n!), INF if n<0
lngamma(x)           Return ln(|gamma(x)|), |x| <= MAXLGM, invalid if x is a non-positive integer.
lngamma1p(x)         Return ln(|gamma(1+x)|) with increased accuracy for x near 0
lngammas(x,s)        Return ln(|gamma(x)|), |x| <= MAXLGM, s=-1,1 is the sign of gamma
lngamma_inv(y)       Inverse of lngamma: return x with lngamma(x) = y, y >= -0.12142, x > 1.4616
pentagamma(x)        Return the pentagamma function psi'''(x), INF if x is a negative integer
poch1(a,x)           Return (pochhammer(a,x)-1)/x, psi(a) if x=0; accurate even for small |x|
pochhammer(a,x)      Return the Pochhammer symbol (a)_x = gamma(a+x)/gamma(a)
polygamma(n,x)       Return the polygamma function: n'th derivative of psi; n>=0, x>0 for n>NXPGMAX
psi(x)               Return the psi (digamma) function of x, INF if x is a non-positive integer
psi_inv(y)           Inverse of psi, return x with psi(x)=y
rgamma(x)            Return the reciprocal gamma function rgamma = 1/gamma(x)
signgamma(x)         Return sign(gamma(x)), useless for 0 or negative integer
tetragamma(x)        Return the tetragamma function psi''(x), NAN/RTE if x is a negative integer
trigamma(x)          Return the trigamma function of x, INF if x is a negative integer
Zeta functions, polylogarithms, and related
cl2(x)               Return the Clausen function: integral(-ln(|2sin(t/2)|),t=0..x) = Im(Li_2(exp(ix)))
dilog(x)             Return dilog(x) = Re(Li_2(x)), Li_2(x) = -integral(ln(1-t)/t, t=0..x)
DirichletBeta(s)     Return the Dirichlet beta function sum((-1)^n/(2n+1)^s, n=0..INF)
DirichletLambda(s)   Return the Dirichlet lambda function sum(1/(2n+1)^s, n=0..INF), s<>1
eta(s)               Return the Dirichlet eta function
etaint(n)            Return the Dirichlet function eta(n) for integer arguments
etam1(s)             Return Dirichlet eta(s)-1
fermi_dirac(n,x)     Return the integer order Fermi-Dirac integral F_n(x) = 1/n!*integral(t^n/(exp(t-x)+1), t=0..INF)
fermi_dirac_m05(x)   Return the complete Fermi-Dirac integral F(-1/2,x)
fermi_dirac_p05(x)   Return the complete Fermi-Dirac integral F(1/2,x)
fermi_dirac_p15(x)   Return the complete Fermi-Dirac integral F(3/2,x)
fermi_dirac_p25(x)   Return the complete Fermi-Dirac integral F(5/2,x)
harmonic(x)          Return the harmonic number function H(x) = psi(x+1) + EulerGamma
harmonic2(x,r)       Return the generalized harmonic function H(x,r) = zeta(r) - zetah(r,x+1); x >= -1
LegendreChi(s,x)     Return Legendre's Chi-function chi(s,x); s>=0, |x|<=1, x<>1 if s<=1
LerchPhi(z,s,a)      Return the Lerch transcendent Phi(z,s,a) = sum(z^n/(n+a)^s, n=0..INF), |z|<=1, s>=0, a>0; s>1 if z=1.
lobachevsky_c(x)     Return the Lobachevski function L(x) = integral(-ln(|cos(t)|), t=0..x)
lobachevsky_s(x)     Return the Lobachevski function Lambda(x) = integral(-ln(|2sin(t)|), t=0..x)
polylog(n,x)         Return the polylogarithm Li_n(x) of integer order; x<1 for n >= 0.
polylogr(s,x)        Return the polylogarithm Li_s(x) of real order; s>=0, |x|<=1, x<>1 if s<=1
primezeta(x)         Return the prime zeta function P(x) = sum(1/p^x, p prime), x > 1/5; for x<1 the real part of P(x) is returned.
ti2(x)               Return the inverse tangent integral, Ti_2(x) = integral(arctan(t)/t, t=0..x)
trilog(x)            Return the trilogarithm function trilog(x) = Re(Li_3(x))
zeta(s)              Return the Riemann zeta function at s, s<>1
zeta1p(x)            Return the Riemann zeta function at 1+x, x<>0
zetah(s,a)           Return the Hurwitz zeta function zetah(s,a) = sum(1/(i+a)^s, i=0..INF), s<>1, a>0
zetaint(n)           Return zeta(n) for integer arguments, n<>1
zetam1(s)            Return Riemann zeta(s)-1, s<>1
Orthogonal polynomials, Legendre functions, and related
chebyshev_t(n,x)     Return Tn(x), the Chebyshev polynomial of the first kind, degree n
chebyshev_u(n,x)     Return Un(x), the Chebyshev polynomial of the second kind, degree n
chebyshev_v(n,x)     Return Vn(x), the Chebyshev polynomial of the third kind, degree n >= 0
chebyshev_w(n,x)     Return Wn(x), the Chebyshev polynomial of the fourth kind, degree n >= 0
gegenbauer_c(n,a,x)  Return Cn(a,x), the nth Gegenbauer (ultraspherical) polynomial with parameter a > -0.5
hermite_h(n,x)       Return Hn(n,x), the nth Hermite polynomial, degree n >= 0
jacobi_p(n,a,b,x)    Return Pn(a,b,x), the nth Jacobi polynomial with parameters a,b
laguerre(n,a,x)      Return Ln(a,x), the nth generalized Laguerre polynomial with parameter a; degree n must be >= 0. x >=0 and a > -1 are the standard ranges.
laguerre_ass(n,m,x)  Return the associated Laguerre polynomial Ln(m,x); n,m >= 0
laguerre_l(n,x)      Return the nth Laguerre polynomial Ln(0,x); n >= 0
legendre_p(l,x)      Return P_l(x), the Legendre polynomial/function P_l, degree l
legendre_plm(l,m,x)  Return the associated Legendre polynomial P_lm(x)
legendre_q(l,x)      Return Q_l(x), the Legendre function of the second kind, degree l >= 0, |x| <> 1
legendre_qlm(l,m,x)  Return Q(l,m,x), the associated Legendre function of the second kind; l >= 0, l+m >= 0, |x|<>1
toroidal_plm(l,m,x)  Return the toroidal harmonic function P(l-0.5,m,x); l,m=0,1; x >= 1
toroidal_qlm(l,m,x)  Return the toroidal harmonic function Q(l-0.5,m,x); l=0,1; x > 1
spherical_harmonic   Return Re and Im of the spherical harmonic function Y_lm(theta,phi)
zernike_r(n,m,r)     Return the Zernike radial polynomial Rnm(r), r >= 0, n >= m >= 0, n-m even
Hypergeometric functions and related
CylinderD(v,x)       Return Whittaker's parabolic cylinder function D_v(x)
CylinderU(a,x)       Return the parabolic cylinder function U(a,x)
CylinderV(a,x)       Return the parabolic cylinder function V(a,x) with 2a integer >= 0
HermiteH(v,x)        Return the Hermite function H_v(x) of degree v
hyperg_0F1(b,x)      Return the confluent hypergeometric limit function 0F1(;b;x)
hyperg_0F1r(b,x)     Return the regularized confluent hypergeometric limit function 0F1(;b;x)/Gamma(b)
hyperg_1F1(a,b,x)    Return the confluent hypergeometric function 1F1(a,b,x); Kummer's function M(a,b,x)
hyperg_1F1r(a,b,x)   Return the regularized Kummer hypergeometric function 1F1(a,b,x)/Gamma(b)
hyperg_2F1(a,b,c,x)  Return the Gauss hypergeometric function 2F1(a,b;c;x)
hyperg_2F1r(a,b,c,x) Return the regularized Gauss hypergeometric function 2F1(a,b,c,x)/Gamma(c)
hyperg_u(a,b,x)      Return Tricomi's confluent hypergeometric function U(a,b,x). If x<0, then a must be an integer and a<0 or 1+a-b an integer < 0.
WhittakerM(k,m,x)    Return the Whittaker M function = exp(-x/2)*x^(0.5+m) * 1F1(m-k-0.5,2m+1,x)
WhittakerW(k,m,x)    Return the Whittaker W function = exp(-x/2)*x^(0.5+m) * U(m-k-0.5,2m+1,x)
Statistical distributions
beta_cdf(a,b,x)      Return the cumulative beta distribution function; a>0, b>0
beta_inv(a,b,y)      Return the functional inverse of the beta distribution function; a > 0, b > 0, 0 <= y <= 1
beta_pdf(a,b,x)      Return the probability density function of the beta distribution with parameters a and b
binomial_cdf(p,n,k)  Return the cumulative binomial distribution function with number of trials n >= 0 and success probability 0 <= p <= 1
binomial_pmf(p,n,k)  Return the binomial distribution probability mass function with number of trials n >= 0 and success probability 0 <= p <= 1
cauchy_cdf(a,b,x)    Return the cumulative Cauchy distribution function with location a and scale b > 0
cauchy_inv(a,b,y)    Return the functional inverse of Cauchy distribution function with location a and scale b > 0
cauchy_pdf(a,b,x)    Return the Cauchy probability density function with location a and scale b > 0
chi2_cdf(nu,x)       Return the cumulative chi-square distribution with nu>0 degrees of freedom, x >= 0
chi2_inv(nu,p)       Return the functional inverse of the chi-square distribution, nu>0, 0 <= p < 1
chi2_pdf(nu,x)       Return the probability density function of the chi-square distribution, nu>0
evt1_pdf(a,b,x)      Return the probability density function of the Extreme Value Type I distribution with location a and scale b > 0, result = exp(-(x-a)/b)/b * exp(-exp(-(x-a)/b))
evt1_cdf(a,b,x)      Return the cumulative Extreme Value Type I distribution function with location a and scale b > 0; result = exp(-exp(-(x-a)/b)).
evt1_inv(a,b,y)      Return the functional inverse of the Extreme Value Type I distribution function with location a and scale b > 0; result = a - b*ln(ln(-y)).
exp_cdf(a,alpha,x)   Return the cumulative exponential distribution function with location a and rate alpha > 0
exp_inv(a,alpha,y)   Return the functional inverse of exponential distribution function with location a and rate alpha > 0
exp_pdf(a,alpha,x)   Return the exponential probability density function with location a and rate alpha > 0
f_cdf(nu1,nu2,x)     Return the cumulative F distribution function; x >= 0, nu1, nu2 > 0
f_inv(nu1,nu2,y)     Return the functional inverse of the F distribution, nu1, nu2 > 0, 0 <= y <= 1
f_pdf(nu1,nu2,x)     Return the probability density function of the F distribution; x >= 0, nu1, nu2 > 0
gamma_cdf(a,b,x)     Return the cumulative gamma distribution function, shape a>0, scale b>0
gamma_inv(a,b,p)     Return the functional inverse of the gamma distribution function, shape a>0, scale b>0
gamma_pdf(a,b,x)     Return the probability density function of a gamma distribution with shape a>0, scale b>0
hypergeo_pmf(n1,n2,n,k)  Return the hypergeometric distribution probability mass function; n,n1,n2 >= 0, n <= n1+n2
hypergeo_cdf(n1,n2,n,k)  Return the cumulative hypergeometric distribution function; n,n1,n2 >= 0, n <= n1+n2
kumaraswamy_pdf(a,b,x)   Return the Kumaraswamy probability density function with shape parameters a,b>0, 0<=x<=1; result = a*b*x^(a-1)*(1-x^a)^(b-1)
kumaraswamy_cdf(a,b,x)   Return the cumulative Kumaraswamy distribution function with shape parameters a,b > 0, 0 <= x <= 1; result = 1-(1-x^a)^b
kumaraswamy_inv(a,b,y)   Return the functional inverse of the Kumaraswamy distribution with shape parameters a,b > 0; result = [1-(1-y)^(1/b)]^(1/a)
levy_pdf(a,b,x)      Return the Levy probability density function with location a and scale parameter b > 0
levy_cdf(a,b,x)      Return the cumulative Levy distribution function with location a and scale parameter b > 0
levy_inv(a,b,y)      Return the functional inverse of the Levy distribution with location a and scale parameter b > 0
laplace_cdf(a,b,x)   Return the cumulative Laplace distribution function with location a and scale b > 0
laplace_inv(a,b,y)   Return the functional inverse of the Laplace distribution with location a and scale b > 0
laplace_pdf(a,b,x)   Return the Laplace probability density function with location a and scale b > 0, result = exp(-abs(x-a)/b) / (2*b)
logistic_cdf(a,b,x)  Return the cumulative logistic distribution function with location a and scale parameter b > 0
logistic_inv(a,b,y)  Return the functional inverse of the logistic distribution with location a and scale parameter b > 0
logistic_pdf(a,b,x)  Return the logistic probability density function with location a and scale parameter b > 0, exp(-(x-a)/b)/b/(1+exp(-(x-a)/b))^2
lognormal_cdf(a,b,x) Return the cumulative log-normal distribution function with location a and scale parameter b > 0, zero for x <= 0.
lognormal_inv(a,b,y) Return the functional inverse of the log-normal distribution  with location a and scale parameter b > 0, 0 < y < 1.
lognormal_pdf(a,b,x) Return the log-normal probability density function with location a and scale parameter b > 0, zero for x <= 0.
maxwell_pdf(b,x)     Return the Maxwell probability density function with scale b > 0, x >= 0
maxwell_cdf(b,x)     Return the cumulative Maxwell distribution function with scale b > 0, x >= 0
maxwell_inv(b,y)     Return the functional inverse of the Maxwell distribution with scale b > 0
moyal_pdf(a,b x)     Return the Moyal probability density function with location a and scale parameter b > 0
moyal_cdf(a,b,x)     Return the cumulative Moyal distribution function with location a and scale parameter b > 0
moyal_inv(a,b,y)     Return the functional inverse of the Moyal distribution with location a and scale parameter b > 0
negbinom_cdf(p,r,k)  Return the cumulative negative binomial distribution function with r > 0 and success probabiility 0 <= p <= 1
negbinom_pmf(p,r,k)  Return the negative binomial distribution probability mass function with r > 0 and success probabiility 0 <= p <= 1
normal_cdf(mu,sd,x)  Return the normal (Gaussian) distribution density function with mean mu and standard deviation sd > 0
normal_inv(mu,sd,y)  Return the functional inverse of the normal (Gaussian) distribution with mean mu and standard deviation sd > 0, 0 < y < 1
normal_pdf(mu,sd,x)  Return the normal (Gaussian) probability density function with mean mu and standard deviation sd > 0
normstd_cdf(x)       Return the standard normal distribution function
normstd_inv(y)       Return the inverse standard normal distribution function, 0 < y < 1
normstd_pdf(x)       Return the std. normal probability density function exp(-x^2/2)/sqrt(2*Pi)
pareto_cdf(k,a,x)    Return the cumulative Pareto distribution function minimum value k > 0 and shape a, x >= a > 0, result = 1-(k/x)^a
pareto_inv(k,a,y)    Return the functional inverse of the Pareto distribution with minimum value k > 0 and shape a, x >= a > 0
pareto_pdf(k,a,x)    Return the Pareto probability density function with minimum value k > 0 and shape a, x >= a > 0, result = (a/x)*(k/x)^a
poisson_cdf(mu,k)    Return the cumulative Poisson distribution function with mean mu >= 0
poisson_pmf(mu,k)    Return the Poisson distribution probability mass function with mean mu >= 0
rayleigh_pdf(b,x)    Return the Rayleigh probability density function with scale b > 0, x >= 0; result = x*exp(-0.5*(x/b)^2)/b^2
rayleigh_cdf(b,x)    Return the cumulative Rayleigh distribution function with scale b > 0, x >= 0
rayleigh_inv(b,y)    Return the functional inverse of the Rayleigh distribution with scale b > 0
triangular_cdf(a,b,c,x)  Return the cumulative triangular distribution function with lower limit a, upper limit b, mode c;  a < b, a <= c <= b
triangular_inv(a,b,c,y)  Return the functional inverse of the triangular distribution with lower limit a, upper limit b, mode c; a < b, a <= c <= b, 0 <= y <= 1
triangular_pdf(a,b,c,x)  Return the triangular probability density function with lower limit a, upper limit b, mode c;  a < b, a <= c <= b
t_cdf(nu,t)          Return the cumulative Student t distribution with nu>0 degrees of freedom
t_inv(nu,p)          Return the functional inverse of Student's t distribution, nu>0, 0 <= p <= 1
t_pdf(nu,x)          Return the probability density function of Student's t distribution, nu>0
uniform_cdf(a,b,x)   Return the cumulative uniform distribution function on [a,b], a < b
uniform_inv(a,b,y)   Return the functional inverse of the uniform distribution on [a,b], a < b
uniform_pdf(a,b,x)   Return the uniform probability density function on [a,b], a < b
weibull_cdf(a,b,x)   Return the cumulative Weibull distribution function with shape parameter a > 0 and scale parameter b > 0
weibull_inv(a,b,y)   Return the functional inverse of the Weibull distribution shape parameter a > 0 and scale parameter b > 0
weibull_pdf(a,b,x)   Return the Weibull probability density function with shape a > 0 and scale b > 0, result = a*x^(a-1)*exp(-(x/b)^a)/ b^a; x > 0
zipf_pmf(r,k)        Return the Zipf distribution probability mass function k^(-(r+1))/zeta(r+1), r>0, k>0
zipf_cdf(r,k)        Return the cumulative Zipf distribution function H(k,r+1)/zeta(r+1), r>0, k>0
Other special functions
agm(x,y)             Return the arithmetic-geometric mean of |x| and |y|
bernoulli(n)         Return the nth Bernoulli number, 0 if n<0 or odd n >= 3
bernpoly(n,x)        Return the Bernoulli polynomial B_n(x), 0 <= n <= MaxBernoulli
catalan(x)           Return the Catalan function C(x) = binomial(2x,x)/(x+1)
cosint(n,x)          Return cosint(n,x) = integral(cos(t)^n, t=0..x), n >= 0
debye(n,x)           Return the Debye function D(n,x) = n/x^n*integral(t^n/(exp(t)-1),t=0..x) of order n>0, x>=0
fibpoly(n,x)         Return the Fibonacci polynomial F_n(x)
LambertW(x)          Return the Lambert W function (principal branch), x >= -1/e
LambertW1(x)         Return the Lambert W function (-1 branch), -1/e <= x < 0
li_inv(x)            Return the functional inverse of li(x), li(li_inv(x))=x
lucpoly(n,x)         Return the Lucas polynomial L_n(x)
RiemannR(x)          Return the Riemann prime counting function R(x), x >= 1/16
sinint(n,x)          Return sinint(n,x) = integral(sin(t)^n, t=0..x), n >= 0

AMTools and DAMTools functions

Function minimization
localmin(f,a,b,eps,t,x,fx,ic): Brent's algorithm (with guaranteed convergence) for finding a local minimum of the function f in the interval (a,b). x is the approximate minimum abscissa, fx=f(x). eps and t define a tolerance tol = eps*|x|+t. f is never evaluated for two points closer together than tol. eps shall not be < 2*eps_x, preferably not smaller than sqrt(eps_x). ic is the iteration count, -1 if a=b, 0 if max. count = 5000 exceeded. The algorithm combines golden section search and successive parabolic interpolation using only function (not derivative) evaluations.

mbrent(f,a,b,t,x,fx,ic): Find a local minimum of the function f in the interval (a,b). x is the approximate minimum abscissa, fx = f(x). Simplified version of procedure localmin with fixed eps = 0.5*sqrt(eps_x). ic is the iteration count, -1 if a=b, 0 if max. = 5000 exceeded.

fmin(f,a,b,tol): Find a local minimum of the function f in the interval (a,b) and return the approximate minimum abscissa. fmin is a simple shell version for the procedure localmin using eps = 0.5*sqrt(eps_x) and t = tol/3.

Zero finding
zbrent(f,a,b,t,ic,err): Brent/Dekker algorithm with guaranteed convergence for finding a zero of a function: Return a zero x of the function f in the interval [a,b] to within a tolerance 6*eps_x*|x|+2*t, where t is a positive tolerance; assumes that f(a) and f(b) have different signs. ic is the iteration count; err is an error code (0: no error, -1: if f(a) and f(b) have the same sign, -2: max. iteration count exceeded). The algorithm is based on a combination of successive interpolations and bisection.

zeroin(f,a,b,t): Return a zero x of the function f in the interval [a, b] to within a tolerance 6*eps_x*|x| + 2*t, where t is a positive tolerance, assumes that f(a) and f(b) have different signs. Simplified version of zbrent.

Numerical integration
quanc8(fun,a,b,abserr,relerr, result, ..): Estimate the integral of fun(x) from a to b to a user provided tolerance using an automatic adaptive routine based on the 8-panel Newton-Cotes rule; used for smooth functions on finite intervals.

quagk(f,a,b,epsabs,result,abserr,ier): General global adaptive quadrature of f over (a,b) based on Gauss-Kronrod rules for the subintervals, with acceleration by Wynn's epsilon algorithm. Simple shell for qags and qagi, a or b may be infinite.

qags(f,a,b, .. , result, .. ,ier): Global adaptive quadrature of f over (a,b) based on 21-point Gauss-Kronrod rule for the subintervals, with acceleration by Wynn's epsilon algorithm. Simplified user interface to procedure qagse.

qagi(f,bound,inf, .. ,result, .. ,ier): Global adaptive quadrature of f over an infinite interval based on transformed 15-point Gauss-Kronrod rule for the subintervals, with acceleration by Wynn's epsilon algorithm. Simplified user interface to procedure qagie.

qawc(f,a,b,c, .. , result, .. ,ier): Adaptive quadrature of the function f(x)/(x-c) over the finite interval (a,b) with the singularity at c. The routine calculates an approximation to the Cauchy principal value. Simplified user interface to procedure qawce.

qk21(f,a,b,result, ..): Integrate function f over (a,b) with 21-point Gauss-Kronrod rule.

qagse(f,a,b, .. ,result, .. ,ier, ..): Global adaptive quadrature of f over (a,b) based on 21-point Gauss-Kronrod rule for the subintervals, with acceleration by Wynn's epsilon algorithm.

qagie(f,bound,inf, .. ,result, .. ,ier, ..): Global adaptive quadrature of f over an infinite interval based on transformed 15-point Gauss-Kronrod rule for the subintervals, with acceleration by Wynn's epsilon algorithm.

qawce(f,a,b,c, .. , result, .. ,ier, ..): Adaptive quadrature of the function f(x)/(x-c) over the finite interval (a,b) with the singularity at c with c<>a and c<>b. The routine calculates an approximation to the Cauchy principal value.

intde(f,a,b,eps,result, .. ,ier): Automatic quadrature of f(x) over the finite interval (a,b) using Double Exponential transformation.

intdei(f,a,eps,result, .. ,ier): Automatic quadrature of f(x) over (a,INF) using Double Exponential transformation when f(x) has no oscillatory factor.

intdeo(f,a,omega,eps,result, .. ,ier): Automatic quadrature of f(x) over (a,INF) using Double Exponential transformation when f(x) has an oscillatory factor, i.e. f has the form f(x) = g(x)*sin(omega*x + theta) as x -> INF.

Convergence acceleration
levinu1(an,n,nmax,qnum,qden,extsum,sn,ierr): Perform one step of the Levin u-transformation for sum(an, n=0..). Note: The driver routine has to analyze the convergence of the process.

wynneps1(an,n,nmax,sum,e,extlim,sn,ierr): Perform one step of Wynn's epsilon algorithm for the sequence an, n=0.. or the sum(an, n=0..) Note: The driver routine has to analyze the convergence of the whole process.

Solving quadratic and cubic equations
squad(a,b,c,x1,y1,x2,y2): 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,x1,y1,x2,y2): 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.

cubsolve(a,b,c,d,x,x1,y1,x2,y2): Solve the cubic equation ax^3 + bx^2 + cx + d = 0: compute a real root x (may be INF if a~0) and two complex zeros x1 + i*y1, x2 + i*y2 where y2 = -y1 may be zero, i.e. there are three reel roots.

AMath and DAMath complex functions

Complex arithmetic and basic functions
cabs(z)                Return the complex absolute value |z| = sqrt(z.re^2 + z.im^2)
cadd(x,y; z)           Return the complex sum z = x + y
carg(z)                Return the principle value of the argument or phase angle arg(z) = arctan2(z.im, z.re)
ccis(x; z)             Return z = exp(i*x) = cos(x) + i*sin(x)
cconj(z; w)            Return the complex conjugate w = z.re - i*z.im
cdiv(x,y; z)           Return the quotient z = x/y
cinv(z; w)             Return the complex inverse w = 1/z
cmul(x,y; z)           Return the complex product z = x*y
cneg(z; w)             Return the negative w = -z
cpolar(z; r,theta)     Return the polar form z = r*exp(i*theta) with r = |z|, theta = arg z
cpoly(z,a,n; w)        Evaluate polynomial; return a[0] + a[1]*z + ... + a[n-1]*z^(n-1)
cpolyr(z,a,n; w)       Evaluate polynomial; return a[0] + a[1]*z + ... + a[n-1]*z^(n-1) with real a[]
cset(z; x,y)           Set real and imaginary part of z = x+iy
csgn(z)                Return the sign of z. Result = isign(z.re) if z.re<>0, isign(z.im) otherwise
csqr(z; w)             Return the square w = z^2
csqrt(z; w)            Return the complex principal square root w = sqrt(z)
csqrt1mz2(z; w)        Return the complex principal square root w = sqrt(1-z^2)
csub(x,y; z)           Return the complex difference z = x - y
rdivc(x,y; z)          Return the quotient z = x/y for real x
Complex transcendental functions
cagm(x,y; w)           Return the 'optimal' arithmetic-geometric mean w = AGM(x,y)
cagm1(z; w)            Return the 'optimal' arithmetic-geometric mean w = AGM(1,z)
carccos(z; w)          Return the principal value of the complex inverse circular cosine w = arccos(z)
carccosh(z; w)         Return the principal value of the complex inverse hyperbolic cosine w = arccosh(z)
carccot(z; w)          Return the principal value of the complex inverse circular cotangent w = arccot(z) = arctan(1/z)
carccotc(z; w)         Return the principal value of the complex inverse circular cotangent w = arccotc(z) = Pi/2 - arctan(z)
carccoth(z; w)         Return the principal value of the complex inverse hyperbolic cotangent w = arccoth(z) = arctanh(1/z)
carccothc(z; w)        Return the principal value of the complex inverse hyperbolic cotangent w = arccothc(z) = arctanh(z) + i*Pi/2
carccsc(z; w)          Return the principal value of the complex inverse circular cosecant w = arccsc(z) = arcsin(1/z)
carccsch(z; w)         Return the principal value of the complex inverse hyperbolic cosecant w = arccsch(z) = arcsinh(1/z)
carcsec(z; w)          Return the principal value of the complex inverse circular secant w = arcsec(z) = arccos(1/z)
carcsech(z; w)         Return the principal value of the complex inverse hyperbolic secant w = arcsech(z) = arccosh(1/z)
carcsin(z; w)          Return the principal value of the complex inverse circular sine w = arcsin(z)
carcsinh(z; w)         Return the principal value of the complex inverse hyperbolic sine w = arcsinh(z)
carctan(z; w)          Return the principal value of the complex inverse circular tangent w = arctan(z)
carctanh(z; w)         Return the principal value of the complex inverse hyperbolic tangent w = arctanh(z)
ccbrt(z; w)            Return the complex principal cube root w = cbrt(z) = z^(1/3)
ccos(z; w)             Return the complex circular cosine w = cos(z)
ccosh(z; w)            Return the complex hyperbolic cosine w = cosh(z)
ccot(z; w)             Return the complex circular cotangent w = cot(z)
ccoth(z; w)            Return the complex hyperbolic cotangent w = coth(z)
ccsc(z; w)             Return the complex circular cosecant w = csc(z) = 1/sin(z)
ccsch(z; w)            Return the complex hyperbolic cosecant w = csch(z) = 1/sinh(z)
cdilog(z; w)           Return the principal branch of the complex dilogarithm w = -integral(ln(1-t)/t, t=0..z)
cellck(k; w)           Return w = K'(k), the complementary complete elliptic integral of the first kind
cellk(k; w)            Return w = K(k), the complete elliptic integral of the first kind
cexp(z; w)             Return the complex exponential function w = exp(z)
cexpm1(z; w)           Return w = exp(z)-1, accuracy improved for z near 0
cgamma(z; w)           Return the complex Gamma function w = Gamma(z)
cln(z; w)              Return the complex natural logarithm w = ln(z); principal branch ln(|z|) + i*arg(z), accurate near |z|=1
cln1p(z; w)            Return the principal branch of ln(1+z), accuracy improved for z near 0
clngamma(z; w)         Return w = lnGamma(z), the principal branch of the log-Gamma function
clog10(z; w)           Return the principal branch of the base 10 logarithm of z, w = ln(z)/ln(10)
clogbase(b,z; w)       Return the principal branch of the base  b logarithm of z, w = ln(z)/ln(b)
cnroot(z,n; w)         Return the complex principal n'th root w = z^(1/n)
cnroot1(n; z)          Return the principal nth root of unity z = exp(2*Pi*i/n)
cpow(z,a; w)           Return the principal value of the complex power w = z^a = exp(a*ln(z))
cpowx(z,x; w)          Return the principal value w = z^x = |z|^x * exp(i*x*arg(z))
cpsi(z; w)             Return the complex digamma function w = psi(z), z <> 0,-1,-2,...
csec(z; w)             Return the complex circular secant w = sec(z) = 1/cos(z)
csech(z; w)            Return the complex hyperbolic secant w = sech(z) = 1/cosh(z)
csin(z; w)             Return the complex circular sine w = sin(z)
csinh(z; w)            Return the complex hyperbolic sine w = sinh(z)
csurd(z,n; w)          Return the complex n'th root w = z^(1/n) with arg(w) closest to arg(z)
ctan(z; w)             Return the complex circular tangent w = tan(z)
ctanh(z; w)            Return the complex hyperbolic tangent w = tanh(z)

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 http://www.intel.com/products/processor/manuals/
  3. D. Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, 1991; 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. Available from http://www.digizeitschriften.de/en/dms/toc/?PPN=PPN362160546_0007
  11. R. Bulirsch, Numerical Calculation of Elliptic Integrals and Elliptic Functions, part III. Numerische Mathematik 13, 305-315, 1969. Available from http://www.digizeitschriften.de/en/dms/toc/?PPN=PPN362160546_0013
  12. B.C. Carlson, Computing Elliptic Integrals by Duplication, Numerische Mathematik 33, 1-16, 1979. Available from http://www.digizeitschriften.de/en/dms/toc/?PPN=PPN362160546_0033
  13. W.H. Press et al, Numerical Recipes in C, 2nd ed., Cambridge, 1992. http://www.nrbook.com/a/bookcpdf.html
  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. Available from http://oai.dtic.mil/oai/oai?verb=getRecord&metadataPrefix=html&identifier=AD0639052
  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/
  20. Special functions by Wayne Fullerton, http://www.netlib.org/fn/. Almost identical to the FNLIB subset of SLATEC [14].
  21. GNU Scientific Library, GSL-1.14 (March 2010), http://www.gnu.org/software/gsl/
  22. A.J. MacLeod, MISCFUN: A software package to compute uncommon special functions. ACM Trans. on Math. Soft. 22 (1996), pp. 288-301. Fortran source: http://netlib.org/toms/757
  23. R.M. Corless, G.H. Gonnet, D.E.G. Hare, D.J. Jeffrey, D.E. Knuth, On the Lambert W Function, Adv. Comput. Math., 5 (1996), pp. 329-359. http://www.apmaths.uwo.ca/~rcorless/frames/PAPERS/LambertW/LambertW.ps or http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.33.3583
  24. D. Veberic, Having Fun with Lambert W(x) Function, 2010. http://arxiv.org/abs/1003.1628
  25. I. Smith, Examples.xls/txt Version 3.3.4, Personal communication, 2010
  26. N.M. Temme, A Set of Algorithms for the Incomplete Gamma Functions, Probability in the Engineering and Informational Sciences, 8 (1994), pp. 291-307. Available from http://oai.cwi.nl/oai/asset/10080/10080A.pdf
  27. A.R. Didonato, A.H. Morris, Computation of the Incomplete Gamma Function Ratios and their Inverse. ACM TOMS, Vol 12, No 4, Dec 1986, pp. 377-393. Fortran source: ACM TOMS 13 (1987) pp. 318-319; available from http://netlib.org/toms/654
  28. R.P. Brent, Algorithms for Minimization without Derivatives, Englewood Cliffs, 1973. Scanned copy available from the author's site: http://maths.anu.edu.au/~brent/pub/pub011.html
  29. G.E. Forsythe, M.A. Malcolm, C.B. Moler, Computer Methods for Mathematical Computations, Englewood Cliffs, 1977. Fortran code from http://www.netlib.org/fmm/
  30. [NIST]: F.W.J. Olver, D.W. Lozier, R.F. Boisvert, C.W. Clark, NIST Handbook of Mathematical Functions, Cambridge, 2010. Online resource: NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/
  31. R.E. Crandall, Note on fast polylogarithm computation, 2006. For a local copy see the links page.
  32. D.E. Knuth: The Art of computer programming; Volume 1, Fundamental Algorithms, 3rd ed., 1997; Volume 2, Seminumerical Algorithms, 3rd ed., 1998; http://www-cs-faculty.stanford.edu/~knuth/taocp.html
  33. http://functions.wolfram.com/: Formulas and graphics about mathematical functions for the mathematical and scientific community. Also used: http://mathworld.wolfram.com/ ("the web's most extensive mathematical resource") and Wolfram Alpha - Computational Knowledge Engine at http://www.wolframalpha.com/ for online calculation of multi-precision special function reference values.
  34. R. Piessens, E. de Doncker-Kapenga, C.W. Überhuber, D. Kahaner, QUADPACK: A subroutine package for automatic integration (1983). Public domain Fortran source: http://www.netlib.org/quadpack/
  35. L.C. Maximon, The dilogarithm function for complex argument, 2003, Proc. R. Soc. Lond. A, 459, 2807-2819, doi: 10.1098/rspa.2003.1156; http://rspa.royalsocietypublishing.org/content/459/2039/2807.full.pdf
  36. P. Borwein, An Efficient Algorithm for the Riemann Zeta Function, CMS Conference Proc. 27 (2000), pp. 29-34. Available as http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P155.pdf
  37. A.R. DiDonato, A.H. Morris, Algorithm 708: Significant digit computation of the incomplete beta function ratios, ACM TOMS, Vol.18, No.3, 1992, pp.360-373. Fortran source available from http://netlib.org/toms/708
  38. T. Ooura's Fortran and C source code for automatic quadrature using Double Exponential transformation; available from http://www.kurims.kyoto-u.ac.jp/~ooura/intde.html
  39. R: A Language and Environment for Statistical Computing, Version 2.11.1, http://www.r-project.org/
  40. S.V. Aksenov et al., Application of the combined nonlinear-condensation transformation to problems in statistical analysis and theoretical physics. Computer Physics Communications, 150, 1-20, 2003. E-print available from http://arxiv.org/pdf/math/0207086v1
  41. C. Ferreira, J.L. López, Asymptotic expansions of the Hurwitz-Lerch zeta function, J. Math. Anal. Appl. 298 (2004), 210-224. Available from http://dx.doi.org/10.1016/j.jmaa.2004.05.040
  42. Y. L. Luke, Algorithms for the Computation of Mathematical Functions, Academic Press, 1977
  43. J. Pearson, Computation of Hypergeometric Functions, Master's thesis, University of Oxford, 2009. Available as http://people.maths.ox.ac.uk/porterm/research/pearson_final.pdf
  44. K.E. Muller, Computing the confluent hypergeometric function M(a,b,x), Numerische Mathematik 90, 179-196, 2001
  45. N.N. Lebedev, Special Functions and Their Applications, Dover, New York, 1972
  46. S. Graillat, Ph. Langlois, N. Louvet, Compensated Horner Scheme, Research Report No RR2005-04, 2005, Université de Perpignan. Available from http://www-pequan.lip6.fr/~graillat/papers/rr2005-04.pdf or http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.81.2979
  47. M. Goano, Algorithm 745: Computation of the complete and incomplete Fermi-Dirac integral. ACM TOMS, Vol.21, No.3, 1995, pp.221-232. Fortran source available from http://netlib.org/toms/745
  48. D. Dijkstra, A Continued Fraction Expansion for a Generalization of Dawson's Integral, Mathematics of Computation, Vol.31, 503-510 (1977). Available as http://doc.utwente.nl/75001/1/Dijkstra77continued.pdf
  49. W. Gautschi, Evaluation of the Repeated Integrals of the Coerror Function, ACM TOMS, Vol.3, No.3, 1977, pp.240-252. Fortran source available from http://netlib.org/toms/521
  50. A. Erdélyi et al., Higher Transcendental Functions Vol. I-III, California Institute of Technology - Bateman Manuscript Project, 1953-1955, Available via http://en.wikipedia.org/wiki/Bateman_Manuscript_Project
  51. N.M. Temme, On the Numerical Evaluation of the Ordinary Bessel Function of the Second Kind. J. Comput. Phys., 21(3): 343-350 (1976), section 2. Available as http://oai.cwi.nl/oai/asset/10710/10710A.pdf
  52. N.M. Temme, On the Numerical Evaluation of the Modified Bessel Function of the Third Kind, 2nd edition. Preprint, available from http://oai.cwi.nl/oai/asset/7885/7885A.pdf
  53. F.G. Tricomi, Fonctions hypergéométriques confluentes. Mémorial des sciences mathématiques, 140 (1960), p. 1-86. Available from http://www.numdam.org/item?id=MSM_1960__140__1_0
  54. R.C. Forrey, Computing the hypergeometric function, J. Comp. Phys. 137, 79-100 (1997). Available as http://physics.bk.psu.edu/pub/papers/hyper.pdf, Fortran code from http://physics.bk.psu.edu/codes.html.
  55. N.M. Temme, The Numerical Computation of the Confluent Hypergeometric Function U(a,b,z), Numerische Mathematik, 41, p.63-82, 1983. Available from http://www.digizeitschriften.de/en/dms/toc/?PPN=PPN362160546_0041 or http://oai.cwi.nl/oai/asset/10717/10717A.pdf
  56. E.J. Weniger, Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series, 1989, Comput. Phys. Rep. 10, pp.189-371; available as http://arxiv.org/pdf/math/0306302v1.pdf
  57. T. Fessler, W.F. Ford, D.A. Smith, HURRY: An Acceleration Algorithm for Scalar Sequences and Series, ACM TOMS, Vol.9, No.3, 1983, pp.346-354. Fortran source available from http://netlib.org/toms/602
  58. W. Gautschi, On mean convergence of extended Lagrange interpolation, J. Comput. Appl. Math. 43, 1992, pp.19-35. Available as http://www.cs.purdue.edu/homes/wxg/selected_works/section_03/132.pdf
  59. J.C. Mason. Chebyshev polynomials of second, third and fourth kinds. J. Comput. Appl. Math. 49, 1993, 169-178. Available from doi: 10.1016/0377-0427(93)90148-5
  60. W. Kahan, To Solve a Real Cubic Equation (Lecture Notes for a Numerical Analysis Course), 1986. Available as http://www.dtic.mil/dtic/tr/fulltext/u2/a206859.pdf
  61. W. Kahan, Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign Bit, in The State of Art in Numerical Analysis, ed. by A. Iserles and M.J.D. Powell, 1987, pp. 165-211. Available as http://people.freebsd.org/~das/kahan86branch.pdf
  62. PARI/GP: Open Source Number Theory-oriented Computer Algebra System, available from http://pari.math.u-bordeaux.fr/
  63. R.M. Corless, J.H. Davenport, D.J. Jeffrey, and S.M. Watt, "According to Abramowitz and Stegun" or arccoth needn't be uncouth. SIGSAM BULLETIN: Communications on Computer Algebra, 34(2), June 2000. Available as http://www.sigsam.org/bulletin/articles/132/paper12.pdf or from the Ontario Research Centre for Computer Algebra http://www.orcca.on.ca/TechReports/2000/TR-00-17.html
  64. T.J. Dekker, A Floating-point technique for extending the available precision. Numerische Mathematik, 18, 224-242, 1971. Available from http://www.digizeitschriften.de/en/dms/toc/?PPN=PPN362160546_0018
  65. S. Linnainmaa, Software for Doubled-Precision Floating-Point Computations. ACM TOMS 7 (1981), pp. 272-283. doi: 10.1145/355958.355960