#!/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!=
*/