A Simple Gamma Function Approximation
Daniel D. Warner and Robert Simms
June 8, 1999
A First Approximation to the Gamma Function
We can develop an accurate approximation to the Gamma Function following the guidelines in the book, Computer Approximations. We use the best uniform rational approximation to the Gamma Function on the interval [2,3] and use the recurrence relation to reduce other values to this interval.
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.
>
gam23a := proc( xx )
local p, q, x, z;
# Best Uniform Approximation of Degree 7 over Degree 3 on [2, 3] - Index 5234
p := [47.0033628547426492, 25.4191188301489618, 14.1777319483221048, 3.9372323903267183,
1.2660529275231088, .1971127900125424, .0460485036059242, .00320611206970855]:
q := [47.0033628547206298, 5.5468333236447965, -7.5252629999787867, 1]:
x := evalf( xx );
z := (p[1]+x*(p[2]+x*(p[3]+x*(p[4]+x*(p[5]+x*(p[6]+x*(p[7]+x*p[8]))))))) / (q[1]+x*(q[2]+x*(q[3]+x*q[4])));
RETURN( evalf(z) );
end:
>
Gamma1 := proc( xx )
local p, q, x, z;
x := evalf( xx );
if (2 <= x) and (x <= 3) then
z := gam23a(x-2);
elif (0 <= x) and (x < 2) then
z := Gamma1(x+1)/x;
elif (3 < x) and (x <= 12) then
z := (x-1)*Gamma1(x-1);
elif (12 < x) then
z := exp((x-1/2)*ln(x) -x +1/2*ln(2*Pi) +1/(12*x) -1/(360*x^3) +1/(1260*x^5) -1/(1680*x^7));
else # z < 0
z := -Pi/(Gamma1(-x)*x*sin(Pi*x));
fi;
RETURN( evalf( z ) );
end:
Error Calculation
> Digits := 22;
First check several special values.
>
Gamma1(1); Gamma1(2); Gamma1(3); Gamma1(4);
Gamma1(1/2); Gamma1(1/2) - evalf(sqrt(Pi));
Gamma1(3/2); Gamma1(3/2) - evalf(sqrt(Pi)/2);
Gamma1(5); Gamma1(10); Gamma1(15);
> GAMMA(5); GAMMA(10); GAMMA(15);
Plotting the absolute error over the interval [2, 3] confirms the equioscillating behaviour of a best uniform approximation. Moreover, the approximation has 11 degrees of freedom (q[4] can always be taken to be 1), and the error curve has 12 extreme points on the closed interval, which is also predicted by the theory. The relative error decreases as the argument increases since the Gamma Function is increasing on this interval.
> plot((GAMMA - Gamma1), 2..3, title="Absolute Error");
> plot(((GAMMA - Gamma1)/GAMMA), 2..3, title="Relative Error");
Now consider the absolute and relative error on the positive real axis.
> plot((GAMMA - Gamma1), 1..4, title="Asbsolute Error");
> plot(((GAMMA - Gamma1)/GAMMA), 1..4, title="Relative Error");
> plot(((GAMMA - Gamma1)/GAMMA), 1/4..15, title="Relative Error");
A Second Approximation to the Gamma Function
While the first approximation procedure is well behaved in the sense of relative error, it is not quite accurate enough for the FX 2.0. Approximation #5239 with an indicated precision of 15.81 would appear to be more appropriate. Let's redefine gam23 using this approximation and rerun the calculations.
>
gam23b := proc( xx )
local p, q, x, z;
# Best Uniform Approximation of Degree 6 over Degree 6 on [2, 3] - Index 5239
p := [ 3786.0105034825724547, 2077.4597938941873209, 893.58180452374981423, 222.11239616801179483,
48.954346227909938052, 6.1260674503360842987, 0.77807958561330057586]:
q := [ 3786.0105034825719725, 476.79386050368791516, -867.23098753110299445, 83.550058667919769574,
50.788475328895409737, -13.400414785781348262, 1]:
x := evalf( xx );
z := ((((((p[7]*x+p[6])*x+p[5])*x+p[4])*x+p[3])*x+p[2])*x+p[1]) /
((((((q[7]*x+q[6])*x+q[5])*x+q[4])*x+q[3])*x+q[2])*x+q[1]);
RETURN( evalf(z) );
end:
It is also necessary to adjust the breakpoint and the number of terms for the asymptotic expansion.
>
Gamma2 := proc( xx )
local p, q, x, z;
x := evalf( xx );
if (2 <= x) and (x <= 3) then
z := gam23b(x-2);
elif (0 <= x) and (x < 2) then
z := Gamma2(x+1)/x;
elif (3 < x) and (x <= 16) then
z := (x-1)*Gamma2(x-1);
elif (15 < x) then
z := exp((x-1/2)*ln(x)-x+1/2*ln(2*Pi)+1/(12*x)-1/(360*x^3)+1/(1260*x^5)-1/(1680*x^7)+1/(1188*x^9));
else # z < 0
z := -Pi/(Gamma2(-x)*x*sin(Pi*x));
fi;
RETURN( evalf( z ) );
end:
>
Error Calculation
> plot((GAMMA - Gamma2), 2..3, title="Absolute Error");
> plot(((GAMMA - Gamma2)/GAMMA), 1/4..20, title="Relative Error", numpoints=1024);
Clearly this approximation procedure would be adequate for the FX 2.0.
Of course, the recursion steps would need to be replaced by simple loops for efficiency. Also, if this approximation procedure is adopted for the log of the Gamma Function, then the range reduction would need to be carried out as a loop of multiplications followed by a single log of the product rather than a loop of logarithm calculations.
>