gammaln.js

// 定义级数系数
const COEFFICIENTS = [
    0.99999999999999709182, 57.156235665862923517, -59.597960355475491248,
    14.136097974741747174, -0.49191381609762019978, 0.33994649984811888699e-4,
    0.46523628927048575665e-4, -0.98374475304879564677e-4,
    0.15808870322491248884e-3, -0.21026444172410488319e-3,
    0.2174396181152126432e-3, -0.16431810653676389022e-3,
    0.84418223983852743293e-4, -0.2619083840158140867e-4,
    0.36899182659531622704e-5
];

const g = 607 / 128;
const LOGSQRT2PI = Math.log(Math.sqrt(2 * Math.PI));

/**
 * 使用Lanczos近似法计算[伽马函数](https://en.wikipedia.org/wiki/Gamma_function)的对数值。
 * 该函数接受任何大于0的实数n作为输入。
 * 该函数对于普通伽马函数无法处理的大数值n（n > 165）非常有用。
 * 代码基于Lanczos的伽马近似法，定义见[此处](http://my.fit.edu/~gabdo/gamma.txt)。
 *
 * @param {number} n 任何大于零的实数。
 * @returns {number} 输入值的伽马函数的对数值。
 *
 * @example
 * gammaln(500); // 2605.1158503617335
 * gammaln(2.4); // 0.21685932244884043
 */
function gammaln(n) {
    // 如果值不在定义域内，返回无穷大
    if (n <= 0) {
        return Number.POSITIVE_INFINITY;
    }

    // 递减n，因为近似法是为n - 1定义的
    n--;

    // 创建级数近似
    let a = COEFFICIENTS[0];

    for (let i = 1; i < 15; i++) {
        a += COEFFICIENTS[i] / (n + i);
    }

    const tmp = g + 0.5 + n;

    // 返回伽马函数的自然对数
    return LOGSQRT2PI + Math.log(a) - tmp + (n + 0.5) * Math.log(tmp);
}

export default gammaln;