[gnumeric] IMGAMMA: Greatly improve accuracy.

[Morten Welinder <mortenw src gnome org>]

commit 3ea27b2cdaf630ce2b4f634c7a3f99345bc964a2
Author: Morten Welinder <terra gnome org>
Date:   Fri Dec 6 15:53:19 2013 -0500

    IMGAMMA: Greatly improve accuracy.

 plugins/fn-complex/functions.c |    1 +
 src/complex.c                  |   88 ----------------------------
 src/complex.h                  |    3 -
 src/sf-gamma.c                 |  123 ++++++++++++++++++++++++++++++++++++++++
 src/sf-gamma.h                 |    3 +
 5 files changed, 127 insertions(+), 91 deletions(-)
---
diff --git a/plugins/fn-complex/functions.c b/plugins/fn-complex/functions.c
index 93373a1..66c3bac 100644
--- a/plugins/fn-complex/functions.c
+++ b/plugins/fn-complex/functions.c
@@ -29,6 +29,7 @@
 #include <func.h>
 
 #include <complex.h>
+#include <sf-gamma.h>
 #include <parse-util.h>
 #include <cell.h>
 #include <expr.h>
diff --git a/src/complex.c b/src/complex.c
index 96d8791..799c324 100644
--- a/src/complex.c
+++ b/src/complex.c
@@ -178,91 +178,3 @@ complex_invalid_p (complex_t const *src)
 }
 
 /* ------------------------------------------------------------------------- */
-
-void
-complex_gamma (complex_t *dst, complex_t const *src)
-{
-       if (complex_real_p (src)) {
-               complex_init (dst, gnm_gamma (src->re), 0);
-       } else {
-               complex_t z = *src, f;
-
-               complex_fact (&f, src);
-               complex_div (dst, &f, &z);
-       }
-}
-
-/* ------------------------------------------------------------------------- */
-
-void
-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) {
-               /* Fact(z) = -pi / (sin(pi*z) * Gamma(-z)) */
-               complex_t a, b, mz;
-
-               complex_init (&mz, -src->re, -src->im);
-               complex_gamma (&a, &mz);
-
-               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);
-
-               complex_init (&b, -M_PIgnum, 0);
-
-               complex_div (dst, &b, &a);
-       } else {
-               static const gnm_float c[] = {
-                       GNM_const (0.99999999999999709182),
-                       GNM_const (57.156235665862923517),
-                       GNM_const (-59.597960355475491248),
-                       GNM_const (14.136097974741747174),
-                       GNM_const (-0.49191381609762019978),
-                       GNM_const (0.33994649984811888699e-4),
-                       GNM_const (0.46523628927048575665e-4),
-                       GNM_const (-0.98374475304879564677e-4),
-                       GNM_const (0.15808870322491248884e-3),
-                       GNM_const (-0.21026444172410488319e-3),
-                       GNM_const (0.21743961811521264320e-3),
-                       GNM_const (-0.16431810653676389022e-3),
-                       GNM_const (0.84418223983852743293e-4),
-                       GNM_const (-0.26190838401581408670e-4),
-                       GNM_const (0.36899182659531622704e-5)
-               };
-               const gnm_float g = GNM_const(607.0) / 128;
-               const gnm_float sqrt2pi =
-                       GNM_const (2.506628274631000502415765284811045253006986740609938316629923);
-               complex_t zph, zpghde, s, f;
-               int i;
-
-               complex_init (&zph, src->re + 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--) {
-                       complex_t d, q;
-                       complex_init (&d, src->re + i, src->im);
-                       complex_init (&q, c[i], 0);
-                       complex_div (&q, &q, &d);
-                       complex_add (&s, &s, &q);
-               }
-               s.re += c[0];
-
-               complex_init (&f, sqrt2pi * gnm_exp (-g), 0);
-               complex_mul (&s, &s, &f);
-               complex_pow (&f, &zpghde, &zph);
-               complex_mul (&s, &s, &f);
-
-               *dst = s;
-       }
-}
-
-/* ------------------------------------------------------------------------- */
diff --git a/src/complex.h b/src/complex.h
index f06dea4..7efafc4 100644
--- a/src/complex.h
+++ b/src/complex.h
@@ -66,9 +66,6 @@ int complex_from_string (complex_t *dst, char const *src, char *imunit);
 
 int complex_invalid_p (complex_t const *src);
 
