W.Ehrhardt: AMath functions

Last update July 26, 2010

Basic information about the AMath package

Contents

Introduction
AMath functions
QuadSolv functions
Special functions
References

Introduction

The main parts of the AMath 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, the value for power(x,y) is 2.8e-18 (for |x|,|y| < 1000) and 4.8e-18 (for |x|,|y| < 2000).
Unit QuadSolv accurately solves quadratic equations with double coefficients. The functions implement ideas of G.E. Forsythe, W. Kahan, P.H. Sterbenz (high precision calculation of discriminant, scaling by powers of two etc).
The units SpecFun, SpecFunX, SFCommon implement special functions for extended and double precision; the following function groups are available: Bessel functions, elliptic integrals, gamma function and related, error function and related, exponential integrals 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.

AMath functions

Elementary numerical functions

cbrt               Return the cube root of x
ceil               Return the smallest integer >= x; |x| <= MaxLongint
ceilx              Return the smallest integer >= x
floor              Return the largest integer <= x; |x| <= MaxLongint
floorx             Return the largest integer <= x
fmod               Return x mod y, y<>0, sign(result) = sign(x)
hypot              Return sqrt(x*x + y*y)
intpower           Return x^n; via binary exponentiation
modf               Return frac(x) and trunc(x) in ip, |x| <= MaxLongint
nroot              Return the nth root of x; n<>0, x >= 0 if n is even
sqrt1pm1           Return sqrt(1+x)-1, 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
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
cosf1              Return 1-cos(x), accurate even for x near multiples of 2*Pi
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
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
exp3               Return 3^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
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
power              Return x^y; if frac(y)<>0 then x must be > 0
powm1              Return x^y - 1; special code for small x,y
pow1pm1            Return (1+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
sincPi             Return the normalized cardinal sine sincPi(x) = sin(Pi*x)/(Pi*x)
sincos             Return accurate values s=sin(x), c=cos(x)
sincosPi           Return s=sin(Pi*x), c=cos(Pi*x); (s,x)=(0,1) for abs(x) >= 2^64
sinh               Return the hyperbolic sine of x, 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

Floating point functions

copysign           Return abs(x)*sign(y)
copysignd          Return abs(x)*sign(y)
frexp              Return the mantissa m and exponent e of x with x = m*2^e, 0.5 <= abs(m) < 1; if x is 0, +-INF, NaN or denormal, return m=x, e=0
frexpd             Return the mantissa m and exponent e of d with d = m*2^e, 0.5 <= abs(m) < 1; if d is 0, +-INF, NaN or denormal, return m=d, e=0
ilogb              Return base 2 exponent of x, -MaxLongint for x=0/denormal, MaxLongint if INF or Nan
IsInf              Return true if x is +INF or -INF
IsInfD             Return true if d is +INF or -INF
IsNaN              Return true if x is a NaN
IsNaND             Return true if d is a NaN
IsNaNorInf         Return true if x is a NaN or infinite
IsNaNorInfD        Return true if d is a NaN or infinite
ldexp              Return x*2^e
ldexpd             Return d*2^e
predd              Return next representable double after d in the direction -Inf
predx              Return next representable extended after x in the direction -Inf
scalbn             Return x*2^e
succd              Return next representable double after d in the direction +Inf
succx              Return next representable extended after x in the direction +Inf

FPU control functions

Get8087CW          Return the FPC control word
GetExceptionMask   Return the current exception mask
GetPrecisionMode   Return the current precision control mode
GetRoundMode       Return the current rounding mode
Set8087CW          Set new FPU control word
SetExceptionMask   Set new exception mask
SetPrecisionMode   Set new precision control mode and return the old precision
SetRoundMode       Set new rounding mode and return the old mode

Polynomial, Vector, Statistic Operations

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

Dbl2Hex            Return d as a big-endian hex string
Ext2Hex            Return x as a big-endian hex string
fisEQd             Return true if x and y are bit-identical
fisEQx             Return true if x and y are bit-identical
fisNEd             Return true if x and y are not bit-identical
fisNEx             Return true if x and y are not bit-identical
maxd               Return the maximum of two doubles; x,y <> NAN
maxx               Return the maximum of two extendeds; x,y <> NAN
mind               Return the minimum of two doubles; x,y <> NAN
minx               Return the minimum of two extendeds; x,y <> NAN

QuadSolv functions

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

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

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

Special functions

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

Bessel functions

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; not suitable for large n or x.
bessel_jn(n,x)     Return J_n(x), the Bessel function of the 1st kind, order n; not suitable for large n or x.
bessel_kn(n,x)     Return K_n(x), the modified Bessel function of the 2nd kind, order n, x>0, not suitable for large n.
bessel_yn(n,x)     Return Y_n(x), the Bessel function of the 2nd kind, order n, x>0, not suitable for large n or x.

Elliptic integrals

ell_rc(x,y)        Return Carlson's degenerate elliptic integral RC; x>=0, y<>0
ell_rf(x,y,z)      Return Carlson's elliptic integral of the 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 the complementary complete elliptic integral of the 1st kind, kc<>0
cel2(kc,a,b)       Return the general complete elliptic integral of the 2nd kind, kc<>0
cel(kc,p,a,b)      Return the general complete elliptic integral, kc<>0, Cauchy principle value if p<0
el1(x,kc)          Return the complementary complete elliptic integral of the 1st kind
el2(x,kc,a,b)      Return the general complete elliptic integral of the 2nd kind, kc<>0
el3(x,kc,p)        Return the incomplete elliptic integral of the 3rd kind, 1+p*x^2<>0

EllipticF(z,k)     Return the incomplete elliptic integral of the 1st kind; |z|<=1, |k*z|<1
EllipticK(k)       Return the complete elliptic integral of the 1st kind, |k| < 1
EllipticCK(kc)     Return the complementary complete elliptic integral of the 1st kind, kc<>0
EllipticE(z,k)     Return the incomplete elliptic integrals of the 2nd kind, |z|<=1, |k|<=1
EllipticEC(k)      Return the complete elliptic integral of the 2nd kind, |k| < 1
EllipticCE(kc)     Return the complementary complete elliptic integral of the 2nd kind
EllipticPi(z,nu,k) Return the incomplete elliptic integral of the 3rd kind, |z|<=1, |k*z|<1
EllipticPiC(nu,k)  Return the complete elliptic integral of the 3rd kind, |k|<1
EllipticCPi(nu,k)  Return the complementary complete elliptic integral of the 3rd kind, k<>0

Error function and 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)
erf_inv(x)         Return the inverse function of erf, erf(erf_inv(x)) = x, -1 < x < 1
erfc(x)            Return complementary error function erfc(x) = 1-erf(x)
erfc_inv(x)        Return the inverse function of erfc, erfc(erfc_inv(x)) = x, 0 < x < 2
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)

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 modified cosine integral, cin(x) = integral((1-cos(t))/t, t=0..x)
cinh(x)            Return the modified 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)
en(n,x)            Return the exponential integral E_n(x) = integral(exp(-x*t)/t^n, t=1..Inf), x > 0
li(x)              Return the logarithmic integral li(x)  = PV-integral(1/ln(t), t=0..x), x > 1
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

