/* 
Author Mihai Preda, 2006. 
The author disclaims copyright to this source code.

The method lgamma() is adapted from FDLIBM 5.3 (http://www.netlib.org/fdlibm/), 
which comes with this copyright notice:
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================

The Lanczos and Stirling approximations are based on:
http://en.wikipedia.org/wiki/Lanczos_approximation
http://en.wikipedia.org/wiki/Stirling%27s_approximation
http://www.gnu.org/software/gsl/
http://jakarta.apache.org/commons/math/
http://my.fit.edu/~gabdo/gamma.txt
*/

class Gamma {
    private static String ulps(double v, double ref) {
        double ulp = ref == 0 ? Math.ulp(.1) : Math.ulp(ref);
        int ulps = (int)Math.floor((v - ref)/ulp + .5);
        //return ulps != 0 ? ""+ulps : "";
        return ""+ulps;
    }
    
    public static void main(String argv[]) {
        double a, b, c, d, e, ulp;        
        for (int i = 1; i < 171; ++i) {
            a = Math.log(factorial(i));
            b = lgamma(i+1);
            e = lanczosLGamma15(i);
            c = f(i);
            d = stirlingLGamma(i);

            System.out.printf("%3d | %6s | %6s | %6s | %6s |\n", i, ulps(b, a), ulps(e, a), ulps(c, a), i >= 10 ? ulps(d, a): "-1000+");
            //System.out.printf("<tr><td>%d</td><td>%s</td><td>%s</td><td>%s</td><td>%s</td></tr>\n", i, ulps(b, a), ulps(e, a), ulps(c, a), i >= 10 ? ulps(d, a): "-1000+");
        }
        
        boolean doBenchmark = true;
        if (doBenchmark) {
            final int N = 20000;
            long t1, t2;
            
            t1 = System.currentTimeMillis();
            for (int r = 0; r < N; ++r) {
                for (int i = 1; i < 171; ++i) {
                    b = lgamma(i+1);
                }
            }
            t2 = System.currentTimeMillis();
            System.out.println("fdlibm's  : " + (t2 - t1));
            
            t1 = System.currentTimeMillis();
            for (int r = 0; r < N; ++r) {
                for (int i = 1; i < 171; ++i) {
                    e = lanczosLGamma15(i);
                }
            }
            t2 = System.currentTimeMillis();
            System.out.println("Lanczos 15: " + (t2 - t1));
            
            t1 = System.currentTimeMillis();
            for (int r = 0; r < N; ++r) {
                for (int i = 1; i < 171; ++i) {
                    c = f(i);
                }
            }
            t2 = System.currentTimeMillis();
            System.out.println("f : " + (t2 - t1));

            /*
            t1 = System.currentTimeMillis();
            for (int r = 0; r < N; ++r) {
                for (int i = 1; i < 171; ++i) {
                    c = lanczosLGamma9(i);
                }
            }
            t2 = System.currentTimeMillis();
            System.out.println("Lanczos 8 : " + (t2 - t1));
            */
            
            t1 = System.currentTimeMillis();
            for (int r = 0; r < N; ++r) {
                for (int i = 1; i < 171; ++i) {
                    d = stirlingLGamma(i);                
                }
            }
            t2 = System.currentTimeMillis();
            System.out.println("Stirling  :  " + (t2 - t1));
        }
    }

    private static final double
        zero = 0.0,
        one  = 1.0,
        two  = 2.0,
        half = .5,
        SQRT2PI    = 2.50662827463100024157,
        LN_SQRT2PI = 0.9189385332046727418;
    
    private static final int HI(double x) {
        return (int)(Double.doubleToLongBits(x) >> 32);
    }
  
    private static final int LO(double x) {
        return (int)Double.doubleToLongBits(x);
    }

