Function logfactorial

numpy/random/src/distributions/logfactorial.c:142–158

 */

double logfactorial(int64_t k)
{
    const double halfln2pi = 0.9189385332046728;

    if (k < (int64_t) (sizeof(logfact)/sizeof(logfact[0]))) {
        /* Use the lookup table. */
        return logfact[k];
    }

    /*
     *  Use the Stirling series, truncated at the 1/k**3 term.
     *  (In a Python implementation of this approximation, the result
     *  was within 2 ULP of the best 64 bit floating point value for
     *  k up to 10000000.)
     */
    return (k + 0.5)*log(k) - k + (halfln2pi + (1.0/k)*(1/12.0 - 1/(360.0*k*k)));
}