//
// normaldistribution.cpp
//

#if defined(_MSC_VER)
#define NOMINMAX // to suppress MSVC's definitions of min() and max()
// These four pragmas are the equivalent to /fp:fast.
#pragma float_control(except, off)
#pragma float_control(precise, off)
#pragma fp_contract(on)
#pragma fenv_access(off)
#endif

#include "normaldistribution.h"
#include <float.h>

namespace {
// The asymptotic expansion  Φ(z) = φ(z)/|z|·[1-1/z^2+...],  Abramowitz & Stegun (26.2.12), suffices for Φ(z) to have
// relative accuracy of 1.64E-16 for z<=-10 with 17 terms inside the square brackets (not counting the leading 1).
// This translates to a maximum of about 9 iterations below, which is competitive with a call to erfc() and never
// less accurate when z<=-10. Note that, as mentioned in section 4 (and discussion of figures 2 and 3) of George
// Marsaglia's article "Evaluating the Normal Distribution" (available at http://www.jstatsoft.org/v11/a05/paper),
// for values of x approaching -8 and below, the error of any cumulative normal function is actually dominated by
// the hardware (or compiler implementation) accuracy of exp(-x²/2) which is not reliably more than 14 digits when
// x becomes large. Still, we should switch to the asymptotic only when it is beneficial to do so.
const double norm_cdf_asymptotic_expansion_first_threshold = -10.0;
const double norm_cdf_asymptotic_expansion_second_threshold = -1 / sqrt(DBL_EPSILON);
} // namespace

double norm_cdf(double z) {
    if (z <= norm_cdf_asymptotic_expansion_first_threshold) {
        // Asymptotic expansion for very negative z following (26.2.12) on page 408
        // in M. Abramowitz and A. Stegun, Pocketbook of Mathematical Functions, ISBN 3-87144818-4.
        double sum = 1;
        if (z >= norm_cdf_asymptotic_expansion_second_threshold) {
            double zsqr = z * z, i = 1, g = 1, x, y, a = DBL_MAX, lasta;
            do {
                lasta = a;
                x = (4 * i - 3) / zsqr;
                y = x * ((4 * i - 1) / zsqr);
                a = g * (x - y);
                sum -= a;
                g *= y;
                ++i;
                a = fabs(a);
            } while (lasta > a && a >= fabs(sum * DBL_EPSILON));
        }
        return -norm_pdf(z) * sum / z;
    }
    return 0.5 * erfc_cody(-z * ONE_OVER_SQRT_TWO);
}

double inverse_norm_cdf(double u) {
    //
    // ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3
    //
    // Produces the normal deviate Z corresponding to a given lower
    // tail area of u; Z is accurate to about 1 part in 10**16.
    // see http://lib.stat.cmu.edu/apstat/241
    //
    const double split1 = 0.425;
    const double split2 = 5.0;
    const double const1 = 0.180625;
    const double const2 = 1.6;

    // Coefficients for P close to 0.5
    const double A0 = 3.3871328727963666080E0;
    const double A1 = 1.3314166789178437745E+2;
    const double A2 = 1.9715909503065514427E+3;
    const double A3 = 1.3731693765509461125E+4;
    const double A4 = 4.5921953931549871457E+4;
    const double A5 = 6.7265770927008700853E+4;
    const double A6 = 3.3430575583588128105E+4;
    const double A7 = 2.5090809287301226727E+3;
    const double B1 = 4.2313330701600911252E+1;
    const double B2 = 6.8718700749205790830E+2;
    const double B3 = 5.3941960214247511077E+3;
    const double B4 = 2.1213794301586595867E+4;
    const double B5 = 3.9307895800092710610E+4;
    const double B6 = 2.8729085735721942674E+4;
    const double B7 = 5.2264952788528545610E+3;
    // Coefficients for P not close to 0, 0.5 or 1.
    const double C0 = 1.42343711074968357734E0;
    const double C1 = 4.63033784615654529590E0;
    const double C2 = 5.76949722146069140550E0;
    const double C3 = 3.64784832476320460504E0;
    const double C4 = 1.27045825245236838258E0;
    const double C5 = 2.41780725177450611770E-1;
    const double C6 = 2.27238449892691845833E-2;
    const double C7 = 7.74545014278341407640E-4;
    const double D1 = 2.05319162663775882187E0;
    const double D2 = 1.67638483018380384940E0;
    const double D3 = 6.89767334985100004550E-1;
    const double D4 = 1.48103976427480074590E-1;
    const double D5 = 1.51986665636164571966E-2;
    const double D6 = 5.47593808499534494600E-4;
    const double D7 = 1.05075007164441684324E-9;
    // Coefficients for P very close to 0 or 1
    const double E0 = 6.65790464350110377720E0;
    const double E1 = 5.46378491116411436990E0;
    const double E2 = 1.78482653991729133580E0;
    const double E3 = 2.96560571828504891230E-1;
    const double E4 = 2.65321895265761230930E-2;
    const double E5 = 1.24266094738807843860E-3;
    const double E6 = 2.71155556874348757815E-5;
    const double E7 = 2.01033439929228813265E-7;
    const double F1 = 5.99832206555887937690E-1;
    const double F2 = 1.36929880922735805310E-1;
    const double F3 = 1.48753612908506148525E-2;
    const double F4 = 7.86869131145613259100E-4;
    const double F5 = 1.84631831751005468180E-5;
    const double F6 = 1.42151175831644588870E-7;
    const double F7 = 2.04426310338993978564E-15;

    if (u <= 0)
        return log(u);
    if (u >= 1)
        return log(1 - u);

    const double q = u - 0.5;
    if (fabs(q) <= split1) {
        const double r = const1 - q * q;
        return q * (((((((A7 * r + A6) * r + A5) * r + A4) * r + A3) * r + A2) * r + A1) * r + A0) /
               (((((((B7 * r + B6) * r + B5) * r + B4) * r + B3) * r + B2) * r + B1) * r + 1.0);
    } else {
        double r = q < 0.0 ? u : 1.0 - u;
        r = sqrt(-log(r));
        double ret;
        if (r < split2) {
            r = r - const2;
            ret = (((((((C7 * r + C6) * r + C5) * r + C4) * r + C3) * r + C2) * r + C1) * r + C0) /
                  (((((((D7 * r + D6) * r + D5) * r + D4) * r + D3) * r + D2) * r + D1) * r + 1.0);
        } else {
            r = r - split2;
            ret = (((((((E7 * r + E6) * r + E5) * r + E4) * r + E3) * r + E2) * r + E1) * r + E0) /
                  (((((((F7 * r + F6) * r + F5) * r + F4) * r + F3) * r + F2) * r + F1) * r + 1.0);
        }
        return q < 0.0 ? -ret : ret;
    }
}