It works!

src/math.rs

@@ -1,5 +1,130 @@
 use statrs::distribution::{Normal, Univariate, Continuous};
 
+const S2PI: f64 = 2.50662827463100050242E0;
+
+const P0: [f64; 5] = [
+    -5.99633501014107895267E1,
+    9.80010754185999661536E1,
+    -5.66762857469070293439E1,
+    1.39312609387279679503E1,
+    -1.23916583867381258016E0,
+];
+
+const Q0: [f64; 8] = [
+    1.95448858338141759834E0,
+    4.67627912898881538453E0,
+    8.63602421390890590575E1,
+    -2.25462687854119370527E2,
+    2.00260212380060660359E2,
+    -8.20372256168333339912E1,
+    1.59056225126211695515E1,
+    -1.18331621121330003142E0,
+];
+
+const P1: [f64; 9] = [
+    4.05544892305962419923E0,
+    3.15251094599893866154E1,
+    5.71628192246421288162E1,
+    4.40805073893200834700E1,
+    1.46849561928858024014E1,
+    2.18663306850790267539E0,
+    -1.40256079171354495875E-1,
+    -3.50424626827848203418E-2,
+    -8.57456785154685413611E-4,
+];
+
+const Q1: [f64; 8] = [
+    1.57799883256466749731E1,
+    4.53907635128879210584E1,
+    4.13172038254672030440E1,
+    1.50425385692907503408E1,
+    2.50464946208309415979E0,
+    -1.42182922854787788574E-1,
+    -3.80806407691578277194E-2,
+    -9.33259480895457427372E-4,
+];
+
+const P2: [f64; 9] = [
+    3.23774891776946035970E0,
+    6.91522889068984211695E0,
+    3.93881025292474443415E0,
+    1.33303460815807542389E0,
+    2.01485389549179081538E-1,
+    1.23716634817820021358E-2,
+    3.01581553508235416007E-4,
+    2.65806974686737550832E-6,
+    6.23974539184983293730E-9,
+];
+
+const Q2: [f64; 8] = [
+    6.02427039364742014255E0,
+    3.67983563856160859403E0,
+    1.37702099489081330271E0,
+    2.16236993594496635890E-1,
+    1.34204006088543189037E-2,
+    3.28014464682127739104E-4,
+    2.89247864745380683936E-6,
+    6.79019408009981274425E-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 code = 1;
+    let mut y = y0;
+
+    if y > (1.0 - 0.13533528323661269189) {
+        y = 1.0 - y;
+        code = 0;
+    }
+
+    if y > 0.13533528323661269189 {
+        y = y - 0.5;
+        let y2 = y * y;
+        let x = y + y * (y2 * polevl(y2, &P0, 4) / p1evl(y2, &Q0, 8));
+        let x = x * S2PI;
+
+        return x;
+    }
+
+    let x = (-2.0 * y.ln()).sqrt();
+    let x0 = x - x.ln() / x;
+
+    let z = 1.0 / 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)
+    };
+
+    let mut x = x0 - x1;
+
+    if code != 0 {
+        x = -x;
+    }
+
+    x
+}
+
 pub fn cdf(x: f64) -> f64 {
     Normal::new(0.0, 1.0).unwrap().cdf(x)
 }
@@ -8,6 +133,10 @@ 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 {