use statrs::distribution::{Continuous, Normal, Univariate};

const S2PI: f64 = 2.50662827463100050242_e0;

const P0: [f64; 5] = [
    -5.99633501014107895267_e1,
    9.80010754185999661536_e1,
    -5.66762857469070293439_e1,
    1.39312609387279679503_e1,
    -1.23916583867381258016_e0,
];

const Q0: [f64; 8] = [
    1.95448858338141759834_e0,
    4.67627912898881538453_e0,
    8.63602421390890590575_e1,
    -2.25462687854119370527_e2,
    2.00260212380060660359_e2,
    -8.20372256168333339912_e1,
    1.59056225126211695515_e1,
    -1.18331621121330003142_e0,
];

const P1: [f64; 9] = [
    4.05544892305962419923_e0,
    3.15251094599893866154_e1,
    5.71628192246421288162_e1,
    4.40805073893200834700_e1,
    1.46849561928858024014_e1,
    2.18663306850790267539_e0,
    -1.40256079171354495875_e-1,
    -3.50424626827848203418_e-2,
    -8.57456785154685413611_e-4,
];

const Q1: [f64; 8] = [
    1.57799883256466749731_e1,
    4.53907635128879210584_e1,
    4.13172038254672030440_e1,
    1.50425385692907503408_e1,
    2.50464946208309415979_e0,
    -1.42182922854787788574_e-1,
    -3.80806407691578277194_e-2,
    -9.33259480895457427372_e-4,
];

const P2: [f64; 9] = [
    3.23774891776946035970_e0,
    6.91522889068984211695_e0,
    3.93881025292474443415_e0,
    1.33303460815807542389_e0,
    2.01485389549179081538_e-1,
    1.23716634817820021358_e-2,
    3.01581553508235416007_e-4,
    2.65806974686737550832_e-6,
    6.23974539184983293730_e-9,
];

const Q2: [f64; 8] = [
    6.02427039364742014255_e0,
    3.67983563856160859403_e0,
    1.37702099489081330271_e0,
    2.16236993594496635890_e-1,
    1.34204006088543189037_e-2,
    3.28014464682127739104_e-4,
    2.89247864745380683936_e-6,
    6.79019408009981274425_e-9,
];

fn polevl(x: f64, coef: &[f64], n: usize) -> f64 {
    let mut ans = coef[0];

    for i in 1..n + 1 {
        ans = ans * x + coef[i];
    }

    ans
}

fn p1evl(x: f64, coef: &[f64], n: usize) -> f64 {
    let mut ans = x + coef[0];

    for i in 1..n {
        ans = ans * x + coef[i];
    }

    ans
}

fn ndtri(y0: f64) -> f64 {
    let mut y = y0;

    let code = if y > (1.0 - 0.13533528323661269189) {
        y = 1.0 - y;

        false
    } else {
        true
    };

    if y > 0.13533528323661269189 {
        y = y - 0.5;
        let y2 = y * y;

        return (y + y * (y2 * polevl(y2, &P0, 4) / p1evl(y2, &Q0, 8))) * S2PI;
    }

    let x = (-2.0 * y.ln()).sqrt();
    let z = 1.0 / x;

    let x0 = x - x.ln() / x;
    let x1 = if x < 8.0 {
        z * polevl(z, &P1, 8) / p1evl(z, &Q1, 8)
    } else {
        z * polevl(z, &P2, 8) / p1evl(z, &Q2, 8)
    };

    if code {
        x1 - x0
    } else {
        x0 - x1
    }
}

pub fn cdf(x: f64) -> f64 {
    Normal::new(0.0, 1.0).unwrap().cdf(x)
}

pub fn pdf(x: f64) -> f64 {
    Normal::new(0.0, 1.0).unwrap().pdf(x)
}

pub fn icdf(x: f64) -> f64 {
    ndtri(x)
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_cdf() {
        assert_eq!(cdf(0.5), 0.6914624612740131);
    }
}