[gnumeric] IMIGAMMA: improve coverage

[Morten Welinder <mortenw src gnome org>]

commit d439922c95c733905f56b571164e96d7aa888af3
Author: Morten Welinder <terra gnome org>
Date:   Tue Feb 2 21:34:33 2016 -0500

    IMIGAMMA: improve coverage
    
    Use the asymptotic expansion when |z|>|a| provided we get enough
    precision before the series diverges.

 ChangeLog      |    7 +++
 NEWS           |    2 +-
 src/sf-gamma.c |  127 +++++++++++++++++++++++++++++++++++++++++++++++++++++--
 3 files changed, 130 insertions(+), 6 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index edfd717..088909a 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+2016-02-02  Morten Welinder  <terra gnome org>
+
+       * src/sf-gamma.c (complex_igamma): Try asymptotic expansion.
+       (gamma_error_factor): Extend to all positive numbers.
+       (pochhammer_small_n): Allow any x > 1.
+       (qbetaf): Use pochhammer_small_n as long as a > 1 and |b| < 1.
+
 2016-02-01  Morten Welinder  <terra gnome org>
 
        * configure.ac (yacc, lex): Fail if the required program isn't
diff --git a/NEWS b/NEWS
index cb11dd3..43136b3 100644
--- a/NEWS
+++ b/NEWS
@@ -27,7 +27,7 @@ Morten:
        * Fix problem with database functions.  [#761305]
        * Fix font problem for ssconvert to pdf.  [#761296]
        * Fix bison check.  [#761398]
-
+       * Improve IMIGAMMA coverage.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.26
diff --git a/src/sf-gamma.c b/src/sf-gamma.c
index 93178d0..02f14b8 100644
--- a/src/sf-gamma.c
+++ b/src/sf-gamma.c
@@ -522,7 +522,7 @@ qbetaf (gnm_float a, gnm_float b, GnmQuad *mant, int *exp2)
                b = s;
        }
 
-       if (a > 20 && gnm_abs (b) < 1) {
+       if (a > 1 && gnm_abs (b) < 1) {
                void *state;
                if (qgammaf (b, &mb, &eb))
                        return 1;
@@ -621,7 +621,7 @@ gnm_lbeta3 (gnm_float a, gnm_float b, int *sign)
  *
  *   Gamma(x) = sqrt(2Pi) * x^(x-1/2) * exp(-x) * E(x)
  *
- * x should be >20 and the result is, roughly, 1+1/(12x).
+ * x should be >0 and the result is, roughly, 1+1/(12x).
  */
 static void
 gamma_error_factor (GnmQuad *res, const GnmQuad *x)
@@ -648,21 +648,65 @@ gamma_error_factor (GnmQuad *res, const GnmQuad *x)
                GNM_const(86684309913600.),
                GNM_const(514904800886784000.)
        };
-       GnmQuad zn;
+       GnmQuad zn, xpn;
        int i;
+       gnm_float xval = gnm_quad_value (x);
+       int n;
+
+       g_return_if_fail (xval > 0);
+
+       // We want x >= 20 for the asymptotic expansion
+       n = (xval < 20) ? (int)gnm_floor (21 - xval) : 0;
+       gnm_quad_init (&xpn, n);
+       gnm_quad_add (&xpn, &xpn, x);
 
        gnm_quad_init (&zn, 1);
        *res = zn;
 
        for (i = 0; i < (int)G_N_ELEMENTS (num); i++) {
                GnmQuad t, c;
-               gnm_quad_mul (&zn, &zn, x);
+               gnm_quad_mul (&zn, &zn, &xpn);
                gnm_quad_init (&c, den[i]);
                gnm_quad_mul (&t, &zn, &c);
                gnm_quad_init (&c, num[i]);
                gnm_quad_div (&t, &c, &t);
                gnm_quad_add (res, res, &t);
        }
+
+       if (n > 0) {
+               int i;
+               GnmQuad en, xxn, xph;
+
+               // Gamma(x) = sqrt(2Pi) * x^(x-1/2) * exp(-x) * E(x)
+               // Gamma(x+n) = sqrt(2Pi) * (x+n)^(x+n-1/2) * exp(-x-n) * E(x+n)
+
+               // E(x+n) / E(x) =
+               // Gamma(x+n)/Gamma(x) * (x^(x-1/2) * exp(-x)) / ((x+n)^(x+n-1/2) * exp(-x-n)) =
+               // (x*(x+1)*...*(x+n-1)) * exp(n) * (x/(x+n))^(x-1/2) / (x+n)^n =
+               // ((x+1)*...*(x+n-1)) * exp(n) * (x/(x+n))^(x+1/2) / (x+n)^(n-1) =
+               // ((x+1)/(x+n)*...*(x+n-1)/(x+n)) * exp(n) * (x/(x+n))^(x+1/2)
+
+               for (i = 1; i < n; i++) {
+                       // *= (x+i)/(x+n)
+                       GnmQuad xpi;
+                       gnm_quad_init (&xpi, i);
+                       gnm_quad_add (&xpi, &xpi, x);
+                       gnm_quad_div (res, res, &xpi);
+                       gnm_quad_mul (res, res, &xpn);
+               }
+
+               // /= exp(n)
+               gnm_quad_init (&en, n);
+               gnm_quad_exp (&en, NULL, &en);
+               gnm_quad_div (res, res, &en);
+
+               // /= (x/(x+n))^(x+1/2)
+               gnm_quad_init (&xph, 0.5);
+               gnm_quad_add (&xph, &xph, x);
+               gnm_quad_div (&xxn, x, &xpn);
+               gnm_quad_pow (&xxn, NULL, &xxn, &xph);
+               gnm_quad_div (res, res, &xxn);
+       }
 }
 
 /* ------------------------------------------------------------------------- */
@@ -674,7 +718,7 @@ pochhammer_small_n (gnm_float x, gnm_float n, GnmQuad *res)
        gnm_float r;
        gboolean debug = FALSE;
 
-       g_return_if_fail (x >= 20);
+       g_return_if_fail (x >= 1);
        g_return_if_fail (gnm_abs (n) <= 1);
 
        /*
@@ -1277,6 +1321,23 @@ complex_fact (gnm_complex *dst, gnm_complex const *src)
 
 /* ------------------------------------------------------------------------- */
 
+// D(a,z) := z^a * exp(-z) / Gamma (a + 1)
+static void
+complex_temme_D (gnm_complex *res, gnm_complex const *a, gnm_complex const *z)
+{
+       gnm_complex t, ez, fa;
+
+       // The idea here is to control intermediate sizes and to avoid
+       // accuracy problems caused by exp and pow.  For now, do neither.
+
+       gnm_complex_pow (&t, z, a);
+       gnm_complex_exp (&ez, z);
+       gnm_complex_div (&t, &t, &ez);
+       complex_gamma (&fa, a);
+       gnm_complex_div (res, &t, &fa);
+}
+
+
 typedef void (*GnmComplexContinuedFraction) (gnm_complex *ai, gnm_complex *bi,
                                             size_t i, const gnm_complex *args);
 
@@ -1397,6 +1458,56 @@ igamma_upper_cf (gnm_complex *dst, const gnm_complex *a, const gnm_complex *z)
        gnm_complex_mul (dst, &res, &f);
 }
 
+static gboolean
+igamma_asymp (gnm_complex *dst, const gnm_complex *a, const gnm_complex *z)
+{
+       gnm_float am = gnm_complex_mod (a);
+       gnm_float zm = gnm_complex_mod (z);
+       gnm_float n0;
+       gnm_complex s, t, am1;
+       gboolean debug = TRUE;
+       size_t i;
+
+       if (am >= zm)
+               return FALSE;
+
+       // Things start to diverge here
+       n0 = a->re + gnm_sqrt (zm * zm - a->im * a->im);
+       (void)n0;
+
+       // FIXME: Verify this condition for whether we have enough
+       // precision at term n0
+       if (2 * zm < GNM_MANT_DIG * M_LN2gnum)
+               return FALSE;
+
+       gnm_complex_real (&s, 0);
+
+       gnm_complex_init (&am1, a->re - 1, a->im);
+       complex_temme_D (&t, &am1, z);
+
+       for (i = 0; i < 100; i++) {
+               gnm_complex api;
+
+               gnm_complex_add (&s, &s, &t);
+
+               if (gnm_complex_mod (&t) <= gnm_complex_mod (&s) * (16 * GNM_EPSILON)) {
+                       if (debug)
+                               g_printerr ("igamma_asymp converged.\n");
+                       *dst = s;
+                       return TRUE;
+               }
+
+               gnm_complex_div (&t, &t, z);
+               gnm_complex_init (&api, a->re + i, a->im);
+               gnm_complex_mul (&t, &t, &api);
+       }
+
+       if (debug)
+               g_printerr ("igamma_asymp failed to converge.\n");
+
+       return FALSE;
+}
+
 
 void
 complex_igamma (gnm_complex *dst, const gnm_complex *a, const gnm_complex *z,
@@ -1421,6 +1532,12 @@ complex_igamma (gnm_complex *dst, const gnm_complex *a, const gnm_complex *z,
                goto fixup;
        }
 
+       if (igamma_asymp (&res, a, z)) {
+               have_lower = TRUE;
+               have_regularized = TRUE;
+               goto fixup;
+       }
+
        igamma_upper_cf (&res, a, z);
        have_lower = FALSE;
        have_regularized = FALSE;