[gnumeric] dlnorm: fix and go crazy with accuracy.

[Morten Welinder <mortenw src gnome org>]

commit e6e608d8378ffb7455d36e6eee737ff3980f8953
Author: Morten Welinder <terra gnome org>
Date:   Wed Jan 1 11:41:44 2014 -0500

    dlnorm: fix and go crazy with accuracy.

 ChangeLog     |    4 ++++
 src/numbers.h |    2 ++
 src/sf-dpq.c  |   37 +++++++++++++++++++++++++++++++++----
 3 files changed, 39 insertions(+), 4 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 6e8bb8d..84b60b2 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2014-01-01  Morten Welinder  <terra gnome org>
+
+       * src/sf-dpq.c (dlnorm): Go crazy with accuracy.
+
 2013-12-31  Morten Welinder  <terra gnome org>
 
        * src/sf-dpq.c (pfuncinverter, discpfuncinverter): Extract from
diff --git a/src/numbers.h b/src/numbers.h
index cc0643f..cc9fb11 100644
--- a/src/numbers.h
+++ b/src/numbers.h
@@ -115,6 +115,7 @@ typedef long double gnm_float;
 #define gnm_quad_exp go_quad_expl
 #define gnm_quad_expm1 go_quad_expm1l
 #define gnm_quad_init go_quad_initl
+#define gnm_quad_log go_quad_logl
 #define gnm_quad_mul go_quad_mull
 #define gnm_quad_mul12 go_quad_mul12l
 #define gnm_quad_one go_quad_onel
@@ -224,6 +225,7 @@ typedef double gnm_float;
 #define gnm_quad_exp go_quad_exp
 #define gnm_quad_expm1 go_quad_expm1
 #define gnm_quad_init go_quad_init
+#define gnm_quad_log go_quad_log
 #define gnm_quad_mul go_quad_mul
 #define gnm_quad_mul12 go_quad_mul12
 #define gnm_quad_one go_quad_one
diff --git a/src/sf-dpq.c b/src/sf-dpq.c
index 1359fc2..9dcdba3 100644
--- a/src/sf-dpq.c
+++ b/src/sf-dpq.c
@@ -400,15 +400,44 @@ pnorm2 (gnm_float x1, gnm_float x2)
 gnm_float
 dlnorm (gnm_float x, gnm_float logmean, gnm_float logsd, gboolean give_log)
 {
+       void *state;
+       GnmQuad qx, qlx, qs, qy, qt;
+       static GnmQuad qsqrt2pi;
+       gnm_float r;
+
        if (gnm_isnan (x) || gnm_isnan (logmean) || gnm_isnan (logsd))
                return x + logmean + logsd;
 
        if (logsd <= 0)
                return gnm_nan;
-    
-       return (x > 0)
-               ? dnorm (gnm_log (x), logmean, logsd, give_log)
-               : R_D__0;
+
+       if (x <= 0)
+               return R_D__0;
+
+       state = gnm_quad_start ();
+       if (qsqrt2pi.h == 0)
+               gnm_quad_sqrt (&qsqrt2pi, &gnm_quad_2pi);
+       gnm_quad_init (&qx, x);
+       gnm_quad_log (&qlx, &qx);
+       gnm_quad_init (&qt, logmean);
+       gnm_quad_sub (&qy, &qlx, &qt);
+       gnm_quad_init (&qs, logsd);
+       gnm_quad_div (&qy, &qy, &qs);
+       gnm_quad_mul (&qy, &qy, &qy);
+       qy.h *= -0.5; qy.l *= -0.5;
+       gnm_quad_mul (&qt, &qs, &qx);
+       gnm_quad_mul (&qt, &qt, &qsqrt2pi);
+       if (give_log) {
+               gnm_quad_log (&qt, &qt);
+               gnm_quad_sub (&qy, &qy, &qt);
+       } else {
+               gnm_quad_exp (&qy, NULL, &qy);
+               gnm_quad_div (&qy, &qy, &qt);
+       }
+       r = gnm_quad_value (&qy);
+       gnm_quad_end (state);
+
+       return r;
 }
 
 /* ------------------------------------------------------------------------ */