数学

JavaScriptでガンマ関数

コメントする

ここの下の方にあるランチョス近似のコード（Python）を JavaScript にしてみた。

/*
 * Γ(x) = (x - 1)! = Gamma(x)
 */
var Gamma = (function()
{
    var p = [
        676.5203681218851,
        -1259.1392167224028,
        771.32342877765313,
        -176.61502916214059,
        12.507343278686905,
        -0.13857109526572012,
        9.9843695780195716e-6,
        1.5056327351493116e-7
    ];
    var g = 7;

    function Gamma( z )
    {
        if( z < 0.5 ){
            return Math.PI/( Math.sin(Math.PI*z)*Gamma(1-z) );
        }

        z --;
        var x = 0.99999999999980993;

        for( var i=p.length; i--; ){
            x += p[i] / (z + i + 1);
        }

        var t = z + g + 0.5;
        var G = Math.sqrt(Math.PI*2) *Math.exp( (z+0.5)*Math.log(t)-t ) *x;

        if( Math.floor(z) === z ){
            return Math.round( G );
        } else {
            return G;
        }
    }

    return Gamma;
})();

少し書き換えて対数ガンマ関数を計算できるようにしたものも置いておきます。

/*
 * ln Γ(x) = GammaLN(x)
 */
var GammaLN = (function()
{
    var p = [
        676.5203681218851,
        -1259.1392167224028,
        771.32342877765313,
        -176.61502916214059,
        12.507343278686905,
        -0.13857109526572012,
        9.9843695780195716e-6,
        1.5056327351493116e-7
    ];
    var g     = 7;
    var LN_PI = Math.log( Math.PI );

    function GammaLN( z )
    {
        if( z < 0.5 ){
            return LN_PI -Math.log(Math.sin( Math.PI*z )) -GammaLN(1-z);
        }

        z --;
        var x = 0.99999999999980993;

        for( var i=p.length; i--; ){
            x += p[i] / (z + i + 1);
        }

        var t = z + g + 0.5;
        return ( LN_PI+Math.LN2 )/2 +Math.log(x) +(z+0.5)*Math.log(t) -t;
    }

    return GammaLN;
})();

ついでに逆数を計算できるようにしたものも。

/*
 * 1/Γ(x) = GammaInv(x)
 */
var GammaInv = (function()
{
    var p = [
        676.5203681218851,
        -1259.1392167224028,
        771.32342877765313,
        -176.61502916214059,
        12.507343278686905,
        -0.13857109526572012,
        9.9843695780195716e-6,
        1.5056327351493116e-7
    ];
    var g = 7;

    function GammaInv( z )
    {
        if( z < 0.5 ){
            return Math.sin(Math.PI*z) /( Math.PI*GammaInv(1-z) );
        }

        z --;
        var x = 0.99999999999980993;

        for( var i=p.length; i--; ){
            x += p[i] / (z + i + 1);
        }

        var t = z + g + 0.5;
        return Math.exp( t-(z+0.5)*Math.log(t) ) /Math.sqrt(Math.PI*2) /x;
    }

    return GammaInv;
})();