[gnumeric] gamma: improve precision for large arguments.

[Morten Welinder <mortenw src gnome org>]

commit 6c4e82ceb7ce8542bb2ea88c5c5e5986a9f04ec4
Author: Morten Welinder <terra gnome org>
Date:   Wed Aug 28 12:53:09 2013 -0400

    gamma: improve precision for large arguments.

 ChangeLog      |    5 +++++
 src/mathfunc.c |   52 ++++++++++++++++++++++++++++------------------------
 2 files changed, 33 insertions(+), 24 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 3fdd567..1d6c5ce 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2013-08-28  Morten Welinder  <terra gnome org>
+
+       * src/mathfunc.c (gnm_gamma): Improve precision for large
+       arguments.
+
 2013-08-27  Morten Welinder <terra gnome org>
 
        * configure.ac: Post-release bump.
diff --git a/src/mathfunc.c b/src/mathfunc.c
index d9db981..001fbc8 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -126,30 +126,6 @@ gnm_acoth (gnm_float x)
        return gnm_atanh (1 / x);
 }
 
-gnm_float
-gnm_gamma (gnm_float x)
-{
-       if (gnm_isnan (x))
-               return x;
-
-       if (x > 0) {
-               if (x >= G_MAXINT)
-                       return gnm_pinf;
-
-               if (x == gnm_floor (x))
-                       return fact ((int)x - 1);
-
-               /* We can probably do better than this. */
-               return gnm_exp (gnm_lgamma (x));
-       } else if (x == gnm_floor (x))
-               return gnm_nan;
-       else {
-               return M_PIgnum /
-                       (gnm_sin (x * M_PIgnum) * gnm_gamma (1 - x));
-       }
-}
-
-
 /* ------------------------------------------------------------------------- */
 /* --- BEGIN MAGIC R SOURCE MARKER --- */
 
@@ -8866,3 +8842,31 @@ pochhammer (gnm_float x, gnm_float n, gboolean give_log)
 }
 
 /* ------------------------------------------------------------------------- */
+
+gnm_float
+gnm_gamma (gnm_float x)
+{
+       if (gnm_isnan (x))
+               return x;
+
+       if (x > 0) {
+               if (x == gnm_floor (x) && x <= 100)
+                       return fact ((int)x - 1);
+
+               if (x > 10) {
+                       return gnm_pow (x / M_Egnum, x) /
+                               gnm_sqrt (x / M_2PIgnum)
+                               * gnm_exp (lgammacor (x));
+               }
+
+               /* We can probably do better than this. */
+               return gnm_exp (gnm_lgamma (x));
+       } else if (x == gnm_floor (x))
+               return gnm_nan;
+       else {
+               return M_PIgnum /
+                       (gnm_sin (x * M_PIgnum) * gnm_gamma (1 - x));
+       }
+}
+
+/* ------------------------------------------------------------------------- */