4 files changed, 45 insertions, 0 deletions
diff --git a/arith.c b/arith.c
index abf9f1be..917973c8 100644
--- a/arith.c
+++ b/arith.c
@@ -2725,6 +2725,30 @@ val cum_norm_dist(val arg)
   }
 }
 
+/*
+ * Source:
+ * Odeh and Evans approximation.
+ * 1974
+ */
+val inv_cum_norm(val arg)
+{
+  val arg_flo = to_float(lit("inv-cum-norm"), arg);
+  double p = c_flo(arg_flo);
+  int is_upper_half = (p >= 0.5);
+  double r = is_upper_half ? 1 - p : p;
+  if (r < 1E-20) {
+    return flo(is_upper_half ? 10 : -10);
+  } else {
+    double y = sqrt(-2*log(r));
+    double a = ((((4.53642210148e-5*y + 0.0204231210245)*y +
+                  0.342242088547)*y + 1)*y + 0.322232431088);
+    double b = ((((0.0038560700634*y + 0.10353775285)*y +
+                  0.531103462366)*y + 0.588581570495)*y + 0.099348462606);
+    double z = y - a / b;
+    return flo(is_upper_half ? z : -z);
+  }
+}
+
 static val rising_product(val m, val n)
 {
   val acc;
diff --git a/arith.h b/arith.h
index 235d4ba1..7b9ecc6b 100644
--- a/arith.h
+++ b/arith.h
@@ -37,6 +37,7 @@ val in_int_ptr_range(val bignum);
 ucnum c_unum(val num);
 val unum(ucnum u);
 val cum_norm_dist(val x);
+val inv_cum_norm(val p);
 val n_choose_k(val n, val k);
 val n_perm_k(val n, val k);
 val divides(val d, val n);
diff --git a/eval.c b/eval.c
index a614335d..1aa56866 100644
--- a/eval.c
+++ b/eval.c
@@ -5903,6 +5903,7 @@ void eval_init(void)
   reg_fun(intern(lit("exp"), user_package), func_n1(expo));
   reg_fun(intern(lit("sqrt"), user_package), func_n1(sqroot));
   reg_fun(intern(lit("cum-norm-dist"), user_package), func_n1(cum_norm_dist));
+  reg_fun(intern(lit("inv-cum-norm"), user_package), func_n1(inv_cum_norm));
   reg_fun(intern(lit("n-choose-k"), user_package), func_n2(n_choose_k));
   reg_fun(intern(lit("n-perm-k"), user_package), func_n2(n_perm_k));
   reg_fun(intern(lit("fixnump"), user_package), func_n1(fixnump));
diff --git a/txr.1 b/txr.1
index 5c7aecef..72ad9949 100644
--- a/txr.1
+++ b/txr.1
@@ -34099,6 +34099,25 @@ distribution function: the integral, of the normal distribution function, from
 negative infinity to the
 .metn argument .
 
+.coNP Function @ inv-cum-norm
+.synb
+.mets (inv-cum-norm << argument )
+.syne
+.desc
+The
+.code inv-cum-norm
+function calculates an approximate to the inverse of the cumulative
+normal distribution function.  The argument, a value expected to lie
+in the range [0, 1], represents the integral of the normal distribution
+function from negative infinity to some domain point
+.IR p .
+The function calculates the approximate value of
+.IR p .
+The minimum value returned is -10, and the maximum value returned is 10,
+regardless of how closely the argument approaches, respectively,
+the 0 or 1 integral endpoints.  For values less than zero, or exceeding 1, the
+values returned, respectively, are -10 and 10.
+
 .coNP Functions @ n-choose-k and @ n-perm-k
 .synb
 .mets (n-choose-k < n << k )