# Gamma function

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 (${\displaystyle \Gamma }$) function (in the real field only).

If your language has the function as built-in or you know a library which has it, compare your implementation's results with the results of the built-in/library function.

The Gamma function can be defined as:

${\displaystyle \Gamma (x)=\displaystyle \int _{0}^{\infty }t^{x-1}e^{-t}dt}$

This suggests a straightforward (but inefficient) way of computing the ${\displaystyle \Gamma }$ through numerical integration.

Better suggested methods:

## 360 Assembly

For maximum compatibility, this program uses only the basic instruction set.

GAMMAT   CSECT         USING GAMMAT,R13SAVEAR   B     STM-SAVEAR(R15)         DC    17F'0'         DC    CL8'GAMMAT'STM      STM   R14,R12,12(R13)         ST    R13,4(R15)         ST    R15,8(R13)         LR    R13,R15*        ----  CODE         LE    F4,=E'0'         LH    R2,NILOOPI    EQU   *         AE    F4,=E'1'         xi=xi+1         LER   F0,F4         DE    F0,=E'10'        x=xi/10         STE   F0,X         LE    F6,X         SE    F6,=E'1'         xx=x-1.0         LH    R4,NT         BCTR  R4,0         SLA   R4,2         LE    F0,T(R4)                  STE   F0,SUM           sum=t(nt)         LH    R3,NT         BCTR  R3,0         SH    R4,=H'4'LOOPJ    CH    R3,=H'1'         for j=nt-1 downto 1         BL    ENDLOOPJ         LE    F0,SUM         MER   F0,F6            sum*xx         LE    F2,T(R4)         t(j)         AER   F0,F2         STE   F0,SUM           sum=sum*xx+t(j)         BCTR  R3,0         SH    R4,=H'4'         B     LOOPJENDLOOPJ EQU   *         LE    F0,=E'1'         DE    F0,SUM         STE   F0,GAMMA         gamma=1/sum         LE    F0,X         BAL   R14,CONVERT         MVC   BUF(8),CONVERTM         LE    F0,GAMMA         BAL   R14,CONVERT         MVC   BUF+9(13),CONVERTM         WTO   MF=(E,WTOMSG)		           BCT   R2,LOOPI*        ----  END CODE         CNOP  0,4         L     R13,4(0,R13)         LM    R14,R12,12(R13)         XR    R15,R15         BR    R14*        ----  DATANI       DC    H'30'NT       DC    AL2((TEND-T)/4)T        DC    E'1.00000000000000000000'         DC    E'0.57721566490153286061'         DC    E'-0.65587807152025388108'         DC    E'-0.04200263503409523553'         DC    E'0.16653861138229148950'         DC    E'-0.04219773455554433675'         DC    E'-0.00962197152787697356'         DC    E'0.00721894324666309954'         DC    E'-0.00116516759185906511'         DC    E'-0.00021524167411495097'         DC    E'0.00012805028238811619'         DC    E'-0.00002013485478078824'         DC    E'-0.00000125049348214267'         DC    E'0.00000113302723198170'         DC    E'-0.00000020563384169776'         DC    E'0.00000000611609510448'         DC    E'0.00000000500200764447'         DC    E'-0.00000000118127457049'         DC    E'0.00000000010434267117'         DC    E'0.00000000000778226344'         DC    E'-0.00000000000369680562'         DC    E'0.00000000000051003703'         DC    E'-0.00000000000002058326'         DC    E'-0.00000000000000534812'         DC    E'0.00000000000000122678'         DC    E'-0.00000000000000011813'         DC    E'0.00000000000000000119'         DC    E'0.00000000000000000141'         DC    E'-0.00000000000000000023'         DC    E'0.00000000000000000002'TEND     DS    0E X        DS    ESUM      DS    EGAMMA    DS    EWTOMSG   DS    0F         DC    AL2(L'BUF),XL2'0000'BUF      DC    CL80' '*        Subroutine             Convertion Float->DisplayCONVERT  CNOP  0,4         ME    F0,CONVERTC         STE   F0,CONVERTF         MVI   CONVERTS,X'00'         L     R9,CONVERTF         CH    R9,=H'0'         BZ    CONVERT7         BNL   CONVERT1         is negative?         MVI   CONVERTS,X'80'         N     R9,=X'7FFFFFFF'  sign bitCONVERT1 LR    R8,R9         N     R8,=X'00FFFFFF'         BNZ   CONVERT2         SR    R9,R9         B     CONVERT7CONVERT2 LR    R8,R9         N     R8,=X'FF000000'  characteristic         SRL   R8,24         CH    R8,=H'64'         BH    CONVERT3         SR    R9,R9         B     CONVERT7CONVERT3 CH    R8,=H'72'        2**32         BNH   CONVERT4         L     R9,=X'7FFFFFFF'  biggest         B     CONVERT6CONVERT4 SR    R8,R8         SLDL  R8,8         CH    R8,=H'72'         BL    CONVERT5         CH    R9,=H'0'         BP    CONVERT5         L     R9,=X'7FFFFFFF'         B     CONVERT6CONVERT5 SH    R8,=H'72'         LCR   R8,R8         SLL   R8,2         SRL   R9,0(R8)CONVERT6 SR    R8,R8         IC    R8,CONVERTS         CH    R8,=H'0'         sign bit set?         BZ    CONVERT7         LCR   R9,R9CONVERT7 ST    R9,CONVERTB         CVD   R9,CONVERTP         MVC   CONVERTD,=X'402020202120202020202020'          ED    CONVERTD,CONVERTP+2         MVC   CONVERTM(6),CONVERTD          MVI   CONVERTM+6,C'.'         MVC   CONVERTM+7(6),CONVERTD+6         BR    R14*CONVERTC DC    E'1E6'           X'45F42400'CONVERTF DS    FCONVERTB DS    F                CONVERTS DS    XCONVERTM DS    CL13       CONVERTD DS    CL12CONVERTP DS    PL8*         EQUREGS         EQUREGS REGS=FPR         END   GAMMAT
Output:
     0.1      9.513504
0.2      4.590844
0.3      2.991569
0.4      2.218160
0.5      1.772453
0.6      1.489192
0.7      1.298056
0.8      1.164229
0.9      1.068628
1.0      1.000000
1.1      0.951350
1.2      0.918168
1.3      0.897470
1.4      0.887263
1.5      0.886227
1.6      0.893515
1.7      0.908638
1.8      0.931383
1.9      0.961766
2.0      1.000000
2.1      1.046486
2.2      1.101803
2.3      1.166712
2.4      1.242169
2.5      1.329341
2.6      1.429626
2.7      1.544686
2.8      1.676492
2.9      1.827354
3.0      1.999999



The implementation uses Taylor series coefficients of Γ(x+1)-1, |x| < ∞. The coefficients are taken from Mathematical functions and their approximations by Yudell L. Luke.

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;

Test program:

with Ada.Text_IO;  use Ada.Text_IO;with Gamma; procedure Test_Gamma isbegin   for I in 1..10 loop      Put_Line (Long_Float'Image (Gamma (Long_Float (I) / 3.0)));   end loop;end Test_Gamma;
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
# 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/sumEND # 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))  FIOD;FOR x FROM 10 BY 10 TO 70 DO  printf((sample exp fmt, x, sample(x)))OD
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


## ANSI Standard BASIC

Translation of: BBC Basic
- Lanczos method.
100 DECLARE EXTERNAL FUNCTION FNlngamma110 120 DEF FNgamma(z) = EXP(FNlngamma(z))130 140 FOR x = 0.1 TO 2.05 STEP 0.1150    PRINT USING$("#.#",x), USING$("##.############", FNgamma(x))160 NEXT x170 END180 190 EXTERNAL FUNCTION FNlngamma(z)200 DIM lz(0 TO 6)210 RESTORE220 MAT READ lz230 DATA 1.000000000190015, 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.0012086509738662, -0.000005395239385240 IF z < 0.5 THEN 250    LET FNlngamma = LOG(PI / SIN(PI * z)) - FNlngamma(1.0 - z)260    EXIT FUNCTION270 END IF280 LET z = z - 1.0290 LET b = z + 5.5300 LET a = lz(0)310 FOR i = 1 TO 6320    LET a  = a + lz(i) / (z + i)330 NEXT i340 LET FNlngamma = (LOG(SQR(2*PI)) + LOG(a) - b) + LOG(b) * (z+0.5)350 END FUNCTION

## AutoHotkey

Search autohotkey.com: function
Source: AutoHotkey forum by Laszlo

/*Here is the upper incomplete Gamma function. Omitting or settingthe 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*/

## BBC BASIC

Uses the Lanczos approximation.

#### Stirling's Approximation

def gamma_by_stirling:  def pow(x): if x == 0 then 1 elif x == 1 then . else x * log | exp end;  ((1|atan) * 8) as $twopi | . as$x  | (($twopi/$x) | sqrt) * ( ($x / (1|exp)) | pow($x));

#### Examples

Stirling's method produces poor results, so to save space, the examples contrast the Taylor series and Lanczos methods with built-in tgamma:

def pad(n): tostring | . + (n - length) * " "; "                 i:      gamma                lanczos              tgamma",(range(1;11) | . / 3.0 | "\(pad(18)): \(gamma|pad(18)) : \(gamma_by_lanczos|pad(18)) : \(tgamma)")
Output:

## Modula-3

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


## МК-61/52

П9	9	П0	ИП9	ИП9	1	+	*	Вx	L0
05	1	+	П9	^	ln	1	-	*	ИП9
1	2	*	1/x	+	e^x	<->	/	2	пи
*	ИП9	/	КвКор	*	^	ВП	3	+	Вx
-	С/П


## Nim

const 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 ] proc gamma(x): float =  let y = x.float - 1.0  result = a[a.high]  for n in countdown(high(a) - 1, low(a)):    result = result * y + a[n]  result = 1.0 / result for i in 1..10:  echo gamma(i.float / 3.0)
Output:
2.678938534707748
1.3541179394264
1.0
0.8929795115692493
0.9027452929509336
1.0
1.190639348758999
1.50457548825154
1.999999999993968
2.778158479338573

## 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 0end 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
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

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+function/double.html">double(i));    endfor    w = x + function/double.html">double(cg) + 0.5;    g = sqrt(2.0*pi) * w**(x+0.5) * exp(-w) * t;  endifendfunction  for i = 1:10  printf("%f %f\n", gamma(i/3.0), lacz_gamma(i/3.0));endfor
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.

## Oforth

import: math [    676.5203681218851,  -1259.1392167224028, 771.32342877765313,   -176.61502916214059, 12.507343278686905, -0.13857109526572012,    9.9843695780195716e-6, 1.5056327351493116e-7] const: Gamma.Lanczos : gamma(z)| i t |   z 0.5 < ifTrue: [ Pi dup z * sin 1.0 z - gamma * / return ]   z 1.0 - ->z   0.99999999999980993 Gamma.Lanczos size loop: i [ i Gamma.Lanczos at z i + / + ]   z Gamma.Lanczos size + 0.5 - ->t   2 Pi * sqrt *    t z 0.5 + powf *   t neg exp * ;
Output:
>20 seq apply(#[ 10.0 / dup . gamma .cr ])
0.1 9.51350769866874
0.2 4.5908437119988
0.3 2.99156898768759
0.4 2.21815954375769
0.5 1.77245385090552
0.6 1.48919224881282
0.7 1.29805533264756
0.8 1.1642297137253
0.9 1.06862870211932
1 1
1.1 0.951350769866874
1.2 0.918168742399761
1.3 0.897470696306277
1.4 0.887263817503076
1.5 0.886226925452759
1.6 0.893515349287691
1.7 0.908638732853292
1.8 0.931383770980243
1.9 0.961765831907388
2 1

## PARI/GP

### Built-in

gamma(x)

### Double-exponential integration

[[+oo],k] means that the function approaches ${\displaystyle +\infty }$ as ${\displaystyle \exp(-kx).}$

Gamma(x)=intnum(t=0,[+oo,1],t^(x-1)/exp(t))

### Romberg integration

Gamma(x)=intnumromb(t=0,9,t^(x-1)/exp(t),0)+intnumromb(t=9,max(x,99)^9,t^(x-1)/exp(t),2)

### Stirling approximation

Stirling(x)=x--;sqrt(2*Pi*x)*(x/exp(1))^x

## Perl

use strict;use warnings;use constant pi => 4*atan2(1, 1);use constant e  => exp(1); # Normally would be:  use Math::MPFR# but this will use it if it's installed and ignore otherwisemy $have_MPFR = eval { require Math::MPFR; Math::MPFR->import(); 1; }; sub Gamma { my$z = shift;    my $method = shift // 'lanczos'; if ($method eq 'lanczos') {        use constant g => 9;        $z < .5 ? pi / sin(pi *$z) / Gamma(1 - $z,$method) :        sqrt(2* pi) *        ($z + g - .5)**($z - .5) *        exp(-($z + g - .5)) * do { my @coeff = qw{ 1.000000000000000174663 5716.400188274341379136 -14815.30426768413909044 14291.49277657478554025 -6348.160217641458813289 1301.608286058321874105 -108.1767053514369634679 2.605696505611755827729 -0.7423452510201416151527e-2 0.5384136432509564062961e-7 -0.4023533141268236372067e-8 }; my ($sum, $i) = (shift(@coeff), 0);$sum += $_ / ($z + $i++) for @coeff;$sum;        }    } elsif ($method eq 'taylor') {$z <  .5 ? Gamma($z+1,$method)/$z :$z > 1.5 ? ($z-1)*Gamma($z-1, $method) : do { my$s = 0; ($s *=$z-1) += $_ for qw{ 0.00000000000000000002 -0.00000000000000000023 0.00000000000000000141 0.00000000000000000119 -0.00000000000000011813 0.00000000000000122678 -0.00000000000000534812 -0.00000000000002058326 0.00000000000051003703 -0.00000000000369680562 0.00000000000778226344 0.00000000010434267117 -0.00000000118127457049 0.00000000500200764447 0.00000000611609510448 -0.00000020563384169776 0.00000113302723198170 -0.00000125049348214267 -0.00002013485478078824 0.00012805028238811619 -0.00021524167411495097 -0.00116516759185906511 0.00721894324666309954 -0.00962197152787697356 -0.04219773455554433675 0.16653861138229148950 -0.04200263503409523553 -0.65587807152025388108 0.57721566490153286061 1.00000000000000000000 }; 1/$s;	}    } elsif ($method eq 'stirling') { no warnings qw(recursion);$z < 100 ? Gamma($z + 1,$method)/$z : sqrt(2*pi*$z)*($z/e + 1/(12*e*$z))**$z /$z;    } elsif ($method eq 'MPFR') { my$result = Math::MPFR->new();        Math::MPFR::Rmpfr_gamma($result, Math::MPFR->new($z), 0);        $result; } else { die "unknown method '$method'" }} for my $method (qw(MPFR lanczos taylor stirling)) { next if$method eq 'MPFR' && !$have_MPFR; printf "%10s: ",$method;    print join(' ', map { sprintf "%.12f", Gamma($_/3,$method) } 1 .. 10);    print "\n";}
Output:
      MPFR: 2.678938534708 1.354117939426 1.000000000000 0.892979511569 0.902745292951 1.000000000000 1.190639348759 1.504575488252 2.000000000000 2.778158480438
lanczos: 2.678938534708 1.354117939426 1.000000000000 0.892979511569 0.902745292951 1.000000000000 1.190639348759 1.504575488252 2.000000000000 2.778158480438
taylor: 2.678938534708 1.354117939426 1.000000000000 0.892979511569 0.902745292951 1.000000000000 1.190639348759 1.504575488252 2.000000000000 2.778158480438
stirling: 2.678938532866 1.354117938504 0.999999999306 0.892979510955 0.902745292336 0.999999999306 1.190639347940 1.504575487227 1.999999998611 2.778158478527

sub Γ(\z) {    constant g = 9;    z < .5 ?? pi/ sin(pi * z) / Γ(1 - z) !!    sqrt(2*pi) *    (z + g - 1/2)**(z - 1/2) *    exp(-(z + g - 1/2)) *    [+] <        1.000000000000000174663     5716.400188274341379136   -14815.30426768413909044    14291.49277657478554025    -6348.160217641458813289     1301.608286058321874105     -108.1767053514369634679        2.605696505611755827729       -0.7423452510201416151527e-2        0.5384136432509564062961e-7       -0.4023533141268236372067e-8    > Z* 1, |map 1/(z + *), 0..*} say Γ($_) for 1/3, 2/3 ... 10/3; Output: 2.67893853470775 1.3541179394264 1 0.892979511569248 0.902745292950934 1 1.190639348759 1.50457548825155 2 2.77815848043766 ## Phix Translation of: C sequence c = repeat(0,12) function gamma(atom z) atom accm = c[1] if accm=0 then accm = sqrt(2*PI) c[1] = accm atom k1_factrl = 1 -- (k - 1)!*(-1)^k with 0!==1 for k=2 to 12 do c[k] = exp(13-k)*power(13-k,k-1.5)/k1_factrl k1_factrl *= -(k-1) end for end if for k=2 to 12 do accm += c[k]/(z+k-1) end for accm *= exp(-(z+12))*power(z+12,z+0.5) -- Gamma(z+1) return accm/zend function procedure sq(atom x, atom mul)atom p = x*mul printf(1,"%18.16g,%18.15g\n",{x,p*p})end procedure procedure si(atom x) printf(1,"%18.15f\n",{x})end procedure sq(gamma(-3/2),3/4)sq(gamma(-1/2),-1/2)sq(gamma(1/2),1)si(gamma(1))sq(gamma(3/2),2)si(gamma(2))sq(gamma(5/2),4/3)si(gamma(3))sq(gamma(7/2),8/15)si(gamma(4)) Output:  2.363271801207354, 3.14159265358979 -3.544907701811032, 3.14159265358979 1.772453850905515, 3.14159265358979 1.000000000000001 0.8862269254527643, 3.14159265358984 1.000000000000010 1.329340388179146, 3.14159265358984 2.000000000000024 3.323350970447942, 3.14159265358998 6.000000000000175  ## PicoLisp Translation of: Ada (scl 28) (de *A ~(flip (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 ) ) ) (de gamma (X) (let (Y (- X 1.0) Sum (car *A)) (for A (cdr *A) (setq Sum (+ A (*/ Sum Y 1.0))) ) (*/ 1.0 1.0 Sum) ) ) Output: : (for I (range 1 10) (prinl (round (gamma (*/ I 1.0 3)) 14)) ) 2.67893853470775 1.35411793942640 1.00000000000000 0.89297951156925 0.90274529295093 1.00000000000000 1.19063934875900 1.50457548825154 1.99999999999397 2.77815847933858 ## PL/I /* From Rosetta Fortran */test: procedure options (main); declare i fixed binary; on underflow ; put skip list ('Lanczos', 'Builtin' ); do i = 1 to 10; put skip list (lanczos_gamma(i/3.0q0), gamma(i/3.0q0) ); end; lanczos_gamma: procedure (a) returns (float (18)) recursive; declare a float (18); declare pi float (18) value (3.14159265358979324E0); declare cg fixed binary initial ( 7 ); /* these precomputed values are taken by the sample code in Wikipedia, */ /* and the sample itself takes them from the GNU Scientific Library */ declare p(0:8) float (18) static initial ( 0.99999999999980993e0, 676.5203681218851e0, -1259.1392167224028e0, 771.32342877765313e0, -176.61502916214059e0, 12.507343278686905e0, -0.13857109526572012e0, 9.9843695780195716e-6, 1.5056327351493116e-7 ); declare ( t, w, x ) float (18); declare i fixed binary; x = a; if x < 0.5 then return ( pi / ( sin(pi*x) * lanczos_gamma(1.0-x) ) ); else do; x = x - 1.0; t = p(0); do i = 1 to cg+2; t = t + p(i)/(x+i); end; w = x + float(cg) + 0.5; return ( sqrt(2*pi) * w**(x+0.5) * exp(-w) * t ); end; end lanczos_gamma; end test; Output: Lanczos Builtin 2.67893853470774706E+0000 2.678938534707747630E+0000 1.35411793942640071E+0000 1.354117939426400420E+0000 1.00000000000000021E+0000 1.000000000000000000E+0000 8.92979511569249470E-0001 8.929795115692492110E-0001 9.02745292950933961E-0001 9.027452929509336110E-0001 1.00000000000000048E+0000 1.000000000000000000E+0000 1.19063934875899964E+0000 1.190639348758998950E+0000 1.50457548825155704E+0000 1.504575488251556020E+0000 2.00000000000000154E+0000 2.000000000000000000E+0000 2.77815848043766660E+0000 2.778158480437664210E+0000  ## PowerShell I would download the Math.NET Numerics dll(s). Documentation and download at: http://cyber-defense.sans.org/blog/2015/06/27/powershell-for-math-net-numerics/comment-page-1/  Add-Type -Path "C:\Program Files (x86)\Math\MathNet.Numerics.3.12.0\lib\net40\MathNet.Numerics.dll" 1..20 | ForEach-Object {[MathNet.Numerics.SpecialFunctions]::Gamma($_ / 10)} 
Output:
9.51350769866874
4.5908437119988
2.99156898768759
2.21815954375769
1.77245385090552
1.48919224881282
1.29805533264756
1.1642297137253
1.06862870211932
1
0.951350769866874
0.918168742399759
0.897470696306277
0.887263817503075
0.88622692545276
0.89351534928769
0.908638732853289
0.931383770980245
0.961765831907388
1


## PureBasic

Below is PureBasic code for:

• Complete Gamma function
• Natural Logarithm of the Complete Gamma function
• Factorial function
Procedure.d Gamma(x.d) ; AKJ  01-May-10; Complete Gamma function for x>0 and x<2 (approx); Extended outside this range via: Gamma(x+1) = x*Gamma(x); Based on http://rosettacode.org/wiki/Gamma_function  [Ada]Protected Dim A.d(28)A(0) = 1.0A(1) = 0.5772156649015328606A(2) =-0.6558780715202538811A(3) =-0.0420026350340952355A(4) = 0.1665386113822914895A(5) =-0.0421977345555443368 ; was ...33675A(6) =-0.0096219715278769736A(7) = 0.0072189432466630995A(8) =-0.0011651675918590651A(9) =-0.0002152416741149510A(10) = 0.0001280502823881162A(11) =-0.0000201348547807882A(12) =-0.0000012504934821427A(13) = 0.0000011330272319817A(14) =-0.0000002056338416978A(15) = 0.0000000061160951045A(16) = 0.0000000050020076445A(17) =-0.0000000011812745705A(18) = 0.0000000001043426712A(19) = 0.0000000000077822634A(20) =-0.0000000000036968056A(21) = 0.0000000000005100370A(22) =-0.0000000000000205833A(23) =-0.0000000000000053481A(24) = 0.0000000000000012268A(25) =-0.0000000000000001181A(26) = 0.0000000000000000012A(27) = 0.0000000000000000014A(28) =-0.0000000000000000002;A(29) = 0.00000000000000000002Protected y.d, Prod.d, Sum.d, NIf x<=0.0: ProcedureReturn 0.0: EndIf ; Errory = x-1.0: Prod = 1.0While y>=1.0  Prod*y: y-1.0 ; Recurse using Gamma(x+1) = x*Gamma(x)WendSum= A(28)For N = 27 To 0 Step -1  Sum*y+A(N)Next NIf Sum=0.0: ProcedureReturn Infinity(): EndIfProcedureReturn Prod / SumEndProcedure Procedure.d GammLn(x.d) ; AKJ  01-May-10; Returns Ln(Gamma()) or 0 on error; Based on Numerical Recipes gamma.hProtected j, tmp.d, y.d, ser.dProtected Dim cof.d(13)cof(0) = 57.1562356658629235cof(1) = -59.5979603554754912cof(2) = 14.1360979747417471cof(3) = -0.491913816097620199cof(4) = 0.339946499848118887e-4cof(5) = 0.465236289270485756e-4cof(6) = -0.983744753048795646e-4cof(7) = 0.158088703224912494e-3cof(8) = -0.210264441724104883e-3cof(9) = 0.217439618115212643e-3cof(10) = -0.164318106536763890e-3cof(11) = 0.844182239838527433e-4cof(12) = -0.261908384015814087e-4cof(13) = 0.368991826595316234e-5If x<=0: ProcedureReturn 0: EndIf ; Bad argumenty = xtmp = x+5.2421875tmp = (x+0.5)*Log(tmp)-tmpser = 0.999999999999997092For j=0 To 13  y + 1: ser + cof(j)/yNext jProcedureReturn tmp+Log(2.5066282746310005*ser/x)EndProcedure Procedure Factorial(x) ; AKJ  01-May-10  ProcedureReturn Gamma(x+1)EndProcedure
Examples
Debug "Gamma()"For i = 12 To 15  Debug StrD(i/3.0, 3)+"   "+StrD(Gamma(i/3.0))Next iDebug ""Debug "Ln(Gamma(5.0)) = "+StrD(GammLn(5.0), 16) ; Ln(24)Debug ""Debug "Factorial 6 = "+StrD(Factorial(6), 0) ; 72
Output:
[Debug] Gamma():
[Debug] 4.000   6.0000000000
[Debug] 4.333   9.2605282681
[Debug] 4.667   14.7114047740
[Debug] 5.000   24.0000000000
[Debug]
[Debug] Ln(Gamma(5.0)) = 3.1780538303479458
[Debug]
[Debug] Factorial 6 = 720

## 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) 
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 loosely converted from the Octave code.

Translation of: Octave
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)))   }} # Checksz <- (1:10)/3all.equal(gamma(z), stirling(z))             # Mean relative difference: 0.07181942all.equal(gamma(z), nemes(z))                # Mean relative difference: 0.003460549all.equal(as.complex(gamma(z)), lanczos(z))  # TRUEall.equal(gamma(z), spouge(z))               # TRUEdata.frame(z=z, stirling=stirling(z), nemes=nemes(z), lanczos=lanczos(z), spouge=spouge(z), builtin=gamma(z))
Output:
          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


