[gnumeric] Stat: improve various q-functions.

[Morten Welinder <mortenw src gnome org>]

commit c0230971a8fa604d54b442c1b5983b4d01c381d0
Author: Morten Welinder <terra gnome org>
Date:   Thu Oct 31 15:01:00 2013 -0400

    Stat: improve various q-functions.

 ChangeLog              |    3 +++
 NEWS                   |    1 +
 plugins/fn-r/ChangeLog |    5 +++++
 plugins/fn-r/extra.c   |   12 ++++++++++++
 samples/amath.gnumeric |  Bin 39204 -> 39480 bytes
 src/mathfunc.c         |   14 +++++++++++++-
 6 files changed, 34 insertions(+), 1 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 67b3525..1a52928 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -4,6 +4,9 @@
        avoids accuracy problems since any cdf is completely flat in terms
        of floating point values when going sufficiently far into the
        tail.
+       (qtukey): Ditto.
+       (qbeta): Ditto.
+       (J_bessel): Extend xlrg_BESS_IJ to match current R.
 
 2013-10-22  Morten Welinder  <terra gnome org>
 
diff --git a/NEWS b/NEWS
index 1153f2a..a0e41b9 100644
--- a/NEWS
+++ b/NEWS
@@ -6,6 +6,7 @@ Morten:
        * Improve accuracy of R.QCAUCHY.
        * Acquire more special function test cases.
        * Improve accuracy of R.QGAMMA and thus R.QCHISQ.
+       * Improve accuracy of R.QBETA, R.QF, R.QTUKEY, R.QSNORM, and R.QST.
 
 Xabier Rodríguez Calvar:
        * Fix dialog button order. [#710378]
diff --git a/plugins/fn-r/ChangeLog b/plugins/fn-r/ChangeLog
index 8bc072a..2b37a01 100644
--- a/plugins/fn-r/ChangeLog
+++ b/plugins/fn-r/ChangeLog
@@ -1,3 +1,8 @@
+2013-10-31  Morten Welinder  <terra gnome org>
+
+       * extra.c (qsnorm, qst): Flip tail if p>0.9.  Improves accuracy
+       when we get really close to 1.
+
 2013-10-30  Morten Welinder  <terra gnome org>
 
        * extra.c (qcauchy): Flip tail if p>0.5.  Improves accuracy when
diff --git a/plugins/fn-r/extra.c b/plugins/fn-r/extra.c
index 1aef7a9..0b800e3 100644
--- a/plugins/fn-r/extra.c
+++ b/plugins/fn-r/extra.c
@@ -136,6 +136,12 @@ qsnorm (gnm_float p, gnm_float shape, gnm_float location, gnm_float scale,
        if (shape == 0.)
                return qnorm (p, location, scale, lower_tail, log_p);
 
+       if (!log_p && p > 0.9) {
+               /* We're far into the tail.  Flip.  */
+               p = 1 - p;
+               lower_tail = !lower_tail;
+       }
+
        x0 = 0.0;
        params[0] = shape;
        params[1] = location;
@@ -293,6 +299,12 @@ qst (gnm_float p, gnm_float n, gnm_float shape,
        if (shape == 0.)
                return qt (p, n, lower_tail, log_p);
 
+       if (!log_p && p > 0.9) {
+               /* We're far into the tail.  Flip.  */
+               p = 1 - p;
+               lower_tail = !lower_tail;
+       }
+
        x0 = 0.0;
        params[0] = n;
        params[1] = shape;
diff --git a/samples/amath.gnumeric b/samples/amath.gnumeric
index d4c756c..7bd084a 100644
Binary files a/samples/amath.gnumeric and b/samples/amath.gnumeric differ
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 5e26bce..1273ddf 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -3814,7 +3814,7 @@ And routine specific :
 #define enten_BESS     1e308
 
 #define exparg_BESS    709.
-#define xlrg_BESS_IJ   1e4
+#define xlrg_BESS_IJ   1e5
 #define xlrg_BESS_Y    1e8
 #define thresh_BESS_Y  16.
 
@@ -6644,6 +6644,12 @@ qtukey (gnm_float p, gnm_float nmeans, gnm_float df, gnm_float nranges, gboolean
 {
        gnm_float x0, shape[3];
 
+       if (!log_p && p > 0.9) {
+               /* We're far into the tail.  Flip.  */
+               p = 1 - p;
+               lower_tail = !lower_tail;
+       }
+
        /* This is accurate for nmeans==2 and nranges==1.   */
        x0 = M_SQRT2gnum * qt ((1 + p) / 2, df, lower_tail, log_p);
 
@@ -7832,6 +7838,12 @@ qbeta (gnm_float p, gnm_float pin, gnm_float qin, gboolean lower_tail, gboolean
 
        if (pin < 0. || qin < 0.) ML_ERR_return_NAN;
 
+       if (!log_p && p > 0.9) {
+               /* We're far into the tail.  Flip.  */
+               p = 1 - p;
+               lower_tail = !lower_tail;
+       }
+
        if (pin >= 1 && qin >= 1)
                x0 = abramowitz_stegun_26_5_22 (p, pin, qin, lower_tail, log_p);
        else {