[gnumeric] R.QHYPER: Fix problem at left edge of support.

[Morten Welinder <mortenw src gnome org>]

commit 0746e83fb186bb09ef14ac467a5d1a637becc25b
Author: Morten Welinder <terra gnome org>
Date:   Thu Dec 3 10:38:14 2015 -0500

    R.QHYPER: Fix problem at left edge of support.

 ChangeLog    |    5 +++
 NEWS         |    1 +
 src/sf-dpq.c |   24 +++++++++++--
 src/sstest.c |  103 ++++++++++++++++++++++++++++++++++++++--------------------
 4 files changed, 94 insertions(+), 39 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 507d828..3b5ed07 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2015-12-03  Morten Welinder  <terra gnome org>
+
+       * src/sf-dpq.c (discpfuncinverter): Fix problem at left edge of
+       support.  Fixes R.QHYPER(0.1,3,99,13)
+
 2015-10-19  Morten Welinder  <terra gnome org>
 
        * src/func-builtin.c (gnumeric_table): Make sure to invalidate
diff --git a/NEWS b/NEWS
index 63b3e5f..c738b86 100644
--- a/NEWS
+++ b/NEWS
@@ -6,6 +6,7 @@ Morten:
        * Add R.DGUMBEL, R.PGUMBEL, R.QGUMBEL.
        * Fix conditional format problem.  [#750271]
        * Test also RANDWEIBULL and RANDLOGNORM.
+       * Fix problem with R.QHYPER
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.24
diff --git a/src/sf-dpq.c b/src/sf-dpq.c
index 6bf9206..adf3f89 100644
--- a/src/sf-dpq.c
+++ b/src/sf-dpq.c
@@ -221,6 +221,7 @@ discpfuncinverter (gnm_float p, const gnm_float shape[],
 {
        gboolean have_xlow = gnm_finite (xlow);
        gboolean have_xhigh = gnm_finite (xhigh);
+       gboolean check_left = TRUE;
        gnm_float step;
        int i;
 
@@ -253,16 +254,31 @@ discpfuncinverter (gnm_float p, const gnm_float shape[],
                g_printerr ("x=%.20g  e=%.20g\n", x0, ex0);
 #endif
                if (!lower_tail) ex0 = -ex0;
-               if (ex0 <= 0)
-                       xlow = x0, have_xlow = TRUE;
-               if (ex0 >= 0)
+
+               if (ex0 == 0)
+                       return x0;
+               else if (ex0 < 0)
+                       xlow = x0, have_xlow = TRUE, check_left = FALSE;
+               else if (ex0 > 0)
                        xhigh = x0, have_xhigh = TRUE, step = -gnm_abs (step);
 
                if (i > 1 && have_xlow && have_xhigh) {
                        gnm_float xmid = gnm_floor ((xlow + xhigh) / 2);
                        if (xmid - xlow < 0.5 ||
-                           xmid - xlow < gnm_abs (xlow) * GNM_EPSILON)
+                           xmid - xlow < gnm_abs (xlow) * GNM_EPSILON) {
+                               if (check_left) {
+                                       /*
+                                        * The lower edge of the support might
+                                        * have a probability higher than what
+                                        * we are looking for.
+                                        */
+                                       gnm_float e = pfunc (xlow, shape, lower_tail, log_p) - p;
+                                       if (!lower_tail) e = -e;
+                                       if (e >= 0)
+                                               return xhigh = xlow;
+                               }
                                return xhigh;
+                       }
                        x0 = xmid;
                } else {
                        gnm_float x1 = x0 + step;
diff --git a/src/sstest.c b/src/sstest.c
index ddd9ed5..fd8b198 100644
--- a/src/sstest.c
+++ b/src/sstest.c
@@ -563,6 +563,22 @@ rand_fractile_test (gnm_float const *vals, int N, int nf,
                g_printerr (".\n");
        }
 
+       for (i = 1; i < nf - 1; i++) {
+               if (!(fractiles[i] <= fractiles[i + 1])) {
+                       g_printerr ("Severe fractile ordering problem.\n");
+                       return FALSE;
+               }
+
+               if (probs && !(probs[i] <= probs[i + 1])) {
+                       g_printerr ("Severe cumulative probabilities ordering problem.\n");
+                       return FALSE;
+               }
+       }
+       if (probs && (probs[1] < 0 || probs[nf - 1] > 1)) {
+               g_printerr ("Severe cumulative probabilities range problem.\n");
+               return FALSE;
+       }
+
        for (i = 0; i <= nf; i++)
                fractilecount[i] = 0;
 
@@ -1002,6 +1018,8 @@ test_random_randnorm (int N)
        char *expr;
        gnm_float T;
        int i;
+       gnm_float fractiles[10];
+       const int nf = G_N_ELEMENTS (fractiles);
 
        expr = g_strdup_printf ("=RANDNORM(%.10" GNM_FORMAT_g ",%.10" GNM_FORMAT_g ")",
                                mean_target, var_target);
@@ -1029,6 +1047,12 @@ test_random_randnorm (int N)
                ok = FALSE;
        }
 
+       /* Fractile test */
+       for (i = 1; i < nf; i++)
+               fractiles[i] = qnorm (i / (double)nf, mean_target, gnm_sqrt (var_target), TRUE, FALSE);
+       if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
+               ok = FALSE;
+
        if (adtest < 0.05) {
                g_printerr ("Anderson Darling Test rejected [%.10" GNM_FORMAT_g "]\n", adtest);
                ok = FALSE;
@@ -1816,6 +1840,8 @@ test_random_randhyperg (int N)
        gnm_float param_nb = gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 4)));
        gnm_float s = param_nr + param_nb;
        gnm_float param_n = gnm_floor (random_01 () * (s + 1));
+       param_nr = 3, param_nb = 99, param_n = 13;
+
        gnm_float mean_target = param_n * param_nr / s;
        gnm_float var_target = s > 1
                ? mean_target * (param_nb / s) * (s - param_n) / (s - 1)
@@ -2361,47 +2387,54 @@ test_random (void)
        const char *test_name = "test_random";
        const int N = 20000;
        const int High_N = 200000;
+       const char *single = g_getenv ("SSTEST_RANDOM");
 
        mark_test_start (test_name);
 
-       test_random_rand (N);
-       test_random_randuniform (N);
-       test_random_randnorm (High_N);
-       test_random_randsnorm (High_N);
-       test_random_randexp (N);
-       test_random_randgamma (N);
-       test_random_randbeta (N);
-       test_random_randtdist (N);
-       test_random_randfdist (N);
-       test_random_randchisq (N);
-       test_random_randcauchy (N);
-
-       test_random_randbernoulli (N);
-       test_random_randdiscrete (N);
-       test_random_randbinom (N);
-       test_random_randnegbinom (High_N);
-       test_random_randhyperg (High_N);
-       test_random_randbetween (N);
-       test_random_randpoisson (High_N);
-       test_random_randgeom (High_N);
-       test_random_randlog (N);
-       test_random_randweibull (N);
-       test_random_randlognorm (N);
-
+#define CHECK1(NAME,C) \
+       do { if (!single || strcmp(single,#NAME) == 0) test_random_ ## NAME (C); } while (0)
+
+       /* Continuous */
+       CHECK1 (rand, N);
+       CHECK1 (randuniform, N);
+       CHECK1 (randbeta, N);
+       CHECK1 (randcauchy, N);
+       CHECK1 (randchisq, N);
+       CHECK1 (randexp, N);
+       CHECK1 (randfdist, N);
+       CHECK1 (randgamma, N);
+       CHECK1 (randlog, N);
+       CHECK1 (randlognorm, N);
+       CHECK1 (randnorm, High_N);
+       CHECK1 (randsnorm, High_N);
+       CHECK1 (randtdist, N);
+       CHECK1 (randweibull, N);
 #if 0
-       test_random_randexppow (N);
-       test_random_randnormtail (N);
-       test_random_randgumbel (N);
-       test_random_randlandau (N);
-       test_random_randlaplace (N);
-       test_random_randlevy (N);
-       test_random_randlogistic (N);
-       test_random_randpareto (N);
-       test_random_randrayleigh (N);
-       test_random_randrayleightail (N);
-       test_random_randstdist (N);
+       CHECK1 (randexppow, N);
+       CHECK1 (randgumbel, N);
+       CHECK1 (randlandau, N);
+       CHECK1 (randlaplace, N);
+       CHECK1 (randlevy, N);
+       CHECK1 (randlogistic, N);
+       CHECK1 (randnormtail, N);
+       CHECK1 (randpareto, N);
+       CHECK1 (randrayleigh, N);
+       CHECK1 (randrayleightail, N);
+       CHECK1 (randstdist, N);
 #endif
 
+       /* Discrete */
+       CHECK1 (randbernoulli, N);
+       CHECK1 (randbetween, N);
+       CHECK1 (randbinom, N);
+       CHECK1 (randdiscrete, N);
+       CHECK1 (randgeom, High_N);
+       CHECK1 (randhyperg, High_N);
+       CHECK1 (randnegbinom, High_N);
+       CHECK1 (randpoisson, High_N);
+
+#undef CHECK1
+
        mark_test_end (test_name);
 }