## Racket

#lang racket(define (gamma number)  (if (> 1/2 number)      (/ pi (* (sin (* pi number))               (gamma (- 1.0 number))))      (let ((n (sub1 number))            (c '(0.99999999999980993 676.5203681218851 -1259.1392167224028                 771.32342877765313 -176.61502916214059 12.507343278686905	             -0.13857109526572012 9.9843695780195716e-6 1.5056327351493116e-7)))        (* (sqrt (* pi 2))           (expt (+ n 7 0.5) (+ n 0.5))           (exp (- (+ n 7 0.5)))           (+ (car c)              (apply +                (for/list ((i (in-range (length (cdr c)))) (x (in-list (cdr c))))                  (/ x (+ 1 n i))))))))) (map gamma '(0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0));->;'(9.513507698668736;  4.590843711998802;  2.9915689876875904;  2.218159543757687;  1.7724538509055159;  1.489192248812818;  1.2980553326475577;  1.1642297137253037;  1.068628702119319;  1.0)

## REXX

This version uses a Taylor series 80-digits coefficients with much more accuracy.
As a result, the gamma value for   ½   is now   25   digits more accurate than the previous version
(which only used   20   digit coefficients).

Note:   The Taylor series isn't much good above values of   .

/*REXX program calculates GAMMA using the Taylor series coefficients; ≈80 decimal digits*/                            /*The GAMMA function symbol is the Greek capital letter:  Γ */numeric digits 90                                /*be able to handle extended precision.*/parse arg LO HI .                                /*allow specification of gamma arg/args*/                                                 /* [↓]  either show a range or a ···   */        do j=word(LO 1, 1)  to word(HI LO 9, 1)  /*              ··· single gamma value.*/        say 'gamma('j") ="  gamma(j)             /*compute gamma of J and display value.*/        end   /*j*/                              /* [↑]  default LO is one;  HI is nine.*/exit                                             /*stick a fork in it,  we're all done. *//*──────────────────────────────────────────────────────────────────────────────────────*/gamma:  procedure;   parse arg x;    xm=x-1;    sum=0                                /*coefficients thanks to: Arne Fransén & Staffan Wrigge.*/ #.1 =  1                       /* [↓]    #.2   is the   Euler-Mascheroni  constant.    */ #.2 =  0.57721566490153286060651209008240243104215933593992359880576723488486772677766467 #.3 = -0.65587807152025388107701951514539048127976638047858434729236244568387083835372210 #.4 = -0.04200263503409523552900393487542981871139450040110609352206581297618009687597599 #.5 =  0.16653861138229148950170079510210523571778150224717434057046890317899386605647425 #.6 = -0.04219773455554433674820830128918739130165268418982248637691887327545901118558900 #.7 = -0.00962197152787697356211492167234819897536294225211300210513886262731167351446074 #.8 =  0.00721894324666309954239501034044657270990480088023831800109478117362259497415854 #.9 = -0.00116516759185906511211397108401838866680933379538405744340750527562002584816653#.10 = -0.00021524167411495097281572996305364780647824192337833875035026748908563946371678#.11 =  0.00012805028238811618615319862632816432339489209969367721490054583804120355204347#.12 = -0.00002013485478078823865568939142102181838229483329797911526116267090822918618897#.13 = -0.00000125049348214267065734535947383309224232265562115395981534992315749121245561#.14 =  0.00000113302723198169588237412962033074494332400483862107565429550539546040842730#.15 = -0.00000020563384169776071034501541300205728365125790262933794534683172533245680371#.16 =  0.00000000611609510448141581786249868285534286727586571971232086732402927723507435#.17 =  0.00000000500200764446922293005566504805999130304461274249448171895337887737472132#.18 = -0.00000000118127457048702014458812656543650557773875950493258759096189263169643391#.19 =  0.00000000010434267116911005104915403323122501914007098231258121210871073927347588#.20 =  0.00000000000778226343990507125404993731136077722606808618139293881943550732692987#.21 = -0.00000000000369680561864220570818781587808576623657096345136099513648454655443000#.22 =  0.00000000000051003702874544759790154813228632318027268860697076321173501048565735#.23 = -0.00000000000002058326053566506783222429544855237419746091080810147188058196444349#.24 = -0.00000000000000534812253942301798237001731872793994898971547812068211168095493211#.25 =  0.00000000000000122677862823826079015889384662242242816545575045632136601135999606#.26 = -0.00000000000000011812593016974587695137645868422978312115572918048478798375081233#.27 =  0.00000000000000000118669225475160033257977724292867407108849407966482711074006109#.28 =  0.00000000000000000141238065531803178155580394756670903708635075033452562564122263#.29 = -0.00000000000000000022987456844353702065924785806336992602845059314190367014889830#.30 =  0.00000000000000000001714406321927337433383963370267257066812656062517433174649858#.31 =  0.00000000000000000000013373517304936931148647813951222680228750594717618947898583#.32 = -0.00000000000000000000020542335517666727893250253513557337960820379352387364127301#.33 =  0.00000000000000000000002736030048607999844831509904330982014865311695836363370165#.34 = -0.00000000000000000000000173235644591051663905742845156477979906974910879499841377#.35 = -0.00000000000000000000000002360619024499287287343450735427531007926413552145370486#.36 =  0.00000000000000000000000001864982941717294430718413161878666898945868429073668232#.37 = -0.00000000000000000000000000221809562420719720439971691362686037973177950067567580#.38 =  0.00000000000000000000000000012977819749479936688244144863305941656194998646391332#.39 =  0.00000000000000000000000000000118069747496652840622274541550997151855968463784158#.40 = -0.00000000000000000000000000000112458434927708809029365467426143951211941179558301#.41 =  0.00000000000000000000000000000012770851751408662039902066777511246477487720656005#.42 = -0.00000000000000000000000000000000739145116961514082346128933010855282371056899245#.43 =  0.00000000000000000000000000000000001134750257554215760954165259469306393008612196#.44 =  0.00000000000000000000000000000000004639134641058722029944804907952228463057968680#.45 = -0.00000000000000000000000000000000000534733681843919887507741819670989332090488591#.46 =  0.00000000000000000000000000000000000032079959236133526228612372790827943910901464#.47 = -0.00000000000000000000000000000000000000444582973655075688210159035212464363740144#.48 = -0.00000000000000000000000000000000000000131117451888198871290105849438992219023663#.49 =  0.00000000000000000000000000000000000000016470333525438138868182593279063941453996#.50 = -0.00000000000000000000000000000000000000001056233178503581218600561071538285049997#.51 =  0.00000000000000000000000000000000000000000026784429826430494783549630718908519485#.52 =  0.00000000000000000000000000000000000000000002424715494851782689673032938370921241#=52;                           do k=#  by -1  for #                                sum=sum*xm  +  #.k                                end   /*k*/return 1/sum
output   when using the input of:   0.5
gamma(0.5) = 1.77245385090551602729816748334114518279754945612238712821380509003635917689651032047826593


