#include <stdint.h>
#include <string.h>
#include <float.h>
#include <math.h>
#define POLY4(x, x2, c3, c2, c1, c0) ((c0 + x*c1) + x2*(c2 + x*c3))
#define POLY8(x, x2, x4, c7, c6, c5, c4, c3, c2, c1, c0) ((x4)*POLY4(x, x2, c7, c6, c5, c4) + POLY4(x, x2, c3, c2, c1, c0))
#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;
  double x2 = x*x, x4 = x2*x2;

  /* coefficients for gamma=7, kmax=8  Lanczos method */
  //denom = ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6)*(x+7)*(x+8));
  denom = POLY8(x, x2, x4, 36, 546, 4536, 22449, 67284, 118124, 109584, 40320) + x4*x4;

  //Ag  = 0.99999999999980993;
  //Ag += (((((((24.45833333334874451493*x + 1053.042534722144725)*x + 19629.51066503033455781)*x + 205786.215499718336944265)*x + 1312212.26587868165589257)*x + 5092745.01776084957)*x
  //   + 11137345.4846599424858181)*x + 10579290.0990759986696)/denom;

  /*Ag = (((((((24.45833333334874451493*x + 1053.042534722144725)*x + 19629.51066503033455781)*x + 205786.215499718336944265)*x + 1312212.26587868165589257)*x + 5092745.01776084957)*x
     + 11137345.4846599424858181)*x + 10579290.0990759986696);*/
  Ag = x4*POLY4(x, x2, 24.45833333334874451493, 1053.042534722144725, 19629.51066503033455781, 205786.215499718336944265) +
          POLY4(x, x2, 1312212.26587868165589257, 5092745.01776084957, 11137345.4846599424858181, 10579290.0990759986696);
  Ag = 0.99999999999980993 + (Ag/denom);

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

double lngamma_lanczos(double x)
{
  double Ag, term2;

  /* coefficients for gamma=7, kmax=8  Lanczos method */
  Ag = 0.99999999999980993;
  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);

  term2 = LogRootTwoPi_ + log(Ag);
  return term2;
}

#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;
  }
}

static inline uint64_t RDTSC(void)
{
    uint32_t a, d;
    __asm__ volatile("lfence \n\t rdtsc" : "=a" (a), "=d" (d));
    return ((uint64_t)d << 32) + a;
}

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);

#if 1
    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-52) {
                nerrs++;
                if (diff > fabs(maxdiff)) maxdiff = diff;
                if (!(i & 0xffff) && (j > 1))
                    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, lngamma_lanczos: %u, largest error: %.7e\n", nerrs, maxdiff); nerrs = 0;
#endif

    {
        float sampledata[64];
        uint64_t time_start, duration, totalduration = 0, a, b;
        struct timeval tv;
        unsigned int it, nb_its = 16384;

        printf("Speed test...\n");
        /* we measure during about 1 seconds */
        for (i = 0; i < 16384; i++) {
            gettimeofday(&tv, NULL);
            time_start = 1000000U*tv.tv_sec + tv.tv_usec;
            a = RDTSC();
            for (it = 0; it < 16384; it++) {
                double ret = lngamma_lanczos_rational((double)it * 0x1p-16);
                __asm__ __volatile__("":"+x"(ret));
            }
            b = RDTSC();
            b -= a;
            gettimeofday(&tv, NULL);
            duration = 1000000U*tv.tv_sec + tv.tv_usec - time_start;
            //sampledata[i] = duration / 16384.0f;
            totalduration += duration;
        }
        printf("Time: %0.1f µs/call [total time=%0.5f s / %lu instructions, its=%u]\n",
               (float)totalduration / nb_its, (float)totalduration / 1000000.0f, b, nb_its);
        //stats(sampledata, 64);

        /* we measure during about 1 seconds */
        for (i = 0; i < 16384; i++) {
            gettimeofday(&tv, NULL);
            time_start = 1000000U*tv.tv_sec + tv.tv_usec;
            a = RDTSC();
            for (it = 0; it < 16384; it++) {
                double ret = lngamma_lanczos((double)it * 0x1p-16);
                __asm__ __volatile__("":"+x"(ret));
            }
            b = RDTSC();
            b -= a;
            gettimeofday(&tv, NULL);
            duration = 1000000U*tv.tv_sec + tv.tv_usec - time_start;
            //sampledata[i] = duration / 16384.0f;
            totalduration += duration;
        }
        printf("Time: %0.1f µs/call [total time=%0.5f s / %lu instructions, its=%u]\n",
               (float)totalduration / nb_its, (float)totalduration / 1000000.0f, b, nb_its);
        //stats(sampledata, 64);
    }

    return 0;
}

#endif