implement binary splitting acb_hypgeom summation

This commit is contained in:
Fredrik Johansson 2015-07-13 15:47:04 +02:00
parent d7519e76ff
commit 57b1dd98d1
7 changed files with 458 additions and 3 deletions

View file

@ -45,12 +45,21 @@ void acb_hypgeom_pfq_sum_forward(acb_t s, acb_t t, acb_srcptr a, long p, acb_src
void acb_hypgeom_pfq_sum_rs(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q,
const acb_t z, long n, long prec);
void acb_hypgeom_pfq_sum_bs(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q,
const acb_t z, long n, long prec);
void acb_hypgeom_pfq_sum_fme(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q,
const acb_t z, long n, long prec);
void acb_hypgeom_pfq_sum(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q,
const acb_t z, long n, long prec);
void acb_hypgeom_pfq_sum_bs_invz(acb_t s, acb_t t,
acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec);
void acb_hypgeom_pfq_sum_invz(acb_t s, acb_t t, acb_srcptr a, long p,
acb_srcptr b, long q, const acb_t z, const acb_t zinv, long n, long prec);
void acb_hypgeom_pfq_direct(acb_t res, acb_srcptr a, long p, acb_srcptr b, long q,
const acb_t z, long n, long prec);

View file

@ -32,7 +32,32 @@ acb_hypgeom_pfq_sum(acb_t s, acb_t t, acb_srcptr a, long p,
if (n > 4 && prec >= 128
&& _acb_vec_bits(a, p) * p + _acb_vec_bits(b, q) * q + 10 < prec / 2)
{
acb_hypgeom_pfq_sum_rs(s, t, a, p, b, q, z, n, prec);
if (prec >= 256 && acb_bits(z) < prec * 0.01)
acb_hypgeom_pfq_sum_bs(s, t, a, p, b, q, z, n, prec);
else
acb_hypgeom_pfq_sum_rs(s, t, a, p, b, q, z, n, prec);
}
else if (prec >= 1500 && n >= 30 + 100000 / (prec - 1000))
{
acb_hypgeom_pfq_sum_fme(s, t, a, p, b, q, z, n, prec);
}
else
{
acb_hypgeom_pfq_sum_forward(s, t, a, p, b, q, z, n, prec);
}
}
void
acb_hypgeom_pfq_sum_invz(acb_t s, acb_t t, acb_srcptr a, long p,
acb_srcptr b, long q, const acb_t z, const acb_t zinv, long n, long prec)
{
if (n > 4 && prec >= 128
&& _acb_vec_bits(a, p) * p + _acb_vec_bits(b, q) * q + 10 < prec / 2)
{
if (prec >= 256 && acb_bits(zinv) < prec * 0.01)
acb_hypgeom_pfq_sum_bs_invz(s, t, a, p, b, q, zinv, n, prec);
else
acb_hypgeom_pfq_sum_rs(s, t, a, p, b, q, z, n, prec);
}
else if (prec >= 1500 && n >= 30 + 100000 / (prec - 1000))
{

202
acb_hypgeom/pfq_sum_bs.c Normal file
View file

@ -0,0 +1,202 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2015 Fredrik Johansson
******************************************************************************/
#include "acb_hypgeom.h"
/*
[S(k+1)] = [ R(k) 0 ] [S(k)]
[T(k+1)] [ 1 1 ] [T(k)]
[S(k+1)] = [ P(k) / Q(k) 0 ] [S(k)]
[T(k+1)] [ 1 1 ] [T(k)]
1 [ P(k) ]
---- [ ]
Q(k) [ Q(k) Q(k) ]
[[A2 0] [B2 C2]] . [[A1 0] [B1 C1]] = [[A1 A2 0] [A1 B2 + B1 C2 C1 C2]
A1 B2 + B1 B2 = B2 (A1 + B1) -- use to save time?
*/
static void
bsplit(acb_t A1, acb_t B1, acb_t C1,
acb_srcptr a, long p,
acb_srcptr b, long q,
const acb_t z,
long aa,
long bb,
long prec,
int invz)
{
if (bb - aa == 1)
{
long 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);
}
/* acb_set(B1, C1); but we skip this */
}
else
{
long m;
acb_t A2, B2, C2;
acb_init(A2);
acb_init(B2);
acb_init(C2);
m = aa + (bb - aa) / 2;
bsplit(A1, B1, C1, a, p, b, q, z, aa, m, prec, invz);
bsplit(A2, B2, C2, a, p, b, q, z, m, bb, prec, invz);
if (bb - m == 1) /* B2 = C2 */
{
if (m - aa == 1)
acb_add(B2, A1, C1, prec);
else
acb_add(B2, A1, B1, prec);
acb_mul(B1, B2, C2, prec);
}
else
{
if (m - aa == 1)
acb_mul(B1, C1, C2, prec);
else
acb_mul(B1, B1, C2, prec);
acb_addmul(B1, A1, B2, prec);
}
acb_mul(A1, A1, A2, prec);
acb_mul(C1, C1, C2, prec);
acb_clear(A2);
acb_clear(B2);
acb_clear(C2);
}
}
void
acb_hypgeom_pfq_sum_bs(acb_t s, acb_t t,
acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
{
acb_t u, v, w;
if (n < 4)
{
acb_hypgeom_pfq_sum_forward(s, t, a, p, b, q, z, n, prec);
return;
}
acb_init(u);
acb_init(v);
acb_init(w);
bsplit(u, v, w, a, p, b, q, z, 0, n, prec, 0);
acb_div(t, u, w, prec);
acb_div(s, v, w, prec);
acb_clear(u);
acb_clear(v);
acb_clear(w);
}
void
acb_hypgeom_pfq_sum_bs_invz(acb_t s, acb_t t,
acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
{
acb_t u, v, w;
if (n < 4)
{
acb_init(u);
acb_inv(u, z, prec);
acb_hypgeom_pfq_sum_forward(s, t, a, p, b, q, u, n, prec);
acb_clear(u);
return;
}
acb_init(u);
acb_init(v);
acb_init(w);
bsplit(u, v, w, a, p, b, q, z, 0, n, prec, 1);
acb_div(t, u, w, prec);
acb_div(s, v, w, prec);
acb_clear(u);
acb_clear(v);
acb_clear(w);
}

View file

@ -0,0 +1,98 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "acb_hypgeom.h"
int main()
{
long iter;
flint_rand_t state;
printf("pfq_sum_bs....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000; iter++)
{
acb_ptr a, b;
acb_t z, s1, s2, t1, t2;
long i, p, q, n, prec1, prec2;
p = n_randint(state, 5);
q = n_randint(state, 5);
n = n_randint(state, 300);
prec1 = 2 + n_randint(state, 500);
prec2 = 2 + n_randint(state, 500);
acb_init(z);
acb_init(s1);
acb_init(s2);
acb_init(t1);
acb_init(t2);
acb_randtest_special(z, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
acb_randtest_special(s1, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
acb_randtest_special(t1, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
acb_randtest_special(s2, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
acb_randtest_special(t2, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
a = _acb_vec_init(p);
b = _acb_vec_init(q);
for (i = 0; i < p; i++)
acb_randtest(a + i, state, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
for (i = 0; i < q; i++)
acb_randtest(b + i, state, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
acb_hypgeom_pfq_sum_forward(s1, t1, a, p, b, q, z, n, prec1);
acb_hypgeom_pfq_sum_bs(s2, t2, a, p, b, q, z, n, prec2);
if (!acb_overlaps(s1, s2) || !acb_overlaps(t1, t2))
{
printf("FAIL: overlap\n\n");
printf("z = "); acb_print(a); printf("\n\n");
printf("s1 = "); acb_print(s1); printf("\n\n");
printf("s2 = "); acb_print(s2); printf("\n\n");
printf("t1 = "); acb_print(t1); printf("\n\n");
printf("t2 = "); acb_print(t2); printf("\n\n");
abort();
}
_acb_vec_clear(a, p);
_acb_vec_clear(b, q);
acb_clear(z);
acb_clear(s1);
acb_clear(s2);
acb_clear(t1);
acb_clear(t2);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

View file

@ -0,0 +1,107 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2014 Fredrik Johansson
******************************************************************************/
#include "acb_hypgeom.h"
int main()
{
long iter;
flint_rand_t state;
printf("pfq_sum_invz....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 2000; iter++)
{
acb_ptr a, b;
acb_t z, zinv, s1, s2, t1, t2;
long i, p, q, n, prec1, prec2;
p = n_randint(state, 5);
q = n_randint(state, 5);
n = n_randint(state, 300);
prec1 = 2 + n_randint(state, 500);
prec2 = 2 + n_randint(state, 500);
acb_init(z);
acb_init(zinv);
acb_init(s1);
acb_init(s2);
acb_init(t1);
acb_init(t2);
acb_randtest_special(z, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
acb_randtest_special(s1, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
acb_randtest_special(t1, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
acb_randtest_special(s2, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
acb_randtest_special(t2, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
acb_inv(zinv, z, 1 + n_randint(state, 500));
a = _acb_vec_init(p);
b = _acb_vec_init(q);
for (i = 0; i < p; i++)
acb_randtest(a + i, state, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
for (i = 0; i < q; i++)
acb_randtest(b + i, state, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
if (n_randint(state, 2))
acb_hypgeom_pfq_sum_invz(s1, t1, a, p, b, q, zinv, z, n, prec1);
else
acb_hypgeom_pfq_sum_bs_invz(s1, t1, a, p, b, q, z, n, prec1);
acb_hypgeom_pfq_sum(s2, t2, a, p, b, q, zinv, n, prec2);
if (!acb_overlaps(s1, s2) || !acb_overlaps(t1, t2))
{
printf("FAIL: overlap\n\n");
printf("n = %ld\n\n", n);
printf("z = "); acb_printd(z, 15); printf("\n\n");
printf("s1 = "); acb_printd(s1, 15); printf("\n\n");
printf("s2 = "); acb_printd(s2, 15); printf("\n\n");
printf("t1 = "); acb_printd(t1, 15); printf("\n\n");
printf("t2 = "); acb_printd(t2, 15); printf("\n\n");
abort();
}
_acb_vec_clear(a, p);
_acb_vec_clear(b, q);
acb_clear(z);
acb_clear(zinv);
acb_clear(s1);
acb_clear(s2);
acb_clear(t1);
acb_clear(t2);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

View file

@ -179,7 +179,7 @@ void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b,
{
mag_t C1, Cn, alpha, nu, sigma, rho, zinv, tmp, err;
acb_struct aa[3];
acb_t s, t, w;
acb_t s, t, w, winv;
int R, p, q, is_real;
if (!acb_is_finite(a) || !acb_is_finite(b) || !acb_is_finite(z))
@ -204,6 +204,7 @@ void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b,
acb_init(s);
acb_init(t);
acb_init(w);
acb_init(winv);
/* special case, for incomplete gamma
[todo: also when they happen to be exact and with difference 1...] */
@ -228,6 +229,7 @@ void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b,
acb_neg(w, z);
acb_inv(w, w, prec);
acb_neg(winv, z);
if (n < 0)
n = acb_hypgeom_pfq_choose_n(aa, p, aa + p, q, w, prec);
@ -243,7 +245,7 @@ void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b,
}
else
{
acb_hypgeom_pfq_sum(s, t, aa, p, aa + p, q, w, n, prec);
acb_hypgeom_pfq_sum_invz(s, t, aa, p, aa + p, q, w, winv, n, prec);
if (R == 1)
{
@ -307,5 +309,6 @@ void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b,
acb_clear(s);
acb_clear(t);
acb_clear(w);
acb_clear(winv);
}

View file

@ -86,6 +86,8 @@ or remove a 1 from the `a_i` parameters if there is one.
.. function:: void acb_hypgeom_pfq_sum_rs(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
.. function:: void acb_hypgeom_pfq_sum_bs(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
.. function:: void acb_hypgeom_pfq_sum_fme(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
.. function:: void acb_hypgeom_pfq_sum(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
@ -97,6 +99,8 @@ or remove a 1 from the `a_i` parameters if there is one.
The *forward* version computes the sum using forward
recurrence.
The *bs* version computes the sum using binary splitting.
The *rs* version computes the sum in reverse order
using rectangular splitting. It only computes a
magnitude bound for the value of *t*.
@ -106,6 +110,13 @@ or remove a 1 from the `a_i` parameters if there is one.
The default version automatically chooses an algorithm
depending on the inputs.
.. function:: void acb_hypgeom_pfq_sum_bs_invz(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t w, long n, long prec)
.. function:: void acb_hypgeom_pfq_sum_invz(acb_t s, acb_t t, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, const acb_t w, long n, long prec)
Like :func:`acb_hypgeom_pfq_sum`, but taking advantage of
`w = 1/z` possibly having few bits.
.. function:: void acb_hypgeom_pfq_direct(acb_t res, acb_srcptr a, long p, acb_srcptr b, long q, const acb_t z, long n, long prec)
Computes