Commit 7c1a9029 authored by Morten Welinder's avatar Morten Welinder

IMIGAMMA: More fixes

parent 4d613e62
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.
......
......@@ -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
"
......
......@@ -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
/* ------------------------------------------------------------------------- */
......
......@@ -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;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment