[gnumeric] IMIGAMMA: improve coverage
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] IMIGAMMA: improve coverage
- Date: Wed, 3 Feb 2016 02:36:51 +0000 (UTC)
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;
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]