Gamma function

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