mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
acb_dirichlet_l_fmpq_afe: handle the singular series
This commit is contained in:
fredrik
parent
11de2dd2f6
commit
b924a3afcd
6 changed files with 277 additions and 36 deletions
|
|
@ -135,13 +135,12 @@ acb_dirichlet_fmpq_sum_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G,
|
|||
arb_t ns, t, u, v, z, z0, z1, x, x2, Ga, Gz1, Gz0, expmz0, z0_prevn, Gz0_prevn, expmz0_prevn;
|
||||
acb_t c;
|
||||
fmpq_t s2;
|
||||
int parity;
|
||||
int parity, gamma_singular;
|
||||
ulong q;
|
||||
|
||||
double abs_tol_mag;
|
||||
double gammainc_mag, gamma_mag, ns_mag;
|
||||
double aa, zz;
|
||||
double cancellation;
|
||||
|
||||
#if VERBOSE
|
||||
double t1, t2, t3;
|
||||
|
|
@ -188,6 +187,10 @@ acb_dirichlet_fmpq_sum_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G,
|
|||
/* s2 = (s+parity)/2 */
|
||||
fmpq_add_ui(s2, s, parity);
|
||||
fmpq_div_2exp(s2, s2, 1);
|
||||
|
||||
gamma_singular = (fmpz_is_one(fmpq_denref(s2)) && fmpz_sgn(fmpq_numref(s2)) <= 0);
|
||||
|
||||
if (!gamma_singular)
|
||||
arb_gamma_fmpq(Ga, s2, gamma_cached_prec);
|
||||
|
||||
for (n = 1; ; n += 1)
|
||||
|
|
@ -222,28 +225,34 @@ acb_dirichlet_fmpq_sum_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G,
|
|||
|
||||
/* Gamma((s+parity)/2, z) (want lower bound, to estimate cancellation) */
|
||||
gammainc_mag = log_gamma_upper_approx(aa, zz) / log(2);
|
||||
/* Gamma((s+parity)/2) */
|
||||
gamma_mag = ARF_EXP(arb_midref(Ga));
|
||||
/* n^-s */
|
||||
ns_mag = -fmpq_get_d(s) * log(n) / log(2);
|
||||
cancellation = FLINT_MAX(gamma_mag - gammainc_mag, 0);
|
||||
cancellation = cancellation; /* suppress warning when VERBOSE == 0 */
|
||||
|
||||
/* Want Gamma(a,z) n^-s with abs_tol --> want Gamma(a,z) with abs_tol * n^s */
|
||||
mag_set_ui_2exp_si(abs_tol_gamma, 1, abs_tol_mag - ns_mag);
|
||||
|
||||
/* Precision needed sans cancellation. */
|
||||
/* wp = Precision needed sans cancellation. */
|
||||
wp = gammainc_mag + ns_mag - abs_tol_mag + 5;
|
||||
wp = FLINT_MAX(wp, 30);
|
||||
|
||||
/* Precision needed with cancellation. */
|
||||
/* wp2 = Precision needed with cancellation. */
|
||||
if (gamma_singular)
|
||||
{
|
||||
/* Max term is roughly n^(-s) * z^((s+parity)/2) * exp(z) */
|
||||
wp2 = -abs_tol_mag + ns_mag + aa * log(zz) + zz / log(2) + 5;
|
||||
wp2 = FLINT_MAX(wp2, 30);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* Estimate of Gamma((s+parity)/2) */
|
||||
gamma_mag = ARF_EXP(arb_midref(Ga));
|
||||
wp2 = FLINT_MAX(gamma_mag, gammainc_mag) + ns_mag - abs_tol_mag + 5;
|
||||
wp2 = FLINT_MAX(wp2, 30);
|
||||
}
|
||||
|
||||
#if VERBOSE
|
||||
printf(" abs_tol_gamma = "); mag_printd(abs_tol_gamma, 5); printf("\n");
|
||||
printf(" gamma(a) = "); arb_printd(Ga, 10); printf("\n");
|
||||
printf(" cancellation = %f\n", cancellation);
|
||||
printf(" wp = %ld wp2 = %ld\n", wp, wp2);
|
||||
#endif
|
||||
|
||||
|
|
@ -286,34 +295,56 @@ acb_dirichlet_fmpq_sum_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G,
|
|||
else
|
||||
{
|
||||
/* Otherwise fallback to the series at 0. */
|
||||
if (gamma_singular)
|
||||
{
|
||||
slong nn;
|
||||
|
||||
nn = *fmpq_numref(s2);
|
||||
|
||||
if (COEFF_IS_MPZ(nn))
|
||||
{
|
||||
arb_indeterminate(Gz0);
|
||||
}
|
||||
else
|
||||
{
|
||||
nn = -nn;
|
||||
|
||||
NN = _arb_hypgeom_gamma_upper_singular_si_choose_N(AE, nn, z0, abs_tol_gamma);
|
||||
_arb_hypgeom_gamma_upper_singular_si_bsplit(Gz0, nn, z0, NN, wp2);
|
||||
arb_add_error_mag(Gz0, AE);
|
||||
|
||||
#if VERBOSE
|
||||
flint_printf(" singular series with N = %wd, z0 = ", NN); arb_printd(z0, 10); printf(" ");
|
||||
mag_printd(AE, 10); printf("\n");
|
||||
#endif
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
NN = _arb_hypgeom_gamma_lower_fmpq_0_choose_N(AE, s2, z0, abs_tol_gamma);
|
||||
_arb_hypgeom_gamma_lower_fmpq_0_bsplit(Gz0, s2, z0, NN, wp2);
|
||||
arb_add_error_mag(Gz0, AE);
|
||||
|
||||
#if VERBOSE
|
||||
flint_printf(" lower series with N = %wd, z0 = ", NN); arb_printd(z0, 10); printf(" ");
|
||||
mag_printd(AE, 10); printf("\n");
|
||||
#endif
|
||||
|
||||
_arb_hypgeom_gamma_lower_fmpq_0_bsplit(Gz0, s2, z0, NN, wp2);
|
||||
arb_add_error_mag(Gz0, AE);
|
||||
|
||||
while (mag_cmp(arb_radref(Ga), abs_tol_gamma) > 0)
|
||||
{
|
||||
gamma_cached_prec *= 2;
|
||||
arb_gamma_fmpq(Ga, s2, gamma_cached_prec);
|
||||
}
|
||||
|
||||
|
||||
#if VERBOSE
|
||||
flint_printf(" lower series with N = %wd: ", NN); arb_printd(Gz0, 10); printf("\n");
|
||||
#endif
|
||||
|
||||
arb_sub(Gz0, Ga, Gz0, wp);
|
||||
|
||||
#if VERBOSE
|
||||
flint_printf(" G(a) - lower series: "); arb_printd(Gz0, 10); printf("\n");
|
||||
#endif
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
_arb_gamma_upper_fmpq_step_bsplit(Gz0, s2, z0_prevn, z0, Gz0_prevn, expmz0_prevn, abs_tol_gamma, wp);
|
||||
|
|
@ -437,17 +468,19 @@ acb_dirichlet_l_fmpq_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G, c
|
|||
q = (G == NULL) ? 1 : G->q;
|
||||
parity = (G == NULL) ? 0 : dirichlet_parity_char(G, chi);
|
||||
|
||||
/* Limit at gamma(-n) is not implemented. */
|
||||
/* Division by gamma((s+parity)/2) at a pole. */
|
||||
if (fmpz_is_one(fmpq_denref(s)))
|
||||
{
|
||||
const fmpz * n = fmpq_numref(s);
|
||||
|
||||
if ((parity == 0 && fmpz_sgn(n) <= 0 && fmpz_is_even(n)) ||
|
||||
(parity == 0 && fmpz_sgn(n) > 0 && fmpz_is_odd(n)) ||
|
||||
(parity == 1 && fmpz_sgn(n) < 0 && fmpz_is_odd(n)) ||
|
||||
(parity == 1 && fmpz_sgn(n) > 0 && fmpz_is_even(n)))
|
||||
(parity == 1 && fmpz_sgn(n) < 0 && fmpz_is_odd(n)))
|
||||
{
|
||||
acb_indeterminate(res);
|
||||
/* Special case for zeta */
|
||||
if (q == 1 && fmpz_is_zero(n))
|
||||
acb_set_d(res, -0.5);
|
||||
else
|
||||
acb_zero(res);
|
||||
return;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -45,7 +45,11 @@ int main()
|
|||
|
||||
if (dirichlet_char_is_primitive(G, chi))
|
||||
{
|
||||
if (n_randint(state, 2))
|
||||
fmpq_set_si(x, -8 + (slong) n_randint(state, 8), 1);
|
||||
else
|
||||
fmpq_randtest(x, state, 2 + n_randint(state, 8));
|
||||
|
||||
acb_set_fmpq(s, x, prec);
|
||||
|
||||
if (n_randint(state, 2))
|
||||
|
|
|
|||
|
|
@ -177,6 +177,8 @@ void _arb_hypgeom_gamma_upper_fmpq_inf_bsplit(arb_t res, const fmpq_t a, const a
|
|||
slong _arb_hypgeom_gamma_lower_fmpq_0_choose_N(mag_t err, const fmpq_t a, const arb_t z, const mag_t abs_tol);
|
||||
void _arb_hypgeom_gamma_lower_fmpq_0_bsplit(arb_t res, const fmpq_t a, const arb_t z, slong N, slong prec);
|
||||
void _arb_gamma_upper_fmpq_step_bsplit(arb_t Gz1, const fmpq_t a, const arb_t z0, const arb_t z1, const arb_t Gz0, const arb_t expmz0, const mag_t abs_tol, slong prec);
|
||||
slong _arb_hypgeom_gamma_upper_singular_si_choose_N(mag_t err, slong n, const arb_t z, const mag_t abs_tol);
|
||||
void _arb_hypgeom_gamma_upper_singular_si_bsplit(arb_t res, slong n, const arb_t z, slong N, slong prec);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
|
|
|
|||
|
|
@ -351,3 +351,137 @@ _arb_hypgeom_gamma_lower_fmpq_0_bsplit(arb_t res, const fmpq_t a, const arb_t z,
|
|||
arb_clear(S);
|
||||
arb_clear(Q);
|
||||
}
|
||||
|
||||
/* bounded by 1/z^n * sum z^k / k! */
|
||||
slong
|
||||
_arb_hypgeom_gamma_upper_singular_si_choose_N(mag_t err, slong n, const arb_t z, const mag_t abs_tol)
|
||||
{
|
||||
slong k;
|
||||
mag_t t, u, zm;
|
||||
|
||||
mag_init(t);
|
||||
mag_init(u);
|
||||
mag_init(zm);
|
||||
|
||||
arb_get_mag(zm, z);
|
||||
|
||||
arb_get_mag_lower(t, z);
|
||||
mag_inv(t, t);
|
||||
mag_pow_ui(t, t, n);
|
||||
|
||||
for (k = 1; ; k++)
|
||||
{
|
||||
mag_mul(t, t, zm);
|
||||
mag_div_ui(t, t, k);
|
||||
|
||||
if (mag_cmp(t, abs_tol) < 0)
|
||||
{
|
||||
mag_div_ui(u, zm, k);
|
||||
mag_geom_series(u, u, 0);
|
||||
mag_mul(u, t, u);
|
||||
|
||||
if (mag_cmp(u, abs_tol) < 0)
|
||||
{
|
||||
mag_swap(err, t);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
mag_clear(t);
|
||||
mag_clear(u);
|
||||
mag_clear(zm);
|
||||
|
||||
return k;
|
||||
}
|
||||
|
||||
static void
|
||||
singular_bsplit(arb_t M, arb_t S, arb_t Q, slong n, const arb_t z, slong na, slong nb, int cont, slong prec)
|
||||
{
|
||||
if (nb - na == 0)
|
||||
{
|
||||
arb_zero(M);
|
||||
arb_zero(S);
|
||||
arb_one(Q);
|
||||
}
|
||||
else if (nb - na == 1)
|
||||
{
|
||||
fmpz_t t;
|
||||
slong k;
|
||||
|
||||
k = na;
|
||||
|
||||
if (k == n)
|
||||
arb_neg(M, z);
|
||||
else
|
||||
arb_mul_si(M, z, n - k, prec);
|
||||
|
||||
arb_set_si(S, (k != n) ? (k + 1) : 0);
|
||||
|
||||
fmpz_init_set_si(t, k + 1);
|
||||
if (k != n)
|
||||
fmpz_mul_si(t, t, k - n);
|
||||
arb_set_fmpz(Q, t);
|
||||
fmpz_clear(t);
|
||||
}
|
||||
else
|
||||
{
|
||||
slong m;
|
||||
arb_t M2, S2, Q2;
|
||||
|
||||
m = na + (nb - na) / 2;
|
||||
|
||||
arb_init(M2);
|
||||
arb_init(S2);
|
||||
arb_init(Q2);
|
||||
|
||||
singular_bsplit(M, S, Q, n, z, na, m, 1, prec);
|
||||
singular_bsplit(M2, S2, Q2, n, z, m, nb, cont, prec);
|
||||
|
||||
arb_mul(S, S, Q2, prec);
|
||||
arb_addmul(S, M, S2, prec);
|
||||
|
||||
if (cont)
|
||||
arb_mul(M, M, M2, prec);
|
||||
|
||||
arb_mul(Q, Q, Q2, prec);
|
||||
|
||||
arb_clear(M2);
|
||||
arb_clear(S2);
|
||||
arb_clear(Q2);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
_arb_hypgeom_gamma_upper_singular_si_bsplit(arb_t res, slong n, const arb_t z, slong N, slong prec)
|
||||
{
|
||||
arb_t M, S, Q;
|
||||
|
||||
arb_init(M);
|
||||
arb_init(S);
|
||||
arb_init(Q);
|
||||
|
||||
N = FLINT_MAX(N, 0);
|
||||
|
||||
singular_bsplit(M, S, Q, n, z, 0, N, 0, prec);
|
||||
|
||||
/* (-1)**n/fac(n) * (digamma(n+1) - ln(z)) - (S/Q)/z**n */
|
||||
arb_pow_ui(M, z, n, prec);
|
||||
arb_mul(Q, Q, M, prec);
|
||||
arb_div(S, S, Q, prec);
|
||||
|
||||
arb_set_ui(M, n + 1);
|
||||
arb_digamma(M, M, prec);
|
||||
arb_log(Q, z, prec);
|
||||
arb_sub(M, M, Q, prec);
|
||||
arb_fac_ui(Q, n, prec);
|
||||
arb_div(M, M, Q, prec);
|
||||
if (n & 1)
|
||||
arb_neg(M, M);
|
||||
|
||||
arb_sub(res, M, S, prec);
|
||||
|
||||
arb_clear(M);
|
||||
arb_clear(S);
|
||||
arb_clear(Q);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -141,6 +141,62 @@ int main()
|
|||
mag_clear(err2);
|
||||
}
|
||||
|
||||
for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
|
||||
{
|
||||
slong na;
|
||||
arb_t z, r1, r2;
|
||||
mag_t tol1, tol2, err1, err2;
|
||||
slong N1, N2, prec1, prec2;
|
||||
|
||||
mag_init(tol1);
|
||||
mag_init(tol2);
|
||||
mag_init(err1);
|
||||
mag_init(err2);
|
||||
|
||||
arb_init(z);
|
||||
arb_init(r1);
|
||||
arb_init(r2);
|
||||
|
||||
prec1 = 2 + n_randint(state, 100);
|
||||
prec2 = 2 + n_randint(state, 100);
|
||||
|
||||
na = n_randint(state, 50);
|
||||
|
||||
arb_set_ui(z, 1 + n_randint(state, 100));
|
||||
arb_div_ui(z, z, 1 + n_randint(state, 100), 200);
|
||||
|
||||
mag_set_ui_2exp_si(tol1, 1, -(slong) n_randint(state, 100));
|
||||
mag_set_ui_2exp_si(tol2, 1, -(slong) n_randint(state, 100));
|
||||
|
||||
N1 = _arb_hypgeom_gamma_upper_singular_si_choose_N(err1, na, z, tol1);
|
||||
N2 = _arb_hypgeom_gamma_upper_singular_si_choose_N(err2, na, z, tol2);
|
||||
|
||||
_arb_hypgeom_gamma_upper_singular_si_bsplit(r1, na, z, N1, prec1);
|
||||
arb_add_error_mag(r1, err1);
|
||||
|
||||
_arb_hypgeom_gamma_upper_singular_si_bsplit(r2, na, z, N2, prec2);
|
||||
arb_add_error_mag(r2, err2);
|
||||
|
||||
if (!arb_overlaps(r1, r2))
|
||||
{
|
||||
flint_printf("FAIL: overlap\n\n");
|
||||
flint_printf("na = %wd", na); flint_printf("\n\n");
|
||||
flint_printf("z = "); arb_printn(z, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("r1 = "); arb_printn(r1, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("r2 = "); arb_printn(r2, 100, 0); flint_printf("\n\n");
|
||||
flint_abort();
|
||||
}
|
||||
|
||||
arb_clear(z);
|
||||
arb_clear(r1);
|
||||
arb_clear(r2);
|
||||
|
||||
mag_clear(tol1);
|
||||
mag_clear(tol2);
|
||||
mag_clear(err1);
|
||||
mag_clear(err2);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
|
|
|
|||
|
|
@ -313,6 +313,18 @@ Internal evaluation functions
|
|||
the Taylor series at zero before term *N*.
|
||||
The truncation error bound has to be added separately.
|
||||
|
||||
.. function:: slong _arb_hypgeom_gamma_upper_singular_si_choose_N(mag_t err, slong n, const arb_t z, const mag_t abs_tol)
|
||||
|
||||
Returns number of terms *N* and sets *err* to the truncation error for evaluating
|
||||
`\Gamma(-n,z)` using the Taylor series at zero, targeting an absolute
|
||||
tolerance of *abs_tol*.
|
||||
|
||||
.. function:: void _arb_hypgeom_gamma_upper_singular_si_bsplit(arb_t res, slong n, const arb_t z, slong N, slong prec)
|
||||
|
||||
Sets *res* to the approximation of `\Gamma(-n,z)` obtained by truncating
|
||||
the Taylor series at zero before term *N*.
|
||||
The truncation error bound has to be added separately.
|
||||
|
||||
.. function:: void _arb_gamma_upper_fmpq_step_bsplit(arb_t Gz1, const fmpq_t a, const arb_t z0, const arb_t z1, const arb_t Gz0, const arb_t expmz0, const mag_t abs_tol, slong prec)
|
||||
|
||||
Given *Gz0* and *expmz0* representing the values `\Gamma(a,z_0)` and `\exp(-z_0)`,
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue