mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
pfq_series_sum: implement binary splitting and rectangular splitting -> greatly speeds up Y_n(z), K_n(z) and other functions at high precision, and speeds up high-order parameter derivatives
This commit is contained in:
Fredrik Johansson
parent
0153a518f4
commit
5d77920f43
9 changed files with 1017 additions and 145 deletions
|
|
@ -53,6 +53,34 @@ slong acb_hypgeom_pfq_series_choose_n(const acb_poly_struct * a, slong p,
|
|||
const acb_poly_struct * b, slong q,
|
||||
const acb_poly_t z, slong len, slong prec);
|
||||
|
||||
void
|
||||
acb_hypgeom_pfq_series_sum_forward(acb_poly_t s, acb_poly_t t,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_struct * b, slong q,
|
||||
const acb_poly_t z, int regularized,
|
||||
slong n, slong len, slong prec);
|
||||
|
||||
void
|
||||
acb_hypgeom_pfq_series_sum_bs(acb_poly_t s, acb_poly_t t,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_struct * b, slong q,
|
||||
const acb_poly_t z, int regularized,
|
||||
slong n, slong len, slong prec);
|
||||
|
||||
void
|
||||
acb_hypgeom_pfq_series_sum_rs(acb_poly_t s, acb_poly_t t,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_struct * b, slong q,
|
||||
const acb_poly_t z, int regularized,
|
||||
slong n, slong len, slong prec);
|
||||
|
||||
void
|
||||
acb_hypgeom_pfq_series_sum(acb_poly_t s, acb_poly_t t,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_struct * b, slong q,
|
||||
const acb_poly_t z, int regularized,
|
||||
slong n, slong len, slong prec);
|
||||
|
||||
void acb_hypgeom_pfq_series_direct(acb_poly_t res,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_struct * b, slong q,
|
||||
|
|
|
|||
|
|
@ -146,150 +146,6 @@ acb_hypgeom_pfq_series_bound_factor(arb_poly_t F,
|
|||
acb_poly_clear(AB);
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_pfq_series_sum_forward(acb_poly_t s, acb_poly_t t,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_struct * b, slong q,
|
||||
const acb_poly_t z, int regularized,
|
||||
slong n, slong len, slong prec)
|
||||
{
|
||||
acb_poly_t u, v;
|
||||
acb_poly_t tmp;
|
||||
slong k, i;
|
||||
|
||||
acb_poly_init(u);
|
||||
acb_poly_init(v);
|
||||
acb_poly_init(tmp);
|
||||
|
||||
if (!regularized)
|
||||
{
|
||||
acb_poly_zero(s);
|
||||
acb_poly_one(t);
|
||||
|
||||
for (k = 0; k < n && acb_poly_length(t) != 0; k++)
|
||||
{
|
||||
acb_poly_add(s, s, t, prec);
|
||||
|
||||
if (p > 0)
|
||||
{
|
||||
acb_poly_add_si(u, a, k, prec);
|
||||
|
||||
for (i = 1; i < p; i++)
|
||||
{
|
||||
acb_poly_add_si(v, a + i, k, prec);
|
||||
acb_poly_mullow(u, u, v, len, prec);
|
||||
}
|
||||
|
||||
acb_poly_mullow(t, t, u, len, prec);
|
||||
}
|
||||
|
||||
if (q > 0)
|
||||
{
|
||||
acb_poly_add_si(u, b, k, prec);
|
||||
|
||||
for (i = 1; i < q; i++)
|
||||
{
|
||||
acb_poly_add_si(v, b + i, k, prec);
|
||||
acb_poly_mullow(u, u, v, len, prec);
|
||||
}
|
||||
|
||||
acb_poly_div_series(t, t, u, len, prec);
|
||||
}
|
||||
|
||||
acb_poly_mullow(t, t, z, len, prec);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_poly_zero(s);
|
||||
|
||||
if (q == 0)
|
||||
acb_poly_one(t);
|
||||
|
||||
for (i = 0; i < q; i++)
|
||||
{
|
||||
if (i == 0)
|
||||
{
|
||||
acb_poly_rgamma_series(t, b + i, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_poly_rgamma_series(u, b + i, len, prec);
|
||||
acb_poly_mullow(tmp, t, u, len, prec);
|
||||
acb_poly_swap(tmp, t);
|
||||
}
|
||||
}
|
||||
|
||||
for (k = 0; k < n; k++)
|
||||
{
|
||||
acb_poly_add(s, s, t, prec);
|
||||
|
||||
if (p > 0)
|
||||
{
|
||||
acb_poly_add_si(u, a, k, prec);
|
||||
|
||||
for (i = 1; i < p; i++)
|
||||
{
|
||||
acb_poly_add_si(v, a + i, k, prec);
|
||||
acb_poly_mullow(tmp, u, v, len, prec);
|
||||
acb_poly_swap(tmp, u);
|
||||
}
|
||||
|
||||
acb_poly_mullow(tmp, t, u, len, prec);
|
||||
acb_poly_swap(tmp, t);
|
||||
}
|
||||
|
||||
if (q > 0)
|
||||
{
|
||||
acb_poly_add_si(u, b, k, prec);
|
||||
|
||||
for (i = 1; i < q; i++)
|
||||
{
|
||||
acb_poly_add_si(v, b + i, k, prec);
|
||||
acb_poly_mullow(tmp, u, v, len, prec);
|
||||
acb_poly_swap(tmp, u);
|
||||
}
|
||||
|
||||
if (acb_poly_length(u) > 0 && !acb_contains_zero(u->coeffs))
|
||||
{
|
||||
acb_poly_div_series(tmp, t, u, len, prec);
|
||||
acb_poly_mullow(t, tmp, z, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* compute term from scratch */
|
||||
acb_poly_one(t);
|
||||
|
||||
for (i = 0; i < p; i++)
|
||||
{
|
||||
acb_poly_rising_ui_series(v, a + i, k + 1, len, prec);
|
||||
acb_poly_mullow(t, t, v, len, prec);
|
||||
}
|
||||
|
||||
for (i = 0; i < q; i++)
|
||||
{
|
||||
acb_poly_add_si(v, b + i, k + 1, prec);
|
||||
acb_poly_rgamma_series(v, v, len, prec);
|
||||
acb_poly_mullow(t, t, v, len, prec);
|
||||
}
|
||||
|
||||
acb_poly_pow_ui_trunc_binexp(v, z, k + 1, len, prec);
|
||||
acb_poly_mullow(t, t, v, len, prec);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_poly_mullow(tmp, t, z, len, prec);
|
||||
acb_poly_swap(tmp, t);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
acb_poly_clear(u);
|
||||
acb_poly_clear(v);
|
||||
acb_poly_clear(tmp);
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_pfq_series_direct(acb_poly_t res,
|
||||
const acb_poly_struct * a, slong p,
|
||||
|
|
@ -362,7 +218,7 @@ acb_hypgeom_pfq_series_direct(acb_poly_t res,
|
|||
arb_poly_init(C);
|
||||
arb_poly_init(T);
|
||||
|
||||
acb_hypgeom_pfq_series_sum_forward(s, t, a, p, b, q, z, regularized, n, len, prec);
|
||||
acb_hypgeom_pfq_series_sum(s, t, a, p, b, q, z, regularized, n, len, prec);
|
||||
|
||||
if (!terminating)
|
||||
{
|
||||
|
|
|
|||
83
acb_hypgeom/pfq_series_sum.c
Normal file
83
acb_hypgeom/pfq_series_sum.c
Normal file
|
|
@ -0,0 +1,83 @@
|
|||
/*
|
||||
Copyright (C) 2016 Fredrik Johansson
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
void
|
||||
acb_hypgeom_pfq_series_sum(acb_poly_t s, acb_poly_t t,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_struct * b, slong q,
|
||||
const acb_poly_t z, int regularized,
|
||||
slong n, slong len, slong prec)
|
||||
{
|
||||
slong abits, zbits, i, j, cb;
|
||||
|
||||
if (n < 4)
|
||||
{
|
||||
acb_hypgeom_pfq_series_sum_forward(s, t, a, p, b, q, z,
|
||||
regularized, n, len, prec);
|
||||
return;
|
||||
}
|
||||
|
||||
abits = 0;
|
||||
zbits = 0;
|
||||
|
||||
for (i = 0; i < p; i++)
|
||||
{
|
||||
for (j = 0; j < FLINT_MIN(acb_poly_length(a + i), n); j++)
|
||||
{
|
||||
cb = acb_bits((a + i)->coeffs + j);
|
||||
abits = FLINT_MAX(abits, cb);
|
||||
}
|
||||
}
|
||||
|
||||
for (i = 0; i < q; i++)
|
||||
{
|
||||
for (j = 0; j < FLINT_MIN(acb_poly_length(b + i), n); j++)
|
||||
{
|
||||
cb = acb_bits((b + i)->coeffs + j);
|
||||
abits = FLINT_MAX(abits, cb);
|
||||
}
|
||||
}
|
||||
|
||||
for (j = 0; j < FLINT_MIN(acb_poly_length(z), n); j++)
|
||||
{
|
||||
cb = acb_bits(z->coeffs + j);
|
||||
zbits = FLINT_MAX(zbits, cb);
|
||||
}
|
||||
|
||||
/* Prefer RS with small coefficients in parameters, large coefficients
|
||||
in z. This only makes sense right now with len <= 2 due to the way
|
||||
acb_poly_div_series works; in the future, it could be used for
|
||||
slightly larger len. */
|
||||
if (len <= 2 && prec > 900 && abits < prec * 0.1 && zbits > prec * 0.1)
|
||||
{
|
||||
acb_hypgeom_pfq_series_sum_rs(s, t, a, p, b, q, z,
|
||||
regularized, n, len, prec);
|
||||
return;
|
||||
}
|
||||
|
||||
/* Prefer BS with small coefficients and high precision, or when
|
||||
computing derivatives of high order. */
|
||||
if ((abits < prec * 0.1 && zbits < prec * 0.1 && prec > 600) || len > 20)
|
||||
{
|
||||
acb_hypgeom_pfq_series_sum_bs(s, t, a, p, b, q, z,
|
||||
regularized, n, len, prec);
|
||||
return;
|
||||
}
|
||||
|
||||
/* TODO: also use bs here when n is large enough, for better
|
||||
numerical stability? */
|
||||
|
||||
acb_hypgeom_pfq_series_sum_forward(s, t, a, p, b, q, z,
|
||||
regularized, n, len, prec);
|
||||
}
|
||||
|
||||
206
acb_hypgeom/pfq_series_sum_bs.c
Normal file
206
acb_hypgeom/pfq_series_sum_bs.c
Normal file
|
|
@ -0,0 +1,206 @@
|
|||
/*
|
||||
Copyright (C) 2016 Fredrik Johansson
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
static void
|
||||
factor(acb_poly_t A, acb_poly_t tmp,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_t z, slong k, slong len, slong prec)
|
||||
{
|
||||
slong i;
|
||||
|
||||
if (p == 0)
|
||||
{
|
||||
if (z == NULL)
|
||||
acb_poly_one(A);
|
||||
else
|
||||
acb_poly_set(A, z);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_poly_add_si(A, a, k, prec);
|
||||
|
||||
for (i = 1; i < p; i++)
|
||||
{
|
||||
acb_poly_add_si(tmp, a + i, k, prec);
|
||||
acb_poly_mullow(A, A, tmp, len, prec);
|
||||
}
|
||||
|
||||
if (z != NULL)
|
||||
{
|
||||
acb_poly_mullow(A, A, z, len, prec);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
static void
|
||||
bsplit(acb_poly_t A1, acb_poly_t B1, acb_poly_t C1,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_struct * b, slong q,
|
||||
const acb_poly_t z,
|
||||
slong aa,
|
||||
slong bb,
|
||||
slong len, slong prec)
|
||||
{
|
||||
if (bb - aa == 1)
|
||||
{
|
||||
factor(A1, B1, a, p, z, aa, len, prec);
|
||||
factor(C1, B1, b, q, NULL, aa, len, prec);
|
||||
/* acb_poly_set(B1, C1); but we skip this */
|
||||
}
|
||||
else
|
||||
{
|
||||
slong m;
|
||||
|
||||
acb_poly_t A2, B2, C2, tmp;
|
||||
|
||||
acb_poly_init(A2);
|
||||
acb_poly_init(B2);
|
||||
acb_poly_init(C2);
|
||||
acb_poly_init(tmp);
|
||||
|
||||
m = aa + (bb - aa) / 2;
|
||||
|
||||
bsplit(A1, B1, C1, a, p, b, q, z, aa, m, len, prec);
|
||||
bsplit(A2, B2, C2, a, p, b, q, z, m, bb, len, prec);
|
||||
|
||||
if (bb - m == 1) /* B2 = C2 */
|
||||
{
|
||||
if (m - aa == 1)
|
||||
acb_poly_add(B2, A1, C1, prec);
|
||||
else
|
||||
acb_poly_add(B2, A1, B1, prec);
|
||||
|
||||
acb_poly_mullow(B1, B2, C2, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
if (m - aa == 1)
|
||||
{
|
||||
acb_poly_mullow(B1, C1, C2, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_poly_mullow(tmp, B1, C2, len, prec);
|
||||
acb_poly_swap(B1, tmp);
|
||||
}
|
||||
|
||||
acb_poly_mullow(tmp, A1, B2, len, prec);
|
||||
acb_poly_add(B1, B1, tmp, prec);
|
||||
}
|
||||
|
||||
acb_poly_mullow(tmp, A1, A2, len, prec);
|
||||
acb_poly_swap(A1, tmp);
|
||||
acb_poly_mullow(tmp, C1, C2, len, prec);
|
||||
acb_poly_swap(C1, tmp);
|
||||
|
||||
acb_poly_clear(A2);
|
||||
acb_poly_clear(B2);
|
||||
acb_poly_clear(C2);
|
||||
acb_poly_clear(tmp);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_pfq_series_sum_bs(acb_poly_t s, acb_poly_t t,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_struct * b, slong q,
|
||||
const acb_poly_t z, int regularized,
|
||||
slong n, slong len, slong prec)
|
||||
{
|
||||
acb_poly_t u, v, w;
|
||||
slong i, start;
|
||||
|
||||
if (n == 0)
|
||||
{
|
||||
acb_hypgeom_pfq_series_sum_forward(s, t, a, p, b, q, z,
|
||||
regularized, n, len, prec);
|
||||
return;
|
||||
}
|
||||
|
||||
start = 0;
|
||||
|
||||
if (regularized)
|
||||
{
|
||||
/*
|
||||
Use basecase algorithm to get past any poles.
|
||||
|
||||
G(x) = 1/Gamma(x)
|
||||
n = 1: s = 0; t = G(b); b = 0 requires skipping
|
||||
n = 2: s = G(b); t = G(b) G(b+1); b = 0, -1 require skipping...
|
||||
|
||||
TODO: when pole is near the end, do the *start* (or middle segment)
|
||||
using fast algorithm; use basecase algorithm for the end.
|
||||
*/
|
||||
|
||||
for (i = 0; i < q; i++)
|
||||
{
|
||||
if (acb_poly_is_zero(b + i))
|
||||
{
|
||||
start = FLINT_MAX(start, 1);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* todo: use a fuzzier test? */
|
||||
if (acb_contains_int((b + i)->coeffs) &&
|
||||
!arb_is_positive(acb_realref((b + i)->coeffs)) &&
|
||||
arf_cmpabs_2exp_si(arb_midref(acb_realref((b + i)->coeffs)),
|
||||
FLINT_BITS - 2) < 0)
|
||||
{
|
||||
slong c = -arf_get_si(arb_midref(acb_realref((b + i)->coeffs)),
|
||||
ARF_RND_NEAR);
|
||||
|
||||
/* if c >= n, terminates earlier, so no problem */
|
||||
if (c < n)
|
||||
{
|
||||
start = FLINT_MAX(start, c + 1);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* We should now have start <= n. */
|
||||
if (start > n) abort();
|
||||
|
||||
acb_hypgeom_pfq_series_sum_forward(s, t, a, p, b, q, z,
|
||||
regularized, start, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_poly_zero(s);
|
||||
acb_poly_one(t);
|
||||
}
|
||||
|
||||
if (start == n)
|
||||
return;
|
||||
|
||||
acb_poly_init(u);
|
||||
acb_poly_init(v);
|
||||
acb_poly_init(w);
|
||||
|
||||
bsplit(u, v, w, a, p, b, q, z, start, n, len, prec);
|
||||
|
||||
if (n - start == 1)
|
||||
acb_poly_set(v, w); /* B1 not set */
|
||||
|
||||
acb_poly_mullow(v, v, t, len, prec);
|
||||
acb_poly_div_series(v, v, w, len, prec);
|
||||
acb_poly_add(s, s, v, prec);
|
||||
|
||||
acb_poly_mullow(t, t, u, len, prec);
|
||||
acb_poly_div_series(t, t, w, len, prec);
|
||||
|
||||
acb_poly_clear(u);
|
||||
acb_poly_clear(v);
|
||||
acb_poly_clear(w);
|
||||
}
|
||||
|
||||
157
acb_hypgeom/pfq_series_sum_forward.c
Normal file
157
acb_hypgeom/pfq_series_sum_forward.c
Normal file
|
|
@ -0,0 +1,157 @@
|
|||
/*
|
||||
Copyright (C) 2015 Fredrik Johansson
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
void
|
||||
acb_hypgeom_pfq_series_sum_forward(acb_poly_t s, acb_poly_t t,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_struct * b, slong q,
|
||||
const acb_poly_t z, int regularized,
|
||||
slong n, slong len, slong prec)
|
||||
{
|
||||
acb_poly_t u, v;
|
||||
acb_poly_t tmp;
|
||||
slong k, i;
|
||||
|
||||
acb_poly_init(u);
|
||||
acb_poly_init(v);
|
||||
acb_poly_init(tmp);
|
||||
|
||||
if (!regularized)
|
||||
{
|
||||
acb_poly_zero(s);
|
||||
acb_poly_one(t);
|
||||
|
||||
for (k = 0; k < n && acb_poly_length(t) != 0; k++)
|
||||
{
|
||||
acb_poly_add(s, s, t, prec);
|
||||
|
||||
if (p > 0)
|
||||
{
|
||||
acb_poly_add_si(u, a, k, prec);
|
||||
|
||||
for (i = 1; i < p; i++)
|
||||
{
|
||||
acb_poly_add_si(v, a + i, k, prec);
|
||||
acb_poly_mullow(u, u, v, len, prec);
|
||||
}
|
||||
|
||||
acb_poly_mullow(t, t, u, len, prec);
|
||||
}
|
||||
|
||||
if (q > 0)
|
||||
{
|
||||
acb_poly_add_si(u, b, k, prec);
|
||||
|
||||
for (i = 1; i < q; i++)
|
||||
{
|
||||
acb_poly_add_si(v, b + i, k, prec);
|
||||
acb_poly_mullow(u, u, v, len, prec);
|
||||
}
|
||||
|
||||
acb_poly_div_series(t, t, u, len, prec);
|
||||
}
|
||||
|
||||
acb_poly_mullow(t, t, z, len, prec);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_poly_zero(s);
|
||||
|
||||
if (q == 0)
|
||||
acb_poly_one(t);
|
||||
|
||||
for (i = 0; i < q; i++)
|
||||
{
|
||||
if (i == 0)
|
||||
{
|
||||
acb_poly_rgamma_series(t, b + i, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_poly_rgamma_series(u, b + i, len, prec);
|
||||
acb_poly_mullow(tmp, t, u, len, prec);
|
||||
acb_poly_swap(tmp, t);
|
||||
}
|
||||
}
|
||||
|
||||
for (k = 0; k < n; k++)
|
||||
{
|
||||
acb_poly_add(s, s, t, prec);
|
||||
|
||||
if (p > 0)
|
||||
{
|
||||
acb_poly_add_si(u, a, k, prec);
|
||||
|
||||
for (i = 1; i < p; i++)
|
||||
{
|
||||
acb_poly_add_si(v, a + i, k, prec);
|
||||
acb_poly_mullow(tmp, u, v, len, prec);
|
||||
acb_poly_swap(tmp, u);
|
||||
}
|
||||
|
||||
acb_poly_mullow(tmp, t, u, len, prec);
|
||||
acb_poly_swap(tmp, t);
|
||||
}
|
||||
|
||||
if (q > 0)
|
||||
{
|
||||
acb_poly_add_si(u, b, k, prec);
|
||||
|
||||
for (i = 1; i < q; i++)
|
||||
{
|
||||
acb_poly_add_si(v, b + i, k, prec);
|
||||
acb_poly_mullow(tmp, u, v, len, prec);
|
||||
acb_poly_swap(tmp, u);
|
||||
}
|
||||
|
||||
if (acb_poly_length(u) > 0 && !acb_contains_zero(u->coeffs))
|
||||
{
|
||||
acb_poly_div_series(tmp, t, u, len, prec);
|
||||
acb_poly_mullow(t, tmp, z, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* compute term from scratch */
|
||||
acb_poly_one(t);
|
||||
|
||||
for (i = 0; i < p; i++)
|
||||
{
|
||||
acb_poly_rising_ui_series(v, a + i, k + 1, len, prec);
|
||||
acb_poly_mullow(t, t, v, len, prec);
|
||||
}
|
||||
|
||||
for (i = 0; i < q; i++)
|
||||
{
|
||||
acb_poly_add_si(v, b + i, k + 1, prec);
|
||||
acb_poly_rgamma_series(v, v, len, prec);
|
||||
acb_poly_mullow(t, t, v, len, prec);
|
||||
}
|
||||
|
||||
acb_poly_pow_ui_trunc_binexp(v, z, k + 1, len, prec);
|
||||
acb_poly_mullow(t, t, v, len, prec);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_poly_mullow(tmp, t, z, len, prec);
|
||||
acb_poly_swap(tmp, t);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
acb_poly_clear(u);
|
||||
acb_poly_clear(v);
|
||||
acb_poly_clear(tmp);
|
||||
}
|
||||
|
||||
289
acb_hypgeom/pfq_series_sum_rs.c
Normal file
289
acb_hypgeom/pfq_series_sum_rs.c
Normal file
|
|
@ -0,0 +1,289 @@
|
|||
/*
|
||||
Copyright (C) 2016 Fredrik Johansson
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
#define FAST_T 0
|
||||
|
||||
void acb_poly_reciprocal_majorant(arb_poly_t res, const acb_poly_t poly, slong prec);
|
||||
|
||||
static void
|
||||
rsplit(acb_poly_t res, acb_poly_t term,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_struct * b, slong q,
|
||||
const acb_poly_t z, slong offset, slong n, slong len, slong prec)
|
||||
{
|
||||
acb_poly_struct * zpow;
|
||||
acb_poly_t s, t, u;
|
||||
slong i, j, k, m, tprec;
|
||||
int is_real;
|
||||
|
||||
#if FAST_T
|
||||
arb_poly_t B, C, D;
|
||||
#else
|
||||
acb_poly_t B, C, D;
|
||||
#endif
|
||||
|
||||
if (n == 0)
|
||||
{
|
||||
acb_poly_zero(res);
|
||||
acb_poly_one(term);
|
||||
return;
|
||||
}
|
||||
|
||||
if (n < 0)
|
||||
abort();
|
||||
|
||||
m = n_sqrt(n);
|
||||
m = FLINT_MIN(m, 150);
|
||||
|
||||
acb_poly_init(s);
|
||||
acb_poly_init(t);
|
||||
acb_poly_init(u);
|
||||
|
||||
#if FAST_T
|
||||
tprec = MAG_BITS;
|
||||
|
||||
arb_poly_init(B);
|
||||
arb_poly_init(C);
|
||||
arb_poly_init(D);
|
||||
|
||||
arb_poly_one(B);
|
||||
arb_poly_one(D);
|
||||
|
||||
is_real = 1;
|
||||
|
||||
for (i = 0; i < p; i++)
|
||||
for (j = 0; j < FLINT_MIN(acb_poly_length(a + i), len); j++)
|
||||
is_real = is_real && acb_is_real((a + i)->coeffs + j);
|
||||
|
||||
for (i = 0; i < q; i++)
|
||||
for (j = 0; j < FLINT_MIN(acb_poly_length(b + i), len); j++)
|
||||
is_real = is_real && acb_is_real((b + i)->coeffs + j);
|
||||
|
||||
for (j = 0; j < FLINT_MIN(acb_poly_length(z), len); j++)
|
||||
is_real = is_real && acb_is_real(z->coeffs + j);
|
||||
#else
|
||||
tprec = 2 * MAG_BITS;
|
||||
is_real = 0; (void) is_real;
|
||||
|
||||
acb_poly_init(B);
|
||||
acb_poly_init(C);
|
||||
acb_poly_init(D);
|
||||
|
||||
acb_poly_one(B);
|
||||
acb_poly_one(D);
|
||||
#endif
|
||||
|
||||
zpow = flint_malloc(sizeof(acb_poly_struct) * (m + 1));
|
||||
|
||||
for (i = 0; i <= m; i++)
|
||||
acb_poly_init(zpow + i);
|
||||
|
||||
for (i = 0; i <= m; i++)
|
||||
{
|
||||
if (i == 0)
|
||||
acb_poly_one(zpow + i);
|
||||
else if (i == 1)
|
||||
acb_poly_set_round(zpow + i, z, prec);
|
||||
else if (i % 2 == 0)
|
||||
acb_poly_mullow(zpow + i, zpow + i / 2, zpow + i / 2, len, prec);
|
||||
else
|
||||
acb_poly_mullow(zpow + i, zpow + i - 1, zpow + 1, len, prec);
|
||||
}
|
||||
|
||||
for (k = n; k >= 0; k--)
|
||||
{
|
||||
j = k % m;
|
||||
|
||||
if (k < n)
|
||||
acb_poly_add(s, s, zpow + j, prec);
|
||||
|
||||
if (k > 0)
|
||||
{
|
||||
if (p > 0)
|
||||
{
|
||||
acb_poly_add_si(u, a, offset + k - 1, prec);
|
||||
|
||||
for (i = 1; i < p; i++)
|
||||
{
|
||||
acb_poly_add_si(t, a + i, offset + k - 1, prec);
|
||||
acb_poly_mullow(u, u, t, len, prec);
|
||||
}
|
||||
|
||||
if (k < n)
|
||||
acb_poly_mullow(s, s, u, len, prec);
|
||||
|
||||
#if FAST_T
|
||||
acb_poly_majorant(C, u, tprec);
|
||||
arb_poly_mullow(B, B, C, len, tprec);
|
||||
#else
|
||||
acb_poly_set_round(C, u, tprec);
|
||||
acb_poly_mullow(B, B, C, len, tprec);
|
||||
|
||||
#endif
|
||||
}
|
||||
|
||||
if (q > 0)
|
||||
{
|
||||
acb_poly_add_si(u, b, offset + k - 1, prec);
|
||||
|
||||
for (i = 1; i < q; i++)
|
||||
{
|
||||
acb_poly_add_si(t, b + i, offset + k - 1, prec);
|
||||
acb_poly_mullow(u, u, t, len, prec);
|
||||
}
|
||||
|
||||
if (k < n)
|
||||
acb_poly_div_series(s, s, u, len, prec);
|
||||
|
||||
#if FAST_T
|
||||
acb_poly_reciprocal_majorant(C, u, tprec);
|
||||
arb_poly_mullow(D, D, C, len, tprec);
|
||||
#else
|
||||
acb_poly_set_round(C, u, tprec);
|
||||
acb_poly_mullow(D, D, C, len, tprec);
|
||||
#endif
|
||||
}
|
||||
|
||||
if (j == 0 && k < n)
|
||||
{
|
||||
acb_poly_mullow(s, s, zpow + m, len, prec);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#if FAST_T
|
||||
arb_poly_div_series(B, B, D, len, tprec);
|
||||
acb_poly_majorant(C, z, tprec);
|
||||
arb_poly_pow_ui_trunc_binexp(C, C, n, len, tprec);
|
||||
arb_poly_mullow(C, B, C, len, tprec);
|
||||
|
||||
acb_poly_fit_length(term, arb_poly_length(C));
|
||||
|
||||
for (i = 0; i < arb_poly_length(C); i++)
|
||||
{
|
||||
arb_get_mag(arb_radref(acb_realref(term->coeffs + i)), C->coeffs + i);
|
||||
arf_zero(arb_midref(acb_realref(term->coeffs + i)));
|
||||
|
||||
if (is_real)
|
||||
arb_zero(acb_imagref(term->coeffs + i));
|
||||
else
|
||||
arb_set(acb_imagref(term->coeffs + i), acb_realref(term->coeffs + i));
|
||||
}
|
||||
|
||||
_acb_poly_set_length(term, arb_poly_length(C));
|
||||
_acb_poly_normalise(term);
|
||||
|
||||
#else
|
||||
acb_poly_div_series(B, B, D, len, tprec);
|
||||
acb_poly_set_round(C, z, tprec);
|
||||
acb_poly_pow_ui_trunc_binexp(C, C, n, len, tprec);
|
||||
acb_poly_mullow(term, B, C, len, tprec);
|
||||
#endif
|
||||
|
||||
acb_poly_set(res, s);
|
||||
|
||||
#if FAST_T
|
||||
arb_poly_clear(B);
|
||||
arb_poly_clear(C);
|
||||
arb_poly_clear(D);
|
||||
#else
|
||||
acb_poly_clear(B);
|
||||
acb_poly_clear(C);
|
||||
acb_poly_clear(D);
|
||||
#endif
|
||||
|
||||
acb_poly_clear(s);
|
||||
acb_poly_clear(t);
|
||||
acb_poly_clear(u);
|
||||
|
||||
for (i = 0; i <= m; i++)
|
||||
acb_poly_clear(zpow + i);
|
||||
flint_free(zpow);
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_pfq_series_sum_rs(acb_poly_t s, acb_poly_t t,
|
||||
const acb_poly_struct * a, slong p,
|
||||
const acb_poly_struct * b, slong q,
|
||||
const acb_poly_t z, int regularized,
|
||||
slong n, slong len, slong prec)
|
||||
{
|
||||
acb_poly_t u, v;
|
||||
slong i, start;
|
||||
|
||||
if (n == 0)
|
||||
{
|
||||
acb_hypgeom_pfq_series_sum_forward(s, t, a, p, b, q, z,
|
||||
regularized, n, len, prec);
|
||||
return;
|
||||
}
|
||||
|
||||
start = 0;
|
||||
|
||||
if (regularized)
|
||||
{
|
||||
for (i = 0; i < q; i++)
|
||||
{
|
||||
if (acb_poly_is_zero(b + i))
|
||||
{
|
||||
start = FLINT_MAX(start, 1);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* todo: use a fuzzier test? */
|
||||
if (acb_contains_int((b + i)->coeffs) &&
|
||||
!arb_is_positive(acb_realref((b + i)->coeffs)) &&
|
||||
arf_cmpabs_2exp_si(arb_midref(acb_realref((b + i)->coeffs)),
|
||||
FLINT_BITS - 2) < 0)
|
||||
{
|
||||
slong c = -arf_get_si(arb_midref(acb_realref((b + i)->coeffs)),
|
||||
ARF_RND_NEAR);
|
||||
|
||||
/* if c >= n, terminates earlier, so no problem */
|
||||
if (c < n)
|
||||
{
|
||||
start = FLINT_MAX(start, c + 1);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* We should now have start <= n. */
|
||||
if (start > n) abort();
|
||||
|
||||
acb_hypgeom_pfq_series_sum_forward(s, t, a, p, b, q, z,
|
||||
regularized, start, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_poly_zero(s);
|
||||
acb_poly_one(t);
|
||||
}
|
||||
|
||||
if (start == n)
|
||||
return;
|
||||
|
||||
acb_poly_init(u);
|
||||
acb_poly_init(v);
|
||||
|
||||
rsplit(u, v, a, p, b, q, z, start, n - start, len, prec);
|
||||
|
||||
acb_poly_mullow(u, u, t, len, prec);
|
||||
acb_poly_add(s, s, u, prec);
|
||||
|
||||
acb_poly_mullow(t, t, v, len, prec);
|
||||
|
||||
acb_poly_clear(u);
|
||||
acb_poly_clear(v);
|
||||
}
|
||||
|
||||
116
acb_hypgeom/test/t-pfq_series_sum_bs.c
Normal file
116
acb_hypgeom/test/t-pfq_series_sum_bs.c
Normal file
|
|
@ -0,0 +1,116 @@
|
|||
/*
|
||||
Copyright (C) 2016 Fredrik Johansson
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
int main()
|
||||
{
|
||||
slong iter;
|
||||
flint_rand_t state;
|
||||
|
||||
flint_printf("pfq_series_sum_bs....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
|
||||
{
|
||||
acb_poly_struct a[4], b[4];
|
||||
acb_poly_t z, s1, s2, t1, t2;
|
||||
slong i, p, q, len1, len2, n, prec1, prec2;
|
||||
int regularized;
|
||||
|
||||
p = n_randint(state, 3);
|
||||
q = n_randint(state, 3);
|
||||
|
||||
len1 = n_randint(state, 8);
|
||||
len2 = n_randint(state, 8);
|
||||
|
||||
prec1 = 2 + n_randint(state, 400);
|
||||
prec2 = 2 + n_randint(state, 400);
|
||||
|
||||
n = n_randint(state, 50);
|
||||
|
||||
regularized = n_randint(state, 2);
|
||||
|
||||
acb_poly_init(z);
|
||||
acb_poly_init(s1);
|
||||
acb_poly_init(s2);
|
||||
acb_poly_init(t1);
|
||||
acb_poly_init(t2);
|
||||
for (i = 0; i < p; i++)
|
||||
acb_poly_init(a + i);
|
||||
for (i = 0; i < q; i++)
|
||||
acb_poly_init(b + i);
|
||||
|
||||
acb_poly_randtest(z, state, 1 + n_randint(state, 10), 1 + n_randint(state, 500), 10);
|
||||
|
||||
for (i = 0; i < p; i++)
|
||||
acb_poly_randtest(a + i, state, 1 + n_randint(state, 10), 1 + n_randint(state, 500), 3);
|
||||
for (i = 0; i < q; i++)
|
||||
acb_poly_randtest(b + i, state, 1 + n_randint(state, 10), 1 + n_randint(state, 500), 3);
|
||||
|
||||
acb_hypgeom_pfq_series_sum_forward(s1, t1, a, p, b, q, z, regularized, n, len1, prec1);
|
||||
acb_hypgeom_pfq_series_sum_bs(s2, t2, a, p, b, q, z, regularized, n, len2, prec2);
|
||||
|
||||
acb_poly_truncate(s1, FLINT_MIN(len1, len2));
|
||||
acb_poly_truncate(s2, FLINT_MIN(len1, len2));
|
||||
acb_poly_truncate(t1, FLINT_MIN(len1, len2));
|
||||
acb_poly_truncate(t2, FLINT_MIN(len1, len2));
|
||||
|
||||
if (!acb_poly_overlaps(s1, s2) || !acb_poly_overlaps(t1, t2))
|
||||
{
|
||||
flint_printf("FAIL: overlap\n\n");
|
||||
flint_printf("iter = %wd\n", iter);
|
||||
flint_printf("n = %wd, len1 = %wd, len2 = %wd, prec1 = %wd, prec2 = %wd\n\n", n, len1, len2, prec1, prec2);
|
||||
flint_printf("p = %wd, q = %wd\n\n", p, q);
|
||||
flint_printf("z = "); acb_poly_printd(z, 15); flint_printf("\n\n");
|
||||
|
||||
for (i = 0; i < p; i++)
|
||||
{
|
||||
flint_printf("a[%wd] = ", i); acb_poly_printd(a + i, 15); flint_printf("\n\n");
|
||||
}
|
||||
|
||||
for (i = 0; i < q; i++)
|
||||
{
|
||||
flint_printf("b[%wd] = ", i); acb_poly_printd(b + i, 15); flint_printf("\n\n");
|
||||
}
|
||||
|
||||
flint_printf("s1 = "); acb_poly_printd(s1, 15); flint_printf("\n\n");
|
||||
flint_printf("s2 = "); acb_poly_printd(s2, 15); flint_printf("\n\n");
|
||||
acb_poly_sub(s1, s1, s2, prec1);
|
||||
flint_printf("diff = "); acb_poly_printd(s1, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("t1 = "); acb_poly_printd(t1, 15); flint_printf("\n\n");
|
||||
flint_printf("t2 = "); acb_poly_printd(t2, 15); flint_printf("\n\n");
|
||||
acb_poly_sub(t1, t1, t2, prec1);
|
||||
flint_printf("diff = "); acb_poly_printd(t1, 15); flint_printf("\n\n");
|
||||
|
||||
abort();
|
||||
}
|
||||
|
||||
acb_poly_clear(z);
|
||||
acb_poly_clear(s1);
|
||||
acb_poly_clear(s2);
|
||||
for (i = 0; i < p; i++)
|
||||
acb_poly_clear(a + i);
|
||||
for (i = 0; i < q; i++)
|
||||
acb_poly_clear(b + i);
|
||||
acb_poly_clear(t1);
|
||||
acb_poly_clear(t2);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
116
acb_hypgeom/test/t-pfq_series_sum_rs.c
Normal file
116
acb_hypgeom/test/t-pfq_series_sum_rs.c
Normal file
|
|
@ -0,0 +1,116 @@
|
|||
/*
|
||||
Copyright (C) 2016 Fredrik Johansson
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
int main()
|
||||
{
|
||||
slong iter;
|
||||
flint_rand_t state;
|
||||
|
||||
flint_printf("pfq_series_sum_rs....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
|
||||
{
|
||||
acb_poly_struct a[4], b[4];
|
||||
acb_poly_t z, s1, s2, t1, t2;
|
||||
slong i, p, q, len1, len2, n, prec1, prec2;
|
||||
int regularized;
|
||||
|
||||
p = n_randint(state, 3);
|
||||
q = n_randint(state, 3);
|
||||
|
||||
len1 = n_randint(state, 8);
|
||||
len2 = n_randint(state, 8);
|
||||
|
||||
prec1 = 2 + n_randint(state, 400);
|
||||
prec2 = 2 + n_randint(state, 400);
|
||||
|
||||
n = n_randint(state, 50);
|
||||
|
||||
regularized = n_randint(state, 2);
|
||||
|
||||
acb_poly_init(z);
|
||||
acb_poly_init(s1);
|
||||
acb_poly_init(s2);
|
||||
acb_poly_init(t1);
|
||||
acb_poly_init(t2);
|
||||
for (i = 0; i < p; i++)
|
||||
acb_poly_init(a + i);
|
||||
for (i = 0; i < q; i++)
|
||||
acb_poly_init(b + i);
|
||||
|
||||
acb_poly_randtest(z, state, 1 + n_randint(state, 10), 1 + n_randint(state, 500), 10);
|
||||
|
||||
for (i = 0; i < p; i++)
|
||||
acb_poly_randtest(a + i, state, 1 + n_randint(state, 10), 1 + n_randint(state, 500), 3);
|
||||
for (i = 0; i < q; i++)
|
||||
acb_poly_randtest(b + i, state, 1 + n_randint(state, 10), 1 + n_randint(state, 500), 3);
|
||||
|
||||
acb_hypgeom_pfq_series_sum_forward(s1, t1, a, p, b, q, z, regularized, n, len1, prec1);
|
||||
acb_hypgeom_pfq_series_sum_rs(s2, t2, a, p, b, q, z, regularized, n, len2, prec2);
|
||||
|
||||
acb_poly_truncate(s1, FLINT_MIN(len1, len2));
|
||||
acb_poly_truncate(s2, FLINT_MIN(len1, len2));
|
||||
acb_poly_truncate(t1, FLINT_MIN(len1, len2));
|
||||
acb_poly_truncate(t2, FLINT_MIN(len1, len2));
|
||||
|
||||
if (!acb_poly_overlaps(s1, s2) || !acb_poly_overlaps(t1, t2))
|
||||
{
|
||||
flint_printf("FAIL: overlap\n\n");
|
||||
flint_printf("iter = %wd\n", iter);
|
||||
flint_printf("n = %wd, len1 = %wd, len2 = %wd, prec1 = %wd, prec2 = %wd\n\n", n, len1, len2, prec1, prec2);
|
||||
flint_printf("p = %wd, q = %wd\n\n", p, q);
|
||||
flint_printf("z = "); acb_poly_printd(z, 15); flint_printf("\n\n");
|
||||
|
||||
for (i = 0; i < p; i++)
|
||||
{
|
||||
flint_printf("a[%wd] = ", i); acb_poly_printd(a + i, 15); flint_printf("\n\n");
|
||||
}
|
||||
|
||||
for (i = 0; i < q; i++)
|
||||
{
|
||||
flint_printf("b[%wd] = ", i); acb_poly_printd(b + i, 15); flint_printf("\n\n");
|
||||
}
|
||||
|
||||
flint_printf("s1 = "); acb_poly_printd(s1, 15); flint_printf("\n\n");
|
||||
flint_printf("s2 = "); acb_poly_printd(s2, 15); flint_printf("\n\n");
|
||||
acb_poly_sub(s1, s1, s2, prec1);
|
||||
flint_printf("diff = "); acb_poly_printd(s1, 15); flint_printf("\n\n");
|
||||
|
||||
flint_printf("t1 = "); acb_poly_printd(t1, 15); flint_printf("\n\n");
|
||||
flint_printf("t2 = "); acb_poly_printd(t2, 15); flint_printf("\n\n");
|
||||
acb_poly_sub(t1, t1, t2, prec1);
|
||||
flint_printf("diff = "); acb_poly_printd(t1, 15); flint_printf("\n\n");
|
||||
|
||||
abort();
|
||||
}
|
||||
|
||||
acb_poly_clear(z);
|
||||
acb_poly_clear(s1);
|
||||
acb_poly_clear(s2);
|
||||
for (i = 0; i < p; i++)
|
||||
acb_poly_clear(a + i);
|
||||
for (i = 0; i < q; i++)
|
||||
acb_poly_clear(b + i);
|
||||
acb_poly_clear(t1);
|
||||
acb_poly_clear(t2);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
|
|
@ -110,11 +110,32 @@ or remove a 1 from the `a_i` parameters if there is one.
|
|||
If `n < 0`, this function chooses a number of terms automatically
|
||||
using :func:`acb_hypgeom_pfq_choose_n`.
|
||||
|
||||
.. function:: void acb_hypgeom_pfq_series_sum_forward(acb_poly_t s, acb_poly_t t, const acb_poly_struct * a, slong p, const acb_poly_struct * b, slong q, const acb_poly_t z, int regularized, slong n, slong len, slong prec)
|
||||
|
||||
.. function:: void acb_hypgeom_pfq_series_sum_bs(acb_poly_t s, acb_poly_t t, const acb_poly_struct * a, slong p, const acb_poly_struct * b, slong q, const acb_poly_t z, int regularized, slong n, slong len, slong prec)
|
||||
|
||||
.. function:: void acb_hypgeom_pfq_series_sum_rs(acb_poly_t s, acb_poly_t t, const acb_poly_struct * a, slong p, const acb_poly_struct * b, slong q, const acb_poly_t z, int regularized, slong n, slong len, slong prec)
|
||||
|
||||
.. function:: void acb_hypgeom_pfq_series_sum(acb_poly_t s, acb_poly_t t, const acb_poly_struct * a, slong p, const acb_poly_struct * b, slong q, const acb_poly_t z, int regularized, slong n, slong len, slong prec)
|
||||
|
||||
Computes `s = \sum_{k=0}^{n-1} T(k)` and `t = T(n)` given parameters
|
||||
and argument that are power series.
|
||||
Does not allow aliasing between input and output variables.
|
||||
We require `n \ge 0` and that *len* is positive.
|
||||
|
||||
If *regularized* is set, the regularized sum is computed, avoiding
|
||||
division by zero at the poles of the gamma function.
|
||||
|
||||
The *forward*, *bs*, *rs* and default versions use forward recurrence,
|
||||
binary splitting, rectangular splitting, and an automatic algorithm
|
||||
choice.
|
||||
|
||||
.. function:: void acb_hypgeom_pfq_series_direct(acb_poly_t res, const acb_poly_struct * a, slong p, const acb_poly_struct * b, slong q, const acb_poly_t z, int regularized, slong n, slong len, slong prec)
|
||||
|
||||
Computes `{}_pf_{q}(z)` directly using the defining series, given
|
||||
parameters and argument that are power series.
|
||||
The result is a power series of length *len*.
|
||||
We require that *len* is positive.
|
||||
|
||||
An error bound is computed automatically as a function of the number
|
||||
of terms *n*. If `n < 0`, the number of terms is chosen
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue