[gnumeric] GAMMA: minor improvements.

[Morten Welinder <mortenw src gnome org>]

commit ec415e26f4600cfa5d85ea364a9fb2caa861e403
Author: Morten Welinder <terra gnome org>
Date:   Tue Dec 3 22:29:12 2013 -0500

    GAMMA: minor improvements.

 ChangeLog      |    8 ++++++++
 NEWS           |    1 +
 src/complex.c  |   15 ++++++++-------
 src/sf-gamma.c |    4 ++--
 4 files changed, 19 insertions(+), 9 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 43acdc1..4027159 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,11 @@
+2013-12-03  welinder  <terra gnome org>
+
+       * src/sf-gamma.c (qgammaf): Avoid losing the least significant bit
+       of the argument for [-1.5;-0.5].
+
+       * src/complex.c (complex_fact): Avoid infinite recursion for 0 <
+       Re z < 1/2.  Avoid some overflow.
+
 2013-12-02  Morten Welinder  <terra gnome org>
 
        * src/sf-bessel.c (bessel_y): Use the J series when possible.
diff --git a/NEWS b/NEWS
index 00ff9c8..921ac60 100644
--- a/NEWS
+++ b/NEWS
@@ -8,6 +8,7 @@ Morten:
        * Improve accuracy of BETA and BETALN.
        * Improve accuracy of BESSELJ and BESSELY.
        * New functions: IMGAMMA, IMFACT.
+       * Avoid some overflows in IMGAMMA.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.9
diff --git a/src/complex.c b/src/complex.c
index ee2c103..96d8791 100644
--- a/src/complex.c
+++ b/src/complex.c
@@ -199,7 +199,7 @@ complex_fact (complex_t *dst, complex_t const *src)
 {
        if (complex_real_p (src)) {
                complex_init (dst, gnm_fact (src->re), 0);
-       } else if (src->re < 0.5) {
+       } else if (src->re < 0) {
                /* Fact(z) = -pi / (sin(pi*z) * Gamma(-z)) */
                complex_t a, b, mz;
 
@@ -209,6 +209,7 @@ complex_fact (complex_t *dst, complex_t const *src)
                complex_init (&b,
                              M_PIgnum * gnm_fmod (src->re, 2),
                              M_PIgnum * src->im);
+               /* Hmm... sin overflows when b.im is large.  */
                complex_sin (&b, &b);
 
                complex_mul (&a, &a, &b);
@@ -237,11 +238,13 @@ complex_fact (complex_t *dst, complex_t const *src)
                const gnm_float g = GNM_const(607.0) / 128;
                const gnm_float sqrt2pi =
                        GNM_const (2.506628274631000502415765284811045253006986740609938316629923);
-               complex_t zph, zpgh, s, f;
+               complex_t zph, zpghde, s, f;
                int i;
 
                complex_init (&zph, src->re + 0.5, src->im);
-               complex_init (&zpgh, src->re + g + 0.5, src->im);
+               complex_init (&zpghde,
+                             (src->re + g + 0.5) / M_Egnum,
+                             src->im / M_Egnum);
                complex_init (&s, 0, 0);
 
                for (i = G_N_ELEMENTS(c) - 1; i >= 1; i--) {
@@ -253,11 +256,9 @@ complex_fact (complex_t *dst, complex_t const *src)
                }
                s.re += c[0];
 
-               complex_init (&f, sqrt2pi, 0);
+               complex_init (&f, sqrt2pi * gnm_exp (-g), 0);
                complex_mul (&s, &s, &f);
-               complex_exp (&f, &zpgh);
-               complex_div (&s, &s, &f);
-               complex_pow (&f, &zpgh, &zph);
+               complex_pow (&f, &zpghde, &zph);
                complex_mul (&s, &s, &f);
 
                *dst = s;
diff --git a/src/sf-gamma.c b/src/sf-gamma.c
index c522247..a3068f6 100644
--- a/src/sf-gamma.c
+++ b/src/sf-gamma.c
@@ -991,7 +991,7 @@ qfactf (gnm_float x, GnmQuad *mant, int *exp2)
 static int
 qgammaf (gnm_float x, GnmQuad *mant, int *exp2)
 {
-       if (gnm_abs (x) > 0.5)
+       if (x < -1.5 || x > 0.5)
                return qfactf (x - 1, mant, exp2);
        else if (gnm_isnan (x) || x == 0)
                return 2;
@@ -1005,7 +1005,7 @@ qgammaf (gnm_float x, GnmQuad *mant, int *exp2)
                rescale_mant_exp (mant, exp2);
                gnm_quad_end (state);
                return 0;
-       }               
+       }
 }
 
 /* ------------------------------------------------------------------------- */