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:
$jq -M -r -n -f Gamma_function_Stirling.jq i: gamma lanczos tgamma0.3333333333333333: 2.6789385347077483 : 2.6789385347077483 : 2.6789385347077480.6666666666666666: 1.3541179394264005 : 1.3541179394263998 : 1.35411793942640051 : 1 : 0.9999999999999998 : 11.3333333333333333: 0.8929795115692493 : 0.8929795115692494 : 0.89297951156924931.6666666666666667: 0.9027452929509336 : 0.9027452929509342 : 0.90274529295093362 : 1 : 1.0000000000000002 : 12.3333333333333335: 1.190639348758999 : 1.1906393487589995 : 1.1906393487589992.6666666666666665: 1.5045754882515399 : 1.5045754882515576 : 1.50457548825155583 : 1.9999999999939684 : 2.0000000000000013 : 23.3333333333333335: 2.778158479338573 : 2.778158480437665 : 2.7781584804376647 Julia Works with: Julia version 0.6 Built-in function: @show gamma(1) By adaptive Gauss-Kronrod integration: using QuadGKgammaquad(t::Float64) = first(quadgk(x -> x ^ (t - 1) * exp(-x), zero(t), Inf, reltol = 100eps(t)))@show gammaquad(1.0) Output: gamma(1) = 1.0 gammaquad(1.0) = 0.9999999999999999 Works with: Julia version 1.0 Library function: using SpecialFunctionsgamma(1/2) - sqrt(pi) Output: 2.220446049250313e-16 Kotlin // version 1.0.6 fun gammaStirling(x: Double): Double = Math.sqrt(2.0 * Math.PI / x) * Math.pow(x / Math.E, x) fun gammaLanczos(x: Double): Double { var xx = x val p = doubleArrayOf( 0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 ) val g = 7 if (xx < 0.5) return Math.PI / (Math.sin(Math.PI * xx) * gammaLanczos(1.0 - xx)) xx-- var a = p[0] val t = xx + g + 0.5 for (i in 1 until p.size) a += p[i] / (xx + i) return Math.sqrt(2.0 * Math.PI) * Math.pow(t, xx + 0.5) * Math.exp(-t) * a} fun main(args: Array<String>) { println(" x\tStirling\t\tLanczos\n") for (i in 1 .. 20) { val d = i / 10.0 print("%4.2f\t".format(d)) print("%17.15f\t".format(gammaStirling(d))) println("%17.15f".format(gammaLanczos(d))) }} Output:  x Stirling Lanczos 0.10 5.697187148977170 9.513507698668736 0.20 3.325998424022393 4.590843711998803 0.30 2.362530036269620 2.991568987687590 0.40 1.841476335936235 2.218159543757687 0.50 1.520346901066281 1.772453850905516 0.60 1.307158857448356 1.489192248812818 0.70 1.159053292113920 1.298055332647558 0.80 1.053370968425609 1.164229713725304 0.90 0.977061507877695 1.068628702119319 1.00 0.922137008895789 1.000000000000000 1.10 0.883489953168704 0.951350769866874 1.20 0.857755335396591 0.918168742399761 1.30 0.842678259448392 0.897470696306278 1.40 0.836744548637082 0.887263817503076 1.50 0.838956552526496 0.886226925452759 1.60 0.848693242152574 0.893515349287691 1.70 0.865621471793884 0.908638732853292 1.80 0.889639635287995 0.931383770980243 1.90 0.920842721894229 0.961765831907388 2.00 0.959502175744492 1.000000000000000  Limbo Translation of: Go A fairly straightforward port of the Go code. (It could almost have been done with sed). A few small differences are in the use of a tuple as a return value for the builtin gamma function, and we import a few functions from the math library so that we don't have to qualify them. implement Lanczos7; include "sys.m"; sys: Sys;include "draw.m";include "math.m"; math: Math; lgamma, exp, pow, sqrt: import math; Lanczos7: module { init: fn(nil: ref Draw->Context, nil: list of string);}; init(nil: ref Draw->Context, nil: list of string){ sys = load Sys Sys->PATH; math = load Math Math->PATH; # We ignore some floating point exceptions: math->FPcontrol(0, Math->OVFL|Math->UNFL); ns : list of real = -0.5 :: 0.1 :: 0.5 :: 1.0 :: 1.5 :: 2.0 :: 3.0 :: 10.0 :: 140.0 :: 170.0 :: nil; sys->print("%5s %24s %24s\n", "x", "math->lgamma", "lanczos7"); while(ns != nil) { x := hd ns; ns = tl ns; # math->lgamma returns a tuple. (i, r) := lgamma(x); g := real i * exp(r); sys->print("%5.1f %24.16g %24.16g\n", x, g, lanczos7(x)); }} lanczos7(z: real): real{ t := z + 6.5; x := 0.99999999999980993 + 676.5203681218851/z - 1259.1392167224028/(z+1.0) + 771.32342877765313/(z+2.0) - 176.61502916214059/(z+3.0) + 12.507343278686905/(z+4.0) - 0.13857109526572012/(z+5.0) + 9.9843695780195716e-6/(z+6.0) + 1.5056327351493116e-7/(z+7.0); return sqrt(2.0) * sqrt(Math->Pi) * pow(t, z - 0.5) * exp(-t) * x;}  Output:  x math->lgamma lanczos7 -0.5 -3.544907701811032 -3.544907701811089 0.1 9.513507698668729 9.51350769866875 0.5 1.772453850905516 1.772453850905516 1.0 1 0.9999999999999999 1.5 0.8862269254527581 0.8862269254527587 2.0 1 1 3.0 2 2.000000000000001 10.0 362880.0000000005 362880.0000000015 140.0 9.615723196940553e+238 9.615723196940235e+238 170.0 4.269068009004526e+304 Infinity  Lua Uses the wp:Reciprocal gamma function to calculate small values. gamma, coeff, quad, qui, set = 0.577215664901, -0.65587807152056, -0.042002635033944, 0.16653861138228, -0.042197734555571function recigamma(z) return z + gamma * z^2 + coeff * z^3 + quad * z^4 + qui * z^5 + set * z^6end function gammafunc(z) if z == 1 then return 1 elseif math.abs(z) <= 0.5 then return 1 / recigamma(z) else return (z - 1) * gammafunc(z-1) endend M2000 Interpreter  Module PrepareLambdaFunctions { Const e = [email protected] Exp= Lambda e (x) -> e^x gammaStirling=lambda e (x As decimal)->Sqrt(2.0 * pi / x) * ((x / e) ^ x) Rad2Deg =Lambda pidivby180=pi/180 (RadAngle)->RadAngle / pidivby180 Dim p(9) p(0)[email protected], [email protected], [email protected], [email protected] p(4)[email protected], [email protected], [email protected], [email protected] p(8)[email protected] gammaLanczos =Lambda p(), Rad2Deg, Exp (x As decimal) -> { Def Decimal a, t If x < 0.5 Then =pi / (Sin(Rad2Deg(pi * x)) *Lambda(1-x)) : Exit x -= [email protected] a=p(0) t = x + [email protected] For i= [email protected] To [email protected] { a += p(i) / (x + i) } = Sqrt(2.0 * pi) * (t ^ (x + 0.5)) * Exp(-t) * a } Push gammaStirling, gammaLanczos}Call PrepareLambdaFunctionsRead gammaLanczos, gammaStirlingFont "Courier New"Form 120, 40document doc$="     χ                        Stirling                     Lanczos"+{}Print $(2,20),"x", "Stirling",@(55),"Lanczos",$(0)PrintFor d = 0.1 To 2 step 0.1      Print   $("0.00"), d, Print$("0.000000000000000"), gammaStirling(d),      Print  $("0.0000000000000000000000000000"), gammaLanczos(d) doc$=format$("{0:-10} {1:-30} {2:-34}",str$(d,"0.00"), str$(gammaStirling(d),"0.000000000000000"), str$(gammaLanczos(d),"0.0000000000000000000000000000"))+{      }Next dPrint $("")clipboard doc$ 
    χ                        Stirling                     Lanczos