    // coefficients for gamma=7, kmax=8  Lanczos method
    private static final double L9[] = {
        0.99999999999980993227684700473478,
        676.520368121885098567009190444019,
        -1259.13921672240287047156078755283,
        771.3234287776530788486528258894,
        -176.61502916214059906584551354,
        12.507343278686904814458936853,
        -0.13857109526572011689554707,
        9.984369578019570859563e-6,
        1.50563273514931155834e-7
    };
    private static final double SQRT2PI_E7 = 0.0022857491179850424; //sqrt(2*pi)/e**7
    static final double lanczosGamma9(double x) {
        if (x <= -1) return Double.NaN;
        double a = L9[0];
        for (int i = 1; i < 9; ++i) {
            a+= L9[i]/(x+i);
        }
        return (SQRT2PI_E7 * a) * Math.pow((x+7.5)/Math.E, x + .5);        
    }
    static final double lanczosLGamma9(double x) {
        if (x <= -1) return Double.NaN;
        double a = L9[0];
        for (int i = 1; i < 9; ++i) {
            a+= L9[i]/(x+i);
        }
        return (LN_SQRT2PI + Math.log(a) - 7.) + (x+.5)*Math.log((x+7.5)/Math.E);
    }

    private static final double[] L15 = {
        0.99999999999999709182,
        57.156235665862923517,
        -59.597960355475491248,
        14.136097974741747174,
        -0.49191381609762019978,
        .33994649984811888699e-4,
        .46523628927048575665e-4,
        -.98374475304879564677e-4,
        .15808870322491248884e-3,
        -.21026444172410488319e-3,
        .21743961811521264320e-3,
        -.16431810653676389022e-3,
        .84418223983852743293e-4,
        -.26190838401581408670e-4,
        .36899182659531622704e-5,
    };
    private static final double G_PLUS_HALF = 607/128. + .5;
 
    static final double lanczosLGamma15(double x) {
        if (x <= -1) return Double.NaN;
        double a = L15[0];
        for (int i = 1; i < 15; ++i) {
            a += L15[i]/(x+i);
        }
        
        double tmp = x + G_PLUS_HALF;
        return (LN_SQRT2PI + Math.log(a)) + (x+.5)*Math.log(tmp) - tmp;
    }

    static final double g(double x) {
        if (x <= -1) return Double.NaN;        
        double tmp = x + 5.2421875;
        return 0.9189385332046727418 + 
            Math.log(
                     0.99999999999999709182 +
                     57.156235665862923517/(x+1) +
                     -59.597960355475491248/(x+2) +
                     14.136097974741747174/(x+3) +
                     -0.49191381609762019978/(x+4) +
                     .33994649984811888699e-4/(x+5) +
                     .46523628927048575665e-4/(x+6) +
                     -.98374475304879564677e-4/(x+7) +
                     .15808870322491248884e-3/(x+8) +
                     -.21026444172410488319e-3/(x+9) +
                     .21743961811521264320e-3/(x+10) +
                     -.16431810653676389022e-3/(x+11) +
                     .84418223983852743293e-4/(x+12) +
                     -.26190838401581408670e-4/(x+13) +
                     .36899182659531622704e-5/(x+14)
                     ) + 
            (x+.5)*Math.log(tmp) 
            - tmp;
    }
    
    static final double f(double x) {
        if (x <= -1) return Double.NaN;        
        final double tmp = x + 5.2421875;
        //final double saveX = x;
        return 0.9189385332046727418 
            + Math.log(
                       0.99999999999999709182 +
                       57.156235665862923517/++x +
                       -59.597960355475491248/++x +
                       14.136097974741747174/++x +
                       -0.49191381609762019978/++x +
                       .33994649984811888699e-4/++x +
                       .46523628927048575665e-4/++x +
                       -.98374475304879564677e-4/++x +
                       .15808870322491248884e-3/++x +
                       -.21026444172410488319e-3/++x +
                       .21743961811521264320e-3/++x +
                       -.16431810653676389022e-3/++x +
                       .84418223983852743293e-4/++x +
                       -.26190838401581408670e-4/++x +
                       .36899182659531622704e-5/++x
                       )
            + (tmp-4.7421875)*Math.log(tmp) - tmp
            //+ (saveX + .5)*Math.log(tmp) + /*Math.sqrt(tmp)*/ - tmp
            ;
    }
             
