Gamma function
You are encouraged to solve this task according to the task description, using any language you may know.
Implement one algorithm (or more) to compute the Gamma (Γ) function (in the real field only).If your language has the function as builtin or you know a library which has it, compare your implementation's results with the results of the builtin/library function. The Gamma function can be defined as
This suggests a straightforward (but inefficient) way of computing the Γ through numerical integration. Better suggested methods:
Ada
The implementation uses Taylor series coefficients of . The coefficients are taken from Mathematical functions and their approximations by Yudell L. Luke. <lang ada> function Gamma (X : Long_Float) return Long_Float is
A : constant array (0..29) of Long_Float := ( 1.00000_00000_00000_00000, 0.57721_56649_01532_86061, -0.65587_80715_20253_88108, -0.04200_26350_34095_23553, 0.16653_86113_82291_48950, -0.04219_77345_55544_33675, -0.00962_19715_27876_97356, 0.00721_89432_46663_09954, -0.00116_51675_91859_06511, -0.00021_52416_74114_95097, 0.00012_80502_82388_11619, -0.00002_01348_54780_78824, -0.00000_12504_93482_14267, 0.00000_11330_27231_98170, -0.00000_02056_33841_69776, 0.00000_00061_16095_10448, 0.00000_00050_02007_64447, -0.00000_00011_81274_57049, 0.00000_00001_04342_67117, 0.00000_00000_07782_26344, -0.00000_00000_03696_80562, 0.00000_00000_00510_03703, -0.00000_00000_00020_58326, -0.00000_00000_00005_34812, 0.00000_00000_00001_22678, -0.00000_00000_00000_11813, 0.00000_00000_00000_00119, 0.00000_00000_00000_00141, -0.00000_00000_00000_00023, 0.00000_00000_00000_00002 ); Y : constant Long_Float := X - 1.0; Sum : Long_Float := A (A'Last);
begin
for N in reverse A'First..A'Last - 1 loop Sum := Sum * Y + A (N); end loop; return 1.0 / Sum;
end Gamma; </lang> Test program: <lang ada> with Ada.Text_IO; use Ada.Text_IO; with Gamma;
procedure Test_Gamma is begin
for I in 1..10 loop Put_Line (Long_Float'Image (Gamma (Long_Float (I) / 3.0))); end loop;
end Test_Gamma; </lang> Sample output:
2.67893853470775E+00 1.35411793942640E+00 1.00000000000000E+00 8.92979511569249E-01 9.02745292950934E-01 1.00000000000000E+00 1.19063934875900E+00 1.50457548825154E+00 1.99999999999397E+00 2.77815847933858E+00
ALGOL 68
- Stirling & Spouge methods.
- Lanczos method.
<lang algol># Coefficients used by the GNU Scientific Library # []LONG REAL p = ( LONG 0.99999 99999 99809 93,
LONG 676.52036 81218 851, -LONG 1259.13921 67224 028, LONG 771.32342 87776 5313, -LONG 176.61502 91621 4059, LONG 12.50734 32786 86905, -LONG 0.13857 10952 65720 12, LONG 9.98436 95780 19571 6e-6, LONG 1.50563 27351 49311 6e-7);
PROC lanczos gamma = (LONG REAL in z)LONG REAL: (
# Reflection formula # LONG REAL z := in z; IF z < LONG 0.5 THEN long pi / (long sin(long pi*z)*lanczos gamma(1-z)) ELSE z -:= 1; LONG REAL x := p[1]; FOR i TO UPB p - 1 DO x +:= p[i+1]/(z+i) OD; LONG REAL t = z + UPB p - LONG 1.5; long sqrt(2*long pi) * t**(z+LONG 0.5) * long exp(-t) * x FI
);
PROC taylor gamma = (LONG REAL x)LONG REAL: BEGIN # good for values between 0 and 1 #
[]LONG REAL a = ( LONG 1.00000 00000 00000 00000, LONG 0.57721 56649 01532 86061, -LONG 0.65587 80715 20253 88108, -LONG 0.04200 26350 34095 23553, LONG 0.16653 86113 82291 48950, -LONG 0.04219 77345 55544 33675, -LONG 0.00962 19715 27876 97356, LONG 0.00721 89432 46663 09954, -LONG 0.00116 51675 91859 06511, -LONG 0.00021 52416 74114 95097, LONG 0.00012 80502 82388 11619, -LONG 0.00002 01348 54780 78824, -LONG 0.00000 12504 93482 14267, LONG 0.00000 11330 27231 98170, -LONG 0.00000 02056 33841 69776, LONG 0.00000 00061 16095 10448, LONG 0.00000 00050 02007 64447, -LONG 0.00000 00011 81274 57049, LONG 0.00000 00001 04342 67117, LONG 0.00000 00000 07782 26344, -LONG 0.00000 00000 03696 80562, LONG 0.00000 00000 00510 03703, -LONG 0.00000 00000 00020 58326, -LONG 0.00000 00000 00005 34812, LONG 0.00000 00000 00001 22678, -LONG 0.00000 00000 00000 11813, LONG 0.00000 00000 00000 00119, LONG 0.00000 00000 00000 00141, -LONG 0.00000 00000 00000 00023, LONG 0.00000 00000 00000 00002 ); LONG REAL y = x - 1; LONG REAL sum := a [UPB a]; FOR n FROM UPB a - 1 DOWNTO LWB a DO sum := sum * y + a [n] OD; 1/sum
END # taylor gamma #;
LONG REAL long e = long exp(1);
PROC sterling gamma = (LONG REAL n)LONG REAL: ( # improves for values much greater then 1 #
long sqrt(2*long pi/n)*(n/long e)**n
);
PROC factorial = (LONG INT n)LONG REAL: (
IF n=0 OR n=1 THEN 1 ELIF n=2 THEN 2 ELSE n*factorial(n-1) FI
);
REF[]LONG REAL fm := NIL;
PROC sponge gamma = (LONG REAL x)LONG REAL: (
INT a = 12; # alter to get required precision # REF []LONG REAL fm := NIL; LONG REAL res; IF fm :=: REF[]LONG REAL(NIL) THEN fm := HEAP[0:a-1]LONG REAL; fm[0] := long sqrt(LONG 2*long pi); FOR k TO a-1 DO fm[k] := (((k-1) MOD 2=0) | 1 | -1) * long exp(a-k) *
(a-k) **(k-LONG 0.5) / factorial(k-1)
OD FI; res := fm[0]; FOR k TO a-1 DO res +:= fm[k] / ( x + k ) OD; res *:= long exp(-(x+a)) * (x+a)**(x + LONG 0.5); res/x
);
FORMAT real fmt = $g(-real width, real width - 2)$; FORMAT long real fmt16 = $g(-17, 17 - 2)$; # accurate to about 16 decimal places #
[]STRING methods = ("Genie", "Lanczos", "Sponge", "Taylor","Stirling");
printf(($11xg12xg12xg13xg13xgl$,methods));
FORMAT sample fmt = $"gamma("g(-3,1)")="f(real fmt)n(UPB methods-1)(", "f(long real fmt16))l$; FORMAT sqr sample fmt = $"gamma("g(-3,1)")**2="f(real fmt)n(UPB methods-1)(", "f(long real fmt16))l$; FORMAT sample exp fmt = $"gamma("g(-3)")="g(-15,11,0)n(UPB methods-1)(","g(-18,14,0))l$;
PROC sample = (LONG REAL x)[]LONG REAL:
(gamma(SHORTEN x), lanczos gamma(x), sponge gamma(x), taylor gamma(x), sterling gamma(x));
FOR i FROM 1 TO 20 DO
LONG REAL x = i / LONG 10; printf((sample fmt, x, sample(x))); IF i = 5 THEN # insert special case of a half # printf((sqr sample fmt, x, gamma(SHORTEN x)**2, lanczos gamma(x)**2, sponge gamma(x)**2, taylor gamma(x)**2, sterling gamma(x)**2)) FI
OD; FOR x FROM 10 BY 10 TO 70 DO
printf((sample exp fmt, x, sample(x)))
OD</lang> Output:
Genie Lanczos Sponge Taylor Stirling gamma(0.1)=9.5135076986687, 9.513507698668730, 9.513507698668731, 9.513509522249043, 5.697187148977169 gamma(0.2)=4.5908437119988, 4.590843711998802, 4.590843711998803, 4.590843743037192, 3.325998424022393 gamma(0.3)=2.9915689876876, 2.991568987687590, 2.991568987687590, 2.991568988322729, 2.362530036269620 gamma(0.4)=2.2181595437577, 2.218159543757688, 2.218159543757688, 2.218159543764845, 1.841476335936235 gamma(0.5)=1.7724538509055, 1.772453850905517, 1.772453850905516, 1.772453850905353, 1.520346901066281 gamma(0.5)**2=3.1415926535898, 3.141592653589795, 3.141592653589793, 3.141592653589216, 2.311454699581843 gamma(0.6)=1.4891922488128, 1.489192248812817, 1.489192248812817, 1.489192248812758, 1.307158857448356 gamma(0.7)=1.2980553326476, 1.298055332647558, 1.298055332647558, 1.298055332647558, 1.159053292113920 gamma(0.8)=1.1642297137253, 1.164229713725304, 1.164229713725303, 1.164229713725303, 1.053370968425609 gamma(0.9)=1.0686287021193, 1.068628702119320, 1.068628702119319, 1.068628702119319, 0.977061507877695 gamma(1.0)=1.0000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 0.922137008895789 gamma(1.1)=0.9513507698669, 0.951350769866873, 0.951350769866873, 0.951350769866873, 0.883489953168704 gamma(1.2)=0.9181687423998, 0.918168742399761, 0.918168742399760, 0.918168742399761, 0.857755335396591 gamma(1.3)=0.8974706963063, 0.897470696306277, 0.897470696306277, 0.897470696306277, 0.842678259448392 gamma(1.4)=0.8872638175031, 0.887263817503076, 0.887263817503075, 0.887263817503064, 0.836744548637082 gamma(1.5)=0.8862269254528, 0.886226925452758, 0.886226925452758, 0.886226925452919, 0.838956552526496 gamma(1.6)=0.8935153492877, 0.893515349287691, 0.893515349287690, 0.893515349288799, 0.848693242152574 gamma(1.7)=0.9086387328533, 0.908638732853291, 0.908638732853290, 0.908638732822421, 0.865621471793884 gamma(1.8)=0.9313837709802, 0.931383770980243, 0.931383770980242, 0.931383769950169, 0.889639635287994 gamma(1.9)=0.9617658319074, 0.961765831907388, 0.961765831907387, 0.961765815012982, 0.920842721894229 gamma(2.0)=1.0000000000000, 1.000000000000000, 0.999999999999999, 1.000000010045742, 0.959502175744492 gamma( 10)= 3.6288000000e5, 3.6288000000000e5, 3.6288000000000e5, 4.051218760300e-7, 3.5986956187410e5 gamma( 20)= 1.216451004e17, 1.216451004088e17, 1.216451004088e17, 1.07701514977e-18, 1.211393423381e17 gamma( 30)= 8.841761994e30, 8.841761993740e30, 8.841761993739e30, 7.98891286318e-23, 8.817236530765e30 gamma( 40)= 2.039788208e46, 2.039788208120e46, 2.039788208120e46, 6.97946184592e-25, 2.035543161237e46 gamma( 50)= 6.082818640e62, 6.082818640343e62, 6.082818640342e62, 1.81016585713e-26, 6.072689187876e62 gamma( 60)= 1.386831185e80, 1.386831185457e80, 1.386831185457e80, 9.27306839649e-28, 1.384906385829e80 gamma( 70)= 1.711224524e98, 1.711224524281e98, 1.711224524281e98, 7.57303907062e-29, 1.709188578191e98
C
This implements Stirling's approximation and Spouge's approximation.
<lang c>#include <stdio.h>
- include <stdlib.h>
- include <math.h>
- include <gsl/gsl_sf_gamma.h>
/* very simple approximation */ double st_gamma(double x) {
return sqrt(2.0*M_PI/x)*pow(x/M_E, x);
}
int factorial(int n) {
if ( (n==0) || (n==1) ) return 1; if ( n==2 ) return 2; return n * factorial(n-1);
}
double sp_gamma(double x, double **fm) {
const int a = 6; static double *c = NULL; int k; double res;
if ( c == NULL ) { c = malloc(sizeof(double) * a); *fm = c; c[0] = sqrt(2.0*M_PI); for(k=1; k < a; k++) { c[k] = (((k-1)%2==0)?1:-1) * exp(a-k) *
pow(a-k, k-1.0/2.0) / (double)factorial(k-1);
} } res = c[0]; for(k=1; k < a; k++) { res += c[k] / ( x + k ); } res *= exp(-(x+a)) * pow(x+a, x + 1.0/2.0); return res/x;
}
int main() {
double x, *fm;
printf("%15s%15s%15s\n", "Stirling", "Spouge", "GSL"); for(x=1.0; x <= 10.0; x+=1.0) { printf("%15.8lf%15.8lf%15.8lf\n", st_gamma(x/3.0), sp_gamma(x/3.0, &fm),
gsl_sf_gamma(x/3.0));
} free(fm); return 0;
}</lang>
Output:
Stirling Spouge GSL 2.15697602 2.67893851 2.67893853 1.20285073 1.35411792 1.35411794 0.92213701 0.99999998 1.00000000 0.83974270 0.89297950 0.89297951 0.85919025 0.90274527 0.90274529 0.95950218 0.99999998 1.00000000 1.14910642 1.19063932 1.19063935 1.45849038 1.50457544 1.50457549 1.94540320 1.99999994 2.00000000 2.70976382 2.77815839 2.77815848
Fortran
This code shows two methods: Numerical Integration through Simpson formula, and Lanczos approximation. The results of testing are printed altogether with the values given by the function gamma; this function is defined in the Fortran 2008 standard, and supported by GNU Fortran (and other vendors) as extension; if not present in your compiler, you can remove the last part of the print in order to get it compiled with any Fortran 95 compliant compiler.
<lang fortran>program ComputeGammaInt
implicit none
integer :: i
write(*, "(3A15)") "Simpson", "Lanczos", "Builtin" do i=1, 10 write(*, "(3F15.8)") my_gamma(i/3.0), lacz_gamma(i/3.0), gamma(i/3.0) end do
contains
pure function intfuncgamma(x, y) result(z) real :: z real, intent(in) :: x, y z = x**(y-1.0) * exp(-x) end function intfuncgamma
function my_gamma(a) result(g) real :: g real, intent(in) :: a
real, parameter :: small = 1.0e-4 integer, parameter :: points = 100000
real :: infty, dx, p, sp(2, points), x integer :: i logical :: correction
x = a
correction = .false. ! value with x<1 gives \infty, so we use ! \Gamma(x+1) = x\Gamma(x) ! to avoid the problem if ( x < 1.0 ) then correction = .true. x = x + 1 end if
! find a "reasonable" infinity... ! we compute this integral indeed ! \int_0^M dt t^{x-1} e^{-t} ! where M is such that M^{x-1} e^{-M} ≤ \epsilon infty = 1.0e4 do while ( intfuncgamma(infty, x) > small ) infty = infty * 10.0 end do
! using simpson dx = infty/real(points) sp = 0.0 forall(i=1:points/2-1) sp(1, 2*i) = intfuncgamma(2.0*(i)*dx, x) forall(i=1:points/2) sp(2, 2*i - 1) = intfuncgamma((2.0*(i)-1.0)*dx, x) g = (intfuncgamma(0.0, x) + 2.0*sum(sp(1,:)) + 4.0*sum(sp(2,:)) + & intfuncgamma(infty, x))*dx/3.0
if ( correction ) g = g/a
end function my_gamma
function lacz_gamma(a) result(g) real, intent(in) :: a real :: g
real, parameter :: pi = 3.14159265358979324 integer, parameter :: cg = 7
! these precomputed values are taken by the sample code in Wikipedia, ! and the sample itself takes them from the GNU Scientific Library real, dimension(0:8), parameter :: p = & (/ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, & 771.32342877765313, -176.61502916214059, 12.507343278686905, & -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 /)
real :: t, w, x integer :: i
x = a
if ( x < 0.5 ) then g = pi / ( sin(pi*x) * gamma(1.0-x) ) else x = x - 1.0 t = p(0) do i=1, cg+2 t = t + p(i)/(x+real(i)) end do w = x + real(cg) + 0.5 g = sqrt(2.0*pi) * w**(x+0.5) * exp(-w) * t end if end function lacz_gamma
end program ComputeGammaInt</lang>
Output:
Simpson Lanczos Builtin 2.65968132 2.67893839 2.67893839 1.35269761 1.35411859 1.35411787 1.00000060 1.00000024 1.00000000 0.88656044 0.89297968 0.89297950 0.90179849 0.90274525 0.90274531 0.99999803 1.00000036 1.00000000 1.19070935 1.19063985 1.19063926 1.50460517 1.50457609 1.50457561 2.00000286 2.00000072 2.00000000 2.77815390 2.77816010 2.77815843
Modula-3
<lang modula3>MODULE Gamma EXPORTS Main;
FROM IO IMPORT Put; FROM Fmt IMPORT Extended, Style;
PROCEDURE Taylor(x: EXTENDED): EXTENDED =
CONST a = ARRAY [0..29] OF EXTENDED { 1.00000000000000000000X0, 0.57721566490153286061X0, -0.65587807152025388108X0, -0.04200263503409523553X0, 0.16653861138229148950X0, -0.04219773455554433675X0, -0.00962197152787697356X0, 0.00721894324666309954X0, -0.00116516759185906511X0, -0.00021524167411495097X0, 0.00012805028238811619X0, -0.00002013485478078824X0, -0.00000125049348214267X0, 0.00000113302723198170X0, -0.00000020563384169776X0, 0.00000000611609510448X0, 0.00000000500200764447X0, -0.00000000118127457049X0, 0.00000000010434267117X0, 0.00000000000778226344X0, -0.00000000000369680562X0, 0.00000000000051003703X0, -0.00000000000002058326X0, -0.00000000000000534812X0, 0.00000000000000122678X0, -0.00000000000000011813X0, 0.00000000000000000119X0, 0.00000000000000000141X0, -0.00000000000000000023X0, 0.00000000000000000002X0 }; VAR y := x - 1.0X0; sum := a[LAST(a)];
BEGIN FOR i := LAST(a) - 1 TO FIRST(a) BY -1 DO sum := sum * y + a[i]; END; RETURN 1.0X0 / sum; END Taylor;
BEGIN
FOR i := 1 TO 10 DO Put(Extended(Taylor(FLOAT(i, EXTENDED) / 3.0X0), style := Style.Sci) & "\n"); END;
END Gamma.</lang> Output:
2.6789385347077490e+000 1.3541179394264005e+000 1.0000000000000000e+000 8.9297951156924930e-001 9.0274529295093360e-001 1.0000000000000000e+000 1.1906393487589992e+000 1.5045754882515399e+000 1.9999999999939684e+000 2.7781584793385790e+000
Tcl
<lang tcl>package require math
proc ::tcl::mathfunc::Gamma x {
expr {exp([::math::ln_Gamma $x])}
}</lang> Testing: <lang tcl>% expr Gamma(1) 0.9999999999658638 % expr Gamma(2) 0.9999999999218119 % expr Gamma(0.2) 4.590843711986148 % expr Gamma(30) 8.841761992087971e+30</lang>