From c1535145101cf8758f3716c2e4fd3ddaa722f21c Mon Sep 17 00:00:00 2001
From: Kaz Kylheku <kaz@kylheku.com>
Date: Sat, 11 Jan 2014 20:45:37 -0800
Subject: * arith.c (to_float): Print function name as ~a rather than ~s, so it
 doesn't have quotes around it. (cum_norm_dist): New function.

* arith.h (cum_norm_dist): Declared.

* eval.c: Include arith.h.
(eval_init): Register cum_norm_dist as intrinsic.

* txr.1: Documented cum-norm-dist.
---
 ChangeLog | 13 +++++++++++++
 arith.c   | 56 +++++++++++++++++++++++++++++++++++++++++++++++++++++++-
 arith.h   |  1 +
 eval.c    |  2 ++
 txr.1     | 14 ++++++++++++++
 5 files changed, 85 insertions(+), 1 deletion(-)

diff --git a/ChangeLog b/ChangeLog
index 32199aa9..40c1b7e2 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,16 @@
+2014-01-11  Kaz Kylheku  <kaz@kylheku.com>
+
+	* arith.c (to_float): Print function name as ~a rather than ~s,
+	so it doesn't have quotes around it.
+	(cum_norm_dist): New function.
+
+	* arith.h (cum_norm_dist): Declared.
+
+	* eval.c: Include arith.h.
+	(eval_init): Register cum_norm_dist as intrinsic.
+
+	* txr.1: Documented cum-norm-dist.
+
 2014-01-10  Kaz Kylheku  <kaz@kylheku.com>
 
 	* configure: Detect platforms which don't reveal declarations
diff --git a/arith.c b/arith.c
index 27605034..445cebda 100644
--- a/arith.c
+++ b/arith.c
@@ -927,7 +927,7 @@ static val to_float(val func, val num)
   case FLNUM:
     return num;
   default:
-    uw_throwf(error_s, lit("~s: invalid operand ~s"), func, num, nao);
+    uw_throwf(error_s, lit("~a: invalid operand ~s"), func, num, nao);
   }
 }
 
@@ -1815,6 +1815,60 @@ val maskv(val bits)
   return accum;
 }
 
+/*
+ * Source:
+ * Better Approximations to Cumulative Normal Functions
+ * Graeme West
+ * 2009
+ */
+val cum_norm_dist(val arg)
+{
+  val arg_flo = to_float(lit("cum-norm-dist"), arg);
+  double x = c_flo(arg_flo);
+  double xabs = fabs(x);
+
+  if (xabs > 37.0) {
+    return flo(1.0);
+  } else {
+    double ex = exp(-(xabs * xabs) / 2.0);
+    double retval, accum;
+
+    if (xabs < 7.07106781186547) {
+      accum = 3.52624965998911E-02 * xabs + 0.700383064443688;
+      accum = accum * xabs + 6.37396220353165;
+      accum = accum * xabs + 33.912866078383;
+      accum = accum * xabs + 112.079291497871;
+      accum = accum * xabs + 221.213596169931;
+      accum = accum * xabs + 220.206867912376;
+
+      retval = ex * accum;
+
+      accum = 8.83883476483184E-02 * xabs + 1.75566716318264;
+      accum = accum * xabs + 16.064177579207;
+      accum = accum * xabs + 86.7807322029461;
+      accum = accum * xabs + 296.564248779674;
+      accum = accum * xabs + 637.333633378831;
+      accum = accum * xabs + 793.826512519948;
+      accum = accum * xabs + 440.413735824752;
+
+      retval /= accum;
+    } else {
+      accum = xabs + 0.65;
+      accum = xabs + 4.0 / accum;
+      accum = xabs + 3.0 / accum;
+      accum = xabs + 2.0 / accum;
+      accum = xabs + 1.0 / accum;
+
+      retval = ex / accum / 2.506628274631;
+    }
+
+    if (x > 0)
+      retval = 1.0 - retval;
+
+    return flo(retval);
+  }
+}
+
 void arith_init(void)
 {
   mp_init(&NUM_MAX_MP);
diff --git a/arith.h b/arith.h
index 79dcc327..ef341707 100644
--- a/arith.h
+++ b/arith.h
@@ -30,4 +30,5 @@ val bignum_from_long(long l);
 int highest_bit(int_ptr_t n);
 val normalize(val bignum);
 val in_int_ptr_range(val bignum);
+val cum_norm_dist(val x);
 void arith_init(void);
diff --git a/eval.c b/eval.c
index b8e2155c..12f3ebd0 100644
--- a/eval.c
+++ b/eval.c
@@ -46,6 +46,7 @@
 #endif
 #include "lib.h"
 #include "gc.h"
+#include "arith.h"
 #include "signal.h"
 #include "unwind.h"
 #include "regex.h"
@@ -2351,6 +2352,7 @@ void eval_init(void)
   reg_fun(intern(lit("log"), user_package), func_n1(loga));
   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("fixnump"), user_package), func_n1(fixnump));
   reg_fun(intern(lit("bignump"), user_package), func_n1(bignump));
   reg_fun(intern(lit("floatp"), user_package), func_n1(floatp));
diff --git a/txr.1 b/txr.1
index c52c4161..20d871ef 100644
--- a/txr.1
+++ b/txr.1
@@ -8982,6 +8982,20 @@ must be positive.
 The return value is <base> raised to <exponent>, and reduced to the
 least positive residue modulo <modulus>.
 
+.SS Function cum-norm-dist
+
+.TP
+Syntax:
+
+  (cum-norm-dist <argument>)
+
+.TP
+Description:
+
+The cum-norm-dist function calculates an approximation to the cumulative normal
+distribution function: the integral, of the normal distribution function, from
+negative infinity to the <argument>.
+
 .SS Functions fixnump, bignump, integerp, floatp, numberp
 
 .TP
-- 
cgit v1.2.3