r2rs-nmath 0.1.1

// Translation of nmath's dpq
//
// 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;

pub trait DPQ {
    /// 0
    fn d_0(log: bool) -> Self;
    /// 1
    fn d_1(log: bool) -> Self;
    /// 0
    fn dt_0(lower_tail: bool, log: bool) -> Self;
    /// 1
    fn dt_1(lower_tail: bool, log: bool) -> Self;
    /// 1/2 (lower- or upper tail)
    fn d_half(log: bool) -> Self;
    /// Use 0.5 - self + 0.5 to perhaps gain 1 bit of accuracy
    /// self
    fn d_lval(&self, lower_tail: bool) -> Self;
    /// Use 0.5 - self + 0.5 to perhaps gain 1 bit of accuracy
    /// 1 - self
    fn d_cval(&self, lower_tail: bool) -> Self;
    /// self  in pF(self,..)
    fn d_val(&self, log: bool) -> Self;
    /// self  in qF(self,..)
    fn d_qiv(&self, log: bool) -> Self;
    /// exp(self)
    fn d_exp(&self, log: bool) -> Self;
    /// log(self)
    fn d_log(&self, log: bool) -> Self;
    /// log(1-self)
    fn d_clog(&self, log: bool) -> Self;
    /// log(1 - exp(self))  in more stable form than log1p(- R_D_qIv(self))
    fn log1_exp(&self) -> Self;
    /// log(1-exp(self)):  R_D_LExp(self) == (log1p(- R_D_qIv(self))) but even more stable
    fn d_lexp(&self, log: bool) -> Self;
    fn dt_val(&self, lower_tail: bool, log: bool) -> Self;
    fn dt_cval(&self, lower_tail: bool, log: bool) -> Self;
    /// self  in qF
    fn dt_qiv(&self, lower_tail: bool, log: bool) -> Self;
    /// 1 - self in qF
    fn dt_civ(&self, lower_tail: bool, log: bool) -> Self;
    /// exp(self)
    fn dt_exp(&self, lower_tail: bool, log: bool) -> Self;
    /// exp(1 - self)
    fn dt_cexp(&self, lower_tail: bool, log: bool) -> Self;
    /// log(self) in qF
    fn dt_log(&self, lower_tail: bool, log: bool) -> Self;
    /// log(1-self) in qF
    fn dt_clog(&self, lower_tail: bool, log: bool) -> Self;
    fn is_non_integer(&self) -> bool;
    fn is_negative_non_integer(&self) -> bool;
    fn q_p01_check(&self, log: bool) -> bool;
    /// Do the boundaries exactly for q() functions
    /// Often  $left = -\infinity$, and very often $right = \infinity$
    fn q_p01_boundaries(
        &self,
        left: Self,
        right: Self,
        lower_tail: bool,
        log: bool,
    ) -> Option<Self>
    where
        Self: Sized;
    fn p_bounds_01(&self, x_min: Self, x_max: Self, lower_tail: bool, log: bool) -> Option<Self>
    where
        Self: Sized;
    fn p_bounds_inf_01(&self, lower_tail: bool, log: bool) -> Option<Self>
    where
        Self: Sized;
    fn d_fexp(&self, f: Self, log: bool) -> Self;
    fn d_rtxp(&self, x: Self, log: bool) -> Self;
    fn pow_di(&self, x: isize) -> Self;
    /// x * 2^exp
    fn ldexp(&self, exp: Self) -> Self;
}

impl DPQ for f64 {
    fn d_0(log: bool) -> Self {
        if log {
            f64::NEG_INFINITY
        } else {
            0.0
        }
    }

    fn d_1(log: bool) -> Self {
        if log {
            0.0
        } else {
            1.0
        }
    }

    fn dt_0(lower_tail: bool, log: bool) -> Self {
        if lower_tail {
            Self::d_0(log)
        } else {
            Self::d_1(log)
        }
    }

    fn dt_1(lower_tail: bool, log: bool) -> Self {
        if lower_tail {
            Self::d_1(log)
        } else {
            Self::d_0(log)
        }
    }

    fn d_half(log: bool) -> Self {
        if log {
            -std::f64::consts::LN_2
        } else {
            0.5
        }
    }

    fn d_lval(&self, lower_tail: bool) -> Self {
        if lower_tail {
            *self
        } else {
            0.5 - *self + 0.5
        }
    }

    fn d_cval(&self, lower_tail: bool) -> Self {
        if lower_tail {
            0.5 - *self + 0.5
        } else {
            *self
        }
    }

    fn d_val(&self, log: bool) -> Self {
        if log {
            self.ln()
        } else {
            *self
        }
    }

    fn d_qiv(&self, log: bool) -> Self {
        if log {
            self.exp()
        } else {
            *self
        }
    }

    fn d_exp(&self, log: bool) -> Self {
        if log {
            *self
        } else {
            self.exp()
        }
    }

    fn d_log(&self, log: bool) -> Self {
        if log {
            *self
        } else {
            self.ln()
        }
    }

