r2rs-nmath 0.1.1

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

use crate::{
    distribution::{
        beta::{log_pbeta, pbeta},
        norm::{log_pnorm, pnorm},
    },
    func::lbeta,
    traits::DPQ,
};

/// return  $ P[ T <= x ] $ where
/// T ~ t_{n}  (t distrib. with n degrees of freedom).
/// --> ./pnt.c for NON-central
pub fn pt<RE: Into<Real64>, RA: Into<Rational64>>(
    x: RE,
    n: RA,
    mut lower_tail: bool,
) -> Probability64 {
    let x = x.into().unwrap();
    let n = n.into().unwrap();

    let mut val = 0.0;
    let mut nx = 0.0;

    if !x.is_finite() {
        return if x < 0.0 {
            f64::dt_0(lower_tail, false)
        } else {
            f64::dt_1(lower_tail, false)
        }
        .into();
    }
    if !n.is_finite() {
        return pnorm(x, 0.0, 1.0, lower_tail);
    }

    nx = 1.0 + x / n * x;
    /* FIXME: This test is probably losing rather than gaining precision,
     * now that pbeta(*, log = true) is much better.
     * Note however that a version of this test *is* needed for x*x > D_MAX */
    if nx > 1e100 {
        /* <==>  x*x > 1e100 * n  */
        /* Danger of underflow. So use Abramowitz & Stegun 26.5.4
           pbeta(z, a, b) ~ z^a(1-z)^b / aB(a,b) ~ z^a / aB(a,b),
           with z = 1/nx,  a = n/2,  b= 1/2 :
        */
        let mut lval = 0.0_f64;
        lval = -0.5 * n * (2.0 * x.abs().ln() - n.ln())
            - lbeta(0.5 * n, 0.5).unwrap()
            - (0.5 * n).ln();
        val = lval.exp();
    } else {
        val = if n > x * x {
            pbeta(x * x / (n + x * x), 0.5, n / 2.0, false)
        } else {
            pbeta(1.0 / nx, n / 2.0, 0.5, true)
        }
        .unwrap()
    }

    /* Use "1 - v"  if lower_tail  and  x > 0 (but not both):*/
    if x <= 0.0 {
        lower_tail = !lower_tail
    }

    val /= 2.0;
    val = val.d_cval(lower_tail);

    val.into()
}

/// return  $ P[ T <= x ] $ where
/// T ~ t_{n}  (t distrib. with n degrees of freedom).
/// --> ./pnt.c for NON-central
pub fn log_pt<RE: Into<Real64>, RA: Into<Rational64>>(
    x: RE,
    n: RA,
    mut lower_tail: bool,
) -> LogProbability64 {
    let x = x.into().unwrap();
    let n = n.into().unwrap();

    let mut val = 0.0;
    let mut nx = 0.0;

    if !x.is_finite() {
        return if x < 0.0 {
            f64::dt_0(lower_tail, true)
        } else {
            f64::dt_1(lower_tail, true)
        }
        .into();
    }
    if !n.is_finite() {
        return log_pnorm(x, 0.0, 1.0, lower_tail);
    }

    nx = 1.0 + x / n * x;
    /* FIXME: This test is probably losing rather than gaining precision,
     * now that pbeta(*, log = true) is much better.
     * Note however that a version of this test *is* needed for x*x > D_MAX */
    if nx > 1e100 {
        /* <==>  x*x > 1e100 * n  */
        /* Danger of underflow. So use Abramowitz & Stegun 26.5.4
           pbeta(z, a, b) ~ z^a(1-z)^b / aB(a,b) ~ z^a / aB(a,b),
           with z = 1/nx,  a = n/2,  b= 1/2 :
        */
        let mut lval = 0.0;
        lval = -0.5 * n * (2.0 * x.abs().ln() - n.ln())
            - lbeta(0.5 * n, 0.5).unwrap()
            - (0.5 * n).ln();
        val = lval;
    } else {
        val = if n > x * x {
            log_pbeta(x * x / (n + x * x), 0.5, n / 2.0, false)
        } else {
            log_pbeta(1.0 / nx, n / 2.0, 0.5, true)
        }
        .unwrap()
    }

    /* Use "1 - v"  if lower_tail  and  x > 0 (but not both):*/
    if x <= 0.0 {
        lower_tail = !lower_tail
    }

    if lower_tail {
        (-0.5 * val.exp()).ln_1p()
    } else {
        val - std::f64::consts::LN_2
    }
    .into() /* = (.5* pbeta(....)).ln() */
}