W.Ehrhardt: AMath / Special Functions

Last update Feb. 21, 2012

Basic information about the AMath package

Contents

Introduction
AMath functions
Special functions
AMTools functions
References

Introduction

The AMath package contains Pascal/Delphi source for accurate mathematical methods without using multi precision arithmetic. Please note that the high accuracy can only be achieved with the rmNearest rounding mode; it decreases if other modes are used. The main parts of the package are the AMath unit and the Special Functions units:

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 standard one argument elementary transcendental functions have peak relative errors less than 2.2e-19, values for power(x,y) are 2.1e-19 (for |x|,|y| < 1000) and 3.4e-19 (for |x|,|y| < 2000).
The AMTools unit provides accurate and reliable tools for finding zeroes and local minima of functions, numerical integration of one-dimensional functions, and solving quadratic equations:
- The functions localmin, mbrent, and fmin (differing in parameter count / ease of use) use Brent's algorithm with guaranteed convergence for finding a local minimum of a function f in an interval (a,b). The algorithm combines golden section search and successive parabolic interpolation using only function (not derivative) evaluations.
- The functions zbrent and zeroin use the Brent/Dekker algorithm with guaranteed convergence for finding a zero of a continuous function f in the interval [a,b], when f(a) and f(b) have different signs. The algorithm is based on a combination of successive interpolations and bisection.
- The qag* procedures are Pascal translations of Quadpack algorithms by R. Piessens, E. de Doncker-Kapenga, C.W. Überhuber, D. Kahaner. These routines perform global adaptive quadrature of functions over finite or infinite intervals based on Gauss-Kronrod rules for the subintervals and acceleration by Wynn's epsilon algorithm, they can handle rather difficult integrals including integrand functions with local difficulties such as a discontinuities and integrable singularities. quagk is a simple general purpose shell for the qag* routines.
  
  The Quadpack algorithm qawc computes the Cauchy principal value of f(x)/(x-c) using a globally adaptive strategy and modified Clenshaw-Curtis integration on the subintervals containing the point x=c.
  
  The adaptive quanc8 algorithm by G.E. Forsythe, M.A. Malcolm, C.B. Moler estimates the integral of a smooth function over a finite interval using a Newton-Cotes rule.
  
  The procedures intdeo and intdei use the Double Exponential (DE) transformation (developed by M. Mori, T. Ooura, and others) for automatic quadrature of f(x) over the infinite interval (a,+INF) for functions with and without oscillatory factors resp. The DE transformation is very powerful if the abscissas a+x_i can be processed with high accuracy; therefore the two routines are best used with a=0, otherwise there can be severe roundoff problems (a+x_i=a).
- The squad functions accurately solve quadratic equations with double coefficients; they 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 and SpecFunX implement special functions for extended and double precision; the following function groups are available: Bessel functions and related, elliptic integrals/functions and theta functions, gamma function and related, error function and related, zeta functions and polylogarithms, exponential integrals and related, orthogonal polynomials and related, statistical distributions, and other special functions.

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.