gamma(x)           Return gamma(x), |x| <= MAXGAM; invalid if x is non-positive integer
gamma1pm1(x)       Return gamma(1+x)-1, accurate even for x near 0
fac(n)             Return the factorial n!, n+1 <= MAXGAM; 0 if n<0
invgamma(x)        Return 1/gamma(x)
signgamma(x)       Return sign(gamma(x)), useless for 0 or negative integer
lngamma(x)         Return ln(|gamma(x)|), |x| <= MAXLGM, invalid if x is non-positive integer.
psi(x)             Return the psi (digamma) function of x, INF if x is non-positive integer
beta(x,y)          Return beta(x,y)=gamma(x)*gamma(y)/gamma(x+y)
ibeta(a,b,x)       Return the normalized incomplete beta function, a>0, b>0, 0 <= x <= 1
lnbeta(x,y)        Return the logarithm of |beta(x,y)|=|gamma(x)*gamma(y)/gamma(x+y)|
zeta(s)            Return the Riemann zeta function at s, s<>1
zeta1p(x)          Return the Riemann zeta function at 1+x, x<>0

Statistical distributions

beta_pdf(a,b,x)    Return the probability density function of the beta distribution with parameters a and b
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
f_pdf(nu1,nu2,x)   Return the probability density function of F distribution; x >= 0, nu1, nu2 > 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
normstd_pdf(x)     Return the std. normal probability density function exp(-x^2/2)/sqrt(2*Pi)
normstd_cdf(x)     Return the standard normal distribution function
normstd_inv(y)     Return the inverse standard normal distribution function, 0 < y < 1
t_pdf(nu,x)        Return the probability density function of Student's t distribution, nu>0
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

Other special functions

agm(x,y)           Return the arithmetic-geometric mean of |x| and |y|
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)
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
sncndn(x,mc...)    Calculate the Jacobian elliptic functions sn,cn,dn for argument x and complementary parameter mc.
ti2(x)             Return the inverse tangent integral, Ti_2(x) = integral(arctan(t)/t, t=0..x)

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
D. Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, 1991, http://docs.sun.com/source/806-3568/ncg_goldberg.html; extended and edited reprint via http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.22.6768
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