[gnumeric] R.DNORM, RAYLEIGH: Improve accuracy of far-tail cases.

[Morten Welinder <mortenw src gnome org>]

commit a6f1265d5365e7801babc69d091117ec3e4e2e67
Author: Morten Welinder <terra gnome org>
Date:   Sat Apr 11 15:13:10 2015 -0400

    R.DNORM, RAYLEIGH: Improve accuracy of far-tail cases.

 ChangeLog                   |    5 ++++
 NEWS                        |    2 +
 plugins/fn-stat/functions.c |   16 +-----------
 src/sf-dpq.c                |   55 +++++++++++++++++++++++++++++++++++++++---
 src/sf-dpq.h                |    5 ++++
 5 files changed, 64 insertions(+), 19 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 7259b0f..dd5d618 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2015-04-11  Morten Welinder  <terra gnome org>
+
+       * src/sf-dpq.c (dnorm): Improve accuracy in certain far-tail cases.
+       (drayleigh): Import from fn-stat.  Rename.  Improve accuracy.
+
 2015-04-09  Morten Welinder  <terra gnome org>
 
        * src/sheet-filter.c (filter_expr_eval): Fix UMR in the non-match
diff --git a/NEWS b/NEWS
index 1606555..e5fcdd4 100644
--- a/NEWS
+++ b/NEWS
@@ -31,6 +31,8 @@ Morten:
        * Plug leaks.
        * Fix REPLACEB problem.  [#747210]
        * Fix sheet filter problem.
+       * Minor R.DNORM improvement.
+       * Improve RAYLEIGH accuracy.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.21
diff --git a/plugins/fn-stat/functions.c b/plugins/fn-stat/functions.c
index 6b4e10b..4a3c09c 100644
--- a/plugins/fn-stat/functions.c
+++ b/plugins/fn-stat/functions.c
@@ -4507,27 +4507,13 @@ static GnmFuncHelp const help_rayleigh[] = {
        { GNM_FUNC_HELP_END }
 };
 
-static gnm_float
-random_rayleigh_pdf (gnm_float x, gnm_float sigma)
-{
-       if (sigma <= 0)
-               return gnm_nan;
-
-       if (x < 0)
-               return 0;
-       else {
-               gnm_float u = x / sigma;
-               return (u / sigma) * expmx2h (u);
-       }
-}
-
 static GnmValue *
 gnumeric_rayleigh (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
 {
        gnm_float x     = value_get_as_float (argv[0]);
        gnm_float sigma = value_get_as_float (argv[1]);
 
-       return value_new_float (random_rayleigh_pdf (x, sigma));
+       return value_new_float (drayleigh (x, sigma, FALSE));
 }
 
 /***************************************************************************/
diff --git a/src/sf-dpq.c b/src/sf-dpq.c
index 1a8ad35..46177e4 100644
--- a/src/sf-dpq.c
+++ b/src/sf-dpq.c
@@ -8,6 +8,7 @@
 #define R_DT_0 (lower_tail ? R_D__0 : R_D__1)
 #define R_DT_1 (lower_tail ? R_D__1 : R_D__0)
 #define M_1_SQRT_2PI    GNM_const(0.398942280401432677939946059934)  /* 1/sqrt(2pi) */
+#define M_SQRT_2PI GNM_const(2.506628274631000502415765284811045253006986740609938316629923)
 
 /* ------------------------------------------------------------------------- */
 
@@ -305,18 +306,37 @@ discpfuncinverter (gnm_float p, const gnm_float shape[],
 gnm_float
 dnorm (gnm_float x, gnm_float mu, gnm_float sigma, gboolean give_log)
 {
+       gnm_float x0;
+
        if (gnm_isnan (x) || gnm_isnan (mu) || gnm_isnan (sigma))
                return x + mu + sigma;
        if (sigma < 0)
                return gnm_nan;
 
-       x = gnm_abs ((x - mu) / sigma);
+       /* Center.  */
+       x = x0 = gnm_abs (x - mu);
+       x /= sigma;
 
-       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
+       else if (x > 3 + 2 * gnm_sqrt (gnm_log (GNM_MAX)))
+               /* Far into the tail; x > ~100 for long double  */
+               return 0;
+       else if (x > 4) {
+               /*
+                * Split x into xh+xl such that:
+                * 1) xh*xh is exact
+                * 2) 0 <= xl <= 1/65536
+                * 3) 0 <= xl*xh < ~100/65536
+                */
+               gnm_float xh = gnm_floor (x * 65536) / 65536;  /* At most 24 bits */
+               gnm_float xl = (x0 - xh * sigma) / sigma;
+               return M_1_SQRT_2PI *
+                       gnm_exp (-0.5 * (xh * xh)) *
+                       gnm_exp (-xl * (0.5 * xl + xh)) /
+                       sigma;
+       } else
+               /* Near-center case.  */
                return M_1_SQRT_2PI * expmx2h (x) / sigma;
 }
 
@@ -563,3 +583,30 @@ qhyper (gnm_float p, gnm_float NR, gnm_float NB, gnm_float n,
 }
 
 /* ------------------------------------------------------------------------ */
+/**
+ * drayleigh:
+ * @x: observation
+ * @scale: scale parameter
+ * @give_log: if %TRUE, log of the result will be returned instead
+ *
+ * Returns: density of the Rayleigh distribution.
+ */
+
+gnm_float
+drayleigh (gnm_float x, gnm_float scale, gboolean 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);
+               gnm_float f = M_SQRT_2PI * x / scale;
+               return give_log
+                       ? p + gnm_log (f)
+                       : p * f;
+       }
+}
+
+/* ------------------------------------------------------------------------ */
diff --git a/src/sf-dpq.h b/src/sf-dpq.h
index 0a48c08..b9ed62a 100644
--- a/src/sf-dpq.h
+++ b/src/sf-dpq.h
@@ -42,5 +42,10 @@ gnm_float qcauchy (gnm_float p, gnm_float location, gnm_float scale,
 gnm_float qhyper (gnm_float p, gnm_float r, gnm_float b, gnm_float n, gboolean lower_tail, gboolean log_p);
 
 /* ------------------------------------------------------------------------- */
+/* Rayleigh distribution  */
+
+gnm_float drayleigh (gnm_float x, gnm_float scale, gboolean give_log);
+
+/* ------------------------------------------------------------------------- */
 
 #endif