#include "misc.h"
#include <math.h>
#include <stdlib.h> /* for size_t */
#include <gsl/gsl_sf_gamma.h>

static struct {
    int n;
    double f;
} ln_fact_table[LNFACT_MAX + 1] = {
    {0,0},
    {1,0},
    {2,0.69314718055994529},
    {3,1.791759469228055},
    {4,3.1780538303479458},
    {5,4.7874917427820458},
    {6,6.5792512120101012},
    {7,8.5251613610654147},
    {8,10.604602902745251},
    {9,12.801827480081469},
    {10,15.104412573075516},
    {11,17.502307845873887},
    {12,19.987214495661885},
    {13,22.552163853123425},
    {14,25.19122118273868},
    {15,27.89927138384089},
    {16,30.671860106080672},
    {17,33.505073450136891},
    {18,36.395445208033053},
    {19,39.339884187199495},
    {20,42.335616460753485},
    {21,45.380138898476908},
    {22,48.471181351835227},
    {23,51.606675567764377},
    {24,54.784729398112319},
    {25,58.003605222980518},
    {26,61.261701761002001},
    {27,64.557538627006338},
    {28,67.88974313718154},
    {29,71.257038967168015},
    {30,74.658236348830158},
    {31,78.092223553315307},
    {32,81.557959456115043},
    {33,85.054467017581516},
    {34,88.580827542197682},
    {35,92.136175603687093},
    {36,95.719694542143202},
    {37,99.330612454787428},
    {38,102.96819861451381},
    {39,106.63176026064346},
    {40,110.32063971475739},
    {41,114.03421178146171},
    {42,117.77188139974507},
    {43,121.53308151543864},
    {44,125.3172711493569},
    {45,129.12393363912722},
    {46,132.95257503561632},
    {47,136.80272263732635},
    {48,140.67392364823425},
    {49,144.5657439463449},
    {50,148.47776695177302},
    {51,152.40959258449735},
    {52,156.3608363030788},
    {53,160.3311282166309},
    {54,164.32011226319517},
    {55,168.32744544842765},
    {56,172.35279713916279},
    {57,176.39584840699735},
    {58,180.45629141754378},
    {59,184.53382886144948},
    {60,188.6281734236716},
    {61,192.7390472878449},
    {62,196.86618167289001},
    {63,201.00931639928152},
    {64,205.1681994826412},
    {65,209.34258675253685},
    {66,213.53224149456327},
    {67,217.73693411395422},
    {68,221.95644181913033},
    {69,226.1905483237276},
    {70,230.43904356577696},
    {71,234.70172344281826},
    {72,238.97838956183432},
    {73,243.26884900298271},
    {74,247.57291409618688},
    {75,251.89040220972319},
    {76,256.22113555000954},
    {77,260.56494097186322},
    {78,264.92164979855278},
    {79,269.29109765101981},
    {80,273.67312428569369},
    {81,278.06757344036612},
    {82,282.4742926876304},
    {83,286.89313329542699},
    {84,291.32395009427029},
    {85,295.76660135076065},
    {86,300.22094864701415},
    {87,304.68685676566872},
    {88,309.1641935801469},
    {89,313.65282994987905},
    {90,318.1526396202093},
    {91,322.66349912672615},
    {92,327.1852877037752},
    {93,331.71788719692847},
    {94,336.26118197919845},
    {95,340.81505887079902},
    {96,345.37940706226686},
    {97,349.95411804077025},
    {98,354.53908551944079},
    {99,359.1342053695754},
    {100,363.73937555556347},
};

double lnfact(unsigned int n)
{
    /* Returns the logarithm of n!.
     *
     * Uses a lookup table to return results for n < 100. */

    if (n <= LNFACT_MAX)
        return ln_fact_table[n].f;

    return gsl_sf_lnfact(n);
}

double kahan_sum(double *x, size_t n)
{
    /* Returns the sum of the elements of `x` using the Kahan summation algorithm.
     *
     * See https://en.wikipedia.org/wiki/Kahan_summation_algorithm. */
    size_t i;
    double sum, c, y, t;

    sum = 0.0;
    c = 0.0;
    for (i = 0; i < n; i++) {
        y = x[i] - c;
        t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }

    return sum;
}

double interp1d(double x, double *xp, double *yp, size_t n)
{
    /* A fast interpolation routine which assumes that the values in `xp` are
     * evenly spaced.
     *
     * If x < xp[0] returns yp[0] and if x > xp[n-1] returns yp[n-1].  */
    size_t i;

    if (x < xp[0]) return yp[0];
    if (x > xp[n-1]) return yp[n-1];

    i = (x-xp[0])/(xp[1]-xp[0]);

    return yp[i] + (yp[i+1]-yp[i])*(x-xp[i])/(xp[i+1]-xp[i]);
}

int isclose(double a, double b, double rel_tol, double abs_tol)
{
    /* Returns 1 if a and b are "close". This algorithm is taken from Python's
     * math.isclose() function.
     *
     * See https://www.python.org/dev/peps/pep-0485/. */
    return fabs(a-b) <= fmax(rel_tol*fmax(fabs(a),fabs(b)),abs_tol);
}

int allclose(double *a, double *b, size_t n, double rel_tol, double abs_tol)
{
    /* Returns 1 if all the elements of a and b are "close". This algorithm is
     * taken from Python's math.isclose() function.
     *
     * See https://www.python.org/dev/peps/pep-0485/. */
    size_t i;

    for (i = 0; i < n; i++) {
        if (!isclose(a[i],b[i],rel_tol,abs_tol)) return 0;
    }

    return 1;
}

double logsumexp(double *a, size_t n)
{
    /* Returns the log of the sum of the exponentials of the array `a`.
     *
     * This function is designed to reduce underflow when the exponentials of
     * `a` are very small, for example when computing probabilities. */
    size_t i;
    double amax, sum;

    amax = a[0];
    for (i = 0; i < n; i++) {
        if (a[i] > amax) amax = a[i];
    }

    sum = 0.0;

    for (i = 0; i < n; i++) {
        sum += exp(a[i]-amax);
    }

    sum = log(sum);

    return amax + sum;
}

double norm(double x, double mu, double sigma)
{
    /* Returns the PDF for a gaussian random variable with mean `mu` and
     * standard deviation `sigma`. */
    return exp(-pow(x-mu,2)/(2*pow(sigma,2)))/(sqrt(2*M_PI)*sigma);
}

double norm_cdf(double x, double mu, double sigma)
{
    /* Returns the CDF for a gaussian random variable with mean `mu` and
     * standard deviation `sigma`. */
    return erfc(-(x-mu)/(sqrt(2)*sigma))/2.0;
}