Note that:   Γ(½) = ${\displaystyle \pi }$ =

1.77245 38509 05516 02729 81674 83341 14518 27975 49456 12238 71282 13807 78985 29112 84591 03218 13749 50656 73854 46654 16226 82362 +

to 110 digits past the decimal point,   the vinculum (overbar) marks the   difference digit   from the computed value (by this REXX program) of   gamma(½).

## Ring

 decimals(3)gamma = 0.577coeff = -0.655quad = -0.042qui = 0.166set = -0.042 for i=1 to 10    see gammafunc(i / 3.0) + nlnext func recigamma z     return z + gamma * pow(z,2) + coeff * pow(z,3) + quad * pow(z,4) + qui * pow(z,5) + set * pow(z,6) func gammafunc z     if z = 1 return 1     but fabs(z) <= 0.5 return 1 / recigamma(z)     else return (z - 1) * gammafunc(z-1) ok 

## RLaB

RLaB through GSL has the following functions related to the Gamma function, namely, Gamma, GammaRegularizedC, LogGamma, RecGamma, and Pochhammer, where

${\displaystyle Gamma(a)=\Gamma (a)}$, the Gamma function;
${\displaystyle Gamma(a,x)={\frac {1}{\Gamma (a)}}\int _{x}^{\infty }dt\,t^{a-1}\,exp(-t)}$, the regularized Gamma function which is also known as the normalized incomplete Gamma function;
${\displaystyle GammaRegularizedC(a,x)={\frac {1}{\Gamma (a)}}\int _{0}^{x}dt\,t^{a-1}\,exp(-t)}$, which the GSL calls the complementary normalized Gamma function;
${\displaystyle LogGamma(a)=\ln \Gamma (a)}$;
${\displaystyle RecGamma(a)={\frac {1}{\Gamma (a)}}}$;
${\displaystyle Pochhammer(a,x)={\frac {\Gamma (a+x)}{\Gamma (x)}}}$.

