#include <stdint.h>
#include <string.h>
#include <float.h>
#include <math.h>
#define LogRootTwoPi_  0.9189385332046727418
#define SqrtPi_        1.77245385090551602

static double lanczos_term1(double x)
{
    return (x+0.5)*log(x+7.5) - (x+7.5);
}

double lngamma_lanczos_rational(double x)
{
  double Ag, denom, term2;

  /* coefficients for gamma=7, kmax=8  Lanczos method */
  denom = ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6)*(x+7));

  Ag = 0.99999999999980993227684700473478;
  Ag += ((((((24.458333182785471*x + 857.375865044089725)*x + 12770.503696196242686)*x + 103622.185635044379367)*x + 483234.779779163822534)*x + 1226866.777550342082017123)*x + 1322411.262289644971388091)/denom;
  Ag += 1.50563273514931155834e-7/(x+8);

  /* LogRootTwoPi_ + log(Ag(x)) */
  term2 = LogRootTwoPi_ + log(Ag);
  return term2;
}

/* Lanczos method for real x > 0;
 * gamma=7, truncated at 1/(z+8) 
 * [J. SIAM Numer. Anal, Ser. B, 1 (1964) 86]
 */

inline double lngamma_lanczos(double x)
{
  double Ag, val;

  /* coefficients for gamma=7, kmax=8  Lanczos method */
  Ag = 0.99999999999980993227684700473478;
  Ag += 676.520368121885098567009190444019/(x+1);
  Ag -= 1259.13921672240287047156078755283/(x+2);
  Ag += 771.3234287776530788486528258894/(x+3);
  Ag -= 176.61502916214059906584551354/(x+4);
  Ag += 12.507343278686904814458936853/(x+5);
  Ag -= 0.13857109526572011689554707/(x+6);
  Ag += 9.984369578019570859563e-6/(x+7);
  Ag += 1.50563273514931155834e-7/(x+8);

  /* (x+0.5)*log(x+7.5) - (x+7.5) + LogRootTwoPi_ + log(Ag(x)) */
  val = LogRootTwoPi_ + log(Ag);
  return val;
}

#ifdef TEST

#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <sys/time.h>

typedef union _f32 {
    float f;
    uint32_t i;
} _f32;

static double rel_err_diff(double x1, double x2)
{
  if(x1 == 0.0 && x2 == 0.0) {
    return 0.0;
  } else if (x1 == 0.0) {
    return fabs(x2);
  } else if(x1 <= 1e308 && x2 <= 1e308 && (x1 + x2 != 0.0)) {
    return fabs((x1-x2)/(x1+x2));
    //return fabs(x1-x2);
  } else {
    return 1.0;
  }
}

typedef union _d64
{
  double d;
  struct {
    uint32_t lsw;
    uint32_t msw;
  } i;
} _d64;

static const double
ln2_hi  =  6.93147180369123816490e-01,	/* 3fe62e42 fee00000 */
ln2_lo  =  1.90821492927058770002e-10,	/* 3dea39ef 35793c76 */
Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */

int ilogb(double x)
{
    _d64 d;
    int hx;

    d.d = fabs(x);
    hx = (int)(d.i.msw>>20);
	if(hx<0x001) {
		return 0x80000001;	/* ilogb(0) = 0x80000001 */
    } else if (hx<0x7ff) {
        return hx-1023;
	} else {
        return 0x7fffffff;
    }
}

double log(double x)
{
	double val,hfsq,f,s,z,R,w,dk;
	int32_t k,hx,i,j;
    _d64 d;

    if (x < 0.0 || x != x)
        return 0.0/0.0; /* log(-#) = NaN, log(NaN) = NaN */

    d.d = fabs(x);
    hx = d.i.msw;
	k = (hx>>20)-1023;
	hx &= 0x000fffff;
	i = (hx+0x95f64)&0x100000;
	k += (i>>20);
    dk=(double)k;

    d.d= x;
	d.i.msw = (hx|(i^0x3ff00000));	/* normalize x or x/2 */
    x = d.d;
	f = x-1.0;
    if(f==0.0) { return dk*ln2_hi + dk*ln2_lo; }
 	s = f/(2.0+f);
	z = s*s;
	w = z*z;
	R =        w*(Lg2+w*(Lg4+w*Lg6));
	R+= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
	hfsq=0.5*f*f;
	val = (f-(hfsq-s*(hfsq+R)-dk*ln2_lo));
    return dk*ln2_hi + val;
}

static __attribute__((always_inline)) inline void disable_denormals(void)
{
    unsigned int a = 0x9f80;
    __asm__ __volatile__ ("ldmxcsr %0" : "+m" (a));
}

int main(void) {
    _f32 f;
    unsigned int i, j, nerrs = 0;
    double d, g1, g2, diff, maxdiff = 0.0;
    int sgn = 0;

    disable_denormals();
    setvbuf(stdout, NULL, _IONBF, 0);

    printf("lanczos(0): %.7g, term1(0): %.7g, lgamma(1.0): %.7g\n", lngamma_lanczos_rational(0.0), lanczos_term1(0.0), lgamma_r(1.0,&sgn));

    nerrs = 0;
    maxdiff = 0.0;
    for (j = 1; j < 7; j++) {
        f.i = ((j+0x7f)<<23);
        printf("start: %.7f (j: %u)\n", f.f, j);
        for (i = 0; i < (1 << 23U); i++) {
            f.i = ((j+0x7f)<<23)|i;
            d = (double)f.f;
            g1 = lgamma_r(d, &sgn);
            g2 = lngamma_lanczos_rational(d-1.0)+lanczos_term1(d-1.0);
            if ((diff = rel_err_diff(g1, g2)) >=  0x1.0p-51) {
                nerrs++;
                if (diff > fabs(maxdiff)) maxdiff = diff;
                if (!(i & 0xffff))
                    printf("Error: lanczos_rational(%.7g): %.7g != lanczos(%.7g): %.7g (diff: %.10e / %a)\n", f.f, g2, f.f, g1, diff, diff);
            }
        }
        printf("end: %.7f\n", f.f);
    }
    printf("Total errors, cosf: %u, largest error: %.7e\n", nerrs, maxdiff); nerrs = 0;

    return 0;
}

#endif