#include "internal.h"

/*
 * Ahrens, J. H., & Dieter, U. (1982).
 * Computer Generation of Poisson Deviates from Modified Normal Distributions.
 */
static double factorial[10] = {1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880.};
static double coeffs[9] = {
    -.500000000, +.333333333, -.249999856,
    +.200011780, -.166684875, +.142187833,
    -.124196313, +.125005956, -.114265030,
};

static inline
double
log1pmx(double x, double off)
{
    int i;
    double r, t;

    if(-0.25 < x && x < 0.25) {
        r = 0;
        t = 1;
        for(i=0;i<arrlen(coeffs);i++) {
            r += coeffs[i]*t;
            t *= x;
        }

        return x*x*r;
    }
    return log(1+x) - off;
}

static inline
double
procf(double mu, double s, int64 K, double *px, double *py, double *fx, double *fy)
{
    double d, V, X;
    double w, b1, b2, c1, c2, c3, c0, c;

    w  = 0.3989422804014327/s;
    b1 = 0.041666666666666664/mu;
    b2 = 0.3*b1*b1;
    c3 = 0.14285714285714285*b1*b2;
    c2 = b2 - 15.*c3;
    c1 = b1 - 6.*b2 + 45.*c3;
    c0 = 1 - b1 + 3.*b2 - 15.*c3;
    c  = .1069/mu;

    if(K < 10) {
        *px = -mu;
        *py = pow(mu,K) / factorial[K];
    }else{
        d = 0.08333333333333333/K;
        d = d - 4.8*d*d*d;
        V = (mu - K) / K;

        *px = K*log1pmx(V,mu-K) - d;
        *py = 0.3989422804014327/sqrt(K);
    }

    X = (K - mu + 0.5)/s;
    *fx = -0.5*X*X;
    *fy = w*(((c3*X*X + c2)*X*X + c1)*X*X + c0);

    return c;
}

static inline
uint64
bigpoisson(double mu)
{
    int64  L,K;
    double G,s,d,U,E,T;
    double px,py,fx,fy,c;

    s = sqrt(mu);
    d = 6*mu*mu;
    L = floor(mu - 1.1484);

stepN:
    G = mu + s*rng·normal();
    K = floor(G);
    if(K<0)
        goto stepP;
stepI:
    if(K>=L)
        return K;
stepS:
    U = rng·random();
    if(d*U >= (mu-K)*(mu-K)*(mu-K))
        return K;
stepP:
    if(G < 0)
        goto stepE;
stepQ:
    c = procf(mu, s, K, &px, &py, &fx, &fy);
stepE:
    E = rng·exponential(1.0);
    U = rng·random();
    U = U + U - 1;
    T = 1.8 + copysign(E,U);
    if(T < 0.6744)
        goto stepE;
    K = floor(mu + s*T);
    c = procf(mu, s, K, &px, &py, &fx, &fy);
stepH:
    if(c*fabs(U) > (py*exp(px + E) - fy*exp(fx + E)))
        goto stepE;
    return K;
}

uint64
rng·poisson(double mean)
{
    int64  n;
    double z;

    if(mean<10.0) {
        for(n=0, z=rng·exponential(1.0); z<mean; ++n, z+=rng·exponential(1.0))
            ;
        return n;
    }

    return bigpoisson(mean);
}