gamma.js

import factorial from "./factorial.js";

/**
 * 使用Nemes近似法计算值的[伽马函数](https://en.wikipedia.org/wiki/Gamma_function)。
 * n的伽马函数等价于(n-1)!，但与阶乘函数不同的是，伽马函数对所有实数n（除了零和负整数）都有定义（此时返回NaN）。
 * 注意，伽马函数也对复数有定义，但此实现目前不处理复数作为输入值。
 * Nemes近似法在[这里](https://arxiv.org/abs/1003.6020)定义为定理2.2。
 * 负值使用[欧拉反射公式](https://en.wikipedia.org/wiki/Gamma_function#Properties)进行计算。
 *
 * @param {number} n 任何实数，除了零和负整数。
 * @returns {number} 输入值的伽马函数值。
 *
 * @example
 * gamma(11.5); // 11899423.084037038
 * gamma(-11.5); // 2.29575810481609e-8
 * gamma(5); // 24
 */
function gamma(n) {
    if (Number.isInteger(n)) {
        if (n <= 0) {
            // 伽马函数对零或负整数无定义
            return Number.NaN;
        } else {
            // 对整数输入使用阶乘
            return factorial(n - 1);
        }
    }

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

    if (n < 0) {
        // 对负输入使用欧拉反射公式
        // 参见: https://en.wikipedia.org/wiki/Gamma_function#Properties
        return Math.PI / (Math.sin(Math.PI * -n) * gamma(-n));
    } else {
        // Nemes展开近似
        const seriesCoefficient =
            Math.pow(n / Math.E, n) * Math.sqrt(2 * Math.PI * (n + 1 / 6));

        const seriesDenom = n + 1 / 4;

        const seriesExpansion =
            1 +
            1 / 144 / Math.pow(seriesDenom, 2) -
            1 / 12960 / Math.pow(seriesDenom, 3) -
            257 / 207360 / Math.pow(seriesDenom, 4) -
            52 / 2612736 / Math.pow(seriesDenom, 5) +
            5741173 / 9405849600 / Math.pow(seriesDenom, 6) +
            37529 / 18811699200 / Math.pow(seriesDenom, 7);

        return seriesCoefficient * seriesExpansion;
    }
}

export default gamma;