    fn d_clog(&self, log: bool) -> Self {
        if log {
            (-self).ln_1p()
        } else {
            0.5 - *self + 0.5
        }
    }

    fn log1_exp(&self) -> Self {
        if *self > -std::f64::consts::LN_2 {
            (-self.exp_m1()).ln()
        } else {
            (-self.exp()).ln_1p()
        }
    }

    fn d_lexp(&self, log: bool) -> Self {
        if log {
            self.log1_exp()
        } else {
            (-self).ln_1p()
        }
    }

    fn dt_val(&self, lower_tail: bool, log: bool) -> Self {
        if lower_tail {
            self.d_val(log)
        } else {
            self.d_clog(log)
        }
    }

    fn dt_cval(&self, lower_tail: bool, log: bool) -> Self {
        if lower_tail {
            self.d_clog(log)
        } else {
            self.d_val(log)
        }
    }

    fn dt_qiv(&self, lower_tail: bool, log: bool) -> Self {
        if log {
            if lower_tail {
                self.exp()
            } else {
                -self.exp_m1()
            }
        } else {
            self.d_lval(lower_tail)
        }
    }

    fn dt_civ(&self, lower_tail: bool, log: bool) -> Self {
        if log {
            if lower_tail {
                -self.exp_m1()
            } else {
                self.exp()
            }
        } else {
            self.d_cval(lower_tail)
        }
    }

    fn dt_exp(&self, lower_tail: bool, log: bool) -> Self {
        self.d_lval(lower_tail).d_exp(log)
    }

    fn dt_cexp(&self, lower_tail: bool, log: bool) -> Self {
        self.d_cval(lower_tail).d_exp(log)
    }

    fn dt_log(&self, lower_tail: bool, log: bool) -> Self {
        if lower_tail {
            self.d_log(log)
        } else {
            self.d_lexp(log)
        }
    }

    fn dt_clog(&self, lower_tail: bool, log: bool) -> Self {
        if lower_tail {
            self.d_lexp(log)
        } else {
            self.d_log(log)
        }
    }

    fn is_non_integer(&self) -> bool {
        (*self - self.round()).abs() > 1e-7 * self.abs().max(1.0)
    }

    fn is_negative_non_integer(&self) -> bool {
        *self < 0.0 || self.is_non_integer()
    }

    fn q_p01_check(&self, log: bool) -> bool {
        (log && *self > 0.0) || (!log && (*self < 0.0 || *self > 1.0))
    }

    fn q_p01_boundaries(&self, left: Self, right: Self, lower_tail: bool, log: bool) -> Option<Self>
    where
        Self: Sized,
    {
        if log {
            if *self > 0.0 {
                Some(f64::nan())
            } else if *self == 0.0 {
                if lower_tail {
                    Some(right)
                } else {
                    Some(left)
                }
            } else if *self == f64::NEG_INFINITY {
                if lower_tail {
                    Some(left)
                } else {
                    Some(right)
                }
            } else {
                None
            }
        } else if *self < 0.0 || *self > 1.0 {
            Some(f64::nan())
        } else if *self == 0.0 {
            if lower_tail {
                Some(left)
            } else {
                Some(right)
            }
        } else if *self == 1.0 {
            if lower_tail {
                Some(right)
            } else {
                Some(left)
            }
        } else {
            None
        }
    }

    fn p_bounds_01(&self, x_min: Self, x_max: Self, lower_tail: bool, log: bool) -> Option<Self>
    where
        Self: Sized,
    {
        if *self <= x_min {
            Some(Self::dt_0(lower_tail, log))
        } else if *self >= x_max {
            Some(Self::dt_1(lower_tail, log))
        } else {
            None
        }
    }

    fn p_bounds_inf_01(&self, lower_tail: bool, log: bool) -> Option<Self>
    where
        Self: Sized,
    {
        if self.is_infinite() {
            if self.is_sign_positive() {
                Some(Self::dt_1(lower_tail, log))
            } else {
                Some(Self::dt_0(lower_tail, log))
            }
        } else {
            None
        }
    }

    fn d_fexp(&self, x: Self, log: bool) -> Self {
        if log {
            -0.5 * self.ln() + x
        } else {
            x.exp() / self.sqrt()
        }
    }

    fn d_rtxp(&self, x: Self, log: bool) -> Self {
        if log {
            -self.ln() + x
        } else {
            x.exp() / self
        }
    }

    fn pow_di(&self, mut n: isize) -> Self {
        let mut x = *self;
        let mut xn = 1.0;

        if x.is_nan() {
            return f64::nan();
        }

        if n != 0 {
            if !x.is_finite() {
                return x.powf(n as f64);
            }

            let is_neg = n < 0;
            if is_neg {
                n = -n
            };
            loop {
                if n & 1 > 0 {
                    xn *= x;
                }
                n >>= 1;
                if n > 0 {
                    x *= x;
                } else {
                    break;
                }
            }
            if is_neg {
                xn = 1. / xn;
            }
        }
        return xn;
    }

    fn ldexp(&self, exp: Self) -> Self {
        self * exp.exp2()
    }
}