## Ruby

#### Taylor series

$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 * y + an}end (1..10).each {|i| puts format("%.14e", gamma(i/3.0))}
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

#### Built in

(1..10).each{|i| puts Math.gamma(i/3.0)}
Output:
2.678938534707748
1.3541179394264005
1.0
0.8929795115692493
0.9027452929509336
1.0
1.190639348758999
1.5045754882515558
2.0
2.7781584804376647


## Scala

import java.util.Locale._ object Gamma {  def stGamma(x:Double):Double=math.sqrt(2*math.Pi/x)*math.pow((x/math.E), x)   def laGamma(x:Double):Double={    val p=Seq(676.5203681218851, -1259.1392167224028, 771.32342877765313,              -176.61502916214059, 12.507343278686905, -0.13857109526572012,                9.9843695780195716e-6, 1.5056327351493116e-7)     if(x < 0.5) {      math.Pi/(math.sin(math.Pi*x)*laGamma(1-x))    } else {      val x2=x-1      val t=x2+7+0.5      val a=p.zipWithIndex.foldLeft(0.99999999999980993)((r,v) => r+v._1/(x2+v._2+1))      math.sqrt(2*math.Pi)*math.pow(t, x2+0.5)*math.exp(-t)*a    }  }   def main(args: Array[String]): Unit = {    println("Gamma    Stirling             Lanczos")    for(x <- 0.1 to 2.0 by 0.1)      println("%.1f  ->  %.16f   %.16f".formatLocal(ENGLISH, x, stGamma(x), laGamma(x)))  }}
Output:
Gamma    Stirling             Lanczos
0.1  ->  5.6971871489771690   9.5135076986687340
0.2  ->  3.3259984240223925   4.5908437119988030
0.3  ->  2.3625300362696198   2.9915689876875904
0.4  ->  1.8414763359362354   2.2181595437576870
0.5  ->  1.5203469010662807   1.7724538509055159
0.6  ->  1.3071588574483560   1.4891922488128180
0.7  ->  1.1590532921139200   1.2980553326475577
0.8  ->  1.0533709684256085   1.1642297137253035
0.9  ->  0.9770615078776956   1.0686287021193193
1.0  ->  0.9221370088957892   1.0000000000000002
1.1  ->  0.8834899531687038   0.9513507698668728
1.2  ->  0.8577553353965909   0.9181687423997607
1.3  ->  0.8426782594483921   0.8974706963062777
1.4  ->  0.8367445486370817   0.8872638175030760
1.5  ->  0.8389565525264964   0.8862269254527583
1.6  ->  0.8486932421525738   0.8935153492876904
1.7  ->  0.8656214717938840   0.9086387328532912
1.8  ->  0.8896396352879945   0.9313837709802430
1.9  ->  0.9208427218942294   0.9617658319073875
2.0  ->  0.9595021757444918   1.0000000000000010