-void complex_gamma (complex_t *dst, complex_t const *src);
-void complex_fact (complex_t *dst, complex_t const *src);
-
 /* ------------------------------------------------------------------------- */
 
 G_END_DECLS
diff --git a/src/sf-gamma.c b/src/sf-gamma.c
index 661815d..2fac6fa 100644
--- a/src/sf-gamma.c
+++ b/src/sf-gamma.c
@@ -1119,3 +1119,126 @@ lgamma_r (double x, int *signp)
 #endif
 
 /* ------------------------------------------------------------------------- */
+
+
+static const gnm_float lanczos_g = GNM_const (808618867.0) / 134217728;
+
+/*
+ * This Mathematica sniplet computes the Lanczos gamma coefficients:
+ *
+ * Dr[k_]:=DiagonalMatrix[Join[{1},Array[-Binomial[2*#-1,#]*#&,k]]]
+ * c[k_]:= 
Array[If[#1+#2==2,1/2,If[#1>=#2,(-1)^(#1+#2)*4^(#2-1)*(#1-1)*(#1+#2-3)!/(#1-#2)!/(2*#2-2)!,0]]&,{k+1,k+1}]
+ * Dc[k_]:=DiagonalMatrix[Array[2*(2*#-3)!!&,k+1]]
+ * B[k_]:=Array[If[#1==1,1,If[#1<=#2,(-1)^(#2-#1)*Binomial[#1+#2-3,#1+#1-3],0]]&,{k+1,k+1}]
+ * M[k_]:=(Dr[k].B[k]).(c[k].Dc[k])
+ * f[g_,k_]:=Array[Sqrt[2]*(E/(2*(#-1+g)+1))^(#-1/2)&,k+1]
+ * a[g_,k_]:=M[k].f[g,k]*Exp[g]/Sqrt[2*Pi]
+ *
+ * The result of a[g,k] will contain both positive and negative constants.
+ * Most people using the Lanczos series do not understand that a naive
+ * implemetation will suffer significant cancellation errors.  The error
+ * estimates assume the series is computes without loss!
+ *
+ * Following Boost we multiply the entire partial fraction by Pochhammer[z+1,k]
+ * That gives us a polynomium with positive coefficient.  For kicks we toss
+ * the constant factor back in.
+ *
+ * b[g_,k_]:=Sum[a[g,k][[i+1]]/If[i==0,1,(z+i)],{i,0,k}]*Pochhammer[z+1,k]
+ * 
c13b:=Block[{$MaxExtraPrecision=500},FullSimplify[N[b[808618867/134217728,12]*Sqrt[2*Pi]/Exp[808618867/134217728],300]]]
+ *
+ * Finally we recast that in terms of gamma's argument:
+ *
+ * N[CoefficientList[c13b /. z->(zp-1),{zp}],50]
+ */
+static const gnm_float lanczos_num[] = {
+       GNM_const(56906521.913471563880907910335591226868592353221448),
+       GNM_const(103794043.11634454519062710536160702385539539810110),
+       GNM_const(86363131.288138591455469272889778684223419113014358),
+       GNM_const(43338889.324676138347737237405905333160850993321475),
+       GNM_const(14605578.087685068084141699827913592185707234229516),
+       GNM_const(3481712.1549806459088207101896477455646802362321652),
+       GNM_const(601859.61716810987866702265336993523025071425740828),
+       GNM_const(75999.293040145426498753034435989091370919973262979),
+       GNM_const(6955.9996025153761403563101155151989875259157712039),
+       GNM_const(449.94455690631681194468586076509884096232715968614),
+       GNM_const(19.519927882476174828478609662356521362076846583112),
+       GNM_const(0.50984166556566761881251786448046945099926051133936),
+       GNM_const(0.0060618423462489065257837539645559368832224636654970)
+};
+
+/* CoefficientList[Pochhammer[z,12],z] */
+static const guint32 lanczos_denom[G_N_ELEMENTS(lanczos_num)] = {
+       0, 39916800, 120543840, 150917976, 105258076, 45995730,
+       13339535, 2637558, 357423, 32670, 1925, 66, 1
+};
+
+void
+complex_gamma (complex_t *dst, complex_t const *src)
+{
+       if (complex_real_p (src)) {
+               complex_init (dst, gnm_gamma (src->re), 0);
+       } else if (src->re < 0) {
+               /* Gamma(z) = pi / (sin(pi*z) * Gamma(-z+1)) */
+               complex_t a, b, mz;
+
+               complex_init (&mz, -src->re, -src->im);
+               complex_fact (&a, &mz);
+
+               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);
+
+               complex_init (&b, M_PIgnum, 0);
+
+               complex_div (dst, &b, &a);
+       } else {
+               complex_t zmh, zmhd2, zmhpg, f, g, p, q, pq;
+               int i;
+
+               i = G_N_ELEMENTS(lanczos_num) - 1;
+               complex_init (&p, lanczos_num[i], 0);
+               complex_init (&q, lanczos_denom[i], 0);
+               while (--i >= 0) {
+                       complex_mul (&p, &p, src);
+                       p.re += lanczos_num[i];
+                       complex_mul (&q, &q, src);
+                       q.re += lanczos_denom[i];
+               }
+               complex_div (&pq, &p, &q);
+
+               complex_init (&zmh, src->re - 0.5, src->im);
+               complex_init (&zmhpg, zmh.re + lanczos_g, zmh.im);
+               complex_init (&zmhd2, zmh.re * 0.5, zmh.im * 0.5);
+               complex_pow (&f, &zmhpg, &zmhd2);
+
+               complex_exp (&g, &zmh);
+               complex_div (&g, &f, &g);
+               complex_mul (&g, &g, &f);
+
+               complex_mul (dst, &g, &pq);
+       }
+}
+
+/* ------------------------------------------------------------------------- */
+
+void
+complex_fact (complex_t *dst, complex_t const *src)
+{
+       if (complex_real_p (src)) {
+               complex_init (dst, gnm_fact (src->re), 0);
+       } else {
+               /*
+                * This formula is valid for all arguments except zero
+                * which we conveniently handled above.
+                */
+               complex_t gz;
+               complex_gamma (&gz, src);
+               complex_mul (dst, &gz, src);
+       }
+}
+
+/* ------------------------------------------------------------------------- */
diff --git a/src/sf-gamma.h b/src/sf-gamma.h
index b3e2023..a6fe89c 100644
--- a/src/sf-gamma.h
+++ b/src/sf-gamma.h
@@ -2,6 +2,7 @@
 #define GNM_SF_GAMMA_H_
 
 #include <numbers.h>
+#include <complex.h>
 
 gnm_float lgamma1p (gnm_float a);
 gnm_float stirlerr(gnm_float n);
@@ -9,6 +10,8 @@ gnm_float stirlerr(gnm_float n);
 gnm_float gnm_gamma (gnm_float x);
 gnm_float gnm_fact (gnm_float x);
 int       qfactf (gnm_float x, GnmQuad *mant, int *exp2);
+void complex_gamma (complex_t *dst, complex_t const *src);
+void complex_fact (complex_t *dst, complex_t const *src);
 
 gnm_float gnm_lbeta (gnm_float a, gnm_float b);
 gnm_float gnm_beta (gnm_float a, gnm_float b);