A Log Gamma Function Approximation
Daniel D. Warner
June 9, 1999
The Cody-Stoltz Approximation to the Log Gamma Function
The following approximation to the log of the Gamma Function was developed by W. J. Cody and L. Stoltz at the Argonne National Laboratory in 1988. This routine is part of W. J. Cody's specfun package. For background information see the following references.
1. Chebyshev Approximations for the Natural Logarithm of the Gamma Function by W. J. Cody and K. E. Hillstrom, Math. Comp. 21, 1967, pp. 198-203.
2. DGAMMA/DLGAMA, ANL/AMD Program ANLC366S, by K. E. Hillstrom, May, 1969.
3. Computer Approximations by John F. Hart, E. W. Cheney, Charles L. Lawson, Hans J. Maehly, Charles K. Mestenyi, John R. Rice, Henry G. Thacher, Jr., and Chritoph Witzgall. Published by John Wiley & Sons, Inc. 1968.
>
Digits := 22;
XINF := 9.99999999999999*10^99; # lambda for FX 2.0
EPS := 8.0*10^(-14); # Machine epsilon for FX 2.0 is not properly defined.
XBIG := fsolve(lnGAMMA(x) = XINF, x);
FRTBIG := sqrt(sqrt(XBIG));
PNT68 := 0.6796875;
>
lgam1 := proc( xx )
local p, q, d, x, xn, xd, k, z;
# Almost best Uniform Approximation of Degree 7 over Degree 8 on [0.5 - 1.5]
# Note that we assume that xx has already been translated by -1 so that -0.5 <= xx <= 0.5.
p := [ 4.945235359296727046734888, 201.8112620856775083915565, 2290.838373831346393026739,
11319.67205903380828685045, 28557.24635671635335736389, 38484.96228443793359990269,
26377.48787624195437963534, 7225.813979700288197698961 ];
q := [ 67.48212550303777196073036, 1113.332393857199323513008, 7738.757056935398733233834,
27639.87074403340708898585, 54993.10206226157329794414, 61611.22180066002127833352,
36351.27591501940507276287, 8785.536302431013170870835 ];
d := -0.5772156649015328605195174;
x := evalf( xx );
xn := 0;
xd := 1;
for k from 1 to 8 do
xn := xn*x + p[k];
xd := xd*x + q[k];
od;
z := x*(d + x*(xn/xd));
RETURN( evalf(z) );
end:
>
lgam2 := proc( xx )
local p, q, d, x, xn, xd, k, z;
# Almost best Uniform Approximation of Degree 7 over Degree 8 on [1.5 - 4]
# Note that we assume that xx has already been translated by -2 so that -0.5 <= xx <= 2.
p := [ 4.974607845568932035012064, 542.4138599891070494101986, 15506.93864978364947665077,
184793.2904445632425417223, 1088204.769468828767498470, 3338152.967987029735917223,
5106661.678927352456275255, 3074109.054850539556250927 ];
q := [ 183.0328399370592604055942, 7765.049321445005871323047, 133190.3827966074194402448,
1136705.821321969608938755, 5267964.117437946917577538, 13467014.54311101692290052,
17827365.30353274213975932, 9533095.591844353613395747 ];
d := 0.4227843350984671393993777;
x := evalf( xx );
xn := 0;
xd := 1;
for k from 1 to 8 do
xn := xn*x + p[k];
xd := xd*x + q[k];
od;
z := x*(d + x*(xn/xd));
RETURN( evalf(z) );
end:
>
lgam3 := proc( xx )
local p, q, d, x, xn, xd, k, z;
# Almost best Uniform Approximation of Degree 7 over Degree 8 on [4 - 12]
# Note that we assume that xx has already been translated by -4 so that 0 <= xx <= 8.
p := [ 14745.02166059939948905062, 2426813.369486704502836312, 121475557.4045093227939592,
2663432449.630976949898078, 29403789566.34553899906876, 170266573776.5398868392998,
492612579337.7430887588120, 560625185622.3951465078242 ];
q := [ 2690.530175870899333379843, 639388.5654300092398984238, 41355999.30241388052042842,
1120872109.616147941376570, 14886137286.78813811542398, 101680358627.2438228077304,
341747634550.7377132798597, 446315818741.9713286462081 ];
d := 1.791759469228055000094023;
x := evalf( xx );
xn := 0;
xd := -1;
for k from 1 to 8 do
xn := xn*x + p[k];
xd := xd*x + q[k];
od;
z := d + x*(xn/xd);
RETURN( evalf(z) );
end:
>
lgam4 := proc( xx )
local c, x, xs, k, z;
# Assymptotic Expansion with 7 terms on [12 - INF]
# Note that we assume that xx has already been translated by -4 so that 0 <= xx <= 8.
c := [-0.001910444077728, 0.00084171387781295, -0.0005952379913043012,
0.000793650793500350248, -0.002777777777777777777778, 0.083333333333333333333333,
0.0057083835261 ];
x := evalf( xx );
xs := x*x;
z := c[7];
for k from 1 to 6 do
z := z/xs + c[k];
od;
z := z/x;
RETURN( evalf(z) );
end:
>
LGamma1 := proc( xx )
local x, z;
x := evalf( xx );
if (x <= 0) or (XBIG <= x) then
RETURN( XINF );
elif (x <= EPS) then
RETURN( -log(x) );
elif (x < 0.5) then
z := -log(x) + lgam1(x);
elif (x < PNT68) then
z := -log(x) + lgam2(x-1);
elif (x < 1.5) then
z := lgam1(x-1);
elif ( x <= 4) then
z := lgam2(x-2);
elif (x < 12) then
z := lgam3(x-4);
elif (x < FRTBIG) then
z := lgam4(x) + (1/2)*log(2*Pi) + (x - 1/2)*log(x) - x;
else # FRTBIG <= x
z := (1/2)*log(2*Pi) + (x - 1/2)*log(x) - x;
fi;
RETURN( evalf( z ) );
end:
The following approximation was included in the Casio Technical Report on the ZX-938 dated April 23, 1997.
This is the Lanczos algorithm that is presented in the book Numerical Recipes in C.
>
LGamma2 := proc( xx )
local c, x, t, y, k, z;
# loggamma
c := [ 76.18009172947146, -86.50532032941677, 24.01409824083091,
-1.231739572450155, 0.001208650973866179, -0.000005395239384953 ];
x := evalf( xx );
t := (x+5.5) - (x+0.5)*log(x+5.5);
z := 1.000000000190015;
y := x + 1;
for k from 1 to 6 do
z := z + c[k]/y;
y := y + 1;
od;
# z := (x+0.5)*log(x+5.5) - (x+5.5) + log(2.50662827463100005*z/x);
z := (x+0.5)*log(x+5.5) - (x+5.5) + log(sqrt(2*Pi)*z/x);
RETURN( evalf(z) );
end:
> Digits:=20;lnGAMMA(1);lnGAMMA(1.5);lnGAMMA(2);LGamma2(1);LGamma2(1.5);LGamma2(2);
Error Calculation
Since this procedure uses different approximations for several different intervals, it is useful to first look at the absolute error on each interval and then assess the relative error on the entire range. Note that the approximation on the interval from 4 to 12 has the largest relative error. However the maximum relative error is well under the machine precision for the FX 2.0.
> plot((lnGAMMA - LGamma1), 0..0.5, title="Absolute Error");
> plot((lnGAMMA - LGamma1), 0.5..1.5, title="Absolute Error");
> plot((lnGAMMA - LGamma1), 1.5..4, title="Absolute Error");
> plot((lnGAMMA - LGamma1), 4..12, title="Absolute Error");
> plot((lnGAMMA - LGamma1), 12..25, title="Absolute Error");
> plot(((lnGAMMA - LGamma1)/lnGAMMA), 0..25, title="Relative Error");
Now for an examination of the Lanczos approximation, LGamma2(x).
The following results show that the approximation is not quite adequate for the FX 2.0 for small values of x,
and that both the absolute and relative error grow excessively for large values of x.
> plot((lnGAMMA - LGamma2), 0..5, title="Absolute Error");
> plot((lnGAMMA - LGamma2), 0..15, title="Absolute Error");
> plot(((lnGAMMA - LGamma2)/lnGAMMA), 0..25, title="Relative Error");
> plot((lnGAMMA - LGamma2), 0..500, title="Absolute Error");
> plot(((lnGAMMA - LGamma2)/lnGAMMA), 0..500, title="Relative Error");
>