    private static final double 
        SC1 = 0.08333333333333333, 
        SC2 = 0.003472222222222222, 
        SC3 = -0.0026813271604938273,
        SC4 = -2.2947209362139917E-4,
        LC1 = 0.08333333333333333,
        LC2 = -0.002777777777777778,
        LC3 = 7.936507936507937E-4,
        LC4 = -5.952380952380953E-4;
    static final double stirlingGamma(double x) {
        final double 
            r1 = 1./x,
            r2 = r1*r1,
            r4 = r2*r2;
        return SQRT2PI * Math.sqrt(x) * (1 + SC1*r1 + SC2*r2 + SC3*r1*r2 + SC4*r4) * Math.pow(x/Math.E, x);
    }

    static final double stirlingLGamma(double x) {
        final double 
            r1 = 1./x,
            r2 = r1*r1,
            r3 = r1*r2,
            r5 = r2*r3,
            r7 = r3*r3*r1;
        return (x+.5)*Math.log(x) -x + LN_SQRT2PI + LC1*r1 + LC2*r3 + LC3*r5 + LC4*r7;
    }

    static final double FACT[] = {
        1.0,
        40320.0,
        2.0922789888E13,
        6.204484017332394E23,
        2.631308369336935E35,
        8.159152832478977E47,
        1.2413915592536073E61,
        7.109985878048635E74,
        1.2688693218588417E89,
        6.1234458376886085E103,
        7.156945704626381E118,
        1.8548264225739844E134,
        9.916779348709496E149,
        1.0299016745145628E166,
        1.974506857221074E182,
        6.689502913449127E198,
        3.856204823625804E215,
        3.659042881952549E232,
        5.5502938327393044E249,
        1.3113358856834524E267,
        4.7147236359920616E284,
        2.5260757449731984E302,
    };

    static final double factorial(double x) {
        if (x <= -1) {
            return Double.NaN;
        }
        if (x <= 170) {
            if (Math.floor(x) == x) {
                int n = (int)x;
                double extra = x;
                switch (n & 7) {
                case 7: extra *= --x;
                case 6: extra *= --x;
                case 5: extra *= --x;
                case 4: extra *= --x;
                case 3: extra *= --x;
                case 2: extra *= --x;
                case 1: return FACT[n >> 3] * extra;
                case 0: return FACT[n >> 3];
                }
            }
        }
        return Math.exp(lgamma(x+1));
    }

    private static final double
        a0  =  7.72156649015328655494e-02,
        a1  =  3.22467033424113591611e-01,
        a2  =  6.73523010531292681824e-02,
        a3  =  2.05808084325167332806e-02,
        a4  =  7.38555086081402883957e-03,
        a5  =  2.89051383673415629091e-03,
        a6  =  1.19270763183362067845e-03,
        a7  =  5.10069792153511336608e-04,
        a8  =  2.20862790713908385557e-04,
        a9  =  1.08011567247583939954e-04,
        a10 =  2.52144565451257326939e-05,
        a11 =  4.48640949618915160150e-05,
        tc  =  1.46163214496836224576e+00,
        tf  = -1.21486290535849611461e-01,
        tt  = -3.63867699703950536541e-18,
        t0  =  4.83836122723810047042e-01,
        t1  = -1.47587722994593911752e-01,
        t2  =  6.46249402391333854778e-02,
        t3  = -3.27885410759859649565e-02,
        t4  =  1.79706750811820387126e-02,
        t5  = -1.03142241298341437450e-02,
        t6  =  6.10053870246291332635e-03,
        t7  = -3.68452016781138256760e-03,
        t8  =  2.25964780900612472250e-03,
        t9  = -1.40346469989232843813e-03,
        t10 =  8.81081882437654011382e-04,
        t11 = -5.38595305356740546715e-04,
        t12 =  3.15632070903625950361e-04,
        t13 = -3.12754168375120860518e-04,
        t14 =  3.35529192635519073543e-04,
        u0  = -7.72156649015328655494e-02,
        u1  =  6.32827064025093366517e-01,
        u2  =  1.45492250137234768737e+00,
        u3  =  9.77717527963372745603e-01,
        u4  =  2.28963728064692451092e-01,
        u5  =  1.33810918536787660377e-02,
        v1  =  2.45597793713041134822e+00,
        v2  =  2.12848976379893395361e+00,
        v3  =  7.69285150456672783825e-01,
        v4  =  1.04222645593369134254e-01,
        v5  =  3.21709242282423911810e-03,
        s0  = -7.72156649015328655494e-02,
        s1  =  2.14982415960608852501e-01,
        s2  =  3.25778796408930981787e-01,
        s3  =  1.46350472652464452805e-01,
        s4  =  2.66422703033638609560e-02,
        s5  =  1.84028451407337715652e-03,
        s6  =  3.19475326584100867617e-05,
        r1  =  1.39200533467621045958e+00,
        r2  =  7.21935547567138069525e-01,
        r3  =  1.71933865632803078993e-01,
        r4  =  1.86459191715652901344e-02,
        r5  =  7.77942496381893596434e-04,
        r6  =  7.32668430744625636189e-06,
        w0  =  4.18938533204672725052e-01,
        w1  =  8.33333333333329678849e-02,
        w2  = -2.77777777728775536470e-03,
        w3  =  7.93650558643019558500e-04,
        w4  = -5.95187557450339963135e-04,
        w5  =  8.36339918996282139126e-04,
        w6  = -1.63092934096575273989e-03;
    
