sidef/scripts/lanczos_approximation.sf

#!/usr/bin/ruby

#
## Algorithm from Wikipedia: https://en.wikipedia.org/wiki/Lanczos_approximation#Simple_implementation
#

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
    const pi = Num.pi

    if (z.real < 0.5) {
        result = (pi / (sin(pi * z) * gamma(1 - z)))
    }
    else {
        z--
        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.real : result
}

for i in [1/2,4,5,6,30,40,50] {
    say "gamma(%3s) =~ %s"%(i, gamma(i));
}