## Scheme

Translation of: Scala
for Lanczos and Stirling
Translation of: Ruby
for Taylor
 (import (scheme base)        (scheme inexact)        (scheme write)) (define PI 3.14159265358979323846264338327950)  (define e 2.7182818284590452353602875) (define gamma-lanczos  (let ((p '(676.5203681218851 -1259.1392167224028 771.32342877765313              -176.61502916214059 12.507343278686905 -0.13857109526572012             9.9843695780195716e-6 1.5056327351493116e-7)))    (lambda (x)      (if (< x 0.5)        (/ PI (* (sin (* PI x)) (gamma-lanczos (- 1 x))))        (let* ((x2 (- x 1))               (t (+ x2 7 0.5))               (a (do ((ps p (cdr ps))                       (idx 0 (+ 1 idx))                       (res 0.99999999999980993 (+ res                                                    (/ (car ps)                                                      (+ x2 idx 1)))))                    ((null? ps) res))))          (* (sqrt (* 2 PI)) (expt t (+ x2 0.5)) (exp (- t)) a)))))) (define (gamma-stirling x)  (* (sqrt (* 2 (/ PI x))) (expt (/ x e) x))) (define gamma-taylor  (let ((a (reverse             '(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))))    (lambda (x)      (let ((y (- x 1)))        (do ((as a (cdr as))             (res 0 (+ (car as) (* res y))))          ((null? as) (/ 1 res))))))) (do ((i 0.1 (+ i 0.1)))  ((> i 2.01) )  (display (string-append "Gamma ("                          (number->string i)                          "): "                          "\n --- Lanczos : "                          (number->string (gamma-lanczos i))                          "\n --- Stirling: "                          (number->string (gamma-stirling i))                          "\n --- Taylor  : "                          (number->string (gamma-taylor i))                          "\n"))) 
