[gnumeric] dnorm: further minor improvements.

[Morten Welinder <mortenw src gnome org>]

commit bdcba232e6a31007b6fd715e008eb2875fe4f11d
Author: Morten Welinder <terra gnome org>
Date:   Mon Dec 30 13:17:53 2013 -0500

    dnorm: further minor improvements.

 src/mathfunc.c |   14 ++++++++------
 1 files changed, 8 insertions(+), 6 deletions(-)
---
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 9fd35a6..aea76f9 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -236,23 +236,25 @@ gnm_float dnorm(gnm_float x, gnm_float mu, gnm_float sigma, gboolean give_log)
     x = (x - mu) / sigma;
     x = gnm_abs (x);
 
+    if (x >= 2 * gnm_sqrt (GNM_MAX))
+           return R_D__0;
     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.
+            * Split x into two parts, x=x1+x2, such that |x2|<=2^-16.
+            * Assuming that we are using IEEE doubles, that means that
+            * x1*x1 is error free for x<1024 (above which we will underflow
+            * anyway).  If we are not using IEEE doubles then this is
+            * still an improvement over the naive formula.
             */
            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));
+                    gnm_exp((-0.5 * x2 - x1) * x2));
     }
 }