#include <stdint.h>
#include <math.h>
#define LogRootTwoPi_  0.9189385332046727418

static const double LogFactTable[8] = {
    0.0, 0.0, 0.6931471805599453, 1.791759469228055, 3.178053830347946, 4.787491742782046, 6.579251212010101, 8.525161361065415
};

static const double
w1  =  8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */
w2  = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */
w3  =  7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */
w4  = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
w5  =  8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
w6  = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */

double sf_lnfactorial(unsigned int n)
{
    double x = (double)n;
    double xinv = 1.0/x;
    double z = xinv*xinv;
    double w = xinv*(w1+z*(w2+z*(w3+z*(w4+z*(w5+z*w6)))));
    return (x+0.5)*(log(x)-1.0)+w+0.5+LogRootTwoPi_;
}