Output:
Gamma (0.1):
--- Lanczos : 9.513507698668736
--- Stirling: 5.69718714897717
--- Taylor  : 9.513507698668734
Gamma (0.2):
--- Lanczos : 4.590843711998803
--- Stirling: 3.3259984240223925
--- Taylor  : 4.5908437119988035
Gamma (0.30000000000000004):
--- Lanczos : 2.9915689876875904
--- Stirling: 2.3625300362696198
--- Taylor  : 2.991568987687591
Gamma (0.4):
--- Lanczos : 2.218159543757687
--- Stirling: 1.8414763359362354
--- Taylor  : 2.2181595437576886
Gamma (0.5):
--- Lanczos : 1.7724538509055159
--- Stirling: 1.5203469010662807
--- Taylor  : 1.772453850905516
Gamma (0.6):
--- Lanczos : 1.489192248812818
--- Stirling: 1.307158857448356
--- Taylor  : 1.489192248812817
Gamma (0.7):
--- Lanczos : 1.2980553326475577
--- Stirling: 1.15905329211392
--- Taylor  : 1.298055332647558
Gamma (0.7999999999999999):
--- Lanczos : 1.1642297137253035
--- Stirling: 1.0533709684256085
--- Taylor  : 1.1642297137253033
Gamma (0.8999999999999999):
--- Lanczos : 1.0686287021193193
--- Stirling: 0.9770615078776956
--- Taylor  : 1.0686287021193195
Gamma (0.9999999999999999):
--- Lanczos : 1.0000000000000002
--- Stirling: 0.9221370088957892
--- Taylor  : 1.0000000000000002
Gamma (1.0999999999999999):
--- Lanczos : 0.9513507698668728
--- Stirling: 0.8834899531687039
--- Taylor  : 0.9513507698668733
Gamma (1.2):
--- Lanczos : 0.9181687423997607
--- Stirling: 0.8577553353965909
--- Taylor  : 0.9181687423997608
Gamma (1.3):
--- Lanczos : 0.8974706963062777
--- Stirling: 0.842678259448392
--- Taylor  : 0.8974706963062773
Gamma (1.4000000000000001):
--- Lanczos : 0.8872638175030759
--- Stirling: 0.8367445486370818
--- Taylor  : 0.8872638175030753
Gamma (1.5000000000000002):
--- Lanczos : 0.8862269254527583
--- Stirling: 0.8389565525264964
--- Taylor  : 0.886226925452758
Gamma (1.6000000000000003):
--- Lanczos : 0.8935153492876904
--- Stirling: 0.8486932421525738
--- Taylor  : 0.8935153492876904
Gamma (1.7000000000000004):
--- Lanczos : 0.9086387328532912
--- Stirling: 0.865621471793884
--- Taylor  : 0.9086387328532904
Gamma (1.8000000000000005):
--- Lanczos : 0.931383770980243
--- Stirling: 0.8896396352879945
--- Taylor  : 0.9313837709802427
Gamma (1.9000000000000006):
--- Lanczos : 0.9617658319073875
--- Stirling: 0.9208427218942294
--- Taylor  : 0.9617658319073876
Gamma (2.0000000000000004):
--- Lanczos : 1.000000000000001
--- Stirling: 0.9595021757444918
--- Taylor  : 1.0000000000000002


## Scilab

