hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

acb_hypgeom_sum - avoid dividing s by the denominator of the next term, since this could lead to a division by zero when evaluating a hypergeometric polynomial truncated early (e.g. gamma_upper(-5,1))

Browse source
This commit is contained in:
Fredrik Johansson 2016-02-22 16:26:26 +01:00
parent 8d7e82c3f6
commit 89b67ce6db
3 changed files with 114 additions and 52 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

109
acb_hypgeom/pfq_sum_bs.c
View file

@ -44,6 +44,33 @@ A1 B2 + B1 B2 = B2 (A1 + B1) -- use to save time?
*/
static void
factor(acb_t A, acb_t tmp, acb_srcptr a, slong p, const acb_t z, slong k, slong prec)
{
slong i;
if (p == 0)
{
if (z == NULL)
acb_one(A);
else
acb_set(A, z);
}
else
{
acb_add_ui(A, a, k, prec);
for (i = 1; i < p; i++)
{
acb_add_ui(tmp, a + i, k, prec);
acb_mul(A, A, tmp, prec);
}
if (z != NULL)
acb_mul(A, A, z, prec);
}
}
static void
bsplit(acb_t A1, acb_t B1, acb_t C1,
acb_srcptr a, slong p,
@ -56,50 +83,8 @@ bsplit(acb_t A1, acb_t B1, acb_t C1,
{
if (bb - aa == 1)
{
slong i;
if (p == 0)
{
if (invz)
acb_one(A1);
else
acb_set(A1, z);
}
else
{
acb_add_ui(A1, a, aa, prec);
for (i = 1; i < p; i++)
{
acb_add_ui(B1, a + i, aa, prec);
acb_mul(A1, A1, B1, prec);
}
if (!invz)
acb_mul(A1, A1, z, prec);
}
if (q == 0)
{
if (invz)
acb_set(C1, z);
else
acb_one(C1);
}
else
{
acb_add_ui(C1, b, aa, prec);
for (i = 1; i < q; i++)
{
acb_add_ui(B1, b + i, aa, prec);
acb_mul(C1, C1, B1, prec);
}
if (invz)
acb_mul(C1, C1, z, prec);
}
factor(A1, B1, a, p, invz ? NULL : z, aa, prec);
factor(C1, B1, b, q, invz ? z : NULL, aa, prec);
/* acb_set(B1, C1); but we skip this */
}
else
@ -149,7 +134,7 @@ void
acb_hypgeom_pfq_sum_bs(acb_t s, acb_t t,
acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, slong n, slong prec)
{
acb_t u, v, w;
acb_t u, v, w, tmp;
if (n < 4)
{
@ -160,22 +145,34 @@ acb_hypgeom_pfq_sum_bs(acb_t s, acb_t t,
acb_init(u);
acb_init(v);
acb_init(w);
acb_init(tmp);
bsplit(u, v, w, a, p, b, q, z, 0, n, prec, 0);
/* we compute to n-1 instead of n to avoid dividing by 0 in the
denominator when computing a hypergeometric polynomial
that terminates right before a pole */
bsplit(u, v, w, a, p, b, q, z, 0, n - 1, prec, 0);
acb_add(s, u, v, prec); /* s = s + t */
acb_div(s, s, w, prec);
/* split off last factor */
factor(t, tmp, a, p, z, n - 1, prec);
acb_mul(u, u, t, prec);
factor(t, tmp, b, q, NULL, n - 1, prec);
acb_mul(w, w, t, prec);
acb_div(t, u, w, prec);
acb_div(s, v, w, prec);
acb_clear(u);
acb_clear(v);
acb_clear(w);
acb_clear(tmp);
}
void
acb_hypgeom_pfq_sum_bs_invz(acb_t s, acb_t t,
acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, slong n, slong prec)
{
acb_t u, v, w;
acb_t u, v, w, tmp;
if (n < 4)
{
@ -189,14 +186,26 @@ acb_hypgeom_pfq_sum_bs_invz(acb_t s, acb_t t,
acb_init(u);
acb_init(v);
acb_init(w);
acb_init(tmp);
bsplit(u, v, w, a, p, b, q, z, 0, n, prec, 1);
/* we compute to n-1 instead of n to avoid dividing by 0 in the
denominator when computing a hypergeometric polynomial
that terminates right before a pole */
bsplit(u, v, w, a, p, b, q, z, 0, n - 1, prec, 1);
acb_add(s, u, v, prec); /* s = s + t */
acb_div(s, s, w, prec);
/* split off last factor */
factor(t, tmp, a, p, NULL, n - 1, prec);
acb_mul(u, u, t, prec);
factor(t, tmp, b, q, z, n - 1, prec);
acb_mul(w, w, t, prec);
acb_div(t, u, w, prec);
acb_div(s, v, w, prec);
acb_clear(u);
acb_clear(v);
acb_clear(w);
acb_clear(tmp);
}

14
acb_hypgeom/pfq_sum_fme.c
View file

@ -179,8 +179,18 @@ acb_hypgeom_pfq_sum_fme(acb_t s, acb_t t,
acb_ptr * tree;
slong i, k, m, w;
m = n_sqrt(n) / 4; /* tuning parameter */
w = n / FLINT_MAX(m, 1);
/* we compute to n-1 instead of n to avoid dividing by 0 in the
denominator when computing a hypergeometric polynomial
that terminates right before a pole */
if (n > 4)
{
m = n_sqrt(n - 1) / 4; /* tuning parameter */
w = (n - 1) / FLINT_MAX(m, 1);
}
else
{
m = w = 0;
}
if (m < 1 || w < 1 || p > 3 || q > 3)
{

43
acb_hypgeom/test/t-gamma_upper.c
View file

@ -35,6 +35,49 @@ int main()
flint_randinit(state);
/* special accuracy test -- see nemo #38 */
for (iter = 0; iter < 1000; iter++)
{
acb_t a, z, res;
slong prec, goal;
int modified;
acb_init(a);
acb_init(z);
acb_init(res);
acb_set_si(a, n_randint(state, 100) - 50);
do {
acb_set_si(z, n_randint(state, 100) - 50);
} while (acb_is_zero(z));
modified = n_randint(state, 2);
goal = 2 + n_randint(state, 4000);
for (prec = 2 + n_randint(state, 1000); ; prec *= 2)
{
acb_hypgeom_gamma_upper(res, a, z, modified, prec);
if (acb_rel_accuracy_bits(res) > goal)
break;
if (prec > 10000)
{
printf("FAIL (convergence)\n");
flint_printf("a = "); acb_printd(a, 30); flint_printf("\n\n");
flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
flint_printf("res = "); acb_printd(res, 30); flint_printf("\n\n");
abort();
}
}
acb_clear(a);
acb_clear(z);
acb_clear(res);
}
for (iter = 0; iter < 2000; iter++)
{
acb_t a0, a1, b, z, w0, w1, t, u;