r2rs-nmath 0.1.1

// Translation of nmath's ptukey
//
// 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 nonstdfloat::f128;
use num_traits::{Float, ToPrimitive, Zero};
use strafe_type::{FloatConstraint, Positive64, Real64};

use crate::{distribution::norm::pnorm, func::lgamma, traits::DPQ};

/// This function calculates probability integral of Hartley's
/// form of the range.
///
/// * w: value of range
/// * rr: no. of rows or groups
/// * cc: no. of columns or treatments
/// * ir: error flag = 1 if pr_w probability > 1
/// * pr_w: returned probability integral from (0, w)
///
/// program will not terminate if ir is raised.
///
/// * bb: upper limit of legendre integration
/// * iMax: maximum acceptable value of integral
/// * nleg: order of legendre quadrature
/// * ihalf: int ((nleg + 1) / 2)
/// * wlar: value of range above which wincr1 intervals are used to
/// calculate second part of integral, else wincr2 intervals are used.
/// * C1, C2, C3: values which are used as cutoffs for terminating or modifying a calculation.
///
/// looks like this is suboptimal for double precision.
/// (see how C1-C3 are used) \<MM\>
fn wprob(w: f64, rr: f64, cc: f64) -> f64 {
    /* const double iMax  = 1.0; not used if = 1*/
    static C1: f64 = -30.0;
    static C2: f64 = -50.0;
    static C3: f64 = 60.0;
    static bb: f64 = 8.0;
    static wlar: f64 = 3.0;
    static wincr1: f64 = 2.0;
    static wincr2: f64 = 3.0;
    static xleg: [f64; 6] = [
        0.981560634246719250690549090149,
        0.904117256370474856678465866119,
        0.769902674194304687036893833213,
        0.587317954286617447296702418941,
        0.367831498998180193752691536644,
        0.125233408511468915472441369464,
    ];
    static aleg: [f64; 6] = [
        0.047175336386511827194615961485,
        0.106939325995318430960254718194,
        0.160078328543346226334652529543,
        0.203167426723065921749064455810,
        0.233492536538354808760849898925,
        0.249147045813402785000562436043,
    ];
    let mut a = 0.0;
    let mut ac = 0.0;
    let mut pr_w = 0.0;
    let mut b = 0.0;
    let mut binc = 0.0;
    let mut c = 0.0;
    let mut cc1 = 0.0;
    let mut pminus = 0.0;
    let mut pplus = 0.0;
    let mut qexpo = 0.0;
    let mut qsqz = 0.0;
    let mut rinsum = 0.0;
    let mut wi = 0.0;
    let mut wincr = 0.0;
    let mut xx = 0.0;
    let mut blb = f128::zero();
    let mut bub = f128::zero();
    let mut einsum = f128::zero();
    let mut elsum = f128::zero();
    let mut j = 0;
    let mut jj = 0;

    qsqz = w * 0.5;

    /* if w >= 16 then the integral lower bound (occurs for c=20) */
    /* is 0.99999999999995 so return a value of 1.0 */

    if qsqz >= bb {
        return 1.0;
    }

    /* find (f(w/2) - 1) ^ cc */
    /* (first term in integral of hartley's form). */

    pr_w = 2.0 * pnorm(qsqz, 0.0, 1.0, true).unwrap() - 1.0; /* erf(qsqz / M_SQRT2) */
    /* if pr_w ^ cc < 2e-22 then set pr_w = 0 */
    if pr_w >= (C2 / cc).exp() {
        pr_w = pr_w.powf(cc)
    } else {
        pr_w = 0.0
    }

    /* if w is large then the second component of the */
    /* integral is small, so fewer intervals are needed. */

    if w > wlar {
        wincr = wincr1
    } else {
        wincr = wincr2
    }

    /* find the integral of second term of hartley's form */
    /* for the integral of the range for equal-length */
    /* intervals using legendre quadrature.  limits of */
    /* integration are from (w/2, 8).  two or three */
    /* equal-length intervals are used. */

    /* blb and bub are lower and upper limits of integration. */

    blb = f128::new(qsqz);
    binc = (bb - qsqz) / wincr;
    bub = blb + f128::new(binc);
    einsum = f128::new(0.0);

    /* integrate over each interval */

    cc1 = cc - 1.0;
    wi = 1.0;
    while wi <= wincr {
        elsum = f128::new(0.0);
        a = (f128::new(0.5) * (bub + blb)).to_f64().unwrap();

        /* legendre quadrature with order = nleg */

        b = (f128::new(0.5) * (bub - blb)).to_f64().unwrap();

        jj = 1;
        while jj <= 12 {
            if (6) < jj {
                j = 12 - jj + 1;
                xx = xleg[j - 1]
            } else {
                j = jj;
                xx = -xleg[j - 1]
            }
            c = b * xx;
            ac = a + c;

            /* if (-qexpo/2).exp() < 9e-14, */
            /* then doesn't contribute to integral */

            qexpo = ac * ac;
            if qexpo > C3 {
                break;
            }

            pplus = 2.0 * pnorm(ac, 0.0, 1.0, true).unwrap();
            pminus = 2.0 * pnorm(ac, w, 1.0, true).unwrap();

            /* if rinsum ^ (cc-1) < 9e-14, */
            /* then doesn't contribute to integral */

            rinsum = pplus * 0.5 - pminus * 0.5;
            if rinsum >= (C1 / cc1).exp() {
                rinsum = aleg[j - 1] * (-(0.5 * qexpo)).exp() * rinsum.powf(cc1);
                elsum += f128::new(rinsum)
            }
            jj += 1
        }
        elsum *= f128::new(2.0 * b * cc * strafe_consts::_1DSQRT_2TPI);
        einsum += elsum;
        blb = bub;
        bub += f128::new(binc);
        wi += 1.
    }

    /* if pr_w ^ rr < 9e-14, then return 0 */
    pr_w += einsum.to_f64().unwrap();
    if pr_w <= (C1 / rr).exp() {
        return 0.0;
    }

    pr_w = pr_w.powf(rr);
    if pr_w >= 1.0 {
        /* 1 was iMax was eps */
        return 1.0;
    }
    return pr_w;
}

