r2rs-nmath 0.1.1

// Translation of nmath's qnchisq
//
// 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, Rational64, Real64,
};

use crate::{
    distribution::chisq::{log_qchisq, pnchisq_raw, qchisq},
    traits::DPQ,
};

pub fn qnchisq<PR: Into<Probability64>, R: Into<Rational64>, PO: Into<Positive64>>(
    p: PR,
    df: R,
    ncp: PO,
    lower_tail: bool,
) -> Real64 {
    let p = p.into().unwrap();
    qnchisq_inner(p, df, ncp, lower_tail, false)
}

pub fn log_qnchisq<LP: Into<LogProbability64>, R: Into<Rational64>, P: Into<Positive64>>(
    p: LP,
    df: R,
    ncp: P,
    lower_tail: bool,
) -> Real64 {
    let p = p.into().unwrap();
    qnchisq_inner(p, df, ncp, lower_tail, true)
}

fn qnchisq_inner<R: Into<Rational64>, P: Into<Positive64>>(
    mut p: f64,
    df: R,
    ncp: P,
    mut lower_tail: bool,
    log: bool,
) -> Real64 {
    let df = df.into().unwrap();
    let ncp = ncp.into().unwrap();

    static accu: f64 = 1e-13;
    static racc: f64 = 4.0 * 2.2204460492503131e-16;
    /* these two are for the "search" loops, can have less accuracy: */
    static Eps: f64 = 1e-11; /* must be > accu */
    static rEps: f64 = 1e-10; /* relative tolerance ... */

    let mut ux = 0.0;
    let mut lx = 0.0;
    let mut ux0 = 0.0;
    let mut nx = 0.0;
    let mut pp = 0.0;

    if !df.is_finite() {
        return f64::nan().into();
    }

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

    pp = p.d_qiv(log); // exp(p) iff log_p
    if pp > 1.0 - f64::EPSILON {
        let ret = if lower_tail { f64::infinity() } else { 0.0 }; // early under/over flow  iff log_p (FIXME)
        return ret.into();
    }

    /* Invert pnchisq(.) :
     * 1.0 finding an upper and lower bound */
    /* This is Pearson's (1959) approximation,
    which is usually good to 4 figs or so.  */
    let mut b = 0.0;
    let mut c = 0.0;
    let mut ff = 0.0;
    b = ncp * ncp / (df + 3.0 * ncp);
    c = (df + 3.0 * ncp) / (df + 2.0 * ncp);
    ff = (df + 2.0 * ncp) / (c * c);
    let chisq = if log {
        log_qchisq(p, ff, lower_tail).unwrap()
    } else {
        qchisq(p, ff, lower_tail).unwrap()
    };
    ux = b + c * chisq;
    if ux <= 0.0 {
        ux = 1.0
    }
    ux0 = ux;

    if !lower_tail && ncp >= 80.0 {
        /* in this case, pnchisq() works via lower_tail = true */
        if pp < 1e-10 {
            warn!("full precision may not have been achieved in qnchisq");
        }
        if log {
            p = -p.exp_m1();
        } else {
            p = (0.5 - p) + 0.5;
        }
        lower_tail = true
    } else {
        p = pp
    }

    pp = (p * (1.0 + Eps)).min(1.0 - 2.2204460492503131e-16);
    if lower_tail {
        while ux < f64::MAX && pnchisq_raw(ux, df, ncp, Eps, rEps, 10000, true, false) < pp {
            ux *= 2.0
        }
        pp = p * (1.0 - Eps);
        lx = ux0.min(f64::MAX);
        while lx > 2.2250738585072014e-308
            && pnchisq_raw(lx, df, ncp, Eps, rEps, 10000, true, false) > pp
        {
            lx *= 0.5
        }
    } else {
        while ux < f64::MAX && pnchisq_raw(ux, df, ncp, Eps, rEps, 10000, false, false) > pp {
            ux *= 2.0
        }
        pp = p * (1.0 - Eps);
        lx = ux0.min(f64::MAX);
        while lx > 2.2250738585072014e-308
            && pnchisq_raw(lx, df, ncp, Eps, rEps, 10000, false, false) < pp
        {
            lx *= 0.5
        }
    }

    /* 2. interval (lx,ux)  halving : */
    if lower_tail {
        loop {
            nx = 0.5 * (lx + ux);
            if pnchisq_raw(nx, df, ncp, accu, racc, 100000, true, false) > p {
                ux = nx
            } else {
                lx = nx
            }
            if !((ux - lx) / nx > accu) {
                break;
            }
        }
    } else {
        loop {
            nx = 0.5 * (lx + ux);
            if pnchisq_raw(nx, df, ncp, accu, racc, 100000, false, false) < p {
                ux = nx
            } else {
                lx = nx
            }
            if !((ux - lx) / nx > accu) {
                break;
            }
        }
    }
    (0.5 * (ux + lx)).into()
}