Gamma function

From Rosetta Code
Task
Gamma function
You are encouraged to solve this task according to the task description, using any language you may know.
Task

Implement one algorithm (or more) to compute the 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:

This suggests a straightforward (but inefficient) way of computing the through numerical integration.


Better suggested methods:



11l

Translation of: Python
V _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
       ]
F gamma(x)
   V y = x - 1.0
   V sm = :_a.last
   L(n) (:_a.len-2 .. 0).step(-1)
      sm = sm * y + :_a[n]
   R 1.0 / sm

L(i) 1..10
   print(‘#.14’.format(gamma(i / 3.0)))
Output:
2.67893853470775
1.35411793942640
1.00000000000000
0.89297951156925
0.90274529295093
1.00000000000000
1.19063934875900
1.50457548825154
1.99999999999397
2.77815847933857

360 Assembly

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

GAMMAT   CSECT
         USING GAMMAT,R13
SAVEAR   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,NI
LOOPI    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     LOOPJ
ENDLOOPJ 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
*        ----  DATA
NI       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    E
SUM      DS    E
GAMMA    DS    E
WTOMSG   DS    0F
         DC    AL2(L'BUF),XL2'0000'
BUF      DC    CL80' '
*        Subroutine             Convertion Float->Display
CONVERT  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 bit
CONVERT1 LR    R8,R9
         N     R8,=X'00FFFFFF'
         BNZ   CONVERT2
         SR    R9,R9
         B     CONVERT7
CONVERT2 LR    R8,R9
         N     R8,=X'FF000000'  characteristic
         SRL   R8,24
         CH    R8,=H'64'
         BH    CONVERT3
         SR    R9,R9
         B     CONVERT7
CONVERT3 CH    R8,=H'72'        2**32
         BNH   CONVERT4
         L     R9,=X'7FFFFFFF'  biggest
         B     CONVERT6
CONVERT4 SR    R8,R8
         SLDL  R8,8
         CH    R8,=H'72'
         BL    CONVERT5
         CH    R9,=H'0'
         BP    CONVERT5
         L     R9,=X'7FFFFFFF'
         B     CONVERT6
CONVERT5 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,R9
CONVERT7 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    F
CONVERTB DS    F                
CONVERTS DS    X
CONVERTM DS    CL13       
CONVERTD DS    CL12
CONVERTP 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

Ada

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 is
begin
   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/sum
END # taylor gamma #;

LONG REAL long e = long exp(1);

PROC sterling gamma = (LONG REAL n)LONG REAL:
( # improves for values much greater then 1 #
  long sqrt(2*long pi/n)*(n/long e)**n
);
 
PROC factorial = (LONG INT n)LONG REAL:
(
  IF n=0 OR n=1 THEN 1
  ELIF n=2 THEN 2
  ELSE n*factorial(n-1) FI
);
 
REF[]LONG REAL fm := NIL;

PROC sponge gamma = (LONG REAL x)LONG REAL:
(
  INT a = 12; # alter to get required precision #
  REF []LONG REAL fm := NIL;
  LONG REAL res;
 
  IF fm :=: REF[]LONG REAL(NIL) THEN
    fm := HEAP[0:a-1]LONG REAL;
    fm[0] := long sqrt(LONG 2*long pi);
    FOR k TO a-1 DO
      fm[k] := (((k-1) MOD 2=0) | 1 | -1) * long exp(a-k) *
	(a-k) **(k-LONG 0.5) / factorial(k-1)
    OD
  FI;
  res := fm[0];
  FOR k TO a-1 DO
    res +:= fm[k] / ( x + k )
  OD;
  res *:= long exp(-(x+a)) * (x+a)**(x + LONG 0.5);
  res/x
);

FORMAT real fmt = $g(-real width, real width - 2)$;
FORMAT long real fmt16 = $g(-17, 17 - 2)$; # accurate to about 16 decimal places #

[]STRING methods = ("Genie", "Lanczos", "Sponge", "Taylor","Stirling");

printf(($11xg12xg12xg13xg13xgl$,methods));

FORMAT sample fmt = $"gamma("g(-3,1)")="f(real fmt)n(UPB methods-1)(", "f(long real fmt16))l$;
FORMAT sqr sample fmt = $"gamma("g(-3,1)")**2="f(real fmt)n(UPB methods-1)(", "f(long real fmt16))l$;
FORMAT sample exp fmt = $"gamma("g(-3)")="g(-15,11,0)n(UPB methods-1)(","g(-18,14,0))l$;

PROC sample = (LONG REAL x)[]LONG REAL:
  (gamma(SHORTEN x), lanczos gamma(x), sponge gamma(x), taylor gamma(x), sterling gamma(x));

FOR i FROM 1 TO 20 DO
  LONG REAL x = i / LONG 10;
  printf((sample fmt, x, sample(x)));
  IF i = 5 THEN # insert special case of a half #
    printf((sqr sample fmt,
            x, gamma(SHORTEN x)**2,  lanczos gamma(x)**2, sponge gamma(x)**2,
            taylor gamma(x)**2, sterling gamma(x)**2))
  FI
OD;
FOR x FROM 10 BY 10 TO 70 DO
  printf((sample exp fmt, x, sample(x)))
OD
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

Arturo

A: @[
    1.00000000000000000000     0.57721566490153286061 neg 0.65587807152025388108
neg 0.04200263503409523553     0.16653861138229148950 neg 0.04219773455554433675
neg 0.00962197152787697356     0.00721894324666309954 neg 0.00116516759185906511
neg 0.00021524167411495097     0.00012805028238811619 neg 0.00002013485478078824
neg 0.00000125049348214267     0.00000113302723198170 neg 0.00000020563384169776
    0.00000000611609510448     0.00000000500200764447 neg 0.00000000118127457049
    0.00000000010434267117     0.00000000000778226344 neg 0.00000000000369680562
    0.00000000000051003703 neg 0.00000000000002058326 neg 0.00000000000000534812
    0.00000000000000122678 neg 0.00000000000000011813     0.00000000000000000119
    0.00000000000000000141 neg 0.00000000000000000023     0.00000000000000000002 
]

ourGamma: function [x][
    y: x - 1
    result: last A 
    loop ((size A)-1)..0 'n ->
        result: (result*y) + get A n
    result: 1 // result
    return result
]

loop 1..10 'z [
    v1: ourGamma z // 3
    v2: gamma z // 3
    print [
        pad (to :string z)++" =>" 10
        pad (to :string v1)++" ~" 30 
        pad (to :string v2)++" :" 30
        pad (to :string v1-v2) 30
    ]
]
Output:
      1 =>            2.678938534707748 ~            2.678938534707748 :          4.440892098500626e-16 
      2 =>              1.3541179394264 ~              1.3541179394264 :                            0.0 
      3 =>                          1.0 ~                          1.0 :                            0.0 
      4 =>           0.8929795115692493 ~           0.8929795115692493 :                            0.0 
      5 =>           0.9027452929509336 ~           0.9027452929509336 :                            0.0 
      6 =>                          1.0 ~                          1.0 :                            0.0 
      7 =>            1.190639348758999 ~            1.190639348758999 :                            0.0 
      8 =>            1.504575488251335 ~            1.504575488251556 :         -2.204902926905561e-13 
      9 =>            1.999999999908069 ~                          2.0 :         -9.193090733106146e-11 
     10 =>            2.778158462440097 ~            2.778158480437665 :         -1.799756743636749e-08

AutoHotkey

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

/*
Here is the upper incomplete Gamma function. Omitting or setting
the second parameter to 0 we get the (complete) Gamma function.
The code is based on: "Computation of Special Functions" Zhang and Jin,
John Wiley and Sons, 1996
*/

SetFormat FloatFast, 0.9e 

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

GAMMA(a,x=0) {  ; upper incomplete gamma: Integral(t**(a-1)*e**-t, t = x..inf) 
   If (a > 171 || x < 0) 
      Return 2.e308   ; overflow 

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

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

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

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

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

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

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

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

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

AWK

# syntax: GAWK -f GAMMA_FUNCTION.AWK
BEGIN {
    e = (1+1/100000)^100000
    pi = atan2(0,-1)

    print("X    Stirling")
    for (i=1; i<=20; i++) {
      d = i / 10
      printf("%4.2f %9.5f\n",d,gamma_stirling(d))
    }
    exit(0)
}
function gamma_stirling(x) {
    return sqrt(2*pi/x) * pow(x/e,x)
}
function pow(a,b) {
    return exp(b*log(a))
}
Output:
X    Stirling
0.10   5.69719
0.20   3.32600
0.30   2.36253
0.40   1.84148
0.50   1.52035
0.60   1.30716
0.70   1.15906
0.80   1.05338
0.90   0.97707
1.00   0.92214
1.10   0.88349
1.20   0.85776
1.30   0.84268
1.40   0.83675
1.50   0.83896
1.60   0.84870
1.70   0.86563
1.80   0.88965
1.90   0.92085
2.00   0.95951

BASIC

ANSI BASIC

Translation of: BBC BASIC

- Lanczos method.

Works with: Decimal BASIC
100 DECLARE EXTERNAL FUNCTION FNlngamma
110 
120 DEF FNgamma(z) = EXP(FNlngamma(z))
130 
140 FOR x = 0.1 TO 2.05 STEP 0.1
150    PRINT USING$("#.#",x), USING$("##.############", FNgamma(x))
160 NEXT x
170 END
180 
190 EXTERNAL FUNCTION FNlngamma(z)
200 DIM lz(0 TO 6)
210 RESTORE
220 MAT READ lz
230 DATA 1.000000000190015, 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.0012086509738662, -0.000005395239385
240 IF z < 0.5 THEN 
250    LET FNlngamma = LOG(PI / SIN(PI * z)) - FNlngamma(1.0 - z)
260    EXIT FUNCTION
270 END IF
280 LET z = z - 1.0
290 LET b = z + 5.5
300 LET a = lz(0)
310 FOR i = 1 TO 6
320    LET a  = a + lz(i) / (z + i)
330 NEXT i
340 LET FNlngamma = (LOG(SQR(2*PI)) + LOG(a) - b) + LOG(b) * (z+0.5)
350 END FUNCTION
Output:
 .1                      9.513507698670
 .2                      4.590843712000
 .3                      2.991568987689
 .4                      2.218159543760
 .5                      1.772453850902
 .6                      1.489192248811
 .7                      1.298055332647
 .8                      1.164229713725
 .9                      1.068628702119
1.0                      1.000000000000
1.1                       .951350769867
1.2                       .918168742400
1.3                       .897470696306
1.4                       .887263817503
1.5                       .886226925453
1.6                       .893515349288
1.7                       .908638732853
1.8                       .931383770980
1.9                       .961765831907
2.0                      1.000000000000

BASIC256

Translation of: FreeBASIC

- Stirling method.

Translation of: Phix

- Lanczos method.

print " x       Stirling         Lanczos"
print
for i = 1 to 20
    d = i / 10.0
    print d;
    print chr(9); ljust(string(gammaStirling(d)), 13, "0");
    print chr(9); ljust(string(gammaLanczos(d)),  13, "0")
next i
end 

function gammaStirling (x)
    e = exp(1)	# e is not predefined in BASIC256
    return sqr(2.0 * pi / x) * ((x / e) ^ x)
end function

function gammaLanczos (x)
    dim p = {0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7}

    g = 7
    if x < 0.5 then return pi / (sin(pi * x) * gammaLanczos(1-x))
    x -= 1
    a = p[0]
    t = x + g + 0.5

    for i = 1 to 8
	a += p[i] / (x + i)
    next i
    return sqr(2.0 * pi) * (t ^ (x + 0.5)) * exp(-t) * a
end function

BBC BASIC

Uses the Lanczos approximation.

      *FLOAT64
      INSTALL @lib$+"FNUSING"
      
      FOR x = 0.1 TO 2.05 STEP 0.1
        PRINT FNusing("#.#",x), FNusing("##.############", FNgamma(x))
      NEXT
      END
      
      DEF FNgamma(z) = EXP(FNlngamma(z))
      
      DEF FNlngamma(z)
      LOCAL a, b, i%, lz()
      DIM lz(6)
      lz() = 1.000000000190015, 76.18009172947146, -86.50532032941677, \
      \ 24.01409824083091, -1.231739572450155, 0.0012086509738662, -0.000005395239385
      IF z < 0.5 THEN = LN(PI / SIN(PI * z)) - FNlngamma(1.0 - z)
      z -= 1.0
      b = z + 5.5
      a = lz(0)
      FOR i% = 1 TO 6
        a += lz(i%) / (z + i%)
      NEXT
      = (LNSQR(2*PI) + LN(a) - b) + LN(b) * (z+0.5)

Output:

0.1        9.513507698670
0.2        4.590843712000
0.3        2.991568987689
0.4        2.218159543760
0.5        1.772453850902
0.6        1.489192248811
0.7        1.298055332647
0.8        1.164229713725
0.9        1.068628702119
1.0        1.000000000000
1.1        0.951350769867
1.2        0.918168742400
1.3        0.897470696306
1.4        0.887263817503
1.5        0.886226925453
1.6        0.893515349288
1.7        0.908638732853
1.8        0.931383770980
1.9        0.961765831907
2.0        1.000000000000

C

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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_sf_gamma.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/* very simple approximation */
double st_gamma(double x)
{
  return sqrt(2.0*M_PI/x)*pow(x/M_E, x);
}

#define A 12
double sp_gamma(double z)
{
  const int a = A;
  static double c_space[A];
  static double *c = NULL;
  int k;
  double accm;
 
  if ( c == NULL ) {
    double k1_factrl = 1.0; /* (k - 1)!*(-1)^k with 0!==1*/
    c = c_space;
    c[0] = sqrt(2.0*M_PI);
    for(k=1; k < a; k++) {
      c[k] = exp(a-k) * pow(a-k, k-0.5) / k1_factrl;
	  k1_factrl *= -k;
    }
  }
  accm = c[0];
  for(k=1; k < a; k++) {
    accm += c[k] / ( z + k );
  }
  accm *= exp(-(z+a)) * pow(z+a, z+0.5); /* Gamma(z+1) */
  return accm/z;
}

int main()
{
  double x;


  printf("%15s%15s%15s%15s\n", "Stirling", "Spouge", "GSL", "libm");
  for(x=1.0; x <= 10.0; x+=1.0) {
    printf("%15.8lf%15.8lf%15.8lf%15.8lf\n", st_gamma(x/3.0), sp_gamma(x/3.0), 
	   gsl_sf_gamma(x/3.0), tgamma(x/3.0));
  }
  return 0;
}
Output:
       Stirling         Spouge            GSL           libm
     2.15697602     2.67893853     2.67893853     2.67893853
     1.20285073     1.35411794     1.35411794     1.35411794
     0.92213701     1.00000000     1.00000000     1.00000000
     0.83974270     0.89297951     0.89297951     0.89297951
     0.85919025     0.90274529     0.90274529     0.90274529
     0.95950218     1.00000000     1.00000000     1.00000000
     1.14910642     1.19063935     1.19063935     1.19063935
     1.45849038     1.50457549     1.50457549     1.50457549
     1.94540320     2.00000000     2.00000000     2.00000000
     2.70976382     2.77815848     2.77815848     2.77815848

C#

This is just rewritten from the Wikipedia Lanczos article. Works with complex numbers as well as reals.

using System;
using System.Numerics;

static int g = 7;
static double[] p = {0.99999999999980993, 676.5203681218851, -1259.1392167224028,
	     771.32342877765313, -176.61502916214059, 12.507343278686905,
	     -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
		 
Complex Gamma(Complex z)
{
    // Reflection formula
    if (z.Real < 0.5)
	{
        return Math.PI / (Complex.Sin( Math.PI * z) * Gamma(1 - z));
	}
    else
	{
        z -= 1;
        Complex x = p[0];
        for (var i = 1; i < g + 2; i++)
		{
            x += p[i]/(z+i);
		}
        Complex t = z + g + 0.5;
        return Complex.Sqrt(2 * Math.PI) * (Complex.Pow(t, z + 0.5)) * Complex.Exp(-t) * x;
	}
}

C++

#include <math.h>
#include <numbers>
#include <stdio.h>
#include <vector>

// Calculate the coefficients used by Spouge's approximation (based on the C
// implemetation)
std::vector<double> CalculateCoefficients(int numCoeff)
{
    std::vector<double> c(numCoeff);
    double k1_factrl = 1.0;
    c[0] = sqrt(2.0 * std::numbers::pi);
    for(size_t k=1; k < numCoeff; k++)
    {
        c[k] = exp(numCoeff-k) * pow(numCoeff-k, k-0.5) / k1_factrl;
        k1_factrl *= -(double)k;
    }
    return c;
}

// The Spouge approximation
double Gamma(const std::vector<double>& coeffs, double x)
{
        const size_t numCoeff = coeffs.size();
        double accm = coeffs[0];
        for(size_t k=1; k < numCoeff; k++)
        {
            accm += coeffs[k] / ( x + k );
        }
        accm *= exp(-(x+numCoeff)) * pow(x+numCoeff, x+0.5);
        return accm/x;
}

int main()
{
    // estimate the gamma function with 1, 4, and 10 coefficients
    const auto coeff1 = CalculateCoefficients(1);
    const auto coeff4 = CalculateCoefficients(4);
    const auto coeff10 = CalculateCoefficients(10);

    const auto inputs = std::vector<double>{
        0.001, 0.01, 0.1, 0.5, 1.0,
        1.461632145, // minimum of the gamma function
        2, 2.5, 3, 4, 5, 6, 7, 8, 9, 10, 50, 100, 
        150 // causes overflow for this implemetation
        };
    
    printf("%16s%16s%16s%16s%16s\n", "gamma( x ) =", "Spouge 1", "Spouge 4", "Spouge 10", "built-in");
    for(auto x : inputs) 
    {
        printf("gamma(%7.3f) = %16.10g %16.10g %16.10g %16.10g\n", 
            x,
            Gamma(coeff1, x),
            Gamma(coeff4, x), 
            Gamma(coeff10, x), 
            std::tgamma(x)); // built-in gamma function
    }
}
Output:
    gamma( x ) =        Spouge 1        Spouge 4       Spouge 10        built-in
gamma(  0.001) =      921.6767466      999.4237321      999.4237725      999.4237725
gamma(  0.010) =      91.76063453      99.43258106      99.43258512      99.43258512
gamma(  0.100) =      8.834899532      9.513507269      9.513507699      9.513507699
gamma(  0.500) =      1.677913105      1.772453737      1.772453851      1.772453851
gamma(  1.000) =     0.9595021757     0.9999999124                1                1
gamma(  1.462) =     0.8562774501     0.8856030992     0.8856031944     0.8856031944
gamma(  2.000) =     0.9727015986     0.9999998717                1                1
gamma(  2.500) =      1.298145499      1.329340195      1.329340388      1.329340388
gamma(  3.000) =      1.958847928      1.999999681                2                2
gamma(  4.000) =      5.900958398      5.999998905                6                6
gamma(  5.000) =      23.66927282      23.99999523               24               24
gamma(  6.000) =      118.5808531      119.9999749              120              120
gamma(  7.000) =      712.5427759      719.9998442              720              720
gamma(  8.000) =      4993.567678       5039.99889             5040             5040
gamma(  9.000) =      39985.50687      40319.99104            40320            40320
gamma( 10.000) =      360142.0459      362879.9193           362880           362880
gamma( 50.000) =  6.072887637e+62  6.082817933e+62   6.08281864e+62   6.08281864e+62
gamma(100.000) = 9.324924563e+155 9.332620912e+155 9.332621544e+155 9.332621544e+155
gamma(150.000) =              inf              inf              inf 3.808922638e+260

Clojure

(defn gamma 
  "Returns Gamma(z + 1 = number) using Lanczos approximation."
  [number]
  (if (< number 0.5)
       (/ Math/PI (* (Math/sin (* Math/PI number))
	             (gamma (- 1 number))))
       (let [n (dec number)
      	     c [0.99999999999980993 676.5203681218851 -1259.1392167224028
	        771.32342877765313 -176.61502916214059 12.507343278686905
	        -0.13857109526572012 9.9843695780195716e-6 1.5056327351493116e-7]]
         (* (Math/sqrt (* 2 Math/PI))
	    (Math/pow (+ n 7 0.5) (+ n 0.5))
	    (Math/exp (- (+ n 7 0.5)))
            (+ (first c) 
               (apply + (map-indexed #(/ %2 (+ n %1 1)) (next c))))))))
Output:
(map #(printf "%.1f %.4f\n" % (gamma %)) (map #(float (/ % 10)) (range 1 31)))
0.1 9.5135
0.2 4.5908
0.3 2.9916
0.4 2.2182
0.5 1.7725
0.6 1.4892
0.7 1.2981
0.8 1.1642
0.9 1.0686
1.0 1.0000
1.1 0.9514
1.2 0.9182
1.3 0.8975
1.4 0.8873
1.5 0.8862
1.6 0.8935
1.7 0.9086
1.8 0.9314
1.9 0.9618
2.0 1.0000
2.1 1.0465
2.2 1.1018
2.3 1.1667
2.4 1.2422
2.5 1.3293
2.6 1.4296
2.7 1.5447
2.8 1.6765
2.9 1.8274
3.0 2.0000

Common Lisp

; Taylor series coefficients
(defconstant tcoeff
  '( 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))

; number of coefficients
(defconstant numcoeff (length tcoeff))

(defun gamma (x)
  (let ((y (- x 1.0))
        (sum (nth (- numcoeff 1) tcoeff)))
    (loop for i from (- numcoeff 2) downto 0 do 
          (setf sum (+ (* sum y) (nth i tcoeff))))
    (/ 1.0 sum)))

(loop for i from 1 to 10
   do (
     format t "~12,10f~%" (gamma (/ i 3.0))))
Produces:
2.6789380000
1.3541179000
1.0000000000
0.8929794500
0.9027453000
1.0000000000
1.1906393000
1.5045753000
1.9999995000
2.7781580000

Crystal

Taylor Series | Lanczos Method | Builtin Function

# Taylor Series
def 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 ]
end

def taylor_gamma(x)
  y = x.to_f - 1
  1.0 / a.reverse.reduce(0) { |sum, an| sum * y + an }
end

# Lanczos Method
def p
  [ 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 ]
end

def lanczos_gamma(z)
  # Reflection formula
  z = z.to_f
  if z < 0.5
    Math::PI / (Math.sin(Math::PI * z) * lanczos_gamma(1 - z))
  else
    z -= 1
    x = p[0]
    (1..p.size - 1).each { |i| x += p[i] / (z + i) }
    t = z + p.size - 1.5
    Math.sqrt(2 * Math::PI) * t**(z + 0.5) * Math.exp(-t) * x
  end
end

puts "                Taylor Series         Lanczos Method        Builtin Function"
(1..27).each { |i| n = i/3.0; puts "gamma(%.2f) = %.14e  %.14e  %.14e" % [n, taylor_gamma(n), lanczos_gamma(n), Math.gamma(n)] }
Output:
 
                Taylor Series         Lanczos Method        Builtin Function
gamma(0.33) = 2.67893853470775e+00  2.67893853470775e+00  2.67893853470775e+00
gamma(0.67) = 1.35411793942640e+00  1.35411793942640e+00  1.35411793942640e+00
gamma(1.00) = 1.00000000000000e+00  1.00000000000000e+00  1.00000000000000e+00
gamma(1.33) = 8.92979511569249e-01  8.92979511569249e-01  8.92979511569249e-01
gamma(1.67) = 9.02745292950934e-01  9.02745292950935e-01  9.02745292950934e-01
gamma(2.00) = 1.00000000000000e+00  1.00000000000000e+00  1.00000000000000e+00
gamma(2.33) = 1.19063934875900e+00  1.19063934875900e+00  1.19063934875900e+00
gamma(2.67) = 1.50457548825154e+00  1.50457548825156e+00  1.50457548825156e+00
gamma(3.00) = 1.99999999999397e+00  2.00000000000000e+00  2.00000000000000e+00
gamma(3.33) = 2.77815847933857e+00  2.77815848043767e+00  2.77815848043766e+00
gamma(3.67) = 4.01220118377482e+00  4.01220130200415e+00  4.01220130200415e+00
gamma(4.00) = 5.99999141007240e+00  6.00000000000001e+00  6.00000000000000e+00
gamma(4.33) = 9.26006653812473e+00  9.26052826812555e+00  9.26052826812554e+00
gamma(4.67) = 1.46918499266721e+01  1.47114047740152e+01  1.47114047740152e+01
gamma(5.00) = 2.33327665969918e+01  2.40000000000000e+01  2.40000000000000e+01
gamma(5.33) = 2.65211050660964e+01  4.01289558285441e+01  4.01289558285440e+01
gamma(5.67) = 7.70471336505311e+00  6.86532222787379e+01  6.86532222787377e+01
gamma(6.00) = 1.10934146590517e+00  1.20000000000000e+02  1.20000000000000e+02
gamma(6.33) = 1.64621072447163e-01  2.14021097752236e+02  2.14021097752235e+02
gamma(6.67) = 2.72102446536397e-02  3.89034926246181e+02  3.89034926246180e+02
gamma(7.00) = 4.98014348954507e-03  7.20000000000002e+02  7.20000000000000e+02
gamma(7.33) = 9.98845907123850e-04  1.35546695243082e+03  1.35546695243082e+03
gamma(7.67) = 2.17513475446479e-04  2.59356617497454e+03  2.59356617497454e+03
gamma(8.00) = 5.10217006678528e-05  5.04000000000001e+03  5.04000000000000e+03
gamma(8.33) = 1.28035516395359e-05  9.94009098449271e+03  9.94009098449270e+03
gamma(8.67) = 3.41689149138074e-06  1.98840073414715e+04  1.98840073414714e+04
gamma(9.00) = 9.64721467591131e-07  4.03200000000001e+04  4.03200000000000e+04

D

import std.stdio, std.math, std.mathspecial;

real taylorGamma(in real x) pure nothrow @safe @nogc {
    static immutable real[30] table = [
     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];

    immutable real y = x - 1.0L;
    real sm = table[$ - 1];
    foreach_reverse (immutable an; table[0 .. $ - 1])
        sm = sm * y + an;
    return 1.0L / sm;
}

real lanczosGamma(real z) pure nothrow @safe @nogc {
    // Coefficients used by the GNU Scientific Library.
    // http://en.wikipedia.org/wiki/Lanczos_approximation
    enum g = 7;
    static immutable real[9] table =
        [    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.5L) {
        return PI / (sin(PI * z) * lanczosGamma(1 - z));
    } else {
        z -= 1;
        real x = table[0];
        foreach (immutable i; 1 .. g + 2)
            x += table[i] / (z + i);
        immutable real t = z + g + 0.5L;
        return sqrt(2 * PI) * t ^^ (z + 0.5L) * exp(-t) * x;
    }
}

void main() {
    foreach (immutable i; 1 .. 11) {
        immutable real x = i / 3.0L;
        writefln("%f: %20.19e %20.19e %20.19e", x,
                 x.taylorGamma, x.lanczosGamma, x.gamma);
    }
}
Output:
0.333333: 2.6789385347077476335e+00 2.6789385347077470551e+00 2.6789385347077476339e+00
0.666667: 1.3541179394264004169e+00 1.3541179394264007092e+00 1.3541179394264004170e+00
1.000000: 1.0000000000000000000e+00 1.0000000000000002126e+00 1.0000000000000000000e+00
1.333333: 8.9297951156924921124e-01 8.9297951156924947465e-01 8.9297951156924921132e-01
1.666667: 9.0274529295093361132e-01 9.0274529295093396555e-01 9.0274529295093361123e-01
2.000000: 1.0000000000000000000e+00 1.0000000000000004903e+00 1.0000000000000000000e+00
2.333333: 1.1906393487589989474e+00 1.1906393487589996490e+00 1.1906393487589989482e+00
2.666667: 1.5045754882515545787e+00 1.5045754882515570474e+00 1.5045754882515560190e+00
3.000000: 1.9999999999992207405e+00 2.0000000000000015575e+00 2.0000000000000000000e+00
3.333333: 2.7781584802531739378e+00 2.7781584804376666336e+00 2.7781584804376642124e+00

Delphi

Library: System.Math
Translation of: Go
program Gamma_function;

{$APPTYPE CONSOLE}

uses
  System.SysUtils,
  System.Math;

function lanczos7(z: double): Double;
begin
  var t := z + 6.5;
  var x := 0.99999999999980993 + 676.5203681218851 / z - 1259.1392167224028 / (z
    + 1) + 771.32342877765313 / (z + 2) - 176.61502916214059 / (z + 3) +
    12.507343278686905 / (z + 4) - 0.13857109526572012 / (z + 5) +
    9.9843695780195716e-6 / (z + 6) + 1.5056327351493116e-7 / (z + 7);

  Result := Sqrt(2) * Sqrt(pi) * Power(t, z - 0.5) * exp(-t) * x;
end;

begin
  var xs: TArray<double> := [-0.5, 0.1, 0.5, 1, 1.5, 2, 3, 10, 140, 170];
  writeln('    x              Lanczos7');
  for var x in xs do
    writeln(format('%5.1f %24.16g', [x, lanczos7(x)]));
  readln;
end.
Output:
    x              Lanczos7
 -0,5       -3,544907701811089
  0,1        9,513507698668747
  0,5        1,772453850905517
  1,0                        1
  1,5       0,8862269254527583
  2,0                        1
  3,0        2,000000000000002
 10,0        362880,0000000007
140,0    9,615723196940235E238
170,0    4,269068009004271E304

DuckDB

Works with: DuckDB version V1.0

This entry contrasts DuckDB's builtin gamma() function with an implementation based on the Taylor series expansion.

create or replace function gamma_taylor(x) as (
  select gamma
  from ( select [
    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
    ] as a,
    length(a) as n,
    (with recursive cte as (
       select 2 as an, a[n] as acc
       union all
       select an+1 as an,
         (acc * (x-1)) + a[1 + n-an] as acc
       from cte
       where an <= n
     )
     select 1 / last(acc) from cte) as gamma
  )
  limit 1
);

## Comparison
select x, gamma, taylor, 2 * @(gamma-taylor)/(gamma+taylor) as "%"
from (select x, gamma(x) as gamma, gamma_taylor(x) as taylor
      from (select x/3 as x from range(1,11) t(x))) ;
Output:
┌────────────────────┬────────────────────┬────────────────────┬────────────────────────┐
│         x          │       gamma        │       taylor       │           %            │
│       double       │       double       │       double       │         double         │
├────────────────────┼────────────────────┼────────────────────┼────────────────────────┤
│ 0.3333333333333333 │  2.678938534707748 │ 2.6789385347077483 │  1.657705856616488e-16 │
│ 0.6666666666666666 │ 1.3541179394264005 │ 1.3541179394264007 │  1.639773009868613e-16 │
│                1.0 │                1.0 │                1.0 │                    0.0 │
│ 1.3333333333333333 │ 0.8929795115692493 │ 0.8929795115692493 │                    0.0 │
│ 1.6666666666666667 │ 0.9027452929509336 │ 0.9027452929509336 │                    0.0 │
│                2.0 │                1.0 │                1.0 │                    0.0 │
│ 2.3333333333333335 │  1.190639348758999 │ 1.1906393487589988 │  1.864919088693549e-16 │
│ 2.6666666666666665 │ 1.5045754882515558 │ 1.5045754882515459 │  6.641080689967907e-15 │
│                3.0 │                2.0 │ 1.9999999999939666 │ 3.0166980025160257e-12 │
│ 3.3333333333333335 │ 2.7781584804376647 │ 2.7781584787715414 │  5.997221726571757e-10 │
├────────────────────┴────────────────────┴────────────────────┴────────────────────────┤
│ 10 rows                                                                     4 columns │
└───────────────────────────────────────────────────────────────────────────────────────┘

EasyLang

Translation of: AWK
e = 2.718281828459
func stirling x .
   return sqrt (2 * pi / x) * pow (x / e) x
.
print " X    Stirling"
for i to 20
   d = i / 10
   numfmt 2 4
   write d & "    "
   numfmt 3 4
   print stirling d
.

Elixir

Translation of: Ruby
defmodule Gamma do
  @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 ]
     |> Enum.reverse
  def taylor(x) do
    y = x - 1
    1 / Enum.reduce(@a, 0, fn a,sum -> sum * y + a end)
  end
end

Enum.each(Enum.map(1..10, &(&1/3)), fn x ->
  :io.format "~f  ~18.15f~n", [x, Gamma.taylor(x)]
end)
Output:
0.333333   2.678938534707748
0.666667   1.354117939426401
1.000000   1.000000000000000
1.333333   0.892979511569249
1.666667   0.902745292950934
2.000000   1.000000000000000
2.333333   1.190639348758999
2.666667   1.504575488251540
3.000000   1.999999999993968
3.333333   2.778158479338573

F#

Solved using the Lanczos Coefficients described in Numerical Recipes (Press et al)

open System

let gamma z = 
    let lanczosCoefficients = [76.18009172947146;-86.50532032941677;24.01409824083091;-1.231739572450155;0.1208650973866179e-2;-0.5395239384953e-5]
    let rec sumCoefficients acc i coefficients =
        match coefficients with
        | []   -> acc
        | h::t -> sumCoefficients (acc + (h/i)) (i+1.0) t
    let gamma = 5.0
    let x = z - 1.0
    Math.Pow(x + gamma + 0.5, x + 0.5) * Math.Exp( -(x + gamma + 0.5) ) * Math.Sqrt( 2.0 * Math.PI ) * sumCoefficients 1.000000000190015 (x + 1.0) lanczosCoefficients

seq { for i in 1 .. 20 do yield ((double)i/10.0) } |> Seq.iter ( fun v -> System.Console.WriteLine("{0} : {1}", v, gamma v ) )
seq { for i in 1 .. 10 do yield ((double)i*10.0) } |> Seq.iter ( fun v -> System.Console.WriteLine("{0} : {1}", v, gamma v ) )
Output:
0.1 : 9.51350769855015
0.2 : 4.59084371196153
0.3 : 2.99156898767207
0.4 : 2.21815954375051
0.5 : 1.77245385090205
0.6 : 1.48919224881114
0.7 : 1.29805533264677
0.8 : 1.16422971372497
0.9 : 1.06862870211921
1 : 1
1.1 : 0.951350769866919
1.2 : 0.91816874239982
1.3 : 0.897470696306335
1.4 : 0.887263817503124
1.5 : 0.886226925452797
1.6 : 0.893515349287718
1.7 : 0.908638732853309
1.8 : 0.931383770980253
1.9 : 0.961765831907391
2 : 1
10 : 362880.000000085
20 : 1.21645100409886E+17
30 : 8.84176199395902E+30
40 : 2.03978820820436E+46
50 : 6.08281864068541E+62
60 : 1.38683118555266E+80
70 : 1.71122452441801E+98
80 : 8.94618213157899E+116
90 : 1.65079551625067E+136
100 : 9.33262154536104E+155
Translation of: C#

The C# version can be translated to F# to support complex numbers:

open System.Numerics
open System

let rec gamma (z: Complex) =
    let mutable z = z
    let lanczosCoefficients = [| 676.520368121885; -1259.1392167224; 771.323428777653; -176.615029162141; 12.5073432786869; -0.13857109526572; 9.98436957801957E-06; 1.50563273514931E-07 |]

    if z.Real < 0.5 then
        Math.PI / (sin (Math.PI * z) * gamma (1.0 - z))
    else
        let mutable x = Complex.One
        z <- z - 1.0

        for i = 0 to lanczosCoefficients.Length - 1 do
            x <- x + lanczosCoefficients.[i] / (z + Complex(i, 0) + 1.0)

        let t = z + float lanczosCoefficients.Length - 0.5
        sqrt (2.0 * Math.PI) * (t ** (z + 0.5)) * exp (-t) * x

Seq.iter (fun i -> printfn "Gamma(%f) = %A" i (gamma (Complex(i, 0)))) [ 0 .. 100 ]
Seq.iter2 (fun i j -> printfn "Gamma(%f + i%f) = %A" i j (gamma (Complex(i, j)))) [ 0 .. 100 ] [ 0 .. 100 ]
Output:
Gamma(0.000000) = (NaN, NaN)
Gamma(1.000000) = (1,0000000000000049, 0)
Gamma(2.000000) = (1,0000000000000115, 0)
Gamma(3.000000) = (2,0000000000000386, 0)
Gamma(4.000000) = (6,000000000000169, 0)
Gamma(5.000000) = (24,00000000000084, 0)
Gamma(6.000000) = (120,00000000000514, 0)
Gamma(7.000000) = (720,0000000000364, 0)
Gamma(8.000000) = (5040,000000000289, 0)
Gamma(9.000000) = (40320,000000002554, 0)
Gamma(10.000000) = (362880,00000002526, 0)
Gamma(11.000000) = (3628800,0000002664, 0)
Gamma(12.000000) = (39916800,000003114, 0)
Gamma(13.000000) = (479001600,00004023, 0)
Gamma(14.000000) = (6227020800,000546, 0)
Gamma(15.000000) = (87178291200,00801, 0)
Gamma(16.000000) = (1307674368000,1267, 0)
Gamma(17.000000) = (20922789888002,08, 0)
Gamma(18.000000) = (355687428096034,7, 0)
Gamma(19.000000) = (6402373705728658, 0)
Gamma(20.000000) = (1,2164510040884514E+17, 0)
Gamma(21.000000) = (2,4329020081769083E+18, 0)
Gamma(22.000000) = (5,109094217171516E+19, 0)
Gamma(23.000000) = (1,1240007277777331E+21, 0)
Gamma(24.000000) = (2,5852016738887927E+22, 0)
Gamma(25.000000) = (6,20448401733312E+23, 0)
Gamma(26.000000) = (1,5511210043332811E+25, 0)
Gamma(27.000000) = (4,032914611266542E+26, 0)
Gamma(28.000000) = (1,0888869450419632E+28, 0)
Gamma(29.000000) = (3,048883446117507E+29, 0)
Gamma(30.000000) = (8,841761993740793E+30, 0)
Gamma(31.000000) = (2,652528598122235E+32, 0)
Gamma(32.000000) = (8,222838654178949E+33, 0)
Gamma(33.000000) = (2,631308369337265E+35, 0)
Gamma(34.000000) = (8,683317618812963E+36, 0)
Gamma(35.000000) = (2,9523279903964145E+38, 0)
Gamma(36.000000) = (1,0333147966387422E+40, 0)
Gamma(37.000000) = (3,719933267899472E+41, 0)
Gamma(38.000000) = (1,3763753091228064E+43, 0)
Gamma(39.000000) = (5,230226174666675E+44, 0)
Gamma(40.000000) = (2,0397882081200028E+46, 0)
Gamma(41.000000) = (8,159152832480012E+47, 0)
Gamma(42.000000) = (3,3452526613168034E+49, 0)
Gamma(43.000000) = (1,4050061177530564E+51, 0)
Gamma(44.000000) = (6,041526306338149E+52, 0)
Gamma(45.000000) = (2,658271574788788E+54, 0)
Gamma(46.000000) = (1,1962222086549537E+56, 0)
Gamma(47.000000) = (5,502622159812779E+57, 0)
Gamma(48.000000) = (2,586232415112008E+59, 0)
Gamma(49.000000) = (1,2413915592537631E+61, 0)
Gamma(50.000000) = (6,082818640343433E+62, 0)
Gamma(51.000000) = (3,04140932017172E+64, 0)
Gamma(52.000000) = (1,551118753287575E+66, 0)
Gamma(53.000000) = (8,065817517095389E+67, 0)
Gamma(54.000000) = (4,27488328406056E+69, 0)
Gamma(55.000000) = (2,308436973392699E+71, 0)
Gamma(56.000000) = (1,2696403353659833E+73, 0)
Gamma(57.000000) = (7,109985878049497E+74, 0)
Gamma(58.000000) = (4,0526919504882125E+76, 0)
Gamma(59.000000) = (2,350561331283167E+78, 0)
Gamma(60.000000) = (1,386831185457067E+80, 0)
Gamma(61.000000) = (8,3209871127424E+81, 0)
Gamma(62.000000) = (5,075802138772854E+83, 0)
Gamma(63.000000) = (3,1469973260391715E+85, 0)
Gamma(64.000000) = (1,9826083154046777E+87, 0)
Gamma(65.000000) = (1,2688693218589942E+89, 0)
Gamma(66.000000) = (8,247650592083449E+90, 0)
Gamma(67.000000) = (5,443449390775078E+92, 0)
Gamma(68.000000) = (3,647111091819299E+94, 0)
Gamma(69.000000) = (2,48003554243712E+96, 0)
Gamma(70.000000) = (1,711224524281613E+98, 0)
Gamma(71.000000) = (1,1978571669971308E+100, 0)
Gamma(72.000000) = (8,504785885679606E+101, 0)
Gamma(73.000000) = (6,123445837689312E+103, 0)
Gamma(74.000000) = (4,470115461513189E+105, 0)
Gamma(75.000000) = (3,307885441519755E+107, 0)
Gamma(76.000000) = (2,4809140811398187E+109, 0)
Gamma(77.000000) = (1,8854947016662648E+111, 0)
Gamma(78.000000) = (1,451830920283022E+113, 0)
Gamma(79.000000) = (1,1324281178207572E+115, 0)
Gamma(80.000000) = (8,946182130783977E+116, 0)
Gamma(81.000000) = (7,1569457046271725E+118, 0)
Gamma(82.000000) = (5,797126020748008E+120, 0)
Gamma(83.000000) = (4,753643337013366E+122, 0)
Gamma(84.000000) = (3,945523969721089E+124, 0)
Gamma(85.000000) = (3,31424013456571E+126, 0)
Gamma(86.000000) = (2,8171041143808564E+128, 0)
Gamma(87.000000) = (2,4227095383675335E+130, 0)
Gamma(88.000000) = (2,1077572983797526E+132, 0)
Gamma(89.000000) = (1,8548264225741817E+134, 0)
Gamma(90.000000) = (1,650795516091023E+136, 0)
Gamma(91.000000) = (1,4857159644819212E+138, 0)
Gamma(92.000000) = (1,3520015276785438E+140, 0)
Gamma(93.000000) = (1,2438414054642616E+142, 0)
Gamma(94.000000) = (1,1567725070817618E+144, 0)
Gamma(95.000000) = (1,0873661566568553E+146, 0)
Gamma(96.000000) = (1,032997848824012E+148, 0)
Gamma(97.000000) = (9,916779348710516E+149, 0)
Gamma(98.000000) = (9,619275968249195E+151, 0)
Gamma(99.000000) = (9,42689044888421E+153, 0)
Gamma(100.000000) = (9,33262154439535E+155, 0)
Gamma(0.000000 + i0.000000) = (NaN, NaN)
Gamma(1.000000 + i1.000000) = (0,49801566811835923, -0,15494982830180806)
Gamma(2.000000 + i2.000000) = (0,11229424234632254, 0,3236128855019324)
Gamma(3.000000 + i3.000000) = (-0,4401134076370088, -0,0636372431263299)
Gamma(4.000000 + i4.000000) = (0,7058649325913451, -0,49673908399741584)
Gamma(5.000000 + i5.000000) = (-0,9743952418053669, 2,0066898827226805)
Gamma(6.000000 + i6.000000) = (1,0560845455210948, -7,123931816061554)
Gamma(7.000000 + i7.000000) = (-0,26095668519941206, 27,88827411508434)
Gamma(8.000000 + i8.000000) = (1,8442848156317595, -125,96060801752867)
Gamma(9.000000 + i9.000000) = (-94,00399991734474, 643,3621714431141)
Gamma(10.000000 + i10.000000) = (1423,851941789479, -3496,081973308168)
Gamma(11.000000 + i11.000000) = (-16211,00700465313, 18168,810510285286)
Gamma(12.000000 + i12.000000) = (158471,8890918886, -68793,30331463458)
Gamma(13.000000 + i13.000000) = (-1329505,1052081874, -142199,12520863872)
Gamma(14.000000 + i14.000000) = (8576976,67312178, 7218722,503716219)
Gamma(15.000000 + i15.000000) = (-20768001,573587183, -99064686,32101583)
Gamma(16.000000 + i16.000000) = (-490395650,85195, 847486174,9207268)
Gamma(17.000000 + i17.000000) = (9782747798,66319, -2523864726,357996)
Gamma(18.000000 + i18.000000) = (-91408144092,80728, -62548333665,79536)
Gamma(19.000000 + i19.000000) = (80368797570,63837, 1283152922013,1064)
Gamma(20.000000 + i20.000000) = (12322153606702,379, -9813622771583,531)
Gamma(21.000000 + i21.000000) = (-191651224429571,5, -67416801166719,305)
Gamma(22.000000 + i22.000000) = (476610838765573,75, 2709614130691551,5)
Gamma(23.000000 + i23.000000) = (31282423285710508, -23340622982977492)
Gamma(24.000000 + i24.000000) = (-5,062346571412891E+17, -2,8075410996386413E+17)
Gamma(25.000000 + i25.000000) = (-1,1135374386470528E+18, 8,889271476011264E+18)
Gamma(26.000000 + i26.000000) = (1,4103207063357242E+20, -3,1111347801608966E+19)
Gamma(27.000000 + i27.000000) = (-1,20394153792971E+21, -2,100813378567903E+21)
Gamma(28.000000 + i28.000000) = (-2,9772583255514483E+22, 2,9845666640271078E+22)
Gamma(29.000000 + i29.000000) = (6,474322395366405E+23, 4,002110222682179E+23)
Gamma(30.000000 + i30.000000) = (4,982468347052982E+24, -1,3332730971666784E+25)
Gamma(31.000000 + i31.000000) = (-2,701233953257487E+26, -5,3335652308619844E+25)
Gamma(32.000000 + i32.000000) = (-3,5481927258667466E+26, 5,492417944693437E+27)
Gamma(33.000000 + i33.000000) = (1,13499763334942E+29, -3,936060260026878E+27)
Gamma(34.000000 + i34.000000) = (-2,497617380231244E+29, -2,4036706229026682E+30)
Gamma(35.000000 + i35.000000) = (-5,244047476964302E+31, 7,550247174479453E+30)
Gamma(36.000000 + i36.000000) = (1,8326953464053433E+32, 1,1815797667494789E+33)
Gamma(37.000000 + i37.000000) = (2,7496330950708185E+34, -3,7903032739968576E+33)
Gamma(38.000000 + i38.000000) = (-6,153702881069039E+34, -6,593476804299637E+35)
Gamma(39.000000 + i39.000000) = (-1,622188754583588E+37, 3,702408035066972E+35)
Gamma(40.000000 + i40.000000) = (-2,9787072201603815E+37, 4,069596487794187E+38)
Gamma(41.000000 + i41.000000) = (1,0327223226588694E+40, 2,028448840265614E+39)
Gamma(42.000000 + i42.000000) = (9,258591867044077E+40, -2,623834405446389E+41)
Gamma(43.000000 + i43.000000) = (-6,5825696143104E+42, -3,6674407207212435E+42)
Gamma(44.000000 + i44.000000) = (-1,3470021394042014E+44, 1,5970917674147318E+44)
Gamma(45.000000 + i45.000000) = (3,611191294816286E+45, 4,700641025433376E+45)
Gamma(46.000000 + i46.000000) = (1,572050028597055E+47, -6,978410499639601E+46)
Gamma(47.000000 + i47.000000) = (-8,064316457005957E+47, -5,037500759309252E+48)
Gamma(48.000000 + i48.000000) = (-1,5346852327090097E+50, -1,8750204416673013E+49)
Gamma(49.000000 + i49.000000) = (-1,963123049571723E+51, 4,3640542056764855E+51)
Gamma(50.000000 + i50.000000) = (1,112141672863102E+53, 1,0242389193853564E+53)
Gamma(51.000000 + i51.000000) = (4,3124175139534486E+54, -2,2723736598857867E+54)
Gamma(52.000000 + i52.000000) = (-1,9794169465380243E+55, -1,5907071359978385E+56)
Gamma(53.000000 + i53.000000) = (-5,198770735397539E+57, -1,364053254830618E+57)
Gamma(54.000000 + i54.000000) = (-1,120153853109903E+59, 1,4557034579325553E+59)
Gamma(55.000000 + i55.000000) = (3,0502206279504673E+60, 5,6213787579421754E+60)
Gamma(56.000000 + i56.000000) = (2,263901598439912E+62, -1,3869357395249425E+61)
Gamma(57.000000 + i57.000000) = (3,1337184675942048E+63, -7,566798379913671E+63)
Gamma(58.000000 + i58.000000) = (-1,9547779277051327E+65, -2,289059342042218E+65)
Gamma(59.000000 + i59.000000) = (-1,0990589957539216E+67, 2,43674092519427E+66)
Gamma(60.000000 + i60.000000) = (-1,2138821648921022E+68, 4,1070828497360523E+68)
Gamma(61.000000 + i61.000000) = (1,1429060333301634E+70, 1,199617402383565E+70)
Gamma(62.000000 + i62.000000) = (6,356301051435098E+71, -1,4385777048397157E+71)
Gamma(63.000000 + i63.000000) = (8,496122222382356E+72, -2,4629448359526944E+73)
Gamma(64.000000 + i64.000000) = (-6,555666234065589E+74, -8,308825832655404E+74)
Gamma(65.000000 + i65.000000) = (-4,35298630131942E+76, 3,566503620394996E+75)
Gamma(66.000000 + i66.000000) = (-9,093720755985274E+77, 1,5886817012950856E+78)
Gamma(67.000000 + i67.000000) = (3,276704637172729E+79, 7,067540192000104E+79)
Gamma(68.000000 + i68.000000) = (3,3000031716552424E+81, 6,606403155114497E+80)
Gamma(69.000000 + i69.000000) = (1,1027558271573558E+83, -9,8053446601259E+82)
Gamma(70.000000 + i70.000000) = (-4,322638823338879E+83, -6,551055369074648E+84)
Gamma(71.000000 + i71.000000) = (-2,4300476213174805E+86, -1,695901195606269E+86)
Gamma(72.000000 + i72.000000) = (-1,3066738432408389E+88, 3,647419995967721E+87)
Gamma(73.000000 + i73.000000) = (-2,631035116210543E+89, 5,722329713478048E+89)
Gamma(74.000000 + i74.000000) = (1,2168690154152985E+91, 2,7033309001464033E+91)
Gamma(75.000000 + i75.000000) = (1,3482397421660745E+93, 4,28037329878204E+92)
Gamma(76.000000 + i76.000000) = (5,937687180258459E+94, -3,3970662984923895E+94)
Gamma(77.000000 + i77.000000) = (7,877518886977203E+95, -3,258429077036584E+96)
Gamma(78.000000 + i78.000000) = (-8,893580405683355E+97, -1,4068641589865653E+98)
Gamma(79.000000 + i79.000000) = (-8,171120040126758E+99, -1,8182361788970335E+99)
Gamma(80.000000 + i80.000000) = (-3,620672062418723E+101, 2,252383956928435E+101)
Gamma(81.000000 + i81.000000) = (-5,470323981863164E+102, 2,130475015794634E+103)
Gamma(82.000000 + i82.000000) = (5,5003853678128614E+104, 1,0085770500872239E+105)
Gamma(83.000000 + i83.000000) = (5,740274728712029E+106, 1,986132310865518E+106)
Gamma(84.000000 + i84.000000) = (3,002634263016605E+108, -1,245707037705857E+108)
Gamma(85.000000 + i85.000000) = (7,865153067057636E+109, -1,575289481058232E+110)
Gamma(86.000000 + i86.000000) = (-2,2903303039804873E+111, -9,374401780068583E+111)
Gamma(87.000000 + i87.000000) = (-4,291880277570836E+113, -3,196185984412844E+113)
Gamma(88.000000 + i88.000000) = (-2,9995707183640408E+115, 1,1845722184191957E+114)
Gamma(89.000000 + i89.000000) = (-1,2943245318203042E+117, 1,1073023439932073E+117)
Gamma(90.000000 + i90.000000) = (-2,0007723042777158E+118, 9,568042990491017E+118)
Gamma(91.000000 + i91.000000) = (2,4147687416625443E+120, 5,132960890087842E+120)
Gamma(92.000000 + i92.000000) = (2,9270673985442265E+122, 1,5846423875074053E+122)
Gamma(93.000000 + i93.000000) = (1,9583971561168918E+124, -2,5166936071920716E+123)
Gamma(94.000000 + i94.000000) = (8,735240498823332E+125, -7,993043748530299E+125)
Gamma(95.000000 + i95.000000) = (1,6159775690644545E+127, -6,9922194609237035E+127)
Gamma(96.000000 + i96.000000) = (-1,569012926564151E+129, -4,10649302643008E+129)
Gamma(97.000000 + i97.000000) = (-2,2080981691082815E+131, -1,5902953016632018E+131)
Gamma(98.000000 + i98.000000) = (-1,6995881449793942E+133, -8,982249136857376E+131)
Gamma(99.000000 + i99.000000) = (-9,394719018633069E+134, 5,234766160326565E+134)
Gamma(100.000000 + i100.000000) = (-3,3597454530316526E+136, 5,986962556433683E+136)

Factor

! built in
USING: picomath prettyprint ;
0.1 gamma .  ! 9.513507698668723
2.0 gamma .  ! 1.0
10. gamma .  ! 362880.0

Forth

Cristinel Mortici describes this method in Applied Mathematics Letters. "A substantial improvement of the Stirling formula". This algorithm is said to give about 7 good digits, but becomes more inaccurate close to zero. Therefore, a "shift" is performed to move the value returned into the more accurate domain.

8 constant (gamma-shift)

: (mortici)                            ( f1 -- f2)
  -1 s>f f+ 1 s>f
  fover 271828183e-8 f* 12 s>f f* f/
  fover 271828183e-8 f/ f+
  fover f** fswap
  628318530e-8 f* fsqrt f*             \ 2*pi
;

: gamma                                ( f1 -- f2)
  fdup f0< >r fdup f0= r> or abort" Gamma less or equal to zero"
  fdup (gamma-shift) s>f f+ (mortici) fswap
  1 s>f (gamma-shift) 0 do fover i s>f f+ f* loop fswap fdrop f/
;
0.1e gamma f. 9.51348888533932  ok
2e gamma f. 0.999999031674546  ok
10e gamma f. 362879.944850072  ok
70e gamma fe. 171.122444600510E96  ok

This is a word, based on a formula of Ramanujan's famous "lost notebook", which was rediscovered in 1976. His formula contained a constant, which had a value between 1/100 and 1/30. In 2008, E.A. Karatsuba described the function, which determines the value of this constant. Since it contains the gamma function itself, it can't be used in a word calculating the gamma function, so here it is emulated by two symmetrical sigmoidals.

2 constant (gamma-shift)               \ don't change this
                                       \ an approximation of the d(x) function
: ~d(x)                                ( f1 -- f2)
  fdup 10 s>f f<                       \ use first symmetrical sigmoidal
  if                                   \ for range 1-10
    -2705443e-8 fswap 2280802e-6 f/ 1428045e-6 f** 1 s>f f+ f/ 3187831e-8 f+
  else                                 \ use second symmetrical sigmoidal
    -29372563e-9 fswap 1841693e-6 f/ 1052779e-6 f** 1 s>f f+ f/ 3330828e-8 f+
  then 333333333e-10 fover f< if fdrop 1 s>f 30 s>f f/ then
;                                      \ perform some sane clipping to infinity

: (ramanujan)                          ( f1 -- f2)
  fdup fdup f* 4 s>f f*                ( n 4n2)
  fover fover f* fdup f+ f+ fover f+   ( n 8n3+4n2+n)
  fover ~d(x) f+                       ( n 8n3+4n2+n+d[x])
  1 s>f 6 s>f f/ f**                   ( n 8n3+4n2+n+d[x]^1/6)
  fswap fdup 2.7182818284590452353e f/ ( 8n3+4n2+n+d[x]^1/6 n n/e)
  fswap f** f* pi fsqrt f*             ( f)
;

: gamma                                ( f1 -- f2)
  fdup f0< >r fdup f0= r> or abort" Gamma less or equal to zero"
  fdup (gamma-shift) 1- s>f f+ (ramanujan) fswap
  1 s>f (gamma-shift) 0 do fover i s>f f+ f* loop fswap fdrop f/
;
0.1e gamma f. 9.51351721918848  ok
2e gamma f. 0.999999966026125  ok
10e gamma f. 362879.999559333  ok
70e gamma fe. 171.122452428147E96  ok

Fortran

This code shows two methods: Numerical Integration through Simpson formula, and Lanczos approximation. The results of testing are printed altogether with the values given by the function gamma; this function is defined in the Fortran 2008 standard, and supported by GNU Fortran (and other vendors) as extension; if not present in your compiler, you can remove the last part of the print in order to get it compiled with any Fortran 95 compliant compiler.

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

  implicit none

  integer :: i

  write(*, "(3A15)") "Simpson", "Lanczos", "Builtin"
  do i=1, 10
     write(*, "(3F15.8)") my_gamma(i/3.0), lacz_gamma(i/3.0), gamma(i/3.0)
  end do

contains

  pure function intfuncgamma(x, y) result(z)
    real :: z
    real, intent(in) :: x, y
    
    z = x**(y-1.0) * exp(-x)
  end function intfuncgamma


  function my_gamma(a) result(g)
    real :: g
    real, intent(in) :: a

    real, parameter :: small = 1.0e-4
    integer, parameter :: points = 100000

    real :: infty, dx, p, sp(2, points), x
    integer :: i
    logical :: correction

    x = a

    correction = .false.
    ! value with x<1 gives \infty, so we use
    ! \Gamma(x+1) = x\Gamma(x)
    ! to avoid the problem
    if ( x < 1.0 ) then
       correction = .true.
       x = x + 1
    end if

    ! find a "reasonable" infinity...
    ! we compute this integral indeed
    ! \int_0^M dt t^{x-1} e^{-t}
    ! where M is such that M^{x-1} e^{-M} ≤ \epsilon
    infty = 1.0e4
    do while ( intfuncgamma(infty, x) > small )
       infty = infty * 10.0
    end do

    ! using simpson
    dx = infty/real(points)
    sp = 0.0
    forall(i=1:points/2-1) sp(1, 2*i) = intfuncgamma(2.0*(i)*dx, x)
    forall(i=1:points/2) sp(2, 2*i - 1) = intfuncgamma((2.0*(i)-1.0)*dx, x)
    g = (intfuncgamma(0.0, x) + 2.0*sum(sp(1,:)) + 4.0*sum(sp(2,:)) + &
         intfuncgamma(infty, x))*dx/3.0

    if ( correction ) g = g/a

  end function my_gamma

  
  recursive function lacz_gamma(a) result(g)
    real, intent(in) :: a
    real :: g

    real, parameter :: pi = 3.14159265358979324
    integer, parameter :: cg = 7

    ! these precomputed values are taken by the sample code in Wikipedia,
    ! and the sample itself takes them from the GNU Scientific Library
    real, dimension(0:8), parameter :: p = &
         (/ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, &
         771.32342877765313, -176.61502916214059, 12.507343278686905, &
         -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 /)

    real :: t, w, x
    integer :: i

    x = a

    if ( x < 0.5 ) then
       g = pi / ( sin(pi*x) * lacz_gamma(1.0-x) )
    else
       x = x - 1.0
       t = p(0)
       do i=1, cg+2
          t = t + p(i)/(x+real(i))
       end do
       w = x + real(cg) + 0.5
       g = sqrt(2.0*pi) * w**(x+0.5) * exp(-w) * t
    end if
  end function lacz_gamma

end program ComputeGammaInt
Output:
        Simpson        Lanczos        Builtin
     2.65968132     2.67893744     2.67893839
     1.35269761     1.35411859     1.35411787
     1.00000060     1.00000024     1.00000000
     0.88656044     0.89297968     0.89297950
     0.90179849     0.90274525     0.90274531
     0.99999803     1.00000036     1.00000000
     1.19070935     1.19063985     1.19063926
     1.50460517     1.50457609     1.50457561
     2.00000286     2.00000072     2.00000000
     2.77815390     2.77816010     2.77815843

FreeBASIC

Translation of: Java
' FB 1.05.0 Win64

Const pi = 3.1415926535897932
Const e  = 2.7182818284590452

Function gammaStirling (x As Double) As Double
  Return Sqr(2.0 * pi / x) * ((x / e) ^ x)
End Function

Function gammaLanczos (x As Double) As Double
  Dim p(0 To 8) As Double = _ 
  { _
       0.99999999999980993, _ 
     676.5203681218851, _ 
   -1259.1392167224028, _			     	  
     771.32342877765313, _ 
    -176.61502916214059, _ 
      12.507343278686905, _
      -0.13857109526572012, _ 
       9.9843695780195716e-6, _
       1.5056327351493116e-7 _
  }
 
  Dim As Integer g = 7
  If x < 0.5 Then Return pi / (Sin(pi * x) * gammaLanczos(1-x))
  x -= 1
  Dim a As Double = p(0)
  Dim t As Double = x + g + 0.5
  
  For i As Integer = 1 To 8
    a += p(i) / (x + i)
  Next		 
 
  Return Sqr(2.0 * pi) * (t ^ (x + 0.5)) * Exp(-t) * a  
End Function

Print " x", "    Stirling",, "    Lanczos"
Print
For i As Integer = 1 To 20
   Dim As Double d = i / 10.0
   Print   Using "#.##"; d; 
   Print , Using "#.###############"; gammaStirling(d);
   Print , Using "#.###############"; gammaLanczos(d)
Next
Print
Print "Press any key to quit"
Sleep
Output:
 x                Stirling                    Lanczos

0.10          5.697187148977170           9.513507698668738
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.164229713725303
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.908638732853291
1.80          0.889639635287995           0.931383770980243
1.90          0.920842721894229           0.961765831907388
2.00          0.959502175744492           1.000000000000000

Go

package main

import (
    "fmt"
    "math"
)

func main() {
    fmt.Println("    x               math.Gamma                 Lanczos7")
    for _, x := range []float64{-.5, .1, .5, 1, 1.5, 2, 3, 10, 140, 170} {
        fmt.Printf("%5.1f %24.16g %24.16g\n", x, math.Gamma(x), lanczos7(x))
    }
}

func lanczos7(z float64) float64 {
    t := z + 6.5
    x := .99999999999980993 +
        676.5203681218851/z -
        1259.1392167224028/(z+1) +
        771.32342877765313/(z+2) -
        176.61502916214059/(z+3) +
        12.507343278686905/(z+4) -
        .13857109526572012/(z+5) +
        9.9843695780195716e-6/(z+6) +
        1.5056327351493116e-7/(z+7)
    return math.Sqrt2 * math.SqrtPi * math.Pow(t, z-.5) * math.Exp(-t) * x
}
Output:
    x               math.Gamma                 Lanczos7
 -0.5       -3.544907701811032       -3.544907701811087
  0.1        9.513507698668732        9.513507698668752
  0.5        1.772453850905516        1.772453850905517
  1.0                        1                        1
  1.5       0.8862269254527579       0.8862269254527587
  2.0                        1                        1
  3.0                        2                        2
 10.0                   362880        362880.0000000015
140.0    9.61572319694107e+238   9.615723196940201e+238
170.0   4.269068009004746e+304                     +Inf

Groovy

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].reverse()

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

(1..10).each{ printf("%  1.9e\n", gamma(it / 3.0)) }
Output:
  2.678938535e+00
  1.354117939e+00
  1.000000000e+00
  8.929795116e-01
  9.027452930e-01
  1.000000000e+00
  1.190639349e+00
  1.504575488e+00
  2.000000000e+00
  2.778158479e+00

Haskell

Based on HaskellWiki (compatible license):

The Gamma and Beta function as described in 'Numerical Recipes in C++', the approximation is taken from [Lanczos, C. 1964 SIAM Journal on Numerical Analysis, ser. B, vol. 1, pp. 86-96]
cof :: [Double]
cof =
  [ 76.18009172947146
  , -86.50532032941677
  , 24.01409824083091
  , -1.231739572450155
  , 0.001208650973866179
  , -0.000005395239384953
  ]

ser :: Double
ser = 1.000000000190015

gammaln :: Double -> Double
gammaln xx =
  let tmp_ = (xx + 5.5) - (xx + 0.5) * log (xx + 5.5)
      ser_ = ser + sum (zipWith (/) cof [xx + 1 ..])
  in -tmp_ + log (2.5066282746310005 * ser_ / xx)

main :: IO ()
main = mapM_ print $ gammaln <$> [0.1,0.2 .. 1.0]

Or equivalently, as a point-free applicative expression:

import Control.Applicative

cof :: [Double]
cof =
  [ 76.18009172947146
  , -86.50532032941677
  , 24.01409824083091
  , -1.231739572450155
  , 0.001208650973866179
  , -0.000005395239384953
  ]

gammaln :: Double -> Double
gammaln =
  ((+) . negate . (((-) . (5.5 +)) <*> (((*) . (0.5 +)) <*> (log . (5.5 +))))) <*>
  (log .
   ((/) =<<
    (2.5066282746310007 *) .
    (1.000000000190015 +) . sum . zipWith (/) cof . enumFrom . (1 +)))

main :: IO ()
main = mapM_ print $ gammaln <$> [0.1,0.2 .. 1.0]
Output:
2.252712651734255
1.5240638224308496
1.09579799481814
0.7966778177018394
0.572364942924743
0.3982338580692666
0.2608672465316877
0.15205967839984869
6.637623973474716e-2
-4.440892098500626e-16

Icon and Unicon

This works in Unicon. Changing the !10 into (1 to 10) would enable it to work in Icon.

procedure main()
    every write(left(i := !10/10.0,5),gamma(.i))
end

procedure gamma(z)	# Stirling's approximation
    return (2*&pi/z)^0.5 * (z/&e)^z
end
Output:
->gamma
0.1  5.69718714897717
0.2  3.325998424022393
0.3  2.36253003626962
0.4  1.841476335936235
0.5  1.520346901066281
0.6  1.307158857448356
0.7  1.15905329211392
0.8  1.053370968425609
0.9  0.9770615078776954
1.0  0.9221370088957891
->

J

This code shows the built-in method, which works for any value (positive, negative and complex numbers -- but note that negative integer arguments give infinite results).

gamma=: !@<:

<: subtracts one from a number. It's sort of like --lvalue in C, except it always accepts an "rvalue" as an argument (which means it does not modify that argument). And !value finds the factorial of value if value is a positive integer. This illustrates the close relationship between the factorial and gamma functions.

The following direct coding of the task comes from the Stirling's approximation essay on the J wiki:

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

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

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

(Column 1 is the argument, column 2 is the stirling approximation and column 3 uses the builtin support for gamma.)

Java

Implementation of Stirling's approximation and Lanczos approximation.

public class GammaFunction {

	public double st_gamma(double x){
		return Math.sqrt(2*Math.PI/x)*Math.pow((x/Math.E), x);
	}
	
	public double la_gamma(double x){
		double[] p = {0.99999999999980993, 676.5203681218851, -1259.1392167224028,
			     	  771.32342877765313, -176.61502916214059, 12.507343278686905,
			     	  -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
		int g = 7;
		if(x < 0.5) return Math.PI / (Math.sin(Math.PI * x)*la_gamma(1-x));

		x -= 1;
		double a = p[0];
		double t = x+g+0.5;
		for(int i = 1; i < p.length; i++){
			a += p[i]/(x+i);
		}
		
		return Math.sqrt(2*Math.PI)*Math.pow(t, x+0.5)*Math.exp(-t)*a;
	}
	
	public static void main(String[] args) {
		GammaFunction test = new GammaFunction();
		System.out.println("Gamma \t\tStirling \t\tLanczos");
		for(double i = 1; i <= 20; i += 1){
			System.out.println("" + i/10.0 + "\t\t" + test.st_gamma(i/10.0) + "\t" + test.la_gamma(i/10.0));
		}
	}
}
Output:
Gamma 		Stirling 		Lanczos
0.1		5.697187148977169	9.513507698668734
0.2		3.3259984240223925	4.590843711998803
0.3		2.3625300362696198	2.9915689876875904
0.4		1.8414763359362354	2.218159543757687
0.5		1.5203469010662807	1.7724538509055159
0.6		1.307158857448356	1.489192248812818
0.7		1.15905329211392	1.2980553326475577
0.8		1.0533709684256085	1.1642297137253035
0.9		0.9770615078776954	1.0686287021193193
1.0		0.9221370088957891	0.9999999999999998
1.1		0.8834899531687038	0.9513507698668735
1.2		0.8577553353965909	0.9181687423997607
1.3		0.8426782594483921	0.8974706963062777
1.4		0.8367445486370817	0.8872638175030757
1.5		0.8389565525264963	0.8862269254527586
1.6		0.8486932421525738	0.8935153492876909
1.7		0.865621471793884	0.9086387328532916
1.8		0.8896396352879945	0.9313837709802425
1.9		0.9208427218942293	0.9617658319073877
2.0		0.9595021757444916	1.0000000000000002

JavaScript

Implementation of Lanczos approximation.

function gamma(x) {
    var p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
        771.32342877765313, -176.61502916214059, 12.507343278686905,
        -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
    ];

    var g = 7;
    if (x < 0.5) {
        return Math.PI / (Math.sin(Math.PI * x) * gamma(1 - x));
    }

    x -= 1;
    var a = p[0];
    var t = x + g + 0.5;
    for (var i = 1; i < p.length; i++) {
        a += p[i] / (x + i);
    }

    return Math.sqrt(2 * Math.PI) * Math.pow(t, x + 0.5) * Math.exp(-t) * a;
}

jq

Works with: jq version 1.4

Taylor Series

Translation of: Ada
def gamma:
  [
    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
  ] as $a
  | (. - 1) as $y
  | ($a|length) as $n
  | reduce range(2; 1+$n) as $an
      ($a[$n-1]; (. * $y) + $a[$n - $an])
  | 1 / . ;

Lanczos Approximation

# for reals, but easily extended to complex values
def gamma_by_lanczos:
  def pow(x): if x == 0 then 1 elif x == 1 then . else x * log | exp end;
  . as $x
  | ((1|atan) * 4) as $pi
  | if $x < 0.5 then $pi / ((($pi * $x) | sin) * ((1-$x)|gamma_by_lanczos ))
    else
      [   0.99999999999980993, 676.5203681218851,     -1259.1392167224028,
        771.32342877765313,   -176.61502916214059,       12.507343278686905,
         -0.13857109526572012,   9.9843695780195716e-6,   1.5056327351493116e-7] as $p
    | ($x - 1) as $x
    | ($x + 7.5) as $t
    |  reduce range(1; $p|length) as $i
          ($p[0]; . + ($p[$i]/($x + $i) ))
       * ((2*$pi) | sqrt) * ($t | pow($x+0.5)) * ((-$t)|exp) 
    end;

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              tgamma
0.3333333333333333: 2.6789385347077483 : 2.6789385347077483 : 2.678938534707748
0.6666666666666666: 1.3541179394264005 : 1.3541179394263998 : 1.3541179394264005
1                 : 1                  : 0.9999999999999998 : 1
1.3333333333333333: 0.8929795115692493 : 0.8929795115692494 : 0.8929795115692493
1.6666666666666667: 0.9027452929509336 : 0.9027452929509342 : 0.9027452929509336
2                 : 1                  : 1.0000000000000002 : 1
2.3333333333333335: 1.190639348758999  : 1.1906393487589995 : 1.190639348758999
2.6666666666666665: 1.5045754882515399 : 1.5045754882515576 : 1.5045754882515558
3                 : 1.9999999999939684 : 2.0000000000000013 : 2
3.3333333333333335: 2.778158479338573  : 2.778158480437665  : 2.7781584804376647

Jsish

Translation of: Javascript
#!/usr/bin/env jsish
/* Gamma function, in Jsish, using the Lanczos approximation */
function gamma(x) {
    var p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
        771.32342877765313, -176.61502916214059, 12.507343278686905,
        -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
    ];
 
    var g = 7;
    if (x < 0.5) {
        return Math.PI / (Math.sin(Math.PI * x) * gamma(1 - x));
    }
 
    x -= 1;
    var a = p[0];
    var t = x + g + 0.5;
    for (var i = 1; i < p.length; i++) {
        a += p[i] / (x + i);
    }
 
    return Math.sqrt(2 * Math.PI) * Math.pow(t, x + 0.5) * Math.exp(-t) * a;
}

if (Interp.conf('unitTest')) {
    for (var i=-5.5; i <= 5.5; i += 0.5) {
        printf('%2.1f %+e\n', i, gamma(i));
    }
}

/*
=!EXPECTSTART!=
-5.5 +1.091265e-02
-5.0 -4.275508e+13
-4.5 -6.001960e-02
-4.0 +2.672193e+14
-3.5 +2.700882e-01
-3.0 -1.425169e+15
-2.5 -9.453087e-01
-2.0 +6.413263e+15
-1.5 +2.363272e+00
-1.0 -2.565305e+16
-0.5 -3.544908e+00
0.0 +inf
0.5 +1.772454e+00
1.0 +1.000000e+00
1.5 +8.862269e-01
2.0 +1.000000e+00
2.5 +1.329340e+00
3.0 +2.000000e+00
3.5 +3.323351e+00
4.0 +6.000000e+00
4.5 +1.163173e+01
5.0 +2.400000e+01
5.5 +5.234278e+01
=!EXPECTEND!=
*/
Output:
prompt$ jsish --U gammaFunction.jsi
-5.5 +1.091265e-02
-5.0 -4.275508e+13
-4.5 -6.001960e-02
-4.0 +2.672193e+14
-3.5 +2.700882e-01
-3.0 -1.425169e+15
-2.5 -9.453087e-01
-2.0 +6.413263e+15
-1.5 +2.363272e+00
-1.0 -2.565305e+16
-0.5 -3.544908e+00
0.0 +inf
0.5 +1.772454e+00
1.0 +1.000000e+00
1.5 +8.862269e-01
2.0 +1.000000e+00
2.5 +1.329340e+00
3.0 +2.000000e+00
3.5 +3.323351e+00
4.0 +6.000000e+00
4.5 +1.163173e+01
5.0 +2.400000e+01
5.5 +5.234278e+01

prompt$ jsish -u gammaFunction.jsi
[PASS] gammaFunction.jsi

Julia

Works with: Julia version 0.6

Built-in function:

@show gamma(1)

By adaptive Gauss-Kronrod integration:

using QuadGK
gammaquad(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 SpecialFunctions
gamma(1/2) - sqrt(pi)
Output:
2.220446049250313e-16

Koka

Based on OCaml implementation

import std/num/float64

fun gamma-lanczos(x)
  val g = 7.0
  // Coefficients used by the GNU Scientific Library
  val c = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
      771.32342877765313, -176.61502916214059, 12.507343278686905,
      -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]
  fun ag(z: float64, d: int)
    if d == 0 then c[0].unjust + ag(z, 1)
    elif d < 8 then c[d].unjust / (z + d.float64) + ag(z, d.inc)
    else c[d].unjust / (z + d.float64)
  fun gamma(z)
    val z' = z - 1.0
    val p = z' + g + 0.5
    sqrt(2.0 * pi) * pow(p, (z' + 0.5)) * exp(0.0 - p) * ag(z', 0)
  gamma(x)

val e = exp(1.0)
fun gamma-stirling(x)
  sqrt(2.0 * pi / x) * pow(x / e, x)

fun gamma-stirling2(x')
  // Extended Stirling method seen in Abramowitz and Stegun 
  val d = [1.0/12.0, 1.0/288.0, -139.0/51840.0, -571.0/2488320.0]
  fun corr(z, x, n)
    if n < d.length - 1 then d[n].unjust / x + corr(z, x*z, n.inc)
    else d[n].unjust / x
  fun gamma(z)
    gamma-stirling(z)*(1.0 + corr(z, z, 0))
  gamma(x')

fun mirror(gma, z)
  if z > 0.5 then gma(z) else pi / sin(pi * z) / gma(1.0 - z)

fun main()
  println("z\tLanczos\t\t\tStirling\t\tStirling2")
  for(1, 20) fn(i)
    val z = i.float64 / 10.0
    println(z.show(1) ++ "\t" ++ mirror(gamma-lanczos, z).show ++ "\t" ++
        mirror(gamma-stirling, z).show ++ "\t" ++ mirror(gamma-stirling2, z).show)
  for(1, 7) fn(i)
    val z = 10.0 * i.float64
    println(z.show ++ "\t" ++ gamma-lanczos(z).show ++ "\t" ++
        gamma-stirling(z).show ++ "\t" ++ gamma-stirling2(z).show)
Output:
z	Lanczos			Stirling		Stirling2
0.1	9.5135076986687359	10.405084329555955	9.5210418318004439
0.2	4.5908437119988017	5.0739927535371763	4.596862295030256
0.3	2.9915689876875904	3.3503395433773222	2.9984402802949961
0.4	2.218159543757686	2.5270578096699556	2.2277588907113128
0.5	1.7724538509055157	2.0663656770612464	1.7883901437260497
0.6	1.4891922488128184	1.3071588574483559	1.4827753636029286
0.7	1.2980553326475577	1.1590532921139201	1.2950806801024195
0.8	1.1642297137253037	1.0533709684256085	1.1627054102439229
0.9	1.068628702119319	0.97706150787769541	1.0677830813298756
1.0	1.0000000000000002	0.92213700889578909	0.99949946853364036
1.1	0.95135076986687361	0.88348995316870382	0.95103799705518899
1.2	0.91816874239976076	0.85775533539659088	0.91796405783487933
1.3	0.89747069630627774	0.84267825944839203	0.8973312868034562
1.4	0.88726381750307537	0.8367445486370817	0.88716548484542823
1.5	0.88622692545275827	0.83895655252649626	0.88615538430170204
1.6	0.89351534928769061	0.8486932421525738	0.89346184003019224
1.7	0.90863873285329122	0.86562147179388405	0.90859770150945562
1.8	0.93138377098024272	0.8896396352879945	0.93135158986107858
1.9	0.96176583190738729	0.92084272189422933	0.96174006762796482
2.0	1.0000000000000002	0.95950217574449159	0.99997898067003532
10	362880.00000000105	359869.56187410367	362879.99717458693
20	1.2164510040883245e+17	1.2113934233805675e+17	1.2164510037907557e+17
30	8.841761993739658e+30	8.8172365307655063e+30	8.8417619934546387e+30
40	2.0397882081197221e+46	2.0355431612365591e+46	2.0397882081041343e+46
50	6.0828186403425409e+62	6.0726891878763362e+62	6.0828186403274418e+62
60	1.3868311854568534e+80	1.3849063858294502e+80	1.3868311854555093e+80
70	1.7112245242813438e+98	1.7091885781910795e+98	1.711224524280615e+98

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

Lambdatalk

Following Javascript, with Lanczos approximation.

{def gamma.p
 {A.new 0.99999999999980993
      676.5203681218851
    -1259.1392167224028   
      771.32342877765313
     -176.61502916214059
       12.507343278686905   
       -0.13857109526572012
        9.9843695780195716e-6
        1.5056327351493116e-7
}}
-> gamma.p

{def gamma.rec
 {lambda {:x :a :i}
  {if {< :i {A.length {gamma.p}}}
   then {gamma.rec :x 
                   {+ :a {/ {A.get :i {gamma.p}} {+ :x :i}} } 
                   {+ :i 1}}
   else :a
}}} 
-> gamma.rec  

{def gamma
 {lambda {:x}
  {if {< :x 0.5}
   then {/ {PI}
           {* {sin {* {PI} :x}}
              {gamma {- 1 :x}}}}
   else {let { {:x {- :x 1}}
               {:t {+ {- :x 1} 7 0.5}}
             } {* {sqrt {* 2 {PI}}}
                  {pow :t {+ :x 0.5}}
                  {exp -:t}
                  {gamma.rec :x {A.first {gamma.p}} 1}} 
}}}}
-> gamma

{S.map {lambda {:i}
               {div} Γ(:i) = {gamma :i}} 
       {S.serie -5.5 5.5 0.5}}

 Γ(-5.5) = 0.010912654781909836 
 Γ(-5) = -42755084646679.17 
 Γ(-4.5) = -0.06001960130050417 
 Γ(-4) = 267219279041745.34 
 Γ(-3.5) = 0.27008820585226917 
 Γ(-3) = -1425169488222640 
 Γ(-2.5) = -0.9453087204829418 
 Γ(-2) = 6413262697001885 
 Γ(-1.5) = 2.363271801207352 
 Γ(-1) = -25653050788007544 
 Γ(-0.5) = -3.5449077018110295 
 Γ(0) = Infinity 
 Γ(0.5) = 1.7724538509055159 
 Γ(1) = 0.9999999999999998 
 Γ(1.5) = 0.8862269254527586 
 Γ(2) = 1.0000000000000002 
 Γ(2.5) = 1.3293403881791384 
 Γ(3) = 2.000000000000001 
 Γ(3.5) = 3.3233509704478426 
 Γ(4) = 6.000000000000007 
 Γ(4.5) = 11.631728396567446 
 Γ(5) = 23.999999999999996 
 Γ(5.5) = 52.34277778455358

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.042197734555571
function recigamma(z)
  return z + gamma * z^2 + coeff * z^3 + quad * z^4 + qui * z^5 + set * z^6
end

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)
  end
end

M2000 Interpreter

Module PrepareLambdaFunctions {
      Const e = 2.7182818284590452@
      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)=0.99999999999980993@, 676.5203681218851@,   -1259.1392167224028@,  771.32342877765313@
      p(4)=-176.61502916214059@,  12.507343278686905@,  -0.13857109526572012@,  0.0000099843695780195716@
      p(8)=0.00000015056327351493116@
      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 -= 1@
            a=p(0)
            t = x + 7.5@
            For i= 1@ To 8@ {
                  a += p(i) / (x + i)
            }
             = Sqrt(2.0 * pi) * (t ^ (x + 0.5)) * Exp(-t) * a  
      }
      Push gammaStirling, gammaLanczos
}
Call PrepareLambdaFunctions
Read gammaLanczos, gammaStirling
Font "Courier New"
Form 120, 40
document doc$="     χ                        Stirling                     Lanczos"+{
}
Print $(2,20),"x", "Stirling",@(55),"Lanczos", $(0)
Print
For 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 d
Print $("")
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 /Wolfram Language

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 */

МК-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
-	С/П

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

Nim

Translation of: Ada

The algorithm is a translation of that from the Ada solution. We have added a comparison with the gamma function provided by the “math” module from Nim standard library (which is, in fact, the C “tgamma” function).

import math, strformat

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): float =
  let y = x - 1
  result = A[^1]
  for n in countdown(A.high - 1, A.low):
    result = result * y + A[n]
  result = 1 / result

echo "Our gamma function     Nim gamma function      Difference"
echo "------------------     ------------------      ----------"
for i in 1..10:
  let val1 = gamma(i.toFloat / 3)
  let val2 = math.gamma(i.toFloat / 3)
  echo &"{val1:18.16f}     {val2:18.16f}     {val1 - val2:11.4e}"
Output:
Our gamma function     Nim gamma function      Difference
------------------     ------------------      ----------
2.6789385347077483     2.6789385347077479      4.4409e-16
1.3541179394264005     1.3541179394264005      0.0000e+00
1.0000000000000000     1.0000000000000000      0.0000e+00
0.8929795115692493     0.8929795115692493      0.0000e+00
0.9027452929509336     0.9027452929509336      0.0000e+00
1.0000000000000000     1.0000000000000000      0.0000e+00
1.1906393487589990     1.1906393487589990      0.0000e+00
1.5045754882515399     1.5045754882515558     -1.5987e-14
1.9999999999939684     2.0000000000000000     -6.0316e-12
2.7781584793385732     2.7781584804376647     -1.0991e-09

OCaml

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

module Lanczos = struct
  (* Lanczos method *)
  (* Coefficients used by the GNU Scientific Library *)
  let g = 7.
  let c = [|0.99999999999980993; 676.5203681218851; -1259.1392167224028;
	    771.32342877765313; -176.61502916214059; 12.507343278686905;
	    -0.13857109526572012; 9.9843695780195716e-6; 1.5056327351493116e-7|]

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

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

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

end

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

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

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

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

let _ =
  Printf.printf "z\t\tLanczos\t\tStirling\tStirling2\n";
  for i = 1 to 20 do
    let z = float i /. 10. in
    Printf.printf "%-10.8g\t%10.8e\t%10.8e\t%10.8e\n"
    		  z 
		  (mirror Lanczos.gamma z)
		  (mirror Stirling.gamma z)
		  (mirror Stirling2.gamma z)
  done;
  for i = 1 to 7 do
    let z = 10. *. float i in
    Printf.printf "%-10.8g\t%10.8e\t%10.8e\t%10.8e\n"
    		  z
		  (Lanczos.gamma z)
		  (Stirling.gamma z)
		  (Stirling2.gamma z)
  done
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+double(i));
    endfor
    w = x + double(cg) + 0.5;
    g = sqrt(2.0*pi) * w**(x+0.5) * exp(-w) * t;
  endif
endfunction


for i = 1:10
  printf("%f %f\n", gamma(i/3.0), lacz_gamma(i/3.0));
endfor
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 as

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

Pascal

A console application in Free Pascal, created with the Lazarus IDE.

Based on the algorithm for ln(Gamma(x)) (x > 0) in Press et al., Numerical Recipes, 3rd edition, pp. 256-7. For x >= 1/2, we simply take the exponential of their value; for x < 1/2 we calculate Gamma(1 - x) and use the reflection formula Gamma(x)*Gamma(1 - x) = pi/sin(pi*x).

program GammaTest;
{$mode objfpc}{$H+}
uses SysUtils;

function Gamma( x : extended) : extended;
const COF : array [0..14] of extended =
(  0.999999999999997092, // may as well include this in the array
  57.1562356658629235,
 -59.5979603554754912,
  14.1360979747417471,
 -0.491913816097620199,
  0.339946499848118887e-4,
  0.465236289270485756e-4,
 -0.983744753048795646e-4,
  0.158088703224912494e-3,
 -0.210264441724104883e-3,
  0.217439618115212643e-3,
 -0.164318106536763890e-3,
  0.844182239838527433e-4,
 -0.261908384015814087e-4,
  0.368991826595316234e-5);
const
  K = 2.5066282746310005;
  PI_OVER_K = PI / K;
var
  j : integer;
  tmp, w, ser : extended;
  reflect : boolean;
begin
  reflect := (x < 0.5);
  if reflect then w := 1.0 - x else w := x;
  tmp := w + 5.2421875;
  tmp := (w + 0.5)*Ln(tmp) - tmp;
  ser := COF[0];
  for j := 1 to 14 do ser := ser + COF[j]/(w + j);
  try
    if reflect then
      result := PI_OVER_K * w * Exp(-tmp) / (Sin(PI*x) * ser)
    else
      result := K * Exp(tmp) * ser / w;
  except
    raise SysUtils.Exception.CreateFmt(
        'Gamma(%g) is undefined or out of floating-point range', [x]);
  end;
end;

// Main routine for testing the Gamma function
var
    x, k : extended;
begin
  WriteLn( 'Is it seamless at x = 1/2 ?');
  x := 0.49999999999999;
  WriteLn( SysUtils.Format( 'Gamma(%g) = %g', [x, Gamma(x)]));
  x := 0.50000000000001;
  WriteLn( SysUtils.Format( 'Gamma(%g) = %g', [x, Gamma(x)]));
  WriteLn( 'Test a few values:');
  WriteLn( SysUtils.Format( 'Gamma(1)   = %g', [Gamma(1)]));
  WriteLn( SysUtils.Format( 'Gamma(2)   = %g', [Gamma(2)]));
  WriteLn( SysUtils.Format( 'Gamma(3)   = %g', [Gamma(3)]));
  WriteLn( SysUtils.Format( 'Gamma(10)  = %g', [Gamma(10)]));
  WriteLn( SysUtils.Format( 'Gamma(101) = %g', [Gamma(101)]));
  WriteLn( '      100! = 9.33262154439442E157');
  WriteLn( SysUtils.Format( 'Gamma(1/2) = %g', [Gamma(0.5)]));
  WriteLn( SysUtils.Format( 'Sqrt(pi)   = %g', [Sqrt(PI)]));
  WriteLn( SysUtils.Format( 'Gamma(-7/2)       =  %g', [Gamma(-3.5)]));
(*
  Note here a bug or misfeature in Lazarus (doesn't occur in Delphi):
  Putting (16.0/105.0)*Sqrt(PI) does not give the required precision.
  We have to explicitly define the integers as extended floating-point.
*)
  k := extended(16.0)/extended(105.0);
  WriteLn( SysUtils.Format( ' (16/105)Sqrt(pi) =  %g', [k*Sqrt(PI)]));
  WriteLn( SysUtils.Format( 'Gamma(-9/2)       = %g', [Gamma(-4.5)]));
  k := extended(32.0)/extended(945.0);
  WriteLn( SysUtils.Format( '-(32/945)Sqrt(pi) = %g', [-k*Sqrt(PI)]));
end.
Output:
Is it seamless at x = 1/2 ?
Gamma(0.49999999999999) = 1.77245385090555
Gamma(0.50000000000001) = 1.77245385090548
Test a few values:
Gamma(1)   = 1
Gamma(2)   = 1
Gamma(3)   = 2
Gamma(10)  = 362880
Gamma(101) = 9.33262154439452E157
      100! = 9.33262154439442E157
Gamma(1/2) = 1.77245385090552
Sqrt(pi)   = 1.77245385090552
Gamma(-7/2)       =  0.270088205852269
 (16/105)Sqrt(pi) =  0.270088205852269
Gamma(-9/2)       = -0.0600196013005043
-(32/945)Sqrt(pi) = -0.0600196013005042


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 otherwise
my $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

Phix

with javascript_semantics
sequence c = repeat(0,12)
function spouge_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/z
end function

function taylor_gamma(atom x)
-- (good for values between 0 and 1, apparently)
    constant t = { 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 }
    atom y = x-1,
         s = t[$]
    for n=length(t)-1 to 1 by -1 do
        s = s*y + t[n]
    end for
    return 1/s
end function

function lanczos_gamma(atom z)
    if z<0.5 then
        return PI / (sin(PI*z)*lanczos_gamma(1-z))
    end if
    -- use a lanczos approximation:
    atom x = 0.99999999999980993,
         t = z + 6.5;
    sequence p = { 676.5203681218851,   
                 -1259.1392167224028, 
                   771.32342877765313,  
                  -176.61502916214059,  
                    12.507343278686905, 
                    -0.13857109526572012,
                     9.9843695780195716e-6,
                     1.5056327351493116e-7 }
    z -= 1
    for i=1 to length(p) do
        x += p[i] / (z + i)
    end for
    return sqrt(2*PI) * power(t,z+0.5) * exp(-t) * x
end function
 
constant sqPI = sqrt(PI)

procedure sq(sequence zm, string fmt="%19.16f")
    atom {z, mul} = zm
    atom e = sqPI/mul
    sequence s = {spouge_gamma(z),
                  taylor_gamma(z),
                  lanczos_gamma(z)},
         error = sq_abs(sq_sub(s,e))
    string t = join(s,",  ",fmt:=fmt)&", "
    integer bdx = smallest(error,return_index:=true)
    atom best = s[bdx],
         p = s[bdx]*mul
    for i=1 to length(s) do -- (potentially mark >1)
        if s[i]=best then
            t[i*22-2..i*22-1] = "*,"
        end if
    end for
    string es = sprintf(fmt,e)
    printf(1,"%5g: %s %s,  %19.16f\n",{z,t,es,p*p})
end procedure

printf(1,"   z    ------ spouge -----   ----- taylor ------   ----- lanczos -----   ---- expected -----  %19.16f\n",PI)
papply({{-3/2,3/4},{-1/2,-1/2},{1/2,1},{1,sqPI},{3/2,2},{2,sqPI},{5/2,4/3},{3,sqPI/2},{7/2,8/15},{4,sqPI/6}},sq)
sq({0.001,sqPI/999.4237725},"%19.15f")
sq({0.01,sqPI/99.43258512},"%19.16f")
sq({0.1,sqPI/9.513507699},"%19.16f")
sq({10,sqPI/362880},"%19.12f")
sq({100,sqPI/9.332621544e155},"%19.13g")
if machine_bits()=64 then
    sq({150,sqPI/3.808922638e260},"%19.13g") -- (fatal power overflow error on 32 bits)
end if
Output:

The closest to the expected result for each z (row) is marked with a trailing asterisk.
The final column is the value of PI (to 16dp) we would get from that best/starred result.

   z    ------ spouge -----   ----- taylor ------   ----- lanczos -----   ---- expected -----   3.1415926535897932
 -1.5:  2.3632718012073547*,  2.3632718095606211,   2.3632718012073532,   2.3632718012073547,   3.1415926535897932
 -0.5: -3.5449077018110320*, -3.5449077018110306,  -3.5449077018110308,  -3.5449077018110321,   3.1415926535897932
  0.5:  1.7724538509055158,   1.7724538509055160*,  1.7724538509055166,   1.7724538509055160,   3.1415926535897932
    1:  0.9999999999999998,   1.0000000000000000*,  1.0000000000000002,   1.0000000000000000,   3.1415926535897932
  1.5:  0.8862269254527577,   0.8862269254527580*,  0.8862269254527583,   0.8862269254527580,   3.1415926535897932
    2:  0.9999999999999994,   1.0000000000000000*,  1.0000000000000005,   1.0000000000000000,   3.1415926535897932
  2.5:  1.3293403881791359,   1.3293403881791365*,  1.3293403881791379,   1.3293403881791370,   3.1415926535897906
    3:  1.9999999999999978,   1.9999999999939679,   2.0000000000000016*,  2.0000000000000000,   3.1415926535897981
  3.5:  3.3233509704478376,   3.3233509583896768,   3.3233509704478456*,  3.3233509704478426,   3.1415926535897990
    4:  5.9999999999999884,   5.9999914100724727,   6.0000000000000063*,  6.0000000000000000,   3.1415926535897998
0.001: 999.423772484595421*, 999.423772484595404,  999.423772484595254,  999.423772500000000,   3.1415926534929476
 0.01: 99.4325851191505990,  99.4325851191506035*, 99.4325851191505828,  99.4325851200000000,   3.1415926535361195
  0.1:  9.5135076986687313,   9.5135076986687318*,  9.5135076986687300,   9.5135076990000000,   3.1415926533710076
   10: 362879.999999996094,       0.000000029163,  362880.000000000725*, 362880.000000000000,   3.1415926535898058
  100: 9.332621544394e+155,   7.510232292979e-39,  9.332621544394e+155*,    9.332621544e+155,   3.1415926538549222
  150:  3.80892263763e+260*,  5.128102530869e-44,   3.80892263763e+260,     3.808922638e+260,   3.1415926529801383

mpfr version

Above translated to mpfr, but spouge only since there's not much point transferring inherent inaccuracies in taylor/lanczos constants, and compared against the builtin.

Library: Phix/mpfr
without javascript_semantics -- (no mpfr_exp(), mpfr_gamma() in pwa/p2js)
constant dp = 30
string fmt = "%5s:  %33s,  %33s,  %32s\n"
requires("1.0.2") -- (mpfr_get_fixed(maxlen), mpfr_gamma)
include mpfr.e
mpfr_set_default_precision(-87) -- 87 decimal places. 
sequence mc = mpfr_inits(40)
 
function mpfr_spouge_gamma(mpfr z)
    mpfr accm = mc[1]
    if mpfr_cmp_si(accm,0)=0 then
        -- mc[1] := sqrt(2*PI)
        mpfr_const_pi(accm)
        mpfr_mul_si(accm,accm,2)
        mpfr_sqrt(accm,accm)
        -- k1_factrl = (k - 1)!*(-1)^k with 0!==1
        mpfr k1_factrl = mpfr_init(1),
             tmk = mpfr_init(),
             p = mpfr_init()
        for k=2 to length(mc) do
            -- mc[k] = exp(13-k)*power(13-k,k-1.5)/k1_factrl
            mpfr_set_si(tmk,length(mc)+1-k)
            mpfr_exp(mc[k],tmk)
            mpfr_set_d(p,k-1.5)
            mpfr_pow(p,tmk,p)
            mpfr_div(p,p,k1_factrl)
            mpfr_mul(mc[k],mc[k],p)
            -- k1_factrl *= -(k-1)
            mpfr_mul_si(k1_factrl,k1_factrl,-(k-1))
        end for
    end if
    accm = mpfr_init_set(accm)
    for k=2 to length(mc) do
        -- accm += mc[k]/(z+k-1)
        mpfr ck = mpfr_init_set(mc[k]),
             zk = mpfr_init_set(z)
        mpfr_add_si(zk,zk,k-1)
        mpfr_div(ck,ck,zk)
        mpfr_add(accm,accm,ck)
    end for
    -- atom zc = z+length(mc)
    -- accm *= exp(-zc)*power(zc,z+0.5) -- Gamma(z+1)
    mpfr p = mpfr_init_set(z),
         ez = mpfr_init(),
         zh = mpfr_init(0.5)
    mpfr_add_si(p,p,length(mc))
    mpfr_neg(ez,p)
    mpfr_exp(ez,ez)
    mpfr_add(zh,zh,z)
    mpfr_pow(p,p,zh)
    mpfr_mul(accm,accm,ez)
    mpfr_mul(accm,accm,p)
    -- return accm/z
    mpfr_div(accm,accm,z)
    return accm
end function

constant mPI = mpfr_init(),
        mqPI = mpfr_init()
mpfr_const_pi(mPI)
string pistr = mpfr_get_fixed(mPI,dp)
mpfr_sqrt(mqPI,mPI)

function makempfr(object x)
    mpfr res
    if string(x) then
        x = split(x,'/')
        res = mpfr_init(x[1])
        mpfr_div_si(res,res,to_integer(x[2]))
    elsif sequence(x) then
        {mpfr x1, object x2} = x
        res = mpfr_init_set(x1)
        if string(x2) then
            mpfr d = mpfr_init(x2)
            mpfr_div(res,res,d)
        else
            mpfr_div_d(res,res,x2)
        end if
    else
        res = mpfr_init(x)
    end if
    return res
end function

procedure mq(sequence zm, integer d=dp)
    mpfr {z, mul} = apply(zm,makempfr)
    mpfr s = mpfr_spouge_gamma(z)
    string t = mpfr_get_fixed(s,d,10,maxlen:=dp+2)
    mpfr e = mpfr_init()
    mpfr_gamma(e,z)
    mpfr p = mpfr_init_set(s)
    mpfr_mul(p,p,mul)
    mpfr_mul(p,p,p)
    string zs = mpfr_get_fixed(z,3),
           es = mpfr_get_fixed(e,d,10,maxlen:=dp+2),
           ps = mpfr_get_fixed(p,dp,10)
    printf(1,fmt,{zs,t,es,ps})
end procedure

printf(1,"   z     %s    %s   %s\n",{pad(" spouge ",dp+2,"BOTH",'-'),pad(" expected ",dp+2,"BOTH",'-'),pistr})
papply({{-1.5,0.75},{-0.5,-0.5},{0.5,1},{1,{mqPI,1}},{1.5,2},{2,{mqPI,1}},{2.5,"4/3"},{3,{mqPI,2}},{3.5,"8/15"},{4,{mqPI,6}}},mq)
mq({"1/1000",{mqPI,"999.4237724845954661149822012996"}},28)
mq({"1/100",{mqPI,"99.43258511915060371353298887051"}},29)
mq({"1/10",{mqPI,"9.513507698668731836292487177265"}},30)
mq({10,{mqPI,362880}},25)
mq({100,{mqPI,"9.332621544394415268169923885627e155"}},0)
mq({150,{mqPI,"3.80892263763056972698595524350735e260"}},0)
Output:
   z     ------------ spouge ------------    ----------- expected -----------   3.141592653589793238462643383279
 -1.5:   2.363271801207354703064223311121,   2.363271801207354703064223311121,  3.141592653589793238462643383279
 -0.5:   -3.544907701811032054596334966e0,   -3.544907701811032054596334966e0,  3.141592653589793238462643383279
  0.5:   1.772453850905516027298167483341,   1.772453850905516027298167483341,  3.141592653589793238462643383279
    1:   1.000000000000000000000000000000,                                  1,  3.141592653589793238462643383279
  1.5:   0.886226925452758013649083741670,   0.886226925452758013649083741670,  3.141592653589793238462643383279
    2:   1.000000000000000000000000000000,                                  1,  3.141592653589793238462643383279
  2.5:   1.329340388179137020473625612505,   1.329340388179137020473625612505,  3.141592653589793238462643383279
    3:   2.000000000000000000000000000000,                                  2,  3.141592653589793238462643383279
  3.5:   3.323350970447842551184064031264,   3.323350970447842551184064031264,  3.141592653589793238462643383279
    4:   6.000000000000000000000000000000,                                  6,  3.141592653589793238462643383279
0.001:   999.4237724845954661149822012996,   999.4237724845954661149822012996,  3.141592653589793238462643383279
0.009:   99.43258511915060371353298887051,   99.43258511915060371353298887051,  3.141592653589793238462643383279
  0.1:   9.513507698668731836292487177265,   9.513507698668731836292487177265,  3.141592653589793238462643383279
   10:   362880.0000000000000000000000000,                             362880,  3.141592653589793238462643383279
  100:   9.33262154439441526816992388e155,   9.33262154439441526816992388e155,  3.141592653589793238462643383279
  150:   3.80892263763056972698595524e260,   3.80892263763056972698595524e260,  3.141592653589793238462643383279

Phixmonti

0.577215664901 var gamma
-0.65587807152056 var coeff
-0.042002635033944 var quad
0.16653861138228 var qui
-0.042197734555571 var theSet

def recigamma var z /# n -- n #/
    z 6 power theSet *
    z 5 power qui *
    z 4 power quad *
    z 3 power coeff *
    z 2 power gamma *
    z + + + + +
enddef

/# without var
def recigamma
    dup 6 power theSet * swap
    dup 5 power qui * swap
    dup 4 power quad * swap
    dup 3 power coeff * swap
    dup 2 power gamma * swap
    + + + + +
enddef
#/

def gammafunc   /# n -- n #/
    dup 1 == if
    else
        dup abs 0.5 <= if
            recigamma 1 swap /
        else
            dup 1 - gammafunc swap 1 - *
        endif
    endif
enddef

0.1 2.1 .1 3 tolist
for
    dup print " = " print gammafunc print nl
endfor

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

Prolog

This version matches Wolfram Alpha to within a few digits at the end, so the last few digits are a bit off. There's an early check to stop evaluating coefficients once the desired accuracy is reached. Removing this check does not improve accuracy vs. Wolfram Alpha.

gamma_coefficients(
     [ 1.00000000000000000000000,  0.57721566490153286060651, -0.65587807152025388107701,
      -0.04200263503409523552900,  0.16653861138229148950170, -0.04219773455554433674820,
      -0.00962197152787697356211,  0.00721894324666309954239, -0.00116516759185906511211,
      -0.00021524167411495097281,  0.00012805028238811618615, -0.00002013485478078823865,
      -0.00000125049348214267065,  0.00000113302723198169588, -0.00000020563384169776071,
       0.00000000611609510448141,  0.00000000500200764446922, -0.00000000118127457048702,
       0.00000000010434267116911,  0.00000000000778226343990, -0.00000000000369680561864,
       0.00000000000051003702874, -0.00000000000002058326053, -0.00000000000000534812253,
       0.00000000000000122677862, -0.00000000000000011812593,  0.00000000000000000118669,
       0.00000000000000000141238, -0.00000000000000000022987,  0.00000000000000000001714
    ]).

tolerance(1e-17).

gamma(X, _) :- X =< 0.0, !, fail.
gamma(X, Y) :-
    X < 1.0, small_gamma(X, Y), !.
gamma(1, 1) :- !.
gamma(1.0, 1) :- !.
gamma(X, Y) :-
    X1 is X - 1,
    gamma(X1, Y1),
    Y is X1 * Y1.
    
small_gamma(X, Y) :-
    gamma_coefficients(Cs),
    recip_gamma(X, 1.0, Cs, 1.0, 0.0, Y0),
    Y is 1 / Y0.

recip_gamma(_, _, [], _, Y, Y) :- !.
recip_gamma(_, _, [], X0, X1, Y) :- tolerance(Tol), abs(X1 - X0) < Tol, Y = X1, !. % early exit
recip_gamma(X, PrevPow, [C|Cs], _, X1, Y) :-
    Power is PrevPow * X,
    X2 is X1 + C*Power,
    recip_gamma(X, Power, Cs, X1, X2, Y).
Output:
% see how close gamma(0.5) is to the square root of pi.
?- gamma(0.5,X), Y is sqrt(pi), Err is abs(X - Y).
X = 1.772453850905516,
Y = 1.7724538509055159,
Err = 2.220446049250313e-16.

?- gamma(1.5,X).
X = 0.886226925452758.

?- gamma(4.9,X).
X = 20.667385961857857.

?- gamma(5,X).
X = 24.

?- gamma(5.01,X).
X = 24.364473447872836.

?- gamma(6.9,X).
X = 597.4941281573107.

?- gamma(6.95,X).
X = 655.7662628554252.

?- gamma(7,X).
X = 720.

% 100!
?- gamma(101.0,X).
X = 9.33262154439441e+157.

% Note when passed integer, gamma(101) returns full big int precision
?- gamma(101,X).
X = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000.

?- gamma(100.98,X).
X = 8.510619261391532e+157.

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.0
A(1) = 0.5772156649015328606
A(2) =-0.6558780715202538811
A(3) =-0.0420026350340952355
A(4) = 0.1665386113822914895
A(5) =-0.0421977345555443368 ; was ...33675
A(6) =-0.0096219715278769736
A(7) = 0.0072189432466630995
A(8) =-0.0011651675918590651
A(9) =-0.0002152416741149510
A(10) = 0.0001280502823881162
A(11) =-0.0000201348547807882
A(12) =-0.0000012504934821427
A(13) = 0.0000011330272319817
A(14) =-0.0000002056338416978
A(15) = 0.0000000061160951045
A(16) = 0.0000000050020076445
A(17) =-0.0000000011812745705
A(18) = 0.0000000001043426712
A(19) = 0.0000000000077822634
A(20) =-0.0000000000036968056
A(21) = 0.0000000000005100370
A(22) =-0.0000000000000205833
A(23) =-0.0000000000000053481
A(24) = 0.0000000000000012268
A(25) =-0.0000000000000001181
A(26) = 0.0000000000000000012
A(27) = 0.0000000000000000014
A(28) =-0.0000000000000000002
;A(29) = 0.00000000000000000002
Protected y.d, Prod.d, Sum.d, N
If x<=0.0: ProcedureReturn 0.0: EndIf ; Error
y = x-1.0: Prod = 1.0
While y>=1.0
  Prod*y: y-1.0 ; Recurse using Gamma(x+1) = x*Gamma(x)
Wend
Sum= A(28)
For N = 27 To 0 Step -1
  Sum*y+A(N)
Next N
If Sum=0.0: ProcedureReturn Infinity(): EndIf
ProcedureReturn Prod / Sum
EndProcedure

Procedure.d GammLn(x.d) ; AKJ  01-May-10
; Returns Ln(Gamma()) or 0 on error
; Based on Numerical Recipes gamma.h
Protected j, tmp.d, y.d, ser.d
Protected Dim cof.d(13)
cof(0) = 57.1562356658629235
cof(1) = -59.5979603554754912
cof(2) = 14.1360979747417471
cof(3) = -0.491913816097620199
cof(4) = 0.339946499848118887e-4
cof(5) = 0.465236289270485756e-4
cof(6) = -0.983744753048795646e-4
cof(7) = 0.158088703224912494e-3
cof(8) = -0.210264441724104883e-3
cof(9) = 0.217439618115212643e-3
cof(10) = -0.164318106536763890e-3
cof(11) = 0.844182239838527433e-4
cof(12) = -0.261908384015814087e-4
cof(13) = 0.368991826595316234e-5
If x<=0: ProcedureReturn 0: EndIf ; Bad argument
y = x
tmp = x+5.2421875
tmp = (x+0.5)*Log(tmp)-tmp
ser = 0.999999999999997092
For j=0 To 13
  y + 1: ser + cof(j)/y
Next j
ProcedureReturn 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 i
Debug ""
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

Procedural

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

Functional

In terms of fold/reduce:

Works with: Python version 3.7
'''Gamma function'''

from functools import reduce


# gamma_ :: [Float] -> Float -> Float
def gamma_(tbl):
    '''Gamma function.'''
    def go(x):
        y = float(x) - 1.0
        return 1.0 / reduce(
            lambda a, x: a * y + x,
            tbl[-2::-1],
            tbl[-1]
        )
    return lambda x: go(x)


# TBL :: [Float]
TBL = [
    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
]


# TEST ----------------------------------------------------
# main :: IO()
def main():
    '''Gamma function over a range of values.'''

    gamma = gamma_(TBL)
    print(
        fTable(' i -> gamma(i/3):\n')(repr)(lambda x: "%0.7e" % x)(
            lambda x: gamma(x / 3.0)
        )(enumFromTo(1)(10))
    )


# GENERIC -------------------------------------------------

# enumFromTo :: (Int, Int) -> [Int]
def enumFromTo(m):
    '''Integer enumeration from m to n.'''
    return lambda n: list(range(m, 1 + n))


# FORMATTING -------------------------------------------------

# fTable :: String -> (a -> String) ->
#                     (b -> String) -> (a -> b) -> [a] -> String
def fTable(s):
    '''Heading -> x display function -> fx display function ->
                     f -> xs -> tabular string.
    '''
    def go(xShow, fxShow, f, xs):
        ys = [xShow(x) for x in xs]
        w = max(map(len, ys))
        return s + '\n' + '\n'.join(map(
            lambda x, y: y.rjust(w, ' ') + ' -> ' + fxShow(f(x)),
            xs, ys
        ))
    return lambda xShow: lambda fxShow: lambda f: lambda xs: go(
        xShow, fxShow, f, xs
    )


# MAIN ---
if __name__ == '__main__':
    main()
Output:
 i -> gamma(i/3):

 1 -> 2.6789385e+00
 2 -> 1.3541179e+00
 3 -> 1.0000000e+00
 4 -> 8.9297951e-01
 5 -> 9.0274529e-01
 6 -> 1.0000000e+00
 7 -> 1.1906393e+00
 8 -> 1.5045755e+00
 9 -> 2.0000000e+00
10 -> 2.7781585e+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]
         for (i in 1:8) {
           x <- x+p[i+1]/(z+i)
         }
         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)))
   }
}

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

Raku

(formerly 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

REXX

Taylor series, 80-digit coefficients

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

Note:   The Taylor series isn't much good above values of   . Already on modest values of x (say > 3) you loose precision. See below for a solution.

/*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:   Γ(½) = =

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(½).

Spouge's approximation, using 87 digit coefficients

Translation of: Phix
Translation of: C

This REXX version is a translation of   Phix   but with more (decimal digits) precision and more steps.

Many of the "normal" high-level mathematical functions aren't available in REXX, so some of them (RYO) are included here.

/*REXX program calculates the gamma function using Spouge's approximation with 87 digits*/
e=2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138
numeric digits length(e)   -  length(.)          /*use the number of decimal digits in E*/
c.=  0
# = 40                                           /*#:  the number of steps in GAMMA func*/
                    call sq gamma(-3/2),  3/4
                    call sq gamma(-1/2), -1/2
                    call sq gamma( 1/2),   1
                    call si gamma(  1 )
                    call sq gamma( 3/2),   2
                    call si gamma(  2 )
                    call sq gamma( 5/2),  4/3
                    call si gamma(  3 )
                    call sq gamma( 7/2),  8/15
                    call si gamma(  4 )
exit                                             /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
gamma: procedure expose c. e #; parse arg z;         #p= # + 1
       accm = c.1
       if accm==0  then do;  accm= sqrt( 2*pi() )
                             c.1 = accm
                             kfact = 1
                                         do k=2  to #
                                         c.k= exp(#p-k) * pow(#p-k, k-1.5) / kfact
                                         kfact = kfact  *  -(k-1)
                                         end   /*k*/
                        end

           do j=2  to #;   accm = accm   +   c.j / (z+j-1)
           end   /*k*/

       return (accm * exp(-(z+#)) * pow(z+#, z+0.5) ) / z
/*──────────────────────────────────────────────────────────────────────────────────────*/
pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348
fmt:    parse arg n,p,a;  _= format(n,p,a);  L= length(_);      return left( strip0(_), L)
isInt:  return datatype(arg(1), 'W')                      /*is the argument an integer? */
sq:     procedure expose #; parse arg x,mu; say fmt(x,9,#)  fmt((x*mu)**2,9,#);   return
si:     procedure expose #; parse arg x;    say fmt(x,9,#);                       return
strip0: procedure; arg _; if pos(., _)\==0  then _= strip(strip(_,'T',0),'T',.);  return _
/*──────────────────────────────────────────────────────────────────────────────────────*/
exp: procedure expose e; arg x; ix= x%1; if abs(x-ix)>.5  then ix=ix+sign(x); x= x-ix; z=1
     _=1;  w=1;    do j=1;  _= _*x/j;    z= (z+_)/1;      if z==w  then leave;         w=z
                   end  /*j*/;           if z\==0  then z= e**ix * z;             return z
/*──────────────────────────────────────────────────────────────────────────────────────*/
ln:     procedure; parse arg x; call e; ig= x>1.5; is= 1-2*(ig\==1); ii= 0; xx= x
          do while ig & xx>1.5 | \ig & xx<.5; _=e
        do k=-1; iz=xx*_**-is; if k>=0&(ig&iz<1|\ig&iz>.5)  then leave; _=_*_; izz=iz; end
        xx= izz; ii= ii+is*2**k;   end   /*while*/;      x= x*e**-ii-1;  z=0;  _= -1;  p=z
          do k=1; _=-_*x;  z=z+_/k;  if z=p  then leave;  p=z; end;  /*k*/;    return z+ii
/*──────────────────────────────────────────────────────────────────────────────────────*/
pow:    procedure; parse arg x,y;  if y=0  then return 1;  if x=0  then return 0
        if isInt(y)  then return x**y;          if isInt(1/y)  then return root(x, 1/y)
        if abs(y//1)=.5  then return sqrt(x)**sign(y)*x**(y%1);     return exp( y*ln(x) )
/*──────────────────────────────────────────────────────────────────────────────────────*/
root:   procedure; parse arg x 1 ox,y 1 oy;     if x=0 | y=1  then return x/1
        if \isInt(y)  then return $pow(x, 1/y)
        if y==2  then return sqrt(x); if y==-2  then return 1/sqrt(x); return rooti(x,y)/1
/*──────────────────────────────────────────────────────────────────────────────────────*/
rooti:  x=abs(x); y=abs(y); a= digits() + 5;  m= y-1;  d= 5
        parse value format(x,2,1,,0) 'E0'  with  ? 'E' _ .;   g= (?/y'E'_ % y) + (x>1)
          do until d==a;   d=min(d+d, a);  numeric digits d;  o=0
            do until o=g;  o=g;  g= format((m*g**y+x)/y/g**m,,d-2);  end;  end
        _= g*sign(ox);  if oy<0  then _= 1/_;                                     return _
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x;  if x=0  then return 0;  d=digits();  numeric digits;  h=d+6
      numeric form; m.=9; parse value format(x,2,1,,0) 'E0' with g "E" _ .; g=g *.5'e'_ %2
            do j=0  while h>9;        m.j=h;                 h=h%2+1;          end  /*j*/
            do k=j+5  to 0  by -1;    numeric digits m.k;    g=(g+x/g)*.5;     end  /*k*/
      numeric digits d;     return g/1
output   when using the default input:
        2.3632718012073547030642233111215269103967         3.1415926535897932384626433832795028841972
       -3.5449077018110320545963349666822903655951         3.1415926535897932384626433832795028841972
        1.7724538509055160272981674833411451827975         3.1415926535897932384626433832795028841972
        1
        0.8862269254527580136490837416705725913988         3.1415926535897932384626433832795028841972
        1
        1.3293403881791370204736256125058588870982         3.1415926535897932384626433832795028841972
        2
        3.3233509704478425511840640312646472177454         3.1415926535897932384626433832795028841972
        6

This program also has some minor issues. The number of steps, as well as e and pi, are hardcoded. This limits the precision to about 80 digits. See below for a solution.

A generic Gamma function

Libraries: How to use
Library: Numbers
Library: Functions
Library: Constants
Library: Settings
Library: Abend

Below program calculates the Gamma function for any (real) element x (except integers < 1) in the specified precision (but over 100 digits it quickly becomes slow).
Improvements over above programs:

Closed formulas added for all half integers and positive integers (much faster than series). Argument reduction (mapping to interval 0.5...1.5) added to the Lanczos solution. Results are now also for larger values of x accurate to about 60 digits. Dynamic determination of digits and iterations added to the Spouge solution. It's now slightly faster on low x values and gives correct results using > 87 digits. Depending on the parameters, the program selects the optimal method for calculating Gamma.

include Settings

say version; say 'Gamma'; say
arg n; if n = '' then n = 100; numeric digits n
say '(Half)integers formulas'
w = '-99.5 -10.5 -5.5 -2.5 -1.5 -0.5 0.5 1 1.5 2 2.5 5 5.5 10 10.5 99 99.5'
numeric digits n
do i = 1 to Words(w)
   x = Word(w,i); call Time('r'); r = Gamma(x); e = Format(Time('e'),,3)
   say 'Formulas' Format(x,4,1) r '('e 'seconds)'
end
say
say 'Lanczos (max 60 decimals) vs Spouge (no limit) vs Stirling (no limit) approximation'
w = '-12.8 -6.4 -3.2 -1.6 -0.8 -0.4 -0.2 -0.1 0.1 0.2 0.4 0.8 1.6 3.2 6.4 12.8'
do i = 1 to Words(w)
   x = Word(w,i)
   numeric digits Min(60,n)
   call Time('r'); r = Gamma(x); e = Format(Time('e'),,3)
   say 'Lanczos ' Format(x,4,1) r '('e 'seconds)'
   numeric digits n
   call Time('r'); r = Gamma(x); e = Format(Time('e'),,3)
   say 'Spouge  ' Format(x,4,1) r '('e 'seconds)'
   if x > 0 then do
      call Time('r'); r = Stirling(x); e = Format(Time('e'),,3)
      say 'Stirling' Format(x,4,1) r '('e 'seconds)'
   end
end
say
say 'Same for a bigger number'
w = '-99.9 99.9'
do i = 1 to Words(w)
   x = Word(w,i)
   numeric digits Min(60,n)
   call Time('r'); r = Gamma(x); e = Format(Time('e'),,3)
   say 'Lanczos ' Format(x,4,1) r '('e 'seconds)'
   numeric digits n
   call Time('r'); r = Gamma(x); e = Format(Time('e'),,3)
   say 'Spouge  ' Format(x,4,1) r '('e 'seconds)'
   if x > 0 then do
      call Time('r'); r = Stirling(x); e = Format(Time('e'),,3)
      say 'Stirling' Format(x,4,1) r '('e 'seconds)'
   end
end
exit

Gamma:
/* Gamma */
procedure expose glob. fact.
arg x
/* Formulas for negative and positive (half)integers */
if x < 0 then do
   if Half(x) then do
      numeric digits Digits()+2
      i = Abs(Floor(x)); y = (-1)**i*2**(2*i)*Fact(i)*Sqrt(Pi())/Fact(2*i)
      numeric digits Digits()-2
      return y+0
   end
end
if x > 0 then do
   if Whole(x) then
      return Fact(x-1)
   if Half(x) then do
      numeric digits Digits()+2
      i = Floor(x); y = Fact(2*i)*Sqrt(Pi())/(2**(2*i)*Fact(i))
      numeric digits Digits()-2
      return y+0
   end
end
p = Digits()
if p < 61 then do
/* Lanczos with predefined coefficients */
/* Map negative x to positive x */
   if x < 0 then
      return Pi()/(Gamma(1-x)*Sin(Pi()*x))
/* Argument reduction to interval (0.5,1.5) */
   numeric digits p+2
   c = Trunc(x); x = x-c
   if x < 0.5 then do
      x = x+1; c = c-1
   end
/* Series coefficients 1/Gamma(x) in 80 digits Fransen & Wrigge */
    c.1 =  1.00000000000000000000000000000000000000000000000000000000000000000000000000000000
    c.2 =  0.57721566490153286060651209008240243104215933593992359880576723488486772677766467
    c.3 = -0.65587807152025388107701951514539048127976638047858434729236244568387083835372210
    c.4 = -0.04200263503409523552900393487542981871139450040110609352206581297618009687597599
    c.5 =  0.16653861138229148950170079510210523571778150224717434057046890317899386605647425
    c.6 = -0.04219773455554433674820830128918739130165268418982248637691887327545901118558900
    c.7 = -0.00962197152787697356211492167234819897536294225211300210513886262731167351446074
    c.8 =  0.00721894324666309954239501034044657270990480088023831800109478117362259497415854
    c.9 = -0.00116516759185906511211397108401838866680933379538405744340750527562002584816653
   c.10 = -0.00021524167411495097281572996305364780647824192337833875035026748908563946371678
   c.11 =  0.00012805028238811618615319862632816432339489209969367721490054583804120355204347
   c.12 = -0.00002013485478078823865568939142102181838229483329797911526116267090822918618897
   c.13 = -0.00000125049348214267065734535947383309224232265562115395981534992315749121245561
   c.14 =  0.00000113302723198169588237412962033074494332400483862107565429550539546040842730
   c.15 = -0.00000020563384169776071034501541300205728365125790262933794534683172533245680371
   c.16 =  0.00000000611609510448141581786249868285534286727586571971232086732402927723507435
   c.17 =  0.00000000500200764446922293005566504805999130304461274249448171895337887737472132
   c.18 = -0.00000000118127457048702014458812656543650557773875950493258759096189263169643391
   c.19 =  0.00000000010434267116911005104915403323122501914007098231258121210871073927347588
   c.20 =  0.00000000000778226343990507125404993731136077722606808618139293881943550732692987
   c.21 = -0.00000000000369680561864220570818781587808576623657096345136099513648454655443000
   c.22 =  0.00000000000051003702874544759790154813228632318027268860697076321173501048565735
   c.23 = -0.00000000000002058326053566506783222429544855237419746091080810147188058196444349
   c.24 = -0.00000000000000534812253942301798237001731872793994898971547812068211168095493211
   c.25 =  0.00000000000000122677862823826079015889384662242242816545575045632136601135999606
   c.26 = -0.00000000000000011812593016974587695137645868422978312115572918048478798375081233
   c.27 =  0.00000000000000000118669225475160033257977724292867407108849407966482711074006109
   c.28 =  0.00000000000000000141238065531803178155580394756670903708635075033452562564122263
   c.29 = -0.00000000000000000022987456844353702065924785806336992602845059314190367014889830
   c.30 =  0.00000000000000000001714406321927337433383963370267257066812656062517433174649858
   c.31 =  0.00000000000000000000013373517304936931148647813951222680228750594717618947898583
   c.32 = -0.00000000000000000000020542335517666727893250253513557337960820379352387364127301
   c.33 =  0.00000000000000000000002736030048607999844831509904330982014865311695836363370165
   c.34 = -0.00000000000000000000000173235644591051663905742845156477979906974910879499841377
   c.35 = -0.00000000000000000000000002360619024499287287343450735427531007926413552145370486
   c.36 =  0.00000000000000000000000001864982941717294430718413161878666898945868429073668232
   c.37 = -0.00000000000000000000000000221809562420719720439971691362686037973177950067567580
   c.38 =  0.00000000000000000000000000012977819749479936688244144863305941656194998646391332
   c.39 =  0.00000000000000000000000000000118069747496652840622274541550997151855968463784158
   c.40 = -0.00000000000000000000000000000112458434927708809029365467426143951211941179558301
   c.41 =  0.00000000000000000000000000000012770851751408662039902066777511246477487720656005
   c.42 = -0.00000000000000000000000000000000739145116961514082346128933010855282371056899245
   c.43 =  0.00000000000000000000000000000000001134750257554215760954165259469306393008612196
   c.44 =  0.00000000000000000000000000000000004639134641058722029944804907952228463057968680
   c.45 = -0.00000000000000000000000000000000000534733681843919887507741819670989332090488591
   c.46 =  0.00000000000000000000000000000000000032079959236133526228612372790827943910901464
   c.47 = -0.00000000000000000000000000000000000000444582973655075688210159035212464363740144
   c.48 = -0.00000000000000000000000000000000000000131117451888198871290105849438992219023663
   c.49 =  0.00000000000000000000000000000000000000016470333525438138868182593279063941453996
   c.50 = -0.00000000000000000000000000000000000000001056233178503581218600561071538285049997
   c.51 =  0.00000000000000000000000000000000000000000026784429826430494783549630718908519485
   c.52 =  0.00000000000000000000000000000000000000000002424715494851782689673032938370921241
/* Series expansion */
   x = x-1; s = 0
   do k = 52 by -1 to 1
      s = s*x+c.k
   end
   y = 1/s
/* Undo reduction */
   if c = -1 then
      y = y/x
   else do
      do i = 1 to c
         y = (x+i)*y
      end
   end
end
else do
   x = x-1
/* Spouge */
/* Estimate digits and iterations */
   q = Floor(p*1.5); a = Floor(p*1.3)
   numeric digits q
/* Series */
   s = 0
   do k = 1 to a-1
      s = s+((-1)**(k-1)*Power(a-k,k-0.5)*Exp(a-k))/(Fact(k-1)*(x+k))
   end
   s = s+Sqrt(2*Pi()); y = Power(x+a,x+0.5)*Exp(-a-x)*s
end
/* Normalize */
numeric digits p
return y+0

Stirling:
/* Sterling */
procedure expose glob. fact.
arg x
return Sqrt(2*Pi()/x) * Power(x/e(),x)

include Constants
include Functions
include Numbers
include Abend
Output
:
REXX-Regina_3.9.6(MT) 5.00 29 Apr 2024
Gamma function in arbitrary precision

(Half)integers formulas (100 digits)
Formulas  -99.5 3.370459273906717035419140191178165368212858243169804822385940657239537725368041400224226278392502618E-157 (0.003 seconds)
Formulas  -10.5 -2.640121820547716316246385325311240439682468432522587656059168154777653141232089673307782521306919765E-7 (0.000 seconds)
Formulas   -5.5 0.01091265478190986298673234429377905644050439299584730891829569009625650045347052040480970101310888490 (0.001 seconds)
Formulas   -2.5 -0.9453087204829418812256893244486107641586930432652731350473641545882193517818838300666403502605571549 (0.000 seconds)
Formulas   -1.5 2.363271801207354703064223311121526910396732608163182837618410386470548379454709575166600875651392887 (0.000 seconds)
Formulas   -0.5 -3.544907701811032054596334966682290365595098912244774256427615579705822569182064362749901313477089331 (0.000 seconds)
Formulas    0.5 1.772453850905516027298167483341145182797549456122387128213807789852911284591032181374950656738544665 (0.000 seconds)
Formulas    1.0 1 (0.000 seconds)
Formulas    1.5 0.8862269254527580136490837416705725913987747280611935641069038949264556422955160906874753283692723327 (0.001 seconds)
Formulas    2.0 1 (0.000 seconds)
Formulas    2.5 1.329340388179137020473625612505858887098162092091790346160355842389683463443274136031212992553908499 (0.000 seconds)
Formulas    5.0 24 (0.000 seconds)
Formulas    5.5 52.34277778455352018114900849241819367949013237611424488006401129409378637307891910622901158181014715 (0.000 seconds)
Formulas   10.0 362880 (0.000 seconds)
Formulas   10.5 1133278.388948785567334574165588892475560298308275159776608723414529483390056004153717630538727607291 (0.000 seconds)
Formulas   99.0 9.426890448883247745626185743057242473809693764078951663494238777294707070023223798882976159207729113E+153 (0.001 seconds)
Formulas   99.5 9.367802114655996591305637999137598430056780942876744878408158310966911153528806387157595583300106562E+154 (0.001 seconds)

Lanczos (60 digits) vs Spouge (61 digits) vs Stirling (61 digits) approximations
Lanczos   -12.8 -1.44241809348812392704183496823417171228485786984071442421769E-9 (0.001 seconds)
Spouge    -12.8 -1.442418093488123927041834968234171712284857869840714424217876E-9 (0.061 seconds)
Lanczos    -6.4 -0.00214311842968855616971089012135323947575979728834356196087007 (0.001 seconds)
Spouge     -6.4 -0.002143118429688556169710890121353239475759797288343561960870050 (0.062 seconds)
Lanczos    -3.2 0.689056412005979742919224040168358601528149721171394266681422 (0.001 seconds)
Spouge     -3.2 0.6890564120059797429192240401683586015281497211713942666814006 (0.061 seconds)
Lanczos    -1.6 2.31058285808092523235318127282049942788600677267861414816871 (0.001 seconds)
Spouge     -1.6 2.310582858080925232353181272820499427886006772678614148168727 (0.063 seconds)
Lanczos    -0.8 -5.73855463999850381650594784491144000429263749786675377223601 (0.001 seconds)
Spouge     -0.8 -5.738554639998503816505947844911440004292637497866753772236066 (0.063 seconds)
Lanczos    -0.4 -3.72298062203204275598583347080335570330149759689981183834669 (0.001 seconds)
Spouge     -0.4 -3.722980622032042755985833470803355703301497596899811838346699 (0.062 seconds)
Lanczos    -0.2 -5.82114856862651686818160469134229346570980884445593876492447 (0.001 seconds)
Spouge     -0.2 -5.821148568626516868181604691342293465709808844455938764924472 (0.062 seconds)
Lanczos    -0.1 -10.6862870211931935489730533569448077816983878506097317904937 (0.001 seconds)
Spouge     -0.1 -10.68628702119319354897305335694480778169838785060973179049371 (0.066 seconds)
Lanczos     0.1 9.51350769866873183629248717726540219255057862608837734305000 (0.000 seconds)
Spouge      0.1 9.513507698668731836292487177265402192550578626088377343050001 (0.063 seconds)
Stirling    0.1 5.697187148977169278607259516105973271247103273209125657453591 (0.003 seconds)
Lanczos     0.2 4.59084371199880305320475827592915200343410999829340301778885 (0.000 seconds)
Spouge      0.2 4.590843711998803053204758275929152003434109998293403017788853 (0.062 seconds)
Stirling    0.2 3.325998424022392556831252338044446631141936051050822874561565 (0.003 seconds)
Lanczos     0.4 2.21815954375768822305905402190767945077056650177146958224198 (0.000 seconds)
Spouge      0.4 2.218159543757688223059054021907679450770566501771469582241978 (0.063 seconds)
Stirling    0.4 1.841476335936235407774504215721792671787348602355998395147546 (0.003 seconds)
Lanczos     0.8 1.16422971372530337363632093826845869314196176889118775298489 (0.000 seconds)
Spouge      0.8 1.164229713725303373636320938268458693141961768891187752984894 (0.059 seconds)
Stirling    0.8 1.053370968425608533997094555812234627247042644411690511716880 (0.003 seconds)
Lanczos     1.6 0.893515349287690261436600032992805368792359423255954841203208 (0.000 seconds)
Spouge      1.6 0.8935153492876902614366000329928053687923594232559548412032077 (0.064 seconds)
Stirling    1.6 0.8486932421525737351870671884548739160413838607490704857265386 (0.003 seconds)
Lanczos     3.2 2.42396547993536801209211236969059225781321007909891679339251 (0.000 seconds)
Spouge      3.2 2.423965479935368012092112369690592257813210079098916793392514 (0.069 seconds)
Stirling    3.2 2.361851203186240411720280147420856619690317833030951977155556 (0.003 seconds)
Lanczos     6.4 240.833779983445938793193921422181728753410453111118999895588 (0.000 seconds)
Spouge      6.4 240.8337799834459387931939214221817287534104531111189998955875 (0.060 seconds)
Stirling    6.4 237.7207525902404787072183649425603774456984237805478743551168 (0.002 seconds)
Lanczos    12.8 289487660.334241583580309266192984495005040185968050409506654 (0.000 seconds)
Spouge     12.8 289487660.3342415835803092661929844950050401859680504095066538 (0.059 seconds)
Stirling   12.8 287609477.0875058075487675889481380625378498188858381901124308 (0.003 seconds)

Same for a bigger number
Lanczos   -99.9 1.72726520939328005075554286550876635176440408500391962571907E-157 (0.000 seconds)
Spouge    -99.9 1.727265209393280050755542865508766351764404085003870332156550733790282578981487471194267631617799447E-157 (0.274 seconds)
Lanczos    99.9 5.89173215164436165675854187059394609460949390669692020855752E+155 (0.000 seconds)
Spouge     99.9 5.891732151644361656758541870593946094609493906696920208557519945524581875362460658552892944423816715E+155 (0.266 seconds)
Stirling   99.9 5.886819525828908144161960074977041612792659661747365945060453973672990129231211036627906746201731112E+155 (0.010 seconds)

Ring

decimals(3)
gamma = 0.577
coeff = -0.655
quad = -0.042
qui = 0.166
set = -0.042

for i=1 to 10
    see gammafunc(i / 3.0) + nl
next

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

, the Gamma function;
, the regularized Gamma function which is also known as the normalized incomplete Gamma function;
, which the GSL calls the complementary normalized Gamma function;
;
;
.

RPL

Translation of: Ada
≪ 1 -
  { 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 }
  → y a
  ≪ a DUP SIZE GET
     a SIZE 1 - 1 FOR n
        y * a n GET + 
     -1 STEP
     INV
≫ ≫ 'GAMMA' STO
.3 GAMMA

The built-in FACT instruction is obviously based on a similar Taylor formula, since it returns same results:

.3 1 - FACT
Output:
2: 2.99156898769
1: 2.99156898769

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

Rust

Stirling

Translation of: AWK
use std::f64::consts;

fn main() {
    let gamma = |x: f64| { assert_ne!(x, 0.0); (2.0*consts::PI/x).sqrt() * (x * (x/consts::E).ln()).exp()};
    (1..=20).for_each(|x| { 
        let x = f64::from(x) / 10.0; 
        println!("{:.02} => {:.10}", x, gamma(x));
    });
}
Output:
0.10 => 5.6971871490
0.20 => 3.3259984240
0.30 => 2.3625300363
0.40 => 1.8414763359
0.50 => 1.5203469011
0.60 => 1.3071588574
0.70 => 1.1590532921
0.80 => 1.0533709684
0.90 => 0.9770615079
1.00 => 0.9221370089
1.10 => 0.8834899532
1.20 => 0.8577553354
1.30 => 0.8426782594
1.40 => 0.8367445486
1.50 => 0.8389565525
1.60 => 0.8486932422
1.70 => 0.8656214718
1.80 => 0.8896396353
1.90 => 0.9208427219
2.00 => 0.9595021757

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))
    end
endfunction

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.e
define τ = 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 Python and mpmath (see below). 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.772156649015328606065121e-1,
-6.558780715202538810770195e-1,
-4.200263503409523552900393e-2,
 1.665386113822914895017008e-1,
-4.219773455554433674820830e-2,
-9.621971527876973562114922e-3,
 7.218943246663099542395010e-3,
-1.165167591859065112113971e-3,
-2.152416741149509728157300e-4,
 1.280502823881161861531986e-4,
-2.013485478078823865568939e-5,
-1.250493482142670657345359e-6,
 1.133027231981695882374130e-6,
-2.056338416977607103450154e-7,
 6.116095104481415817862499e-9,
 5.002007644469222930055665e-9,
-1.181274570487020144588127e-9,
 1.04342671169110051049154e-10,
 7.782263439905071254049937e-12,
-3.696805618642205708187816e-12,
 5.100370287454475979015481e-13,
-2.05832605356650678322243e-14,
-5.348122539423017982370017e-15,
 1.226778628238260790158894e-15,
-1.181259301697458769513765e-16,
 1.186692254751600332579777e-18,
 1.412380655318031781555804e-18,
-2.298745684435370206592479e-19,
 1.714406321927337433383963e-20

function gamma_(x_) {
    external _gamma_coef
    x = x_
    y = 1
    while (x<0.5) y = y/x++
    while (x>1.5) y = --x*y
    z = _gamma_coef[30]
    x--
    for (i=29; 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::1000)/100
u=map(&gamma(),x)
v=map(&gamma_(),x)
max(abs((v-u):/u))
end

Output

9.80341e-15

Here follows the Python program to compute coefficients.

from mpmath import mp

mp.dps = 50

def gamma_coef(n):
    a = [mp.mpf(1), mp.mpf(mp.euler)]
    for k in range(3, n + 1):
        s = sum((-1)**j * mp.zeta(j) * a[k - j - 1] for j in range(2, k))
        a.append((s - a[1] * a[k - 2]) / (1 - k * a[0]))
    return a

def horner(a, x):
    y = 0
    for s in reversed(a):
        y = y * x + s
    return y

gc = gamma_coef(30)

def gamma_approx(x):
    y = mp.mpf(1)
    while x < 0.5:
        y /= x
        x += 1
    while x > 1.5:
        x -= 1
        y *= x
    return y / horner(gc, x - 1)

for x in gc:
    print(mp.nstr(x, 25))

Tcl

Works with: Tcl version 8.5
Library: Tcllib (Package: math)
Library: Tcllib (Package: math::calculus)
package require math
package require math::calculus

# gamma(1) and gamma(1.5)

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

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

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

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

    # method 2 - library function

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

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

    if {$x == entier($x)} {
        puts [format {%4.1f %13.6f %13.6f %13.6f} $x $gamma $libgamma $f]
        set f [expr $f * $x]
    } else {
        puts [format {%4.1f %13.6f %13.6f %13.6f} $x $gamma $libgamma $f2]
        set f2 [expr $f2 * $x]
    }
}
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→X
X^(X-1/2)e^(-X)√(2π)→Y
Disp X,(X-1)!,Y
Pause
End
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→X
e^(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)->Y
Disp X,(X-1)!,Y
Pause
End
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            

TXR

Taylor Series

Translation of: Ada

Separator commas in numeric tokens are supported only as of TXR 283.

(defun gamma (x)
  (/ (rpoly (- x 1.0)
            #( 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))))

(each ((i 1..11))
  (put-line (pic "##.######" (gamma (/ i 3.0)))))
Output:
 2.678939
 1.354118
 1.000000
 0.892980
 0.902745
 1.000000
 1.190639
 1.504575
 2.000000
 2.778158

Stirling

(defun gamma (x)
  (* (sqrt (/ (* 2 %pi%)
              x))
     (expt (/ x %e%) x)))

(each ((i 1..11))
  (put-line (pic "##.######" (gamma (/ i 3.0)))))
Output:
 2.156976
 1.202851
 0.922137
 0.839743
 0.859190
 0.959502
 1.149106
 1.458490
 1.945403
 2.709764

Lanczos

Translation of: Haskell

The Haskell version calculates the natural log of the gamma function, which is why the function is called gammaln; we correct that here by calling exp:

(defun gamma (x)
  (let* ((cof #(76.18009172947146 -86.50532032941677
                24.01409824083091 -1.231739572450155
                0.001208650973866179 -0.000005395239384953))
         (ser0 1.000000000190015)
         (x55 (+ x 5.5))
         (tmp (- x55 (* (+ x 0.5) (log x55))))
         (ser (+ ser0 (sum [mapcar / cof (succ x)]))))
    (exp (- (log (/ (* 2.5066282746310005 ser) x)) tmp))))

(each ((i (rlist 0.1..1.0..0.1 2..10)))
  (put-line (pic "##.# ######.######" i (gamma i))))
Output:
 0.1      9.513508
 0.2      4.590844
 0.3      2.991569
 0.4      2.218160
 0.5      1.772454
 0.6      1.489192
 0.7      1.298055
 0.8      1.164230
 0.9      1.068629
 1.0      1.000000
 2.0      1.000000
 3.0      2.000000
 4.0      6.000000
 5.0     24.000000
 6.0    120.000000
 7.0    720.000000
 8.0   5040.000000
 9.0  40320.000000
10.0 362880.000000

From Wikipedia Python code. Output is identical to above.

(defun gamma (x)
  (if (< x 0.5)
    (/ %pi%
       (* (sin (* %pi% x))
          (gamma (- 1 x))))
    (let* ((cof #(676.5203681218851 -1259.1392167224028
                  771.32342877765313 -176.61502916214059
                  12.507343278686905 -0.13857109526572012
                  9.9843695780195716e-6 1.5056327351493116e-7))
           (ser0 0.99999999999980993)
           (z (pred x))
           (tmp (+ z (len cof) -0.5))
           (ser (+ ser0 (sum [mapcar / cof (succ z)]))))
      (* (sqrt (* 2 %pi%))
         (expt tmp (+ z 0.5))
         (exp (- tmp))
         ser))))

(each ((i (rlist 0.1..1.0..0.1 2..10)))
  (put-line (pic "##.# ######.######" i (gamma i))))

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 lanczos
CLOSE DATABASES ALL
CLEAR
CREATE CURSOR results (ZVal B(1), GamVal B(15))
INDEX ON zval TAG ZVal COLLATE "Machine"
SET ORDER TO 0
o = CREATEOBJECT("lanczos")
FOR i = 1 TO 20
x = i/10
    INSERT INTO results VALUES (x, o.Gamma(x))
ENDFOR
UPDATE results SET GamVal = ROUND(GamVal, 0) WHERE ZVal = INT(ZVal)
*!* This just creates the output text - it is not part of the algorithm
DO 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 Init
DODEFAULT()
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 Double
IF 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	
ENDIF
RETURN lg
ENDFUNC	

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

V (Vlang)

Translation of: go
import math
fn main() {
    println("    x               math.Gamma                 Lanczos7")
    for x in [-.5, .1, .5, 1, 1.5, 2, 3, 10, 140, 170] {
        println("${x:5.1f} ${math.gamma(x):24.16} ${lanczos7(x):24.16}")
    }
}
 
fn lanczos7(z f64) f64 {
    t := z + 6.5
    x := .99999999999980993 +
        676.5203681218851/z -
        1259.1392167224028/(z+1) +
        771.32342877765313/(z+2) -
        176.61502916214059/(z+3) +
        12.507343278686905/(z+4) -
        .13857109526572012/(z+5) +
        9.9843695780195716e-6/(z+6) +
        1.5056327351493116e-7/(z+7)
    return math.sqrt2 * math.sqrt_pi * math.pow(t, z-.5) * math.exp(-t) * x
}
Output:
    x               math.Gamma                 Lanczos7
 -0.5       -3.544907701811032       -3.544907701811087
  0.1        9.513507698668732        9.513507698668752
  0.5        1.772453850905516        1.772453850905517
  1.0                        1                        1
  1.5       0.8862269254527579       0.8862269254527587
  2.0                        1                        1
  3.0                        2                        2
 10.0                   362880        362880.0000000015
140.0    9.61572319694107e+238   9.615723196940201e+238
170.0   4.269068009004746e+304                     +Inf

Wren

Translation of: Kotlin
Library: Wren-fmt
Library: Wren-math

The gamma method in the Math class is based on the Lanczos approximation.

import "./fmt" for Fmt
import "./math" for Math

var stirling = Fn.new { |x| (2 * Num.pi / x).sqrt * (x / Math.e).pow(x) }

System.print(" x\tStirling\t\tLanczos\n")
for (i in 1..20) {
    var d = i / 10
    Fmt.print("$4.2f\t$16.14f\t$16.14f", d, stirling.call(d), Math.gamma(d))
}
Output:
 x	Stirling		Lanczos

0.10	5.69718714897717	9.51350769866875
0.20	3.32599842402239	4.59084371199881
0.30	2.36253003626962	2.99156898768759
0.40	1.84147633593624	2.21815954375769
0.50	1.52034690106628	1.77245385090552
0.60	1.30715885744836	1.48919224881282
0.70	1.15905329211392	1.29805533264756
0.80	1.05337096842561	1.16422971372530
0.90	0.97706150787770	1.06862870211932
1.00	0.92213700889579	1.00000000000000
1.10	0.88348995316870	0.95135076986687
1.20	0.85775533539659	0.91816874239976
1.30	0.84267825944839	0.89747069630628
1.40	0.83674454863708	0.88726381750308
1.50	0.83895655252650	0.88622692545276
1.60	0.84869324215257	0.89351534928769
1.70	0.86562147179388	0.90863873285329
1.80	0.88963963528799	0.93138377098024
1.90	0.92084272189423	0.96176583190739
2.00	0.95950217574449	1.00000000000000

XPL0

Translation of: Ada
function real Gamma (X);
real X, A, Y, Sum;
integer N;
begin
   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   := X - 1.0;
   Sum := A (29);
   for N:= 29-1 downto 0 do
      Sum := Sum * Y + A (N);
   return 1.0 / Sum;
end \Gamma\;

\Test program:
integer I;
begin
   Format(0, 14);
   for I:= 1 to 10 do
      [RlOut(0, Gamma (Float (I) / 3.0));  CrLf(0)];
end
Output:
 2.67893853470775E+000
 1.35411793942640E+000
 1.00000000000000E+000
 8.92979511569249E-001
 9.02745292950934E-001
 1.00000000000000E+000
 1.19063934875900E+000
 1.50457548825154E+000
 1.99999999999397E+000
 2.77815847933857E+000

Yabasic

Translation of: Phix
dim c(12)
 
sub gamma(z)
    local accm, k, k1_factrl
    
    accm = c(1)
    if accm=0 then
        accm = sqrt(2*PI)
        c(1) = accm
        k1_factrl = 1 
        for k=2 to 12
            c(k) = exp(13-k)*(13-k)^(k-1.5)/k1_factrl
            k1_factrl = k1_factrl * -(k-1)
        next
    end if
    for k=2 to 12
        accm = accm + c(k)/(z+k-1)
    next
    accm = accm * exp(-(z+12))*(z+12)^(z+0.5)
    return accm/z
end sub
 
sub si(x)
    print x using "%18.13f"
end sub


for i = 0.1 to 2.1 step .1
    print i, " = "; : si(gamma(i))
next

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