[gnumeric] Math: improve accuracy of dnorm.

[Morten Welinder <mortenw src gnome org>]

commit 42e56b0aceaac31e80985ba42f63d67b75078c86
Author: Morten Welinder <terra gnome org>
Date:   Mon Dec 30 12:30:43 2013 -0500

    Math: improve accuracy of dnorm.
    
    Doing exp(-x*x/2) is not accurate for large x.  Large starts about 5.
    This affects a large number of functions that use dnorm.

 ChangeLog      |    4 ++++
 NEWS           |    1 +
 src/mathfunc.c |   26 ++++++++++++++++++++------
 3 files changed, 25 insertions(+), 6 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index c296ba5..98a200d 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2013-12-30  Morten Welinder  <terra gnome org>
+
+       * src/mathfunc.c (dnorm): Improve accuracy for x>5 (normalized).
+
 2013-12-25  Morten Welinder  <terra gnome org>
 
        * src/item-grid.c (cb_cursor_come_to_rest): Clear ->tip_timer when
diff --git a/NEWS b/NEWS
index dda2157..45c0b5c 100644
--- a/NEWS
+++ b/NEWS
@@ -21,6 +21,7 @@ Morten:
        * Fix timeout related critical in ItemGrid.
        * Fix win32 print crash.  [#719997]
        * Fix solver problem with constraints.
+       * Improve accuracy of r.PNORM for large arguments.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.9
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 6cff799..9fd35a6 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -234,12 +234,26 @@ gnm_float dnorm(gnm_float x, gnm_float mu, gnm_float sigma, gboolean give_log)
        return (x == mu) ? gnm_pinf : R_D__0;
     }
     x = (x - mu) / sigma;
-
-    if(!gnm_finite(x)) return R_D__0;
-    return (give_log ?
-           -(M_LN_SQRT_2PI  +  0.5 * x * x + gnm_log(sigma)) :
-           M_1_SQRT_2PI * gnm_exp(-0.5 * x * x)  /       sigma);
-    /* M_1_SQRT_2PI = 1 / gnm_sqrt(2 * pi) */
+    x = gnm_abs (x);
+
+    if (give_log)
+           return -(M_LN_SQRT_2PI + 0.5 * x * x + gnm_log(sigma));
+    else if (x < 5)
+           return M_1_SQRT_2PI * gnm_exp(-0.5 * x * x) / sigma;
+    else if (x >= 256)
+           return 0;  /* Will underflow anyway. */
+    else {
+           /*
+            * Split x into two parts, x=x1+x2, such that x2 is
+            * small and x1 has less than 26 bits.  That ensures
+            * that x1*x1 is error free.
+            */
+           gnm_float x1 = gnm_floor (x * 65536 + 0.5) / 65536;
+           gnm_float x2 = x - x1;
+           return M_1_SQRT_2PI / sigma *
+                   (gnm_exp(-0.5 * x1 * x1) *
+                    gnm_exp(-(0.5 * x2 + x1) * x2));
+    }
 }
 
 /* ------------------------------------------------------------------------ */