r2rs-nmath 0.1.1

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

use crate::{
    distribution::{chisq::dchisq, pois::dpois_raw},
    traits::DPQ,
};

/// The density of the noncentral chi-squared distribution with "df"
/// degrees of freedom and noncentrality parameter "ncp".
pub fn dnchisq<RE: Into<Real64>, RA: Into<Rational64>, P: Into<Positive64>>(
    x: RE,
    df: RA,
    ncp: P,
    log: bool,
) -> Real64 {
    let x = x.into().unwrap();
    let mut df = df.into().unwrap();
    let ncp = ncp.into().unwrap();

    static eps: f64 = 5e-15;

    let mut i = 0.0;
    let mut ncp2 = 0.0;
    let mut q = 0.0;
    let mut mid = 0.0;
    let mut dfmid = 0.0;
    let mut imax = 0.0;
    let mut sum = f128::zero();
    let mut term = f128::zero();

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

    if x < 0.0 {
        return f64::d_0(log).into();
    }
    if x == 0.0 && df < 2.0 {
        return f64::infinity().into();
    }
    if ncp == 0.0 {
        return dchisq(x, df, log);
    }
    if x == f64::infinity() {
        return f64::d_0(log).into();
    }

    ncp2 = 0.5 * ncp;

    /* find max element of sum */
    imax = ((-(2.0 + df) + ((2.0 - df) * (2.0 - df) + 4.0 * ncp * x).sqrt()) / 4.0).ceil(); /* imax = Inf */
    if imax < 0.0 {
        imax = 0.0
    }
    if imax.is_finite() {
        dfmid = df + 2.0 * imax;
        mid = dpois_raw(imax, ncp2, false) * dchisq(x, dfmid, false).unwrap()
    } else {
        mid = 0.0
    }

    if mid == 0.0 {
        /* underflow to 0 -- maybe numerically correct; maybe can be more accurate,
         * particularly when  give_log = true */
        /* Use  central-chisq approximation formula when appropriate;
         * ((FIXME: the optimal cutoff also depends on (x,df);  use always here? )) */
        return if log || ncp > 1000.0 {
            let nl = df + ncp; /* = "1/(1+b)" Abramowitz & St.*/
            let ic = nl / (nl + ncp);
            dchisq(x * ic, nl * ic, log)
        } else {
            f64::d_0(log).into()
        };
    }

    sum = f128::new(mid);

    /* errorbound := term * q / (1-q)  now subsumed in while() / if() below: */
    /* upper tail */
    term = f128::new(mid);
    df = dfmid;
    i = imax;
    let x2 = x * ncp2;
    loop {
        i += 1.0;
        q = x2 / i / df;
        df += 2.0;
        term *= f128::new(q);
        sum += term;
        if !(q >= 1.0
            || term * f128::new(q) > f128::new((1.0 - q) * eps)
            || term > f128::new(1e-10) * sum)
        {
            break;
        }
    }
    /* lower tail */
    term = f128::new(mid);
    df = dfmid;
    i = imax;
    while i != 0.0 {
        df -= 2.0;
        q = i * df / x2;
        i -= 1.0;
        term *= f128::new(q);
        sum += term;
        if q < 1.0 && term * f128::new(q) <= f128::new((1.0 - q) * eps) {
            break;
        }
    }
    sum.to_f64().unwrap().d_val(log).into()
}