r2rs-nmath 0.1.1

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

use crate::{distribution::tukey::ptukey, traits::DPQ};

/// This function finds percentage point of the studentized range
/// which is used as initial estimate for the secant method.
/// function is adapted from portion of algorithm as 70
/// from applied statistics (1974) ,vol. 23, no. 1
/// by odeh, r. e. and evans, j. o.
///
/// p = percentage point
/// c = no. of columns or treatments
/// v = degrees of freedom
/// qinv = returned initial estimate
///
/// vmax is cutoff above which degrees of freedom
/// is treated as infinity.
fn qinv(p: f64, c: f64, v: f64) -> f64 {
    static p0: f64 = 0.322232421088;
    static q0: f64 = 0.993484626060e-01;
    static p1: f64 = -1.0;
    static q1: f64 = 0.588581570495;
    static p2: f64 = -0.342242088547;
    static q2: f64 = 0.531103462366;
    static p3: f64 = -0.204231210125;
    static q3: f64 = 0.103537752850;
    static p4: f64 = -0.453642210148e-04;
    static q4: f64 = 0.38560700634e-02;
    static c1: f64 = 0.8832;
    static c2: f64 = 0.2368;
    static c3: f64 = 1.214;
    static c4: f64 = 1.208;
    static c5: f64 = std::f64::consts::SQRT_2;
    static vmax: f64 = 120.0;

    let mut ps = 0.0;
    let mut q = 0.0;
    let mut t = 0.0;
    let mut yi = 0.0;

    ps = 0.5 - 0.5 * p;
    yi = ((1.0 / (ps * ps)).ln()).sqrt();
    t = yi
        + ((((yi * p4 + p3) * yi + p2) * yi + p1) * yi + p0)
            / ((((yi * q4 + q3) * yi + q2) * yi + q1) * yi + q0);
    if v < vmax {
        t += (t * t * t + t) / v / 4.0
    }
    q = c1 - c2 * t;
    if v < vmax {
        q += -c3 / v + c4 * t / v
    }
    return t * (q * (c - 1.0).ln() + c5);
}

/// Computes the quantiles of the maximum of rr studentized
/// ranges, each based on cc means and with df degrees of freedom
/// for the standard error, is less than q.
pub fn qtukey<R: Into<Real64>, P1: Into<Positive64>, P2: Into<Positive64>, P3: Into<Positive64>>(
    p: R,
    rr: P1,
    cc: P2,
    df: P3,
    lower_tail: bool,
    log: bool,
) -> Real64 {
    let mut p = p.into().unwrap();
    let rr = rr.into().unwrap();
    let cc = cc.into().unwrap();
    let df = df.into().unwrap();

    static eps: f64 = 0.0001;
    let maxiter = 50;

    let mut ans = 0.0;
    let mut valx0 = 0.0;
    let mut valx1 = 0.0;
    let mut x0 = 0.0;
    let mut x1 = 0.0;
    let mut xabs = 0.0;
    let mut iter = 0;

    if p.is_nan() || rr.is_nan() || cc.is_nan() || df.is_nan() {
        return (p + rr + cc + df).into();
    }

    /* df must be > 1 ; there must be at least two values */
    if df < 2.0 || rr < 1.0 || cc < 2.0 {
        return f64::nan().into();
    }

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

    p = p.dt_qiv(lower_tail, log);

    /* Initial value */

    x0 = qinv(p, cc, df);

    /* Find prob(value < x0) */

    valx0 = ptukey(x0, rr, cc, df, true, false).unwrap() - p;

    /* Find the second iterate and prob(value < x1). */
    /* If the first iterate has probability value */
    /* exceeding p then second iterate is 1 less than */
    /* first iterate; otherwise it is 1 greater. */

    if valx0 > 0.0 {
        x1 = (x0 - 1.0).max(0.0)
    } else {
        x1 = x0 + 1.0
    }
    valx1 = ptukey(x1, rr, cc, df, true, false).unwrap() - p;

    /* Find new iterate */

    iter = 1;
    while iter < maxiter {
        ans = x1 - valx1 * (x1 - x0) / (valx1 - valx0);
        valx0 = valx1;

        /* New iterate must be >= 0 */

        x0 = x1;
        if ans < 0.0 {
            ans = 0.0;
            valx1 = -p
        }

        /* Find prob(value < new iterate) */

        valx1 = ptukey(ans, rr, cc, df, true, false).unwrap() - p;
        x1 = ans;

        /* If the difference between two successive */
        /* iterates is less than eps, stop */

        xabs = (x1 - x0).abs();
        if xabs < eps {
            return ans.into();
        }

        iter += 1
    }

    /* The process did not converge in 'maxiter' iterations */
    warn!("convergence failed in qtukey");
    ans.into()
}