r2rs-nmath 0.1.1

// Translation of nmath's pnorm
//
// Copyright (C) 2020 neek-sss
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <https://www.gnu.org/licenses/>.

use num_traits::Float;
use strafe_type::{FloatConstraint, LogProbability64, Positive64, Probability64, Real64};

use crate::traits::DPQ;

/// The main computation evaluates near-minimax approximations derived
/// from those in "Rational Chebyshev approximations for the error
/// function" by W. J. Cody, Math. Comp., 1969, 631-637.  This
/// transportable program uses rational functions that theoretically
/// approximate the normal distribution function to at least 18
/// significant decimal digits.  The accuracy achieved depends on the
/// arithmetic system, the compiler, the intrinsic functions, and
/// proper selection of the machine-dependent constants.
pub fn pnorm<R1: Into<Real64>, R2: Into<Real64>, P: Into<Positive64>>(
    x: R1,
    mu: R2,
    sigma: P,
    lower_tail: bool,
) -> Probability64 {
    pnorm_inner(x, mu, sigma, lower_tail, false).into()
}

/// The main computation evaluates near-minimax approximations derived
/// from those in "Rational Chebyshev approximations for the error
/// function" by W. J. Cody, Math. Comp., 1969, 631-637.  This
/// transportable program uses rational functions that theoretically
/// approximate the normal distribution function to at least 18
/// significant decimal digits.  The accuracy achieved depends on the
/// arithmetic system, the compiler, the intrinsic functions, and
/// proper selection of the machine-dependent constants.
pub fn log_pnorm<R1: Into<Real64>, R2: Into<Real64>, P: Into<Positive64>>(
    x: R1,
    mu: R2,
    sigma: P,
    lower_tail: bool,
) -> LogProbability64 {
    pnorm_inner(x, mu, sigma, lower_tail, true).into()
}

fn pnorm_inner<R1: Into<Real64>, R2: Into<Real64>, P: Into<Positive64>>(
    x: R1,
    mu: R2,
    sigma: P,
    lower_tail: bool,
    log: bool,
) -> f64 {
    let mut x = x.into().unwrap();
    let mu = mu.into().unwrap();
    let sigma = sigma.into().unwrap();

    let mut p = 0.0;
    let mut cp = 0.0;

    /* Note: The structure of these checks has been carefully thought through.
     * For example, if x == mu and sigma == 0, we get the correct answer 1.
     */
    if !x.is_finite() && mu == x {
        return f64::nan();
    }
    if sigma <= 0.0 {
        if sigma < 0.0 {
            return f64::nan();
        }
        /* sigma = 0 : */
        return if x < mu {
            f64::dt_0(lower_tail, log)
        } else {
            f64::dt_1(lower_tail, log)
        };
    }
    p = (x - mu) / sigma;
    if !p.is_finite() {
        return if x < mu {
            f64::dt_0(lower_tail, log)
        } else {
            f64::dt_1(lower_tail, log)
        };
    }
    x = p;

    pnorm_both(x, &mut p, &mut cp, if lower_tail { 0 } else { 1 }, log);

    if lower_tail {
        p
    } else {
        cp
    }
}

