diff --git a/acb_hypgeom.h b/acb_hypgeom.h index a87a1985..9f784ed8 100644 --- a/acb_hypgeom.h +++ b/acb_hypgeom.h @@ -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); diff --git a/acb_hypgeom/pfq_sum.c b/acb_hypgeom/pfq_sum.c index 81a3ac50..b8169e1a 100644 --- a/acb_hypgeom/pfq_sum.c +++ b/acb_hypgeom/pfq_sum.c @@ -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)) { diff --git a/acb_hypgeom/pfq_sum_bs.c b/acb_hypgeom/pfq_sum_bs.c new file mode 100644 index 00000000..c9ae0443 --- /dev/null +++ b/acb_hypgeom/pfq_sum_bs.c @@ -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); +} + diff --git a/acb_hypgeom/test/t-pfq_sum_bs.c b/acb_hypgeom/test/t-pfq_sum_bs.c new file mode 100644 index 00000000..b2f48651 --- /dev/null +++ b/acb_hypgeom/test/t-pfq_sum_bs.c @@ -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; +} diff --git a/acb_hypgeom/test/t-pfq_sum_invz.c b/acb_hypgeom/test/t-pfq_sum_invz.c new file mode 100644 index 00000000..eb84a3c4 --- /dev/null +++ b/acb_hypgeom/test/t-pfq_sum_invz.c @@ -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; +} diff --git a/acb_hypgeom/u_asymp.c b/acb_hypgeom/u_asymp.c index 9de34094..6a06b8b2 100644 --- a/acb_hypgeom/u_asymp.c +++ b/acb_hypgeom/u_asymp.c @@ -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); } diff --git a/doc/source/acb_hypgeom.rst b/doc/source/acb_hypgeom.rst index 10e68bbe..5519f4cf 100644 --- a/doc/source/acb_hypgeom.rst +++ b/doc/source/acb_hypgeom.rst @@ -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