[gnumeric] R.QGAMMA: Fix accuracy for high p.

[Morten Welinder <mortenw src gnome org>]

commit 660bebbc4f0312d5577e84940bfdcad8c12a5a4a
Author: Morten Welinder <terra gnome org>
Date:   Thu Oct 31 09:58:00 2013 -0400

    R.QGAMMA: Fix accuracy for high p.
    
    This also fixes R.QCHISQ which is defined in terms of R.QGAMMA.
    
    As one goes far into the tail of any cdf, it will become constant (in terms
    of distinct floating point values) for larger and larger stretches.  Anything
    based on inverting the cdf therefore will get accuracy problems.

 ChangeLog      |    7 +++++++
 NEWS           |    1 +
 src/mathfunc.c |    6 ++++++
 3 files changed, 14 insertions(+), 0 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index a4922cd..67b3525 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+2013-10-31  Morten Welinder  <terra gnome org>
+
+       * src/mathfunc.c (qgamma): Flip tail for high values of p.  This
+       avoids accuracy problems since any cdf is completely flat in terms
+       of floating point values when going sufficiently far into the
+       tail.
+
 2013-10-22  Morten Welinder  <terra gnome org>
 
        * src/func.c (function_call_with_exprs): Move flags argument into
diff --git a/NEWS b/NEWS
index 6d6bca1..1153f2a 100644
--- a/NEWS
+++ b/NEWS
@@ -5,6 +5,7 @@ Morten:
        * Fix drop-down sizing (gtk+ regression).  [#710749]
        * Improve accuracy of R.QCAUCHY.
        * Acquire more special function test cases.
+       * Improve accuracy of R.QGAMMA and thus R.QCHISQ.
 
 Xabier Rodríguez Calvar:
        * Fix dialog button order. [#710378]
diff --git a/src/mathfunc.c b/src/mathfunc.c
index ffdf6b6..5e26bce 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -7765,6 +7765,12 @@ qgamma (gnm_float p, gnm_float alpha, gnm_float scale,
        R_Q_P01_check(p);
        if (alpha <= 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;
+       }
+
        /* Make an initial guess, x0, assuming scale==1.  */
        v = 2 * alpha;
        if (v < -1.24 * R_DT_log (p))