Gamma function

From Rosetta Code
Task
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

Translation of: C

- Stirling & Spouge methods.

Translation of: python

- Lanczos method.

Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386

<lang algol68># 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

AutoHotkey

Search autohotkey.com: function
Source: AutoHotkey forum by Laszlo <lang autohotkey> /* Here is the upper incomplete Gamma function. Omitting or setting the second parameter to 0 we get the (complete) Gamma function. The code is based on: "Computation of Special Functions" Zhang and Jin, John Wiley and Sons, 1996

  • /

SetFormat FloatFast, 0.9e

Loop 10 

  MsgBox % GAMMA(A_Index/3) "`n" GAMMA(A_Index*10)

GAMMA(a,x=0) {  ; upper incomplete gamma: Integral(t**(a-1)*e**-t, t = x..inf) 

  If (a > 171 || x < 0) 
     Return 2.e308   ; overflow 

  xam := x > 0 ? -x+a*ln(x) : 0 
  If (xam > 700) 
     Return 2.e308   ; overflow 

  If (x > 1+a) {     ; no need for gamma(a) 
     t0 := 0, k := 60 
     Loop 60 
         t0 := (k-a)/(1+k/(x+t0)), --k 
     Return exp(xam) / (x+t0) 
  } 

  r := 1, ga := 1.0  ; compute ga = gamma(a) ... 
  If (a = round(a))  ; if integer: factorial 
     If (a > 0) 
        Loop % a-1 
           ga *= A_Index 
     Else            ; negative integer 
        ga := 1.7976931348623157e+308 ; Dmax 
  Else {             ; not integer 
     If (abs(a) > 1) { 
        z := abs(a) 
        m := floor(z) 
        Loop %m% 
            r *= (z-A_Index) 
        z -= m 
     } 
     Else 
        z := a 

     gr := (((((((((((((((((((((((       0.14e-14 
         *z - 0.54e-14)             *z - 0.206e-13)          *z + 0.51e-12) 
         *z - 0.36968e-11)          *z + 0.77823e-11)        *z + 0.1043427e-9) 
         *z - 0.11812746e-8)        *z + 0.50020075e-8)      *z + 0.6116095e-8) 
         *z - 0.2056338417e-6)      *z + 0.1133027232e-5)    *z - 0.12504934821e-5) 
         *z - 0.201348547807e-4)    *z + 0.1280502823882e-3) *z - 0.2152416741149e-3) 
         *z - 0.11651675918591e-2)  *z + 0.7218943246663e-2) *z - 0.9621971527877e-2) 
         *z - 0.421977345555443e-1) *z + 0.1665386113822915) *z - 0.420026350340952e-1) 
         *z - 0.6558780715202538)   *z + 0.5772156649015329) *z + 1 

     ga := 1.0/(gr*z) * r 
     If (a < -1) 
        ga := -3.1415926535897931/(a*ga*sin(3.1415926535897931*a)) 
  } 

  If (x = 0)         ; complete gamma requested 
     Return ga 

  s := 1/a           ; here x <= 1+a 
  r := s 
  Loop 60 { 
     r *= x/(a+A_Index) 
     s += r 
     If (abs(r/s) < 1.e-15) 
        break 
  } 
  Return ga - exp(xam)*s

}

/* The 10 results shown: 2.678938535e+000 1.354117939e+000 1.0 8.929795115e-001 9.027452930e-001 3.628800000e+005 1.216451004e+017 8.841761994e+030 2.039788208e+046 6.082818640e+062

1.000000000e+000 1.190639348e+000 1.504575489e+000 2.000000000e+000 2.778158479e+000 1.386831185e+080 1.711224524e+098 8.946182131e+116 1.650795516e+136 9.332621544e+155

  • /

</lang>

C

Library: libgsl

This implements Stirling's approximation and Spouge's approximation.

<lang c>#include <stdio.h>

  1. include <stdlib.h>
  2. include <math.h>
  3. 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.

Works with: Fortran version 2008
Works with: Fortran version 95 with extensions

<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



 recursive 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) * lacz_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.67893744     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

Groovy

Translation of: Ada

<lang groovy>a = [ 1.00000000000000000000, 0.57721566490153286061, -0.65587807152025388108, 

    -0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,
    -0.00962197152787697356, 0.00721894324666309954, -0.00116516759185906511,
    -0.00021524167411495097, 0.00012805028238811619, -0.00002013485478078824,
    -0.00000125049348214267, 0.00000113302723198170, -0.00000020563384169776,
     0.00000000611609510448, 0.00000000500200764447, -0.00000000118127457049,
     0.00000000010434267117, 0.00000000000778226344, -0.00000000000369680562,
     0.00000000000051003703, -0.00000000000002058326, -0.00000000000000534812,
     0.00000000000000122678, -0.00000000000000011813, 0.00000000000000000119,
     0.00000000000000000141, -0.00000000000000000023, 0.00000000000000000002].reverse()

def gamma = { 1.0 / a.inject(0) { sm, a_i -> sm * (it - 1) + a_i } }

(1..10).each{ printf("% 1.9e\n", gamma(it / 3.0)) } </lang> Sample output: 

  2.678938535e+00
  1.354117939e+00
  1.000000000e+00
  8.929795116e-01
  9.027452930e-01
  1.000000000e+00
  1.190639349e+00
  1.504575488e+00
  2.000000000e+00
  2.778158479e+00

Haskell

From HaskellWiki (compatible license):

The Gamma and Beta function as described in 'Numerical Recipes in C++', the approximation is taken from [Lanczos, C. 1964 SIAM Journal on Numerical Analysis, ser. B, vol. 1, pp. 86-96]

<lang haskell>cof :: [Double] cof = [76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.001208650973866179,-0.000005395239384953]

ser :: Double ser = 1.000000000190015

gammaln :: Double -> Double gammaln xx = let tmp' = (xx+5.5) - (xx+0.5)*log(xx+5.5) 

                ser' = ser + sum (zipWith (/) cof [xx+1..])
            in -tmp' + log(2.5066282746310005 * ser' / xx)</lang>

J

This code shows the built-in method, which works for any value (positive, negative and complex numbers).

<lang J> gamma=. !@<: </lang>

The following direct coding of the task comes from http://www.jsoftware.com/jwiki/Essays/Stirling's%20Approximation (copyright issues ignored, since this is simple math):

<lang J>sbase =: %:@(2p1&%) * %&1x1 ^ ] scorr =: 1 1r12 1r288 _139r51840 _571r2488320&p.@% stirlg=: sbase * scorr</lang>

Checking against !@<: we can see that this approximation loses accuracy for small arguments 

   (,. stirlg ,. !@<:) 10 1p1 1x1 1.5 1
     10   362880   362880
3.14159  2.28803  2.28804
2.71828  1.56746  1.56747
    1.5 0.886155 0.886227
      1 0.999499        1

Mathematica

This code shows the built-in method, which works for any value (positive, negative and complex numbers). <lang mathematica> 

Gamma[x]

</lang>

Output integers and half-integers (a space is multiplication in Mathematica): 

1/2	Sqrt[pi]
1	1
3/2	Sqrt[pi]/2
2	1
5/2	(3 Sqrt[pi])/4
3	2
7/2	(15 Sqrt[pi])/8
4	6
9/2	(105 Sqrt[pi])/16
5	24
11/2	(945 Sqrt[pi])/32
6	120
13/2	(10395 Sqrt[pi])/64
7	720
15/2	(135135 Sqrt[pi])/128
8	5040
17/2	(2027025 Sqrt[pi])/256
9	40320
19/2	(34459425 Sqrt[pi])/512
10	362880

Output approximate numbers: 

0.1	9.51351
0.2	4.59084
0.3	2.99157
0.4	2.21816
0.5	1.77245
0.6	1.48919
0.7	1.29806
0.8	1.16423
0.9	1.06863
1.	1.

Output complex numbers: 

I	-0.15495-0.498016 I
2 I	0.00990244-0.075952 I
3 I	0.0112987-0.00643092 I
4 I	0.00173011+0.00157627 I
5 I	-0.000271704+0.000339933 I

Modula-3

Translation of: Ada

<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

OCaml

<lang ocaml>

let e = exp 1. let pi = 4. *. atan 1. let sqrttwopi = sqrt (2. *. pi)

module Lanczos = struct 

 (* Lanczos method *)
 (* Coefficients used by the GNU Scientific Library *)
 let g = 7.
 let c = [|0.99999999999980993; 676.5203681218851; -1259.1392167224028;

771.32342877765313; -176.61502916214059; 12.507343278686905; -0.13857109526572012; 9.9843695780195716e-6; 1.5056327351493116e-7|] 

 let rec ag z d =
   if d = 0 then c.(0) +. ag z 1
   else if d < 8 then c.(d) /. (z +. float d) +. ag z (succ d)
   else c.(d) /. (z +. float d)

 let gamma z =
   let z = z -. 1. in
   let p = z +. g +. 0.5 in
   sqrttwopi *. p ** (z +. 0.5) *. exp (-. p) *. ag z 0

end

module Stirling = struct 

 (* Stirling method *)
 let gamma z =
   sqrttwopi /. sqrt z *. (z /. e) ** z

end

module Stirling2 = struct 

 (* Extended Stirling method seen in Abramowitz and Stegun *)
 let d = [|1./.12.; 1./.288.; -139./.51840.; -571./.2488320.|]

 let rec corr z x n =
   if n < Array.length d - 1 then d.(n) /. x +. corr z (x *. z) (succ n)
   else d.(n) /. x

 let gamma z = Stirling.gamma z *. (1. +. corr z z 0)

end

let mirror gma z = 

 if z > 0.5 then gma z
 else pi /. sin (pi *. z) /. gma (1. -. z)

let _ = 

 Printf.printf "z\t\tLanczos\t\tStirling\tStirling2\n";
 for i = 1 to 20 do
   let z = float i /. 10. in
   Printf.printf "%-10.8g\t%10.8e\t%10.8e\t%10.8e\n"
   		  z

(mirror Lanczos.gamma z) (mirror Stirling.gamma z) (mirror Stirling2.gamma z) 

 done;
 for i = 1 to 7 do
   let z = 10. *. float i in
   Printf.printf "%-10.8g\t%10.8e\t%10.8e\t%10.8e\n"
   		  z

(Lanczos.gamma z) (Stirling.gamma z) (Stirling2.gamma z) 

 done


</lang>

Output 

z               Lanczos         Stirling        Stirling2
0.1             9.51350770e+00  1.04050843e+01  9.52104183e+00
0.2             4.59084371e+00  5.07399275e+00  4.59686230e+00
0.3             2.99156899e+00  3.35033954e+00  2.99844028e+00
0.4             2.21815954e+00  2.52705781e+00  2.22775889e+00
0.5             1.77245385e+00  2.06636568e+00  1.78839014e+00
0.6             1.48919225e+00  1.30715886e+00  1.48277536e+00
0.7             1.29805533e+00  1.15905329e+00  1.29508068e+00
0.8             1.16422971e+00  1.05337097e+00  1.16270541e+00
0.9             1.06862870e+00  9.77061508e-01  1.06778308e+00
1               1.00000000e+00  9.22137009e-01  9.99499469e-01
1.1             9.51350770e-01  8.83489953e-01  9.51037997e-01
1.2             9.18168742e-01  8.57755335e-01  9.17964058e-01
1.3             8.97470696e-01  8.42678259e-01  8.97331287e-01
1.4             8.87263818e-01  8.36744549e-01  8.87165485e-01
1.5             8.86226925e-01  8.38956553e-01  8.86155384e-01
1.6             8.93515349e-01  8.48693242e-01  8.93461840e-01
1.7             9.08638733e-01  8.65621472e-01  9.08597702e-01
1.8             9.31383771e-01  8.89639635e-01  9.31351590e-01
1.9             9.61765832e-01  9.20842722e-01  9.61740068e-01
2               1.00000000e+00  9.59502176e-01  9.99978981e-01
10              3.62880000e+05  3.59869562e+05  3.62879997e+05
20              1.21645100e+17  1.21139342e+17  1.21645100e+17
30              8.84176199e+30  8.81723653e+30  8.84176199e+30
40              2.03978821e+46  2.03554316e+46  2.03978821e+46
50              6.08281864e+62  6.07268919e+62  6.08281864e+62
60              1.38683119e+80  1.38490639e+80  1.38683119e+80
70              1.71122452e+98  1.70918858e+98  1.71122452e+98

Octave

<lang octave>function g = lacz_gamma(a, cg=7) 

 p = [ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, \
       771.32342877765313, -176.61502916214059, 12.507343278686905, \
       -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 ];
 x=a;
 if ( x < 0.5 )
   g = pi / ( sin(pi*x) * lacz_gamma(1.0-x) );
 else
   x = x - 1.0;
   t = p(1);
   for i=1:(cg+1)
     t = t + p(i+1)/(x+double(i));
   endfor
   w = x + double(cg) + 0.5;
   g = sqrt(2.0*pi) * w**(x+0.5) * exp(-w) * t;
 endif

endfunction


for i = 1:10 

 printf("%f %f\n", gamma(i/3.0), lacz_gamma(i/3.0));

endfor</lang>

Output 

2.678939 2.678939
1.354118 1.354118
1.000000 1.000000
0.892980 0.892980
0.902745 0.902745
1.000000 1.000000
1.190639 1.190639
1.504575 1.504575
2.000000 2.000000
2.778158 2.778158

Which suggests that the built-in gamma uses the same approximation.

Python

Translation of: Ada

<lang python>_a = ( 1.00000000000000000000, 0.57721566490153286061, -0.65587807152025388108, 

        -0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,
        -0.00962197152787697356, 0.00721894324666309954, -0.00116516759185906511,
        -0.00021524167411495097, 0.00012805028238811619, -0.00002013485478078824,
        -0.00000125049348214267, 0.00000113302723198170, -0.00000020563384169776,
         0.00000000611609510448, 0.00000000500200764447, -0.00000000118127457049,
         0.00000000010434267117, 0.00000000000778226344, -0.00000000000369680562,
         0.00000000000051003703, -0.00000000000002058326, -0.00000000000000534812,
         0.00000000000000122678, -0.00000000000000011813, 0.00000000000000000119,
         0.00000000000000000141, -0.00000000000000000023, 0.00000000000000000002
      )

def gamma (x): 

  y  = float(x) - 1.0;
  sm = _a[-1];
  for an in _a[-2::-1]:
     sm = sm * y + an;
  return 1.0 / sm;

if __name__ == '__main__': 

   for i in range(1,11):
       print "  %20.14e" % gamma(i/3.0)

</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.77815847933857e+00

R

Lanczos' approximation is a loose

Translation of: Octave

.

<lang r> stirling <- function(z) sqrt(2*pi/z) * (exp(-1)*z)^z

nemes <- function(z) sqrt(2*pi/z) * (exp(-1)*(z + (12*z - (10*z)^-1)^-1))^z

lanczos <- function(z) { 

  if(length(z) > 1)
  {
     sapply(z, lanczos)
  } else
  {
    g <- 7
     p <- c(0.99999999999980993, 676.5203681218851, -1259.1392167224028,
       771.32342877765313, -176.61502916214059, 12.507343278686905,
       -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7)
     z <- as.complex(z) 
     if(Re(z) < 0.5) 
     {
        pi / (sin(pi*z) * lanczos(1-z)) 
     } else
     {
        z <- z - 1
        x <- p[1] + sum(p[-1]/seq.int(z+1, z+g+1))
        tt <- z + g + 0.5
        sqrt(2*pi) * tt^(z+0.5) * exp(-tt) * x
     }
  }

}

spouge <- function(z, a=49) { 

  if(length(z) > 1)
  {
     sapply(z, spouge)
  } else
  {
     z <- z-1
     k <- seq.int(1, a-1)
     ck <- rep(c(1, -1), len=a-1) / factorial(k-1) * (a-k)^(k-0.5) * exp(a-k)
     (z + a)^(z+0.5) * exp(-z - a) * (sqrt(2*pi) + sum(ck/(z+k)))
  }

}

  1. Checks

z <- (1:10)/3 all.equal(gamma(z), stirling(z)) # Mean relative difference: 0.07181942 all.equal(gamma(z), nemes(z)) # Mean relative difference: 0.003460549 all.equal(as.complex(gamma(z)), lanczos(z)) # TRUE all.equal(gamma(z), spouge(z)) # TRUE data.frame(z=z, stirling=stirling(z), nemes=nemes(z), lanczos=lanczos(z), spouge=spouge(z), builtin=gamma(z)) </lang> 

          z  stirling     nemes      lanczos    spouge   builtin
1  0.3333333 2.1569760 2.6290752 2.6789385+0i 2.6789385 2.6789385
2  0.6666667 1.2028507 1.3515736 1.3541179+0i 1.3541179 1.3541179
3  1.0000000 0.9221370 0.9996275 1.0000000+0i 1.0000000 1.0000000
4  1.3333333 0.8397427 0.8928835 0.8929795+0i 0.8929795 0.8929795
5  1.6666667 0.8591902 0.9027098 0.9027453+0i 0.9027453 0.9027453
6  2.0000000 0.9595022 0.9999831 1.0000000+0i 1.0000000 1.0000000
7  2.3333333 1.1491064 1.1906296 1.1906393+0i 1.1906393 1.1906393
8  2.6666667 1.4584904 1.5045690 1.5045755+0i 1.5045755 1.5045755
9  3.0000000 1.9454032 1.9999951 2.0000000+0i 2.0000000 2.0000000
10 3.3333333 2.7097638 2.7781544 2.7781585+0i 2.7781585 2.7781585

Ruby

Translation of: Ada

<lang ruby>$a = [ 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 ]

def gamma(x) 

 y = Float(x) - 1
 1.0 / $a.reverse.inject {|sum, an| sum = sum * y + an}

end

(1..10).each {|i| puts format("%.14e", gamma(i/3.0))}</lang>

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.77815847933857e+00

Tcl

Works with: Tcl version 8.5


Library: tcllib

<lang tcl>package require math package require math::calculus

  1. gamma(1) and gamma(1.5)

set f 1.0 set f2 [expr {sqrt(acos(-1.))/2.}]

for {set x 1.0} {$x <= 10.0} {set x [expr {$x + 0.5}]} { 

   # method 1 - numerical integration, Romberg's method, special
   #            case for an improper integral

   set g1 [math::calculus::romberg  \
               [list apply {{x t} {expr {$t ** ($x-1) * exp(-$t)}}} $x] \
               0 1 -relerror 1e-8]
   set g2 [math::calculus::romberg_infinity \
               [list apply {{x t} {expr {$t ** ($x-1) * exp(-$t)}}} $x] \
               1 Inf -relerror 1e-8]
   set gamma [expr {[lindex $g1 0] + [lindex $g2 0]}]

   # method 2 - library function

   set libgamma [expr {exp([math::ln_Gamma $x])}]

   # method 3 - special forms for integer and half-integer arguments

   if {$x == entier($x)} {
       puts [format {%4.1f %13.6f %13.6f %13.6f} $x $gamma $libgamma $f]
       set f [expr $f * $x]
   } else {
       puts [format {%4.1f %13.6f %13.6f %13.6f} $x $gamma $libgamma $f2]
       set f2 [expr $f2 * $x]
   }

}</lang> Output: 

 1.0      1.000000      1.000000      1.000000
 1.5      0.886228      0.886227      0.886227
 2.0      1.000000      1.000000      1.000000
 2.5      1.329340      1.329340      1.329340
 3.0      2.000000      2.000000      2.000000
 3.5      3.323351      3.323351      3.323351
 4.0      6.000000      6.000000      6.000000
 4.5     11.631731     11.631728     11.631728
 5.0     24.000009     24.000000     24.000000
 5.5     52.342778     52.342778     52.342778
 6.0    120.000000    120.000000    120.000000
 6.5    287.885278    287.885278    287.885278
 7.0    720.000001    720.000000    720.000000
 7.5   1871.254311   1871.254305   1871.254306
 8.0   5040.000032   5039.999999   5040.000000
 8.5  14034.298267  14034.407291  14034.407293
 9.0  40320.000705  40319.999992  40320.000000
 9.5 119292.464880 119292.461971 119292.461995
10.0 362880.010950 362879.999927 362880.000000
Retrieved from "https://rosettacode.org/wiki/Gamma_function?oldid=103893"