r2rs-nmath 0.1.1

// Translation of nmath's qt
//
// 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_type::{FloatConstraint, LogProbability64, Probability64, Rational64, Real64};

use crate::{
    distribution::{
        norm::{log_qnorm, qnorm},
        t::{dt, pt},
    },
    traits::{TrigPI, DPQ},
};

/// The "Student" t distribution quantile function.
///
/// Supplemented by inversion for 0 < ndf < 1.
pub fn qt<P: Into<Probability64>, R: Into<Rational64>>(p: P, ndf: R, lower_tail: bool) -> Real64 {
    let p = p.into().unwrap();
    qt_inner(p, ndf, lower_tail, false)
}

/// The "Student" t distribution quantile function.
///
/// Supplemented by inversion for 0 < ndf < 1.
pub fn log_qt<LP: Into<LogProbability64>, R: Into<Rational64>>(
    p: LP,
    ndf: R,
    lower_tail: bool,
) -> Real64 {
    let p = p.into().unwrap();
    qt_inner(p, ndf, lower_tail, true)
}

fn qt_inner<R: Into<Rational64>>(mut p: f64, ndf: R, lower_tail: bool, log: bool) -> Real64 {
    let ndf = ndf.into().unwrap();
    let eps = 1.0e-12;

    let mut P;
    let mut q;

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

    if ndf < 1.0 {
        /* based on qnt */
        let accu = 1e-13;
        let Eps = 1e-11; /* must be > accu */

        let mut ux;
        let mut lx;
        let mut nx;
        let mut pp;

        let mut iter = 0;

        p = p.dt_qiv(lower_tail, log);

        /* Invert pt(.) :
         * 1. finding an upper and lower bound */
        if p > 1.0 - f64::EPSILON {
            return f64::infinity().into();
        }
        pp = (1.0 - f64::EPSILON).min(p * (1.0 + Eps));
        ux = 1.0;
        while ux < f64::max_value() && pt(ux, ndf, true).unwrap() < pp {
            ux *= 2.0
        }
        pp = p * (1.0 - Eps);
        lx = -1.0;
        while lx > -f64::max_value() && pt(lx, ndf, true).unwrap() > pp {
            lx *= 2.0
        }

        /* 2. interval (lx,ux)  halving
          regula falsi failed on qt(0.1, 0.1)
        */
        loop {
            nx = 0.5 * (lx + ux);
            if pt(nx, ndf, true).unwrap() > p {
                ux = nx;
            } else {
                lx = nx;
            }
            iter += 1;
            if (ux - lx) / nx.abs() < accu || iter > 1000 {
                break;
            }
        }

        // if(iter >= 1000) ML_WARNING(ME_PRECISION, "qt");

        return (0.5 * (lx + ux)).into();
    }

    /* Old comment:
     * FIXME: "This test should depend on  ndf  AND p  !!
     * -----  and in fact should be replaced by
     * something like Abramowitz & Stegun 26.7.5 (p.949)"
     *
     * That would say that if the qnorm value is x then
     * the result is about x + (x^3+x)/4df + (5x^5+16x^3+3x)/96df^2
     * The differences are tiny even if x ~ 1e5, and qnorm is not
     * that accurate in the extreme tails.
     */
    if ndf > 1e20 {
        return if log {
            log_qnorm(p, 0.0, 1.0, lower_tail)
        } else {
            qnorm(p, 0.0, 1.0, lower_tail)
        };
    }

    P = p.d_qiv(log); /* if exp(p) underflows, we fix below */

    let neg = (!lower_tail || P < 0.5) && (lower_tail || P > 0.5);
    let is_neg_lower = lower_tail == neg; /* both TRUE or FALSE == !xor */
    if neg {
        P = 2.0
            * if log {
                if lower_tail {
                    P
                } else {
                    -p.exp_m1()
                }
            } else {
                p.d_lval(lower_tail)
            }
    } else {
        P = 2.0
            * if log {
                if lower_tail {
                    -p.exp_m1()
                } else {
                    P
                }
            } else {
                p.d_cval(lower_tail)
            };
    }
    /* 0 <= P <= 1 ; P = 2*min(P', 1 - P')  in all cases */

    if (ndf - 2.0).abs() < eps {
        /* df ~= 2 */
        if P > f64::min_value() {
            if 3.0 * P < f64::epsilon() {
                /* P ~= 0 */
                q = 1.0 / P.sqrt();
            } else if P > 0.9 {
                /* P ~= 1 */
                q = (1.0 - P) * (2.0 / (P * (2.0 - P))).sqrt();
            } else {
                /* eps/3 <= P <= 0.9 */
                q = (2.0 / (P * (2.0 - P)) - 2.0).sqrt();
            }
        } else {
            /* P << 1, q = 1/sqrt(P) = ... */
            if log {
                q = if is_neg_lower {
                    (-p / 2.0).exp() / f64::SQRT_2()
                } else {
                    1.0 / (-p.exp_m1()).sqrt()
                };
            } else {
                q = f64::infinity();
            }
        }
    } else if ndf < 1.0 + eps {
        /* df ~= 1  (df < 1 excluded above): Cauchy */
        if P == 1.0 {
            q = 0.0; // some versions of tanpi give Inf, some NaN
        } else if P > 0.0 {
            q = 1.0 / (P / 2.0).tan_pi(); /* == - tan((P+1) * M_PI_2) -- suffers for P ~= 0 */
        } else {
            /* P = 0, but maybe = 2*exp(p) ! */
            if log {
                /* 1/tan(e) ~ 1/e */
                q = if is_neg_lower {
                    f64::FRAC_1_PI() * (-p).exp()
                } else {
                    -1.0 / (f64::PI() * p.exp_m1())
                };
            } else {
                q = f64::infinity()
            }
        }
    } else {
        /*-- usual case;  including, e.g.,  df = 1.1 */
        let mut x = 0.0;
        let mut y = f64::nan();
        let mut log_P2 = 0.0; /* -Wall */
        let a = 1.0 / (ndf - 0.5);
        let b = 48.0 / (a * a);
        let mut c = ((20700.0 * a / b - 98.0) * a - 16.0) * a + 96.36;
        let d = ((94.5 / (b + c) - 3.0) / b + 1.0) * (a * f64::FRAC_PI_2()).sqrt() * ndf;

        let P_ok1 = P > f64::min_value() || !log;
        let mut P_ok = P_ok1; // when true (after check below), use "normal scale": log_p=FALSE
        if P_ok1 {
            y = (d * P).powf(2.0 / ndf);
            P_ok = y >= f64::epsilon();
        }
        if !P_ok {
            // log.p && P very.small  ||  (d*P)^(2/df) =: y < eps_c
            log_P2 = if is_neg_lower {
                p.d_log(log)
            } else {
                p.d_lexp(log)
            }; /* == log(P / 2) */
            x = (d.ln() + f64::LN_2() + log_P2) / ndf;
            y = (2.0 * x).exp();
        }

        if (ndf < 2.1 && P > 0.5) || y > 0.05 + a {
            /* P > P0(df) */
            /* Asymptotic inverse expansion about normal */
            if P_ok {
                x = qnorm(0.5 * P, 0.0, 1.0, /*lower_tail*/ true).unwrap();
            } else {
                /* log_p && P underflowed */
                x = log_qnorm(log_P2, 0.0, 1.0, lower_tail).unwrap();
            }

            y = x * x;
            if ndf < 5.0 {
                c += 0.3 * (ndf - 4.5) * (x + 0.6);
            }
            c = (((0.05 * d * x - 5.0) * x - 7.0) * x - 2.0) * x + b + c;
            y = (((((0.4 * y + 6.3) * y + 36.0) * y + 94.5) / c - y - 3.0) / b + 1.0) * x;
            y = (a * y * y).exp_m1();
            q = (ndf * y).sqrt();
        } else if !P_ok && x < -f64::LN_2() * f64::MANTISSA_DIGITS as f64 {
            /* 0.5* log(DBL_EPSILON) */
            /* y above might have underflown */
            q = ndf.sqrt() * (-x).exp();
        } else {
            /* re-use 'y' from above */
            y = ((1.0 / (((ndf + 6.0) / (ndf * y) - 0.089 * d - 0.822) * (ndf + 2.0) * 3.0)
                + 0.5 / (ndf + 4.0))
                * y
                - 1.0)
                * (ndf + 1.0)
                / (ndf + 2.0)
                + 1.0 / y;
            q = (ndf * y).sqrt();
        }

        /* Now apply 2-term Taylor expansion improvement (1-term = Newton):
         * as by Hill (1981) [ref.above] */

        /* FIXME: This can be far from optimal when log_p = TRUE
         *      but is still needed, e.g. for qt(-2, df=1.01, log=TRUE).
         *	Probably also improvable when  lower_tail = FALSE */

        if P_ok1 {
            let M = ((f64::max_value() / 2.0).sqrt() - ndf).abs();
            let mut it = 0;
            loop {
                it += 1;
                y = dt(q, ndf, false).unwrap();
                x = (pt(q, ndf, false).unwrap() - P / 2.0) / y;
                if !(it < 10 && y > 0.0 && x.is_finite() && x.abs() > 1e-14 * q.abs()) {
                    break;
                }
                /* Newton (=Taylor 1 term):
                 *  q += x;
                 * Taylor 2-term : */
                let F = if q.abs() < M {
                    q * (ndf + 1.0) / (2.0 * (q * q + ndf))
                } else {
                    (ndf + 1.0) / (2.0 * (q + ndf / q))
                };
                let del_q = x * (1.0 + x * F);
                if del_q.is_finite() && (q + del_q).is_finite() {
                    q += del_q;
                } else if x.is_finite() && (q + x).is_finite() {
                    q += x;
                } else {
                    // FIXME??  if  q+x = +/-Inf is *better* than q should use it
                    break; // cannot improve  q  with a Newton/Taylor step
                }
            }
        }
    }
    if neg { -q } else { q }.into()
}