[gnumeric] drayleigh: improve accuracy.

[Morten Welinder <mortenw src gnome org>]

commit 686b0624df069505e263a4fd4af9012e2e8ea927
Author: Morten Welinder <terra gnome org>
Date:   Fri Jan 8 13:23:50 2016 -0500

    drayleigh: improve accuracy.

 ChangeLog    |    4 ++++
 src/sf-dpq.c |   16 +++++++++++++++-
 2 files changed, 19 insertions(+), 1 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 714dd47..5c7d5fe 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2016-01-08  Morten Welinder  <terra gnome org>
+
+       * src/sf-dpq.c (drayleigh): Undo last change and improve accuracy.
+
 2016-01-06  Morten Welinder  <terra gnome org>
 
        * src/mathfunc.c (ebd0): Fix problem with overflow.  [#760230]
diff --git a/src/sf-dpq.c b/src/sf-dpq.c
index 6cfe6a9..934e04f 100644
--- a/src/sf-dpq.c
+++ b/src/sf-dpq.c
@@ -639,7 +639,21 @@ qhyper (gnm_float p, gnm_float NR, gnm_float NB, gnm_float n,
 gnm_float
 drayleigh (gnm_float x, gnm_float scale, gboolean give_log)
 {
-       return dweibull (x, 2, M_SQRT2gnum * scale, give_log);
+       // This is tempting, but has lower precision since sqrt(2)
+       // is inexact.
+       //
+       // return dweibull (x, 2, M_SQRT2gnum * scale, give_log);
+
+       if (scale <= 0)
+               return gnm_nan;
+       if (x <= 0)
+               return R_D__0;
+       else {
+               gnm_float p = dnorm (x, 0, scale, give_log);
+               return give_log
+                       ? p + gnm_log (x / scale) + M_LN_SQRT_2PI
+                       : p * x / scale * M_SQRT_2PI;
+       }
 }
 
 /* ------------------------------------------------------------------------ */