The AMath Special Functions Reference Manual with Implementation Notes is available as specialfunctions.pdf.

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)
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
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
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)
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
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
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
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
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
exprel             Return exprel(x) = (exp(x) - 1)/x,  1 for x=0
expx2              Return exp(x*|x|) with damped error amplification in computing exp of the product
gd                 Return the Gudermannian function gd(x)
hav                Return the haversine hav(x) = (1 - cos(x))/2
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
log2               Return base 2 logarithm of x
logbase            Return base b logarithm of x
logN               Delphi alias for logbase
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)
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)
sinh               Return the hyperbolic sine of x, accurate even 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, 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 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
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 or denormal, return m=s, 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
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
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

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)
PolEvalX           Evaluate polynomial; return a[0] + a[1]*x + ... + a[n-1]*x^(n-1)
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.

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

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
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
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 for |k|<=1, |x*k|<=1 for |k|>1
jacobi_arcdn(x,k)    Return the inverse Jacobi elliptic function arcdn(x,k), |x|<=1, x^2 + k^2 > 1
jacobi_arcsn(x,k)    Return the inverse Jacobi elliptic function arcsn(x,k), |x|<=1 for |k|<=1, |x*k|<=1 for |k|>1
jacobi_cn(x,k)       Return the Jacobi elliptic function cn(x,k)
jacobi_dn(x,k)       Return the Jacobi elliptic function dn(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)
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)
erfi(x)              Return the imaginary error function erfi(x) = erf(ix)/i
erfc_inv(x)          Return the inverse function of erfc, erfc(erfc_inv(x)) = x, 0 < x < 2
erfce(x)             Return the exponentially scaled complementary error function erfce(x) = exp(x^2)*erfc(x)
erf_inv(x)           Return the inverse function of erf, erf(erf_inv(x)) = x, -1 < x < 1
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(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

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

beta(x,y)            Return beta(x,y)=gamma(x)*gamma(y)/gamma(x+y)
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, a>0, b>0, 0 <= x <= 1
igammap(a,x)         Return the normalised 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 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
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
pentagamma(x)        Return the pentagamma function psi'''(x), INF if x is a negative integer
pochhammer(x,y)      Return the Pochhammer symbol gamma(x+y)/gamma(x)
polygamma(n,x)       Return the polygamma function: n'th derivative of psi; n>=0, x>0 for n>10.
psi(x)               Return the psi (digamma) function of x, INF if x is a non-positive integer
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 and polylogarithms

cl2(x)               Return the Clausen function: integral(-ln(2*|sin(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
etam1(s)             Return Dirichlet eta(s)-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.
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
ti2(x)               Return the inverse tangent integral, Ti_2(x) = integral(arctan(t)/t, t=0..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>=0, 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 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
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, |x| <= 1
legendre_plm(l,m,x)  Return the associated Legendre polynomial P_lm(x), |x| <= 1
legendre_q(l,x)      Return Q_l(x), the Legendre function of the 2nd kind, degree l >=0, |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

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

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
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
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
RiemannR(x)          Return the Riemann prime counting function R(x), x >= 1/16

AMTools functions

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.

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.

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.

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.

References

[HMF]: M. Abramowitz, I.A. Stegun. Handbook of Mathematical Functions. New York, 1970, http://www.math.sfu.ca/~cbm/aands/
Intel, IA-32 Architecture Software Developer's Manual Volume 2A: Instruction Set Reference, A-M http://www.intel.com/products/processor/manuals/
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
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
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
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
Cephes Mathematical Library, Version 2.8, http://www.moshier.net/#Cephes or http://www.netlib.org/cephes/
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
N.J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., Philadelphia, 2002. http://www.maths.manchester.ac.uk/~higham/asna/
R. Bulirsch, Numerical Calculation of Elliptic Integrals and Elliptic Functions, Numerische Mathematik 7, 78-90, 1965
R. Bulirsch, Numerical Calculation of Elliptic Integrals and Elliptic Functions, part III. Numerische Mathematik 13, 305-315, 1969
B.C. Carlson, Computing Elliptic Integrals by Duplication, Numerische Mathematik 33, 1-16, 1979
W.H. Press et al, Numerical Recipes in C, 2nd ed., Cambridge, 1992. http://www.nrbook.com/a/bookcpdf.html
SLATEC Common Mathematical Library, Version 4.1, July 1993 (general purpose mathematical and statistical routines written in Fortran 77), http://www.netlib.org/slatec/
W. Kahan, On the Cost of Floating-Point Computation Without Extra-Precise Arithmetic, 2004. http://www.eecs.berkeley.edu/~wkahan/Qdrtcs.pdf
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
P.H. Sterbenz, Floating-Point Computation. Englewood Cliffs, 1974. Chap.9.3: Carefully written programs/quadratic equation, p.246ff
W.Y. Sit, Quadratic Programming? 1997. http://www.mmrc.iss.ac.cn/ascm/ascm03/sample.pdf
Boost C++ Libraries, Release 1.42.0, 2010. http://www.boost.org/
Special functions by Wayne Fullerton, http://www.netlib.org/fn/. Almost identical to the FNLIB subset of SLATEC [14].
GNU Scientific Library, GSL-1.14 (March 2010), http://www.gnu.org/software/gsl/
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
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
D. Veberic, Having Fun with Lambert W(x) Function, 2010. http://arxiv.org/abs/1003.1628
I. Smith, Examples.xls/txt Version 3.3.4, Personal communication, 2010
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://repository.cwi.nl/search/fullrecord.php?publnr=10080
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
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
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/
[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/
R.E. Crandall, Note on fast polylogarithm computation, 2006. http://www.reed.edu/~crandall/papers/Polylog.pdf
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
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.
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/
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
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
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
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
R: A Language and Environment for Statistical Computing, Version 2.11.1, http://www.r-project.org/
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
C. Ferreira, J.L. López, Asymptotic expansions of the Hurwitz-Lerch zeta function, J. Math. Anal. Appl. 298 (2004), no. 1, 210-224. Available from http://pcmap.unizar.es/~chelo/investig/publicaciones/chelo/Hurwitz-Lerch/original.pdf