    static final double lgamma(double x) {
        double t,y,z,p,p1,p2,p3,q,r,w;
        int i;

        int hx = HI(x);
        int lx = LO(x);

        /* purge off +-inf, NaN, +-0, and negative arguments */
        int ix = hx&0x7fffffff;
        if (ix >= 0x7ff00000) return Double.POSITIVE_INFINITY;
        if ((ix|lx)==0 || hx < 0) return Double.NaN;
        if (ix<0x3b900000) {	/* |x|<2**-70, return -log(|x|) */
            return -Math.log(x);
        }

        /* purge off 1 and 2 */
        if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;
        /* for x < 2.0 */
        else if(ix<0x40000000) {
            if(ix<=0x3feccccc) { 	/* lgamma(x) = lgamma(x+1)-log(x) */
                r = -Math.log(x);
                if(ix>=0x3FE76944) {y = one-x; i= 0;}
                else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;}
                else {y = x; i=2;}
            } else {
                r = zero;
                if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */
                else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;} /* [1.23,1.73] */
                else {y=x-one;i=2;}
            }
            
            switch(i) {
            case 0:
                z = y*y;
                p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
                p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
                p  = y*p1+p2;
                r  += (p-0.5*y); break;
            case 1:
                z = y*y;
                w = z*y;
                p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12)));	/* parallel comp */
                p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
                p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
                p  = z*p1-(tt-w*(p2+y*p3));
                r += (tf + p); break;
            case 2:	
                p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
                p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
                r += (-0.5*y + p1/p2);
            }
        }
        else if(ix<0x40200000) { 			/* x < 8.0 */
            i = (int)x;
            t = zero;
            y = x-(double)i;
            p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
            q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
            r = half*y+p/q;
            z = one;	/* lgamma(1+s) = log(s) + lgamma(s) */
            switch(i) {
            case 7: z *= (y+6.0);	/* FALLTHRU */
            case 6: z *= (y+5.0);	/* FALLTHRU */
            case 5: z *= (y+4.0);	/* FALLTHRU */
            case 4: z *= (y+3.0);	/* FALLTHRU */
            case 3: z *= (y+2.0);	/* FALLTHRU */
                r += Math.log(z); break;
            }
            /* 8.0 <= x < 2**58 */
        } else if (ix < 0x43900000) {
            t = Math.log(x);
            z = one/x;
            y = z*z;
            w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
            r = (x-half)*(t-one)+w;
        } else 
            /* 2**58 <= x <= inf */
            r =  x*(Math.log(x)-one);
        return r;
    }
}