/// Computes the probability that 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.
///
/// The algorithm is based on that of the reference.
///
/// * q: value of studentized range
/// * rr: no. of rows or groups
/// * cc: no. of columns or treatments
/// * df: degrees of freedom of error term
/// * ir\[0\]: error flag = 1 if wprob probability > 1
/// * ir\[1\]: error flag = 1 if qprob probability > 1
pub fn ptukey<R: Into<Real64>, P1: Into<Positive64>, P2: Into<Positive64>, P3: Into<Positive64>>(
    q: R,
    rr: P1,
    cc: P2,
    df: P3,
    lower_tail: bool,
    log: bool,
) -> Real64 {
    let q = q.into().unwrap();
    let rr = rr.into().unwrap();
    let cc = cc.into().unwrap();
    let df = df.into().unwrap();

    /*  const double eps = 1.0; not used if = 1 */
    static eps1: f64 = -30.0;
    static eps2: f64 = 1.0e-14;
    static dhaf: f64 = 100.0;
    static dquar: f64 = 800.0;
    static deigh: f64 = 5000.0;
    static dlarg: f64 = 25000.0;
    static ulen1: f64 = 1.0;
    static ulen2: f64 = 0.5;
    static ulen3: f64 = 0.25;
    static ulen4: f64 = 0.125;
    static xlegq: [f64; 8] = [
        0.989400934991649932596154173450,
        0.944575023073232576077988415535,
        0.865631202387831743880467897712,
        0.755404408355003033895101194847,
        0.617876244402643748446671764049,
        0.458016777657227386342419442984,
        0.281603550779258913230460501460,
        0.950125098376374401853193354250e-1,
    ];
    static alegq: [f64; 8] = [
        0.271524594117540948517805724560e-1,
        0.622535239386478928628438369944e-1,
        0.951585116824927848099251076022e-1,
        0.124628971255533872052476282192,
        0.149595988816576732081501730547,
        0.169156519395002538189312079030,
        0.182603415044923588866763667969,
        0.189450610455068496285396723208,
    ];
    let mut ans = 0.0;
    let mut f2 = 0.0;
    let mut f21 = 0.0;
    let mut f2lf = 0.0;
    let mut ff4 = 0.0;
    let mut otsum = 0.0;
    let mut qsqz = 0.0;
    let mut rotsum = 0.0;
    let mut t1 = 0.0;
    let mut twa1 = 0.0;
    let mut ulen = 0.0;
    let mut wprb = 0.0;
    let mut i = 0;
    let mut j = 0;
    let mut jj = 0;

    if q.is_nan() || rr.is_nan() || cc.is_nan() || df.is_nan() {
        return f64::nan().into();
    }

    if q <= 0.0 {
        return f64::dt_0(lower_tail, log).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 !q.is_finite() {
        return f64::dt_1(lower_tail, log).into();
    }

    if df > dlarg {
        return wprob(q, rr, cc.dt_val(lower_tail, log)).into();
    }

    /* calculate leading constant */

    f2 = df * 0.5;
    /* lgammafn(u) = gamma(u).ln() */
    f2lf = f2 * df.ln() - df * std::f64::consts::LN_2 - lgamma(f2).unwrap();
    f21 = f2 - 1.0;

    /* integral is divided into unit, half-unit, quarter-unit, or */
    /* eighth-unit length intervals depending on the value of the */
    /* degrees of freedom. */

    ff4 = df * 0.25;
    if df <= dhaf {
        ulen = ulen1
    } else if df <= dquar {
        ulen = ulen2
    } else if df <= deigh {
        ulen = ulen3
    } else {
        ulen = ulen4
    }

    f2lf += ulen.ln();

    /* integrate over each subinterval */

    ans = 0.0;

    i = 1;
    while i <= 50 {
        otsum = 0.0;

        /* legendre quadrature with order = nlegq */
        /* nodes (stored in xlegq) are symmetric around zero. */

        twa1 = (2 * i - 1) as f64 * ulen;

        jj = 1;
        while jj <= 16 {
            if (8) < jj {
                j = jj - 8 - 1;
                t1 = f2lf + f21 * (twa1 + xlegq[j] * ulen).ln() - (xlegq[j] * ulen + twa1) * ff4
            } else {
                j = jj - 1;
                t1 = f2lf + f21 * (twa1 - xlegq[j] * ulen).ln() + (xlegq[j] * ulen - twa1) * ff4
            }

            /* if t1.exp() < 9e-14, then doesn't contribute to integral */
            if t1 >= eps1 {
                if (8) < jj {
                    qsqz = q * ((xlegq[j] * ulen + twa1) * 0.5).sqrt()
                } else {
                    qsqz = q * ((-(xlegq[j] * ulen) + twa1) * 0.5).sqrt()
                }

                /* call wprob to find integral of range portion */

                wprb = wprob(qsqz, rr, cc);
                rotsum = wprb * alegq[j] * t1.exp();
                otsum += rotsum
            }
            /* end legendre integral for interval i */
            /* L200: */
            jj += 1
        }

        /* if integral for interval i < 1e-14, then stop.
         * However, in order to avoid small area under left tail,
         * at least  1 / ulen  intervals are calculated.
         */
        if i as f64 * ulen >= 1.0 && otsum <= eps2 {
            break;
        }

        /* end of interval i */
        /* L330: */

        ans += otsum;
        i += 1
    }

    if otsum > eps2 {
        /* not converged */
        warn!("full precision may not have been achieved in ptukey");
    }

    if ans > 1.0 {
        ans = 1.0
    }
    ans.dt_val(lower_tail, log).into()
}