pub fn pnorm_both(x: f64, cum: &mut f64, ccum: &mut f64, i_tail: i32, log: bool) {
    /* i_tail in {0,1,2} means: "lower", "upper", or "both" :
       if(lower) return  *cum := P[X <= x]
       if(upper) return *ccum := P[X >  x] = 1 - P[X <= x]
    */
    static a: [f64; 5] = [
        2.2352520354606839287,
        161.02823106855587881,
        1067.6894854603709582,
        18154.981253343561249,
        0.065682337918207449113,
    ];
    static b: [f64; 4] = [
        47.20258190468824187,
        976.09855173777669322,
        10260.932208618978205,
        45507.789335026729956,
    ];
    static c: [f64; 9] = [
        0.39894151208813466764,
        8.8831497943883759412,
        93.506656132177855979,
        597.27027639480026226,
        2494.5375852903726711,
        6848.1904505362823326,
        11602.651437647350124,
        9842.7148383839780218,
        1.0765576773720192317e-8,
    ];
    static d: [f64; 8] = [
        22.266688044328115691,
        235.38790178262499861,
        1519.377599407554805,
        6485.558298266760755,
        18615.571640885098091,
        34900.952721145977266,
        38912.003286093271411,
        19685.429676859990727,
    ];
    static p: [f64; 6] = [
        0.21589853405795699,
        0.1274011611602473639,
        0.022235277870649807,
        0.001421619193227893466,
        2.9112874951168792e-5,
        0.02307344176494017303,
    ];
    static q: [f64; 5] = [
        1.28426009614491121,
        0.468238212480865118,
        0.0659881378689285515,
        0.00378239633202758244,
        7.29751555083966205e-5,
    ];

    let mut xden = 0.0;
    let mut xnum = 0.0;
    let mut temp = 0.0;
    let mut del = 0.0;
    let mut eps = 0.0;
    let mut xsq = 0.0;
    let mut y = 0.0;
    let mut i = 0;
    let mut lower = 0;
    let mut upper = 0;

    if x.is_nan() {
        *ccum = x;
        *cum = *ccum;
        return;
    }

    /* Consider changing these : */
    eps = 2.2204460492503131e-16 * 0.5;

    /* i_tail in {0,1,2} =^= {lower, upper, both} */
    lower = !i_tail;
    upper = i_tail;

    let d_2 = |x: f64| x.ldexp(-1.0);

    let mut do_del = |X: f64, xsq: &mut f64, cum: &mut f64, ccum: &mut f64, temp: f64| {
        *xsq = X.ldexp(4.0).trunc().ldexp(-4.0);
        del = (X - *xsq) * (X + *xsq);
        if log {
            *cum = (-*xsq * d_2(*xsq)) - d_2(del) + temp.ln();
            if (lower != 0 && x > 0.0) || (upper != 0 && x <= 0.0) {
                *ccum = (-(-*xsq * d_2(*xsq)).exp() * (-d_2(del)).exp() * temp).ln_1p();
            }
        } else {
            *cum = (-*xsq * d_2(*xsq)).exp() * (-d_2(del)).exp() * temp;
            *ccum = 1.0 - *cum;
        }
    };

    let swap_tail = |cum: &mut f64, ccum: &mut f64, temp: &mut f64| {
        if x > 0.0 {
            *temp = *cum;
            if lower != 0 {
                *cum = *ccum
            }
            *ccum = *temp
        }
    };

    y = x.abs();
    if y <= 0.67448975 {
        /* qnorm(3/4) = .6744.... -- earlier had 0.66291 */
        if y > eps {
            xsq = x * x;
            xnum = a[4] * xsq;
            xden = xsq;
            i = 0;
            while i < 3 {
                xnum = (xnum + a[i]) * xsq;
                xden = (xden + b[i]) * xsq;
                i += 1
            }
        } else {
            xden = 0.0;
            xnum = xden
        }

        temp = x * (xnum + a[3]) / (xden + b[3]);
        if lower != 0 {
            *cum = 0.5 + temp
        }
        if upper != 0 {
            *ccum = 0.5 - temp
        }
        if log {
            if lower != 0 {
                *cum = (*cum).ln()
            }
            if upper != 0 {
                *ccum = (*ccum).ln()
            }
        }
    } else if y <= strafe_consts::SQRT_32 {
        /* Evaluate pnorm for 0.674.. = qnorm(3/4) < |x| <= sqrt(32) ~= 5.657 */

        xnum = c[8] * y;
        xden = y;
        i = 0;
        while i < 7 {
            xnum = (xnum + c[i]) * y;
            xden = (xden + d[i]) * y;
            i += 1
        }
        temp = (xnum + c[7]) / (xden + d[7]);

        do_del(y, &mut xsq, cum, ccum, temp);
        swap_tail(cum, ccum, &mut temp);
    } else if log && y < 1e170
        /*  ^^^^^ MM FIXME: could speedup for log_p and  |x| >> 5.657 !
         * Then, make use of  Abramowitz & Stegun, 26.2.13, p.932,  something like
    
         xsq = x*x;
    
         if(xsq * DBL_EPSILON < 1.)
            del = (1. - (1. - 5./(xsq+6.)) / (xsq+4.)) / (xsq+2.);
         else
            del = 0.;
         *cum = -.5*xsq - M_LN_SQRT_2PI - log(x) + log1p(-del);
         *ccum = log1p(-exp(*cum)); /.* ~ log(1) = 0 *./
    
         swap_tail;
    
         Yes, but xsq might be infinite.;
         well, actually  x = -1.34..e154 = -sqrt(DBL_MAX) already overflows x^2
        */
        || lower != 0 && -37.5193 < x && x < 8.2924
        || upper != 0 && -8.2924 < x && x < 37.5193
    {
        /* else   |x| > sqrt(32) = 5.657 :
        * the next two case differentiations were really for lower=T, log=F
        * Particularly  *not* for  log !

        * Cody had (-37.5193 < x  &&  x < 8.2924) ; R originally had y < 50
        *
        * Note that we do want symmetry(0), lower/upper -> hence use y
        */

        /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
        xsq = 1.0 / (x * x); /* large x such that probs are 0 or 1 */
        xnum = p[5] * xsq;
        xden = xsq;
        i = 0;
        while i < 4 {
            xnum = (xnum + p[i]) * xsq;
            xden = (xden + q[i]) * xsq;
            i += 1
        }
        temp = xsq * (xnum + p[4]) / (xden + q[4]);
        temp = (strafe_consts::_1DSQRT_2TPI - temp) / y;

        do_del(x, &mut xsq, cum, ccum, temp);
        swap_tail(cum, ccum, &mut temp);
    } else if x > 0.0 {
        *cum = f64::d_1(log);
        *ccum = f64::d_0(log)
    } else {
        *cum = f64::d_0(log);
        *ccum = f64::d_1(log)
    };
}
/* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
/* (1./x)*(1./x) might be better */