//! 64 bit float implementation for gamma
//!
//! This implementation is based on the Math.NET library.

// Copyright (c) 2009-2010 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.

use core::f64::NAN;
#[allow(unused_imports)]
use num_traits::Float;

use super::Gamma;

impl Gamma for f64 {
    fn ln_gamma(self) -> Self {
        libm::lgamma(self)
    }
    fn gamma(self) -> Self {
        libm::tgamma(self)
    }
    fn lower_gamma_regularized(self, x: Self) -> Self {
        if self.is_nan() {
            return NAN;
        }
        if x.is_nan() {
            return NAN;
        }
        const EPSILON: f64 = 0.000000000000001;
        const BIG: f64 = 4503599627370496.0;
        const BIG_INV: f64 = 2.22044604925031308085e-16;
        if self < 0.0 {
            return NAN;
        }
        if x < 0.0 {
            return NAN;
        }
        if self < 0.000000000000001 {
            return 1.0;
        }
        if x < 0.000000000000001 {
            return 0.0;
        }
        let ax = self * x.ln() - x - Gamma::ln_gamma(self);
        if ax < -709.78271289338399 {
            return if self < x { 1.0 } else { 0.0 };
        }
        if x <= 1.0 || x <= self {
            let mut r2 = self;
            let mut c2 = 1.0;
            let mut ans2 = 1.0;
            loop {
                r2 += 1.0;
                c2 *= x / r2;
                ans2 += c2;
                if (c2 / ans2) <= EPSILON {
                    break;
                }
            }
            return ax.exp() * ans2 / self;
        }
        let mut c = 0;
        let mut y = 1.0 - self;
        let mut z = x + y + 1.0;
        let mut p3 = 1.0;
        let mut q3 = x;
        let mut p2 = x + 1.0;
        let mut q2 = z * x;
        let mut ans = p2 / q2;
        let mut error;

        loop {
            c += 1;
            y += 1.0;
            z += 2.0;
            let yc = y * (c as f64);
            let p = p2 * z - p3 * yc;
            let q = q2 * z - q3 * yc;
            if q != 0.0 {
                let nextans = p / q;
                error = ((ans - nextans) / nextans).abs();
                ans = nextans;
            } else {
                error = 1.0;
            }
            p3 = p2;
            p2 = p;
            q3 = q2;
            q2 = q;
            if p.abs() > BIG {
                p3 *= BIG_INV;
                p2 *= BIG_INV;
                q3 *= BIG_INV;
                q2 *= BIG_INV;
            }
            if error <= EPSILON {
                break;
            }
        }
        1.0 - ax.exp() * ans
    }
}

#[cfg(test)]
mod tests {
    use crate::gamma::Gamma;
    macro_rules! assert_almost_eq {
        ($a:expr, $b:expr, $prec:expr) => {
            #[allow(trivial_numeric_casts)]
            if (($a - $b) as f64).abs() > $prec {
                panic!(
                    "assertion failed: `abs(left - right) < {:e}`, (left: `{}`, right: `{}`)",
                    $prec, $a, $b
                );
            }
        };
    }

    #[test]
    fn test_gamma_lr() {
        assert!(Gamma::lower_gamma_regularized(f64::NAN, f64::NAN).is_nan());
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(0.1, 1.0),
            0.97587265627367222115949155252812057714751052498477013,
            1e-14
        );
        assert_eq!(
            Gamma::lower_gamma_regularized(0.1, 2.0),
            0.99432617602018847196075251078067514034772764693462125
        );
        assert_eq!(
            Gamma::lower_gamma_regularized(0.1, 8.0),
            0.99999507519205198048686442150578226823401842046310854
        );
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(1.5, 1.0),
            0.42759329552912016600095238564127189392715996802703368,
            1e-13
        );
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(1.5, 2.0),
            0.73853587005088937779717792402407879809718939080920993,
            1e-15
        );
        assert_eq!(
            Gamma::lower_gamma_regularized(1.5, 8.0),
            0.99886601571021467734329986257903021041757398191304284
        );
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(2.5, 1.0),
            0.15085496391539036377410688601371365034788861473418704,
            1e-13
        );
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(2.5, 2.0),
            0.45058404864721976739416885516693969548484517509263197,
            1e-14
        );
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(2.5, 8.0),
            0.99315592607757956900093935107222761316136944145439676,
            1e-15
        );
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(5.5, 1.0),
            0.0015041182825838038421585211353488839717739161316985392,
            1e-15
        );
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(5.5, 2.0),
            0.030082976121226050615171484772387355162056796585883967,
            1e-14
        );
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(5.5, 8.0),
            0.85886911973294184646060071855669224657735916933487681,
            1e-14
        );
        assert_almost_eq!(Gamma::lower_gamma_regularized(100.0, 0.5), 0.0, 1e-188);
        assert_almost_eq!(Gamma::lower_gamma_regularized(100.0, 1.5), 0.0, 1e-141);
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(100.0, 90.0),
            0.1582209891864301681049696996709105316998233457433473,
            1e-13
        );
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(100.0, 100.0),
            0.5132987982791486648573142565640291634709251499279450,
            1e-13
        );
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(100.0, 110.0),
            0.8417213299399129061982996209829688531933500308658222,
            1e-13
        );
        assert_almost_eq!(Gamma::lower_gamma_regularized(100.0, 200.0), 1.0, 1e-14);
        assert_eq!(Gamma::lower_gamma_regularized(500.0, 0.5), 0.0);
        assert_eq!(Gamma::lower_gamma_regularized(500.0, 1.5), 0.0);
        assert_almost_eq!(Gamma::lower_gamma_regularized(500.0, 200.0), 0.0, 1e-70);
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(500.0, 450.0),
            0.0107172380912897415573958770655204965434869949241480,
            1e-14
        );
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(500.0, 500.0),
            0.5059471461707603580470479574412058032802735425634263,
            1e-13
        );
        assert_almost_eq!(
            Gamma::lower_gamma_regularized(500.0, 550.0),
            0.9853855918737048059548470006900844665580616318702748,
            1e-14
        );
        assert_almost_eq!(Gamma::lower_gamma_regularized(500.0, 700.0), 1.0, 1e-15);
        assert_eq!(Gamma::lower_gamma_regularized(1000.0, 10000.0), 1.0);
        assert_eq!(Gamma::lower_gamma_regularized(1e+50, 1e+48), 0.0);
        assert_eq!(Gamma::lower_gamma_regularized(1e+50, 1e+52), 1.0);
    }
}