0.10               5.697187148977170       9.5135076986687024462927178610
0.20               3.325998424022390       4.5908437119987955107204909409
0.30               2.362530036269620       2.9915689876875914865114179656
0.40               1.841476335936240       2.2181595437576816416854441034
0.50               1.520346901066280       1.7724538509055147387430498835
0.60               1.307158857448360       1.4891922488128208508983507496
0.70               1.159053292113920       1.2980553326475564892857625396
0.80               1.053370968425610       1.1642297137253055422419914101
0.90               0.977061507877695       1.0686287021193206646594133376
1.00               0.922137008895789       1.0000000000000007024882980221
1.10               0.883489953168704       0.9513507698668745807357371716
1.20               0.857755335396591       0.9181687423997605348002977483
1.30               0.842678259448392       0.8974706963062785326402091223
1.40               0.836744548637082       0.8872638175030748314253582066
1.50               0.838956552526496       0.8862269254527587632845492097
1.60               0.848693242152574       0.8935153492876912865293624528
1.70               0.865621471793884       0.9086387328532921150064803085
1.80               0.889639635287995       0.9313837709802428420608295699
1.90               0.920842721894229       0.9617658319073891431109375442
2.00               0.959502175744492       1.0000000000000015609456469406


Maple

Built-in method that accepts any value.

GAMMA(17/2);GAMMA(7*I);M := Matrix(2, 3, 'fill' = -3.6);MTM:-gamma(M);
Output:
2027025*sqrt(Pi)*(1/256)
GAMMA(7*I)
Matrix(2, 3, [[.2468571430, .2468571430, .2468571430], [.2468571430, .2468571430, .2468571430]])

Mathematica

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

Gamma[x]

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


Maxima

fpprec: 30$gamma_coeff(n) := block([a: makelist(1, n)], a[2]: bfloat(%gamma), for k from 3 thru n do a[k]: bfloat((sum((-1)^j * zeta(j) * a[k - j], j, 2, k - 1) - a[2] * a[k - 1]) / (1 - k * a[1])), a)$ poleval(a, x) := block([y: 0],   for k from length(a) thru 1 step -1 do      y: y * x + a[k],   y)$gc: gamma_coeff(20)$ gamma_approx(x) := block([y: 1],   while x > 2 do (x: x - 1, y: y * x),   y / (poleval(gc, x - 1)))$gamma_approx(12.3b0) - gamma(12.3b0);/* -9.25224705314470500985141176997b-15 */ Modula-3 Translation of: Ada 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 Translation of: Ada 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 Perl 6 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

(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

 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 Translation of: Ada _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 Translation of: Ada $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  Seed7 Translation of: Ada $ 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
`