[gnumeric] IMIGAMMA: More fixes

[Morten Welinder <mortenw src gnome org>]

commit 7c1a902981f403c11555661ef352569bb1070fa8
Author: Morten Welinder <terra gnome org>
Date:   Sat Feb 6 11:41:56 2016 -0500

    IMIGAMMA: More fixes

 ChangeLog      |    8 ++++++++
 configure.ac   |    2 +-
 src/complex.h  |    4 ++++
 src/sf-gamma.c |   41 +++++++++++++++++++++++++++++++++++------
 4 files changed, 48 insertions(+), 7 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 95702e9..26898e9 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,11 @@
+2016-02-06  Morten Welinder  <terra gnome org>
+
+       * configure.ac (goffice): Require latest for
+       go_complex_from_polar_pi.
+
+       * src/sf-gamma.c (complex_igamma): Apply fixup for upper gamma
+       when x<0 and a is real.
+
 2016-02-04  Morten Welinder  <terra gnome org>
 
        * src/sf-gamma.c (complex_temme_D): Fix factorial computation.
diff --git a/configure.ac b/configure.ac
index e18211e..dd29d63 100644
--- a/configure.ac
+++ b/configure.ac
@@ -161,7 +161,7 @@ PKG_PROG_PKG_CONFIG(0.18)
 
 dnl *****************************
 libspreadsheet_reqs="
-       libgoffice-${GOFFICE_API_VER}   >= 0.10.22
+       libgoffice-${GOFFICE_API_VER}   >= 0.10.27
        libgsf-1                >= 1.14.33
        libxml-2.0              >= 2.4.12
 "
diff --git a/src/complex.h b/src/complex.h
index 55908f5..08a49a9 100644
--- a/src/complex.h
+++ b/src/complex.h
@@ -17,6 +17,7 @@ G_BEGIN_DECLS
 #define gnm_complex_div go_complex_divl
 #define gnm_complex_mod go_complex_modl
 #define gnm_complex_angle go_complex_anglel
+#define gnm_complex_angle_pi go_complex_angle_pil
 #define gnm_complex_real go_complex_reall
 #define gnm_complex_real_p go_complex_real_pl
 #define gnm_complex_zero_p go_complex_zero_pl
@@ -31,6 +32,7 @@ G_BEGIN_DECLS
 #define gnm_complex_scale_real go_complex_scale_reall
 #define gnm_complex_to_polar go_complex_to_polarl
 #define gnm_complex_from_polar go_complex_from_polarl
+#define gnm_complex_from_polar_pi go_complex_from_polar_pil
 #else
 #define gnm_complex GOComplex
 #define gnm_complex_init go_complex_init
@@ -40,6 +42,7 @@ G_BEGIN_DECLS
 #define gnm_complex_div go_complex_div
 #define gnm_complex_mod go_complex_mod
 #define gnm_complex_angle go_complex_angle
+#define gnm_complex_angle_pi go_complex_angle_pi
 #define gnm_complex_real go_complex_real
 #define gnm_complex_real_p go_complex_real_p
 #define gnm_complex_zero_p go_complex_zero_p
@@ -54,6 +57,7 @@ G_BEGIN_DECLS
 #define gnm_complex_scale_real go_complex_scale_real
 #define gnm_complex_to_polar go_complex_to_polar
 #define gnm_complex_from_polar go_complex_from_polar
+#define gnm_complex_from_polar_pi go_complex_from_polar_pi
 #endif
 
 /* ------------------------------------------------------------------------- */
diff --git a/src/sf-gamma.c b/src/sf-gamma.c
index 31df6ff..2081962 100644
--- a/src/sf-gamma.c
+++ b/src/sf-gamma.c
@@ -1468,7 +1468,7 @@ igamma_lower_cf (gnm_complex *dst, const gnm_complex *a, const gnm_complex *z)
 }
 
 static gboolean
-igamma_asymp (gnm_complex *dst, const gnm_complex *a, const gnm_complex *z)
+igamma_upper_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);
@@ -1508,7 +1508,7 @@ igamma_asymp (gnm_complex *dst, const gnm_complex *a, const gnm_complex *z)
 
                if (gnm_complex_mod (&t) <= gnm_complex_mod (&s) * GNM_EPSILON) {
                        if (debug)
-                               g_printerr ("igamma_asymp converged.\n");
+                               g_printerr ("igamma_upper_asymp converged.\n");
                        *dst = s;
                        return TRUE;
                }
@@ -1519,18 +1519,46 @@ igamma_asymp (gnm_complex *dst, const gnm_complex *a, const gnm_complex *z)
        }
 
        if (debug)
-               g_printerr ("igamma_asymp failed to converge.\n");
+               g_printerr ("igamma_upper_asymp failed to converge.\n");
 
        return FALSE;
 }
 
+static void
+fixup_upper_real (gnm_complex *res, gnm_complex const *a, gnm_complex const *z)
+{
+       // This assumes we have an upper gamma regularized result.
+       //
+       // It appears that upper algorithms have trouble with negative real z
+       // (for example, such z being outside the allowed domain) in some cases.
+       
+
+       if (gnm_complex_real_p (z) && z->re < 0 &&
+           gnm_complex_real_p (a) && a->re != gnm_floor (a->re)) {
+               // Everything in the lower power series is real expect z^a
+               // which is not.  So...
+               // 1. Flip to lower gamma
+               // 2. Assume modulus is correct
+               // 3. Use exact angle for lower gamma
+               // 4. Flip back to upper gamma
+               gnm_complex lres = *res;
+               lres.re = 1 - lres.re;
+
+               gnm_complex_from_polar_pi (res,
+                                          gnm_complex_mod (&lres),
+                                          a->re);
+               res->re = 1 - res->re;
+               res->im = 0 - res->im;
+       }
+}
 
 void
 complex_igamma (gnm_complex *dst, const gnm_complex *a, const gnm_complex *z,
                gboolean lower, gboolean regularized)
 {
        gnm_complex res, ga;
-       gboolean have_lower, have_regularized, have_ga = FALSE;
+       gboolean have_lower, have_regularized;
+       gboolean have_ga = FALSE;
 
        if (regularized && gnm_complex_real_p (a) &&
            a->re <= 0 && a->re == gnm_floor (a->re)) {
@@ -1548,13 +1576,14 @@ complex_igamma (gnm_complex *dst, const gnm_complex *a, const gnm_complex *z,
                goto fixup;
        }
 
-       if (igamma_asymp (&res, a, z)) {
+       if (igamma_upper_asymp (&res, a, z)) {
                have_lower = FALSE;
                have_regularized = TRUE;
+               fixup_upper_real (&res, a, z);
                goto fixup;
        }
 
-       if (0 && igamma_lower_cf (&res, a, z)) {
+       if (igamma_lower_cf (&res, a, z)) {
                have_lower = TRUE;
                have_regularized = FALSE;
                goto fixup;