function x=gammal(z)  // Lanczos'    lz=[  1.000000000190015 ..          76.18009172947146  -86.50532032941677      24.01409824083091    ..         -1.231739572450155   1.208650973866179e-3  -5.395239384953129e-6 ]    if z < 0.5 then         x=%pi/sin(%pi*z)-gammal(1-z)    else        z=z-1.0        b=z+5.5        a=lz(1)        for i=1:6            a=a+(lz(i+1)/(z+i))        end        x=exp((log(sqrt(2*%pi))+log(a)-b)+log(b)*(z+0.5))    endendfunction printf("%4s %-9s %-9s\n","x","gamma(x)","gammal(x)")for i=1:30    x=i/10    printf("%4.1f %9f %9f\n",x,gamma(x),gammal(x))end
Output:
   x gamma(x)  gammal(x)
0.1  9.097779  9.097779
0.2  4.180567  4.180567
0.3  2.585167  2.585167
0.4  1.814074  1.814074
0.5  1.772454  1.772454
0.6  1.489192  1.489192
0.7  1.298055  1.298055
0.8  1.164230  1.164230
0.9  1.068629  1.068629
1.0  1.000000  1.000000
1.1  0.951351  0.951351
1.2  0.918169  0.918169
1.3  0.897471  0.897471
1.4  0.887264  0.887264
1.5  0.886227  0.886227
1.6  0.893515  0.893515
1.7  0.908639  0.908639
1.8  0.931384  0.931384
1.9  0.961766  0.961766
2.0  1.000000  1.000000
2.1  1.046486  1.046486
2.2  1.101802  1.101802
2.3  1.166712  1.166712
2.4  1.242169  1.242169
2.5  1.329340  1.329340
2.6  1.429625  1.429625
2.7  1.544686  1.544686
2.8  1.676491  1.676491
2.9  1.827355  1.827355
3.0  2.000000  2.000000


$include "seed7_05.s7i"; include "float.s7i"; const func float: gamma (in float: X) is func result var float: result is 0.0; local const array float: A is [] ( 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); var float: Y is 0.0; var float: Sum is A[maxIdx(A)]; var integer: N is 0; begin Y := X - 1.0; for N range pred(maxIdx(A)) downto minIdx(A) do Sum := Sum * Y + A[N]; end for; result := 1.0 / Sum; end func; const proc: main is func local var integer: I is 0; begin for I range 1 to 10 do writeln((gamma(flt(I) / 3.0)) digits 15); end for; end func; Output: 2.678937911987305 1.354117870330811 1.000000000000000 0.892979443073273 0.902745306491852 1.000000000000000 1.190639257431030 1.504575252532959 1.999999523162842 2.778157949447632 ## Sidef Translation of: Ruby var 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 ] func gamma(x) { var y = (x - 1) 1 / a.reverse.reduce {|sum, an| sum*y + an}} for i in 1..10 { say ("%.14e" % gamma(i/3))} 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  Lanczos approximation: func gamma(z) { var epsilon = 0.0000001 func withinepsilon(x) { abs(x - abs(x)) <= epsilon } var p = [ 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7, ] var result = 0 const pi = Num.pi if (z.re < 0.5) { result = (pi / (sin(pi * z) * gamma(1 - z))) } else { z -= 1 var x = 0.99999999999980993 p.each_kv { |i, v| x += v/(z + i + 1) } var t = (z + p.len - 0.5) result = (sqrt(pi*2) * t**(z+0.5) * exp(-t) * x) } withinepsilon(result.im) ? result.re : result} for i in 1..10 { say ("%.14e" % gamma(i/3))} Output: 2.67893853470774e+00 1.35411793942640e+00 1.00000000000000e+00 8.92979511569252e-01 9.02745292950931e-01 1.00000000000000e+00 1.19063934875900e+00 1.50457548825155e+00 2.00000000000000e+00 2.77815848043767e+00  A simpler implementation: define ℯ = Num.edefine τ = Num.tau func Γ(t) { t < 20 ? (__FUNC__(t + 1) / t) : (sqrt(τ*t) * pow(t/ℯ + 1/(12*ℯ*t), t) / t)} for i in (1..10) { say ("%.14e" % Γ(i/3))} Output: 2.67893831294932e+00 1.35411783267848e+00 9.99999913007168e-01 8.92979437649773e-01 9.02745221785653e-01 9.99999913007168e-01 1.19063925019970e+00 1.50457536964275e+00 1.99999982601434e+00 2.77815825046596e+00  ## Stata This implementation uses the Taylor expansion of 1/gamma(1+x). The coefficients were computed with Maxima (see the Maxima implementation above). The results are compared to Mata's gamma function for each real between 1/100 and 100, by steps of 1/100. mata_gamma_coef = 1.0,5.7721566490153286061e-1,-6.5587807152025388108e-1,-4.2002635034095235529e-2,1.665386113822914895e-1,-4.2197734555544336748e-2,-9.6219715278769735621e-3,7.2189432466630995424e-3,-1.1651675918590651121e-3,-2.1524167411495097282e-4,1.2805028238811618615e-4,-2.0134854780788238656e-5,-1.2504934821426706573e-6,1.1330272319816958824e-6,-2.0563384169776071035e-7,6.1160951044814158179e-9,5.0020076444692229301e-9,-1.1812745704870201446e-9,1.0434267116911005105e-10,7.782263439905071254e-12,-3.6968056186422057082e-12,5.100370287454475979e-13,-2.0583260535665067832e-14,-5.3481225394230179824e-15 function gamma_(x_) { external _gamma_coef x = x_ for (y=1; x>2;) y = --x*y z = _gamma_coef[24] x-- for (i=23; i>=1; i--) z = z*x+_gamma_coef[i] return(y/z)} function map(f,a) { n = rows(a) p = cols(a) b = J(n,p,.) for (i=1; i<=n; i++) { for (j=1; j<=p; j++) { b[i,j] = (*f)(a[i,j]) } } return(b)} x=(1::10000)/100u=map(&gamma(),x)v=map(&gamma_(),x)max(abs((u-v):/u))end Output 1.27664e-13 ## Tcl Works with: Tcl version 8.5 Library: Tcllib (Package: math) Library: Tcllib (Package: math::calculus) package require mathpackage require math::calculus # gamma(1) and gamma(1.5) set f 1.0set 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]    }}
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

## TI-83 BASIC

There is an hidden Gamma function in TI-83. Factorial (!) is implemented in increments of 0.5 :

.5! -> .8862269255


As far as Gamma(n)=(n-1)! , we have everything needed.

### Stirling's approximation

for(I,1,10)I/2→XX^(X-1/2)e^(-X)√(2π)→YDisp X,(X-1)!,YPauseEnd
Output:

The output display for x=0.5 to 5 by 0.5 : x, (x-1)!, Y(x) . Y(x) is Stirling's approximation of Gamma.

        0.5
1.772453851
1.520346901
1
1
.9221370089
1.5
.8862269255
.8389565525
2
1
.9595021757
2.5
1.329340388
1.285978615
3
2
1.945403197
3.5
3.32335097
3.245363748
4
6
5.876543783
4.5
11.6317284
11.41865156
5
24
23.60383359


### Lanczos' approximation

for(I,1,10)I/2→Xe^(ln((1.0+76.18009173/(X+1)-86.50532033/(X+2)+24.01409824/(X+3)-1.231739572/(X+4)+1.208650974E-3/(X+5)-5.395239385E-6/(X+6))√(2π)/X)+(X+.5)ln(X+5.5)-X-5.5)->YDisp X,(X-1)!,YPauseEnd
Output:

