r2rs-nmath 0.1.1

// Translation of nmath's qnorm
//
// 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, FloatConst};
use strafe_consts::_2TPI;
use strafe_type::{FloatConstraint, LogProbability64, Positive64, Probability64, Real64};

use crate::traits::DPQ;

/// Compute the quantile function for the normal distribution.
///
///	The algorithm AS 241 of Wichura is used,
/// and has been improved for the very extreme tail (and log_p=TRUE)
pub fn qnorm<PR: Into<Probability64>, R: Into<Real64>, PO: Into<Positive64>>(
    p: PR,
    mu: R,
    sigma: PO,
    lower_tail: bool,
) -> Real64 {
    let p = p.into().unwrap();
    qnorm_inner(p, mu, sigma, lower_tail, false)
}

/// Compute the quantile function for the normal distribution.
///
///	The algorithm AS 241 of Wichura is used,
/// and has been improved for the very extreme tail (and log_p=TRUE)
pub fn log_qnorm<LP: Into<LogProbability64>, R: Into<Real64>, P: Into<Positive64>>(
    p: LP,
    mu: R,
    sigma: P,
    lower_tail: bool,
) -> Real64 {
    let p = p.into().unwrap();
    qnorm_inner(p, mu, sigma, lower_tail, true)
}

fn qnorm_inner<R: Into<Real64>, P: Into<Positive64>>(
    p: f64,
    mu: R,
    sigma: P,
    lower_tail: bool,
    log: bool,
) -> Real64 {
    let mu = mu.into().unwrap();
    let sigma = sigma.into().unwrap();

    let mut p_ = 0.0; /* real lower_tail prob. p */
    let mut q = 0.0;
    let mut r = 0.0;
    let mut val = 0.0;

    if let Some(ret) = p.q_p01_boundaries(f64::NEG_INFINITY, f64::infinity(), lower_tail, log) {
        return ret.into();
    }

    if sigma == 0.0 {
        return mu.into();
    }

    p_ = p.dt_qiv(lower_tail, log);
    q = p_ - 0.5;

    /*-- use AS 241 --- */
    /* double ppnd16_(double *p, long *ifault)*/
    /*      ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3

            Produces the normal deviate Z corresponding to a given lower
            tail area of P; Z is accurate to about 1 part in 10**16.

            (original fortran code used PARAMETER(..) for the coefficients
             and provided hash codes for checking them...)
    */
    if q.abs() <= 0.425 {
        /* |p~ - 0.5| <= .425  <==> 0.075 <= p~ <= 0.925 */
        r = 0.180625 - q * q; // = .425^2 - q^2  >= 0
        val = q
            * (((((((r * 2509.0809287301226727 + 33430.575583588128105) * r
                + 67265.770927008700853)
                * r
                + 45921.953931549871457)
                * r
                + 13731.693765509461125)
                * r
                + 1971.5909503065514427)
                * r
                + 133.14166789178437745)
                * r
                + 3.387132872796366608)
            / (((((((r * 5226.495278852854561 + 28729.085735721942674) * r
                + 39307.89580009271061)
                * r
                + 21213.794301586595867)
                * r
                + 5394.1960214247511077)
                * r
                + 687.1870074920579083)
                * r
                + 42.313330701600911252)
                * r
                + 1.0)
    } else {
        /* closer than 0.075 from {0,1} boundary :
         *  r := log(p~);  p~ = min(p, 1-p) < 0.075 :  */
        let lp;
        if log && ((lower_tail && q <= 0.0) || (!lower_tail && q > 0.0)) {
            lp = p;
        } else {
            lp = (if q > 0.0 {
                p.dt_civ(lower_tail, log) /* 1-p */
            } else {
                p_ /* = R_DT_Iv(p) ^=  p */
            })
            .ln();
        }
        // r = sqrt( - log(min(p,1-p)) )  <==>  min(p, 1-p) = exp( - r^2 ) :
        r = (-lp).sqrt();

        if r <= 5.0 {
            /* r = sqrt(-log(r))  <==>  min(p, 1-p) = exp( - r^2 ) */
            /* <==> min(p,1-p) >= exp(-25) ~= 1.3888e-11 */
            r += -1.6;
            val = (((((((r * 7.7454501427834140764e-4 + 0.0227238449892691845833) * r
                + 0.24178072517745061177)
                * r
                + 1.27045825245236838258)
                * r
                + 3.64784832476320460504)
                * r
                + 5.7694972214606914055)
                * r
                + 4.6303378461565452959)
                * r
                + 1.42343711074968357734)
                / (((((((r * 1.05075007164441684324e-9 + 5.475938084995344946e-4) * r
                    + 0.0151986665636164571966)
                    * r
                    + 0.14810397642748007459)
                    * r
                    + 0.68976733498510000455)
                    * r
                    + 1.6763848301838038494)
                    * r
                    + 2.05319162663775882187)
                    * r
                    + 1.0)
        } else if r <= 27.0 {
            /* p is very close to  0 or 1: r in (5, 27] :
            	*  r >   5 <==> min(p,1-p)  < exp(-25) = 1.3888..e-11
            	*  r <= 27 <==> min(p,1-p) >= exp(-27^2) = exp(-729) ~= 2.507972e-317
            	* i.e., we are just barely in the range where min(p, 1-p) has not yet underflowed to zero.
            	*/
            // Wichura, p.478: minimax rational approx R_3(t) is for 5 <= t <= 27  (t :== r)
            /* very close to  0 or 1 */
            r += -5.0;
            val = (((((((r * 2.01033439929228813265e-7 + 2.71155556874348757815e-5) * r
                + 0.0012426609473880784386)
                * r
                + 0.026532189526576123093)
                * r
                + 0.29656057182850489123)
                * r
                + 1.7848265399172913358)
                * r
                + 5.4637849111641143699)
                * r
                + 6.6579046435011037772)
                / (((((((r * 2.04426310338993978564e-15 + 1.4215117583164458887e-7) * r
                    + 1.8463183175100546818e-5)
                    * r
                    + 7.868691311456132591e-4)
                    * r
                    + 0.0148753612908506148525)
                    * r
                    + 0.13692988092273580531)
                    * r
                    + 0.59983220655588793769)
                    * r
                    + 1.0)
        } else {
            // r > 27: p is *really* close to 0 or 1 .. practically only when log_p =TRUE
            if r >= 6.4e8 {
                // p is *very extremely* close to 0 or 1
                // Using the asymptotical formula ("0-th order"): qn = sqrt(2*s)
                val = r * f64::SQRT_2();
            } else {
                let s2 = -lp.ldexp(1.0); // = -2*lp = 2s
                let mut x2 = s2 - (_2TPI * s2).ln(); // = xs_1
                                                     // if(r >= 36000.)  # <==> s >= 36000^2   use x2 = xs_1  above
                if r < 36000.0 {
                    x2 = s2 - (_2TPI * x2).ln() - 2.0 / (2.0 + x2); // == xs_2
                    if r < 840.0 {
                        // 27 < r < 840
                        x2 = s2 - (_2TPI * x2).ln()
                            + 2.0 * (-(1.0 - 1.0 / (4.0 + x2)) / (2.0 + x2)).ln_1p(); // == xs_3
                        if r < 109.0 {
                            // 27 < r < 109
                            x2 = s2 - (_2TPI * x2).ln()
                                + 2.0
                                    * (-(1.0 - (1.0 - 5.0 / (6.0 + x2)) / (4.0 + x2)) / (2.0 + x2))
                                        .ln_1p(); // == xs_4
                            if r < 55.0 {
                                // 27 < r < 55
                                x2 = s2 - (_2TPI * x2).ln()
                                    + 2.0
                                        * (-(1.0
                                            - (1.0 - (5.0 - 9.0 / (8.0 + x2)) / (6.0 + x2))
                                                / (4.0 + x2))
                                            / (2.0 + x2))
                                            .ln_1p(); // == xs_5
                            }
                        }
                    }
                }
                val = x2.sqrt();
            }
        }

        if q < 0.0 {
            val = -val
            /* return (q >= 0. r : -r ;*/
        }
    }

    (mu + sigma * val).into()
}