The output display for x=0.5 to 5 by 0.5 : x, (x-1)!, Y(x) . Y(x) is Lanczos's approximation of Gamma.

        0.5
1.772453851
1.772453851
1
1
1
1.5
.8862269255
.8862269254
2
1
1
2.5
1.329340388
1.329340388
3
2
2
3.5
3.32335097
3.32335097
4
6
6
4.5
11.6317284
11.6317284
5
24
24


## Visual FoxPro

Translation of BBC Basic but with OOP extensions. Also some ideas from Numerical Methods (Press et al).

 LOCAL i As Integer, x As Double, o As lanczosCLOSE DATABASES ALLCLEARCREATE CURSOR results (ZVal B(1), GamVal B(15))INDEX ON zval TAG ZVal COLLATE "Machine"SET ORDER TO 0o = CREATEOBJECT("lanczos")FOR i = 1 TO 20x = i/10    INSERT INTO results VALUES (x, o.Gamma(x))ENDFORUPDATE results SET GamVal = ROUND(GamVal, 0) WHERE ZVal = INT(ZVal)*!* This just creates the output text - it is not part of the algorithmDO cursor2txt.prg WITH "results", .T. DEFINE CLASS lanczos As Session#DEFINE FPF 5.5#DEFINE HALF 0.5#DEFINE PY PI()DIMENSION LanCoeff[7]nSize = 0 PROCEDURE InitDODEFAULT()WITH THIS    .LanCoeff[1] = 1.000000000190015    .LanCoeff[2] = 76.18009172947146    .LanCoeff[3] = -86.50532032941677    .LanCoeff[4] = 24.01409824083091    .LanCoeff[5] = -1.231739572450155    .LanCoeff[6] = 0.0012086509738662    .LanCoeff[7] = -0.000005395239385    .nSize = ALEN(.LanCoeff)ENDWITH ENDPROC FUNCTION Gamma(z)RETURN EXP(THIS.LogGamma(z))ENDFUNC FUNCTION LogGamma(z)LOCAL a As Double, b As Double, i As Integer, j As Integer, lg As DoubleIF z < 0.5    lg = LOG(PY/SIN(PY*z)) - THIS.LogGamma(1 - z)ELSE    WITH THIS		z = z - 1 	b = z + FPF	a = .LanCoeff[1]	FOR i = 2 TO .nSize	    j = i - 1	    a = a + .LanCoeff[i]/(z + j)	ENDFOR	lg = (LOG(SQRT(2*PY)) + LOG(a) - b) + LOG(b)*(z + HALF)    ENDWITH	ENDIFRETURN lgENDFUNC	 ENDDEFINE 
Output:
zval	gamval
0.1	9.513507698669704
0.2	4.590843712000122
0.3	2.991568987689402
0.4	2.218159543760185
0.5	1.772453850902053
0.6	1.489192248811141
0.7	1.298055332646772
0.8	1.164229713724969
0.9	1.068628702119210
1.0	1.000000000000000
1.1	0.951350769866919
1.2	0.918168742399821
1.3	0.897470696306335
1.4	0.887263817503125
1.5	0.886226925452796
1.6	0.893515349287718
1.7	0.908638732853309
1.8	0.931383770980253
1.9	0.961765831907391
2.0	1.000000000000000


## zkl

Translation of: D
but without a built in gamma function.
fcn taylorGamma(x){   var table = T(     0x1p+0,                    0x1.2788cfc6fb618f4cp-1,    -0x1.4fcf4026afa2dcecp-1,  -0x1.5815e8fa27047c8cp-5,     0x1.5512320b43fbe5dep-3,  -0x1.59af103c340927bep-5,    -0x1.3b4af28483e214e4p-7,   0x1.d919c527f60b195ap-8,    -0x1.317112ce3a2a7bd2p-10, -0x1.c364fe6f1563ce9cp-13,     0x1.0c8a78cd9f9d1a78p-13, -0x1.51ce8af47eabdfdcp-16,    -0x1.4fad41fc34fbb2p-20,    0x1.302509dbc0de2c82p-20,    -0x1.b9986666c225d1d4p-23,  0x1.a44b7ba22d628acap-28,     0x1.57bc3fc384333fb2p-28, -0x1.44b4cedca388f7c6p-30,     0x1.cae7675c18606c6p-34,   0x1.11d065bfaf06745ap-37,    -0x1.0423bac8ca3faaa4p-38,  0x1.1f20151323cd0392p-41,    -0x1.72cb88ea5ae6e778p-46, -0x1.815f72a05f16f348p-48,     0x1.6198491a83bccbep-50,  -0x1.10613dde57a88bd6p-53,     0x1.5e3fee81de0e9c84p-60,  0x1.a0dc770fb8a499b6p-60,    -0x1.0f635344a29e9f8ep-62,  0x1.43d79a4b90ce8044p-66).reverse();     y  := x.toFloat() - 1.0;    sm := table[1,*].reduce('wrap(sm,an){ sm * y + an },table[0]);     return(1.0 / sm);}
fcn lanczosGamma(z) { z = z.toFloat();    // Coefficients used by the GNU Scientific Library.    // http://en.wikipedia.org/wiki/Lanczos_approximation    const g = 7, PI = (0.0).pi;    exp := (0.0).e.pow;    var table = T(             0.99999_99999_99809_93,           676.52036_81218_851,         -1259.13921_67224_028,           771.32342_87776_5313,          -176.61502_91621_4059,            12.50734_32786_86905,            -0.13857_10952_65720_12,             9.98436_95780_19571_6e-6,             1.50563_27351_49311_6e-7);     // Reflection formula.    if (z < 0.5) {        return(PI / ((PI * z).sin() * lanczosGamma(1.0 - z)));    } else {        z -= 1;        x := table[0];        foreach i in ([1 .. g + 1]){            x += table[i] / (z + i); }        t := z + g + 0.5;        return((2.0 * PI).sqrt() * t.pow(z + 0.5) * exp(-t) * x);    }}
Output:
foreach i in ([1.0 .. 10]) {
x := i / 3.0;
println("%f: %20.19e %20.19e %e".fmt( x,
a:=taylorGamma(x), b:=lanczosGamma(x),(a-b).abs()));
}

0.333333: 2.6789385347077483424e+00 2.6789385347077474542e+00 8.881784e-16
0.666667: 1.3541179394264004632e+00 1.3541179394264002411e+00 2.220446e-16
1.000000: 1.0000000000000000000e+00 1.0000000000000002220e+00 2.220446e-16
1.333333: 8.9297951156924926241e-01 8.9297951156924970650e-01 4.440892e-16
1.666667: 9.0274529295093364212e-01 9.0274529295093353110e-01 1.110223e-16
2.000000: 1.0000000000000000000e+00 1.0000000000000006661e+00 6.661338e-16
2.333333: 1.1906393487589990166e+00 1.1906393487589996827e+00 6.661338e-16
2.666667: 1.5045754882515545159e+00 1.5045754882515582906e+00 3.774758e-15
3.000000: 1.9999999999992210675e+00 2.0000000000000017764e+00 7.807088e-13
3.333333: 2.7781584802531797962e+00 2.7781584804376668885e+00 1.844871e-10
`