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

acb_dirichlet_l_fmpq

Browse source
This commit is contained in:
fredrik 2021-11-25 15:20:28 +01:00
parent e66ba199d6
commit 11de2dd2f6
11 changed files with 1684 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
acb_dirichlet.h
View file

@ -145,6 +145,9 @@ void acb_dirichlet_l_vec_hurwitz(acb_ptr res, const acb_t s, const acb_dirichlet
void acb_dirichlet_l_jet(acb_ptr res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, int deflate, slong len, slong prec);
void acb_dirichlet_l_fmpq_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
void acb_dirichlet_l_fmpq(acb_t res, const fmpq_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
void _acb_dirichlet_l_series(acb_ptr res, acb_srcptr s, slong slen,
const dirichlet_group_t G, const dirichlet_char_t chi,
int deflate, slong len, slong prec);

73
acb_dirichlet/l_fmpq.c Normal file
View file

@ -0,0 +1,73 @@
/*
Copyright (C) 2021 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 "arb_mat.h"
#include "arb_hypgeom.h"
#include "acb_dirichlet.h"
int
acb_dirichlet_l_fmpq_use_afe(ulong q, const fmpq_t s, slong prec)
{
double cutoff;
if ((slong) fmpz_bits(fmpq_numref(s)) - (slong) fmpz_bits(fmpq_denref(s)) > 20)
return 0;
if (fabs(fmpq_get_d(s)) > 10 + prec * 0.01)
return 0;
if (q == 1)
{
if (fmpz_is_one(fmpq_denref(s)))
return 0;
if (fmpz_is_one(fmpq_numref(s)) && fmpz_equal_si(fmpq_denref(s), 2))
return prec > 32000;
return prec > 70000;
}
if (fmpq_is_zero(s))
return 0;
if (fmpz_cmp_ui(fmpq_denref(s), 2) <= 0)
{
if (prec > 30000)
return 1;
if (fmpq_is_one(s))
return q > 1000;
return q > 50;
}
cutoff = 15000.0 / q;
return prec > cutoff;
}
void
acb_dirichlet_l_fmpq(acb_t res, const fmpq_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
{
ulong q;
q = (G == NULL) ? 1 : G->q;
if (acb_dirichlet_l_fmpq_use_afe(q, s, prec))
{
acb_dirichlet_l_fmpq_afe(res, s, G, chi, prec);
if (acb_is_finite(res))
return;
}
acb_set_fmpq(res, s, prec + 10);
acb_dirichlet_l(res, res, G, chi, prec);
}

575
acb_dirichlet/l_fmpq_afe.c Normal file
View file

@ -0,0 +1,575 @@
/*
Copyright (C) 2021 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 "arb_mat.h"
#include "arb_hypgeom.h"
#include "acb_dirichlet.h"
#define VERBOSE 0
static double
log_gamma_upper_approx(double a, double z)
{
if (a < z)
return (a - 1) * log(z) - z;
else
return a * (log(a) - 1);
}
void
acb_dirichlet_root_number2(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
{
acb_dirichlet_root_number(res, G, chi, prec);
if (dirichlet_char_is_real(G, chi))
arb_zero(acb_imagref(res));
}
static void
_arf_trunc(arf_t x)
{
if (arf_sgn(x) < 0)
arf_ceil(x, x);
else
arf_floor(x, x);
}
static void
arb_extract_bits(arb_t t, const arb_t z, slong b)
{
arb_mul_2exp_si(t, z, b);
_arf_trunc(arb_midref(t));
mag_zero(arb_radref(t));
arb_mul_2exp_si(t, t, -b);
}
static void
acb_dirichlet_afe_tail_bound(mag_t res, const fmpq_t sd2, slong N, ulong q, int parity)
{
mag_t pi_n2_q, t, u;
fmpz_t sprime;
mag_init(pi_n2_q);
mag_init(t);
mag_init(u);
fmpz_init(sprime);
/* pi_n2_q = pi * N^2 / q (lower bound) */
mag_const_pi_lower(pi_n2_q);
mag_mul_ui_lower(pi_n2_q, pi_n2_q, N);
mag_mul_ui_lower(pi_n2_q, pi_n2_q, N);
mag_set_ui(t, q);
mag_div_lower(pi_n2_q, pi_n2_q, t);
/* upper bound for sd2 */
fmpz_cdiv_q(sprime, fmpq_numref(sd2), fmpq_denref(sd2));
/* require pi_n2_q > s' */
mag_set_fmpz(t, sprime);
if (fmpz_sgn(sprime) > 0 && mag_cmp(pi_n2_q, t) <= 0)
{
mag_inf(res);
}
else
{
mag_expinv(res, pi_n2_q);
mag_div_ui(res, res, N);
if (!parity)
mag_div_ui(res, res, N);
/* (1 + q/pi) */
mag_set_ui(t, q);
mag_const_pi_lower(u);
mag_div(t, t, u);
mag_add_ui(t, t, 1);
mag_mul(res, res, t);
/* max(1, 2^s') */
if (fmpz_sgn(sprime) > 0)
mag_mul_2exp_fmpz(res, res, sprime);
/* (pi/q)^(s'-1) */
fmpz_sub_ui(sprime, sprime, 1);
if (fmpz_sgn(sprime) >= 0)
{
mag_const_pi(t);
mag_set_ui_lower(u, q);
mag_div(t, t, u);
mag_pow_fmpz(t, t, sprime);
}
else
{
mag_const_pi_lower(t);
mag_set_ui(u, q);
mag_div_lower(t, t, u);
fmpz_neg(sprime, sprime);
mag_pow_fmpz_lower(t, t, sprime);
mag_inv(t, t);
}
mag_mul(res, res, t);
}
mag_clear(pi_n2_q);
mag_clear(t);
mag_clear(u);
fmpz_clear(sprime);
}
void
acb_dirichlet_fmpq_sum_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G, const dirichlet_char_t chi, const mag_t abs_tol, slong prec)
{
slong NN, n, start_bits, bits, wp, wp2, gamma_cached_prec;
mag_t AE, err, abs_tol_gamma;
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;
ulong q;
double abs_tol_mag;
double gammainc_mag, gamma_mag, ns_mag;
double aa, zz;
double cancellation;
#if VERBOSE
double t1, t2, t3;
#endif
mag_init(AE);
mag_init(err);
mag_init(abs_tol_gamma);
arb_init(ns);
arb_init(t);
arb_init(u);
arb_init(v);
arb_init(z);
arb_init(z0);
arb_init(z1);
arb_init(x);
arb_init(x2);
arb_init(Ga);
arb_init(Gz0);
arb_init(Gz1);
arb_init(expmz0);
arb_init(z0_prevn);
arb_init(Gz0_prevn);
arb_init(expmz0_prevn);
acb_init(c);
fmpq_init(s2);
if (G == NULL)
{
parity = 0;
q = 1;
}
else
{
parity = dirichlet_parity_char(G, chi);
q = G->q;
}
acb_zero(res);
/* Initial precision for gamma; may have to be increased later. */
gamma_cached_prec = prec * 1.05 + 30;
/* s2 = (s+parity)/2 */
fmpq_add_ui(s2, s, parity);
fmpq_div_2exp(s2, s2, 1);
arb_gamma_fmpq(Ga, s2, gamma_cached_prec);
for (n = 1; ; n += 1)
{
#if VERBOSE
printf("-----------------------------------------------------------\n");
flint_printf("n = %wd (s+parity)/2 = %f z = %f q = %wu\n", n, fmpq_get_d(s2), 3.1415926535897932 * n * n / q, q);
#endif
acb_dirichlet_afe_tail_bound(err, s2, n, q, parity);
#if VERBOSE
printf(" abs_tol = "); mag_printd(abs_tol, 5); printf("\n");
printf(" truncation error = "); mag_printd(err, 5); printf("\n");
#endif
if (mag_cmp(err, abs_tol) < 0)
{
if (G == NULL || dirichlet_char_is_real(G, chi))
arb_add_error_mag(acb_realref(res), err);
else
acb_add_error_mag(res, err);
break;
}
/* Compute local precision and tolerances. */
abs_tol_mag = mag_get_d_log2_approx(abs_tol);
aa = fmpq_get_d(s2);
zz = 3.1415926535897932385 * n * n / q;
/* 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 = gammainc_mag + ns_mag - abs_tol_mag + 5;
wp = FLINT_MAX(wp, 30);
/* Precision needed with cancellation. */
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
if (G == NULL)
acb_one(c);
else
acb_dirichlet_chi(c, G, chi, n, wp);
if (acb_is_zero(c))
continue;
arb_const_pi(z, wp2);
arb_mul_ui(z, z, n, wp2);
arb_mul_ui(z, z, n, wp2);
arb_div_ui(z, z, q, wp2);
start_bits = 32;
arb_extract_bits(z0, z, start_bits);
/* Can we use the asymptotic series? */
NN = _arb_hypgeom_gamma_upper_fmpq_inf_choose_N(AE, s2, z0, abs_tol_gamma);
#if VERBOSE
t1 = clock();
#endif
/* For each new point, evaluate from scratch or use continuation? The
former seems to be faster. */
if (1)
{
if (NN != -1)
{
/* OK to use the asymptotic series. */
_arb_hypgeom_gamma_upper_fmpq_inf_bsplit(Gz0, s2, z0, NN, wp);
arb_add_error_mag(Gz0, AE);
#if VERBOSE
flint_printf(" asymptotic series with N = %wd: ", NN); arb_printd(Gz0, 10); printf("\n");
#endif
}
else
{
/* Otherwise fallback to the series at 0. */
NN = _arb_hypgeom_gamma_lower_fmpq_0_choose_N(AE, s2, z0, abs_tol_gamma);
#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);
}
#if VERBOSE
printf(" Gz0 = "); arb_printd(Gz0, 10); printf("\n");
#endif
if (n == 1)
{
arb_neg(expmz0, z0);
arb_exp(expmz0, expmz0, wp);
}
else
{
arb_sub(t, z0_prevn, z0, wp);
arb_exp(t, t, wp);
arb_mul(expmz0, expmz0_prevn, t, wp);
}
arb_set(z0_prevn, z0);
arb_set(expmz0_prevn, expmz0);
arb_set(Gz0_prevn, Gz0);
#if VERBOSE
t2 = clock();
#endif
/* Bit-burst steps */
for (bits = start_bits * 2; bits < wp / 8; bits *= 2)
{
arb_extract_bits(z1, z, bits);
_arb_gamma_upper_fmpq_step_bsplit(Gz1, s2, z0, z1, Gz0, expmz0, abs_tol_gamma, wp);
arb_sub(t, z0, z1, wp);
arb_exp(t, t, wp);
arb_mul(expmz0, expmz0, t, wp);
arb_set(Gz0, Gz1);
arb_set(z0, z1);
}
/* Final step, including error bound */
_arb_gamma_upper_fmpq_step_bsplit(Gz1, s2, z0, z, Gz0, expmz0, abs_tol_gamma, wp);
arb_set(Gz0, Gz1);
#if VERBOSE
printf(" Gz0 = "); arb_printd(Gz0, 10); printf(" tol "); mag_printd(abs_tol_gamma, 5); printf("\n");
#endif
/* Multiply by prefactor n^-s */
arb_set_ui(ns, n);
arb_pow_fmpq(ns, ns, s, wp);
arb_div(Gz0, Gz0, ns, wp);
#if VERBOSE
printf(" 1/n^s = "); arb_printn(ns, 5, ARB_STR_NO_RADIUS); printf("\n");
printf(" Gz0 * pre = "); arb_printd(Gz0, 10); printf(" tol "); mag_printd(abs_tol, 5); printf("\n");
#endif
acb_addmul_arb(res, c, Gz0, prec);
#if VERBOSE
printf(" sum = "); acb_printd(res, 10); printf("\n");
t3 = clock();
printf(" time: %f, %f\n", (t2 - t1) / CLOCKS_PER_SEC, (t3 - t2) / CLOCKS_PER_SEC);
#endif
}
mag_clear(AE);
mag_clear(err);
mag_clear(abs_tol_gamma);
arb_clear(ns);
arb_clear(t);
arb_clear(u);
arb_clear(v);
arb_clear(z);
arb_clear(z0);
arb_clear(z1);
arb_clear(x);
arb_clear(x2);
arb_clear(Ga);
arb_clear(Gz0);
arb_clear(Gz1);
arb_clear(expmz0);
arb_clear(z0_prevn);
arb_clear(Gz0_prevn);
arb_clear(expmz0_prevn);
acb_clear(c);
fmpq_clear(s2);
}
#define PI 3.1415926535897932385
#define INV_LOG2 1.4426950408889634074;
/* max(pi/q,s/2)**(s/2-1) * exp(-max(pi/q,s/2))
static double
estimate_sum1_mag(double s, double q)
{
return ((0.5 * s - 1) * log(FLINT_MAX(PI / q, 0.5 * s)) - FLINT_MAX(PI / q, 0.5 * s)) * INV_LOG2;
}
*/
void
acb_dirichlet_l_fmpq_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
{
arb_t t;
acb_t S1, S2, w;
fmpq_t s2;
mag_t tol1, tol2;
double ds, m1, m2, m2pre;
slong origprec, prec1, prec2;
ulong q;
int parity;
/* Todo: implement decomposition for imprimitive characters. */
if (G != NULL && !dirichlet_char_is_primitive(G, chi))
{
acb_indeterminate(res);
return;
}
q = (G == NULL) ? 1 : G->q;
parity = (G == NULL) ? 0 : dirichlet_parity_char(G, chi);
/* Limit at gamma(-n) is not implemented. */
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)))
{
acb_indeterminate(res);
return;
}
}
origprec = prec;
prec = prec * 1.001 + 2 * FLINT_BIT_COUNT(q);
acb_init(S1);
acb_init(S2);
acb_init(w);
arb_init(t);
fmpq_init(s2);
mag_init(tol1);
mag_init(tol2);
ds = fmpq_get_d(s);
m1 = log_gamma_upper_approx(0.5 * (ds + parity), PI / q) * INV_LOG2;
m2 = log_gamma_upper_approx(0.5 * (1.0 - ds + parity), PI / q) * INV_LOG2;
m2pre = (ds - 0.5) * log(PI / q) * INV_LOG2;
mag_one(tol1);
mag_mul_2exp_si(tol1, tol1, FLINT_MAX(m1, m2 + m2pre) - prec);
mag_mul_2exp_si(tol2, tol1, -m2pre);
prec1 = prec - (FLINT_MAX(m1, m2 + m2pre) - m1);
prec1 = FLINT_MAX(prec1, 32);
prec2 = prec - (FLINT_MAX(m1, m2 + m2pre) - (m2 + m2pre));
prec2 = FLINT_MAX(prec2, 32);
#if VERBOSE
printf("mag1 = %ld mag2 = %ld mag2 + pre = %ld prec, prec1, prec2 = %ld, %ld, %ld\n",
(slong) m1, (slong) m2, (slong) (m2 + m2pre), prec, prec1, prec2);
printf("tol1 = %ld tol2 = %ld\n", MAG_EXP(tol1), MAG_EXP(tol2));
#endif
acb_dirichlet_fmpq_sum_afe(S1, s, G, chi, tol1, prec1);
#if VERBOSE
printf("=====================================================\n");
printf("S1 = "); acb_printd(S1, 20); printf(" estimate = "); printf(" %g\n", ldexp(1.0, m1));
printf("=====================================================\n");
#endif
if (q == 1 && fmpz_is_one(fmpq_numref(s)) && fmpz_equal_ui(fmpq_denref(s), 2))
{
acb_mul_2exp_si(res, S1, 1);
}
else
{
/* rootnum (pi/q)^(s-1/2) sum(1-s) */
if (fmpz_is_one(fmpq_numref(s)) && fmpz_equal_ui(fmpq_denref(s), 2))
{
acb_conj(S2, S1);
}
else
{
fmpq_sub_ui(s2, s, 1);
fmpq_neg(s2, s2);
acb_dirichlet_fmpq_sum_afe(S2, s2, G, chi, tol2, prec2);
acb_conj(S2, S2);
}
#if VERBOSE
printf("=====================================================\n");
printf("S1 = "); acb_printd(S1, 20); printf(" estimate = "); printf(" %g\n", ldexp(1.0, m1));
printf("S2 = "); acb_printd(S2, 20); printf(" estimate = "); printf(" %g\n", ldexp(1.0, m2));
printf("=====================================================\n");
#endif
arb_const_pi(t, prec);
arb_div_ui(t, t, q, prec);
fmpq_set_si(s2, 1, 2);
fmpq_sub(s2, s, s2);
arb_pow_fmpq(t, t, s2, prec);
acb_mul_arb(S2, S2, t, prec);
if (q != 1)
{
acb_dirichlet_root_number2(w, G, chi, prec);
acb_mul(S2, S2, w, prec);
}
#if VERBOSE
printf("S2 * prefactor = "); acb_printd(S2, 20); printf(" estimate = "); printf(" %g\n", ldexp(1.0, m2 + m2pre));
#endif
acb_add(res, S1, S2, prec);
}
/* add pi^(s/2) / (s (s-1)) */
if (q == 1)
{
arb_const_pi(t, prec);
fmpq_div_2exp(s2, s, 1);
arb_pow_fmpq(t, t, s2, prec);
fmpq_sub_ui(s2, s, 1);
fmpq_mul(s2, s2, s);
arb_div_fmpz(t, t, fmpq_numref(s2), prec);
arb_mul_fmpz(t, t, fmpq_denref(s2), prec);
acb_add_arb(res, res, t, prec);
}
/* divide by gamma((s+parity)/2) */
fmpq_add_ui(s2, s, parity);
fmpq_div_2exp(s2, s2, 1);
arb_gamma_fmpq(t, s2, prec);
acb_div_arb(res, res, t, prec);
acb_set_round(res, res, origprec);
acb_clear(S1);
acb_clear(S2);
acb_clear(w);
arb_clear(t);
fmpq_clear(s2);
mag_clear(tol1);
mag_clear(tol2);
}

83
acb_dirichlet/test/t-l_fmpq.c Normal file
View 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_dirichlet.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("l_fmpq....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 500 * arb_test_multiplier(); iter++)
{
fmpq_t x;
acb_t s, t, u;
dirichlet_group_t G;
dirichlet_char_t chi;
ulong q, k;
slong prec;
fmpq_init(x);
acb_init(s);
acb_init(t);
acb_init(u);
q = 1 + n_randint(state, 50);
prec = 2 + n_randint(state, 300);
k = n_randint(state, n_euler_phi(q));
dirichlet_group_init(G, q);
dirichlet_char_init(chi, G);
dirichlet_char_index(chi, G, k);
fmpq_randtest(x, state, 2 + n_randint(state, 100));
if (fmpq_get_d(x) < -100.0)
fmpq_neg(x, x);
acb_set_fmpq(s, x, prec);
if (n_randint(state, 2))
acb_dirichlet_l_hurwitz(t, s, NULL, G, chi, prec);
else
acb_dirichlet_l(t, s, G, chi, prec);
acb_dirichlet_l_fmpq(u, x, G, chi, prec);
if (!acb_overlaps(t, u))
{
flint_printf("FAIL: overlap\n\n");
flint_printf("iter = %wd q = %wu k = %wu prec = %wd\n\n", iter, q, k, prec);
flint_printf("x = "); fmpq_print(x); flint_printf("\n\n");
flint_printf("s = "); acb_printn(s, 100, 0); flint_printf("\n\n");
flint_printf("t = "); acb_printn(t, 100, 0); flint_printf("\n\n");
flint_printf("u = "); acb_printn(u, 100, 0); flint_printf("\n\n");
acb_sub(t, t, u, prec);
flint_printf("t - u = "); acb_printd(t, 30); flint_printf("\n\n");
flint_abort();
}
dirichlet_char_clear(chi);
dirichlet_group_clear(G);
fmpq_clear(x);
acb_clear(s);
acb_clear(t);
acb_clear(u);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

84
acb_dirichlet/test/t-l_fmpq_afe.c Normal file
View file

@ -0,0 +1,84 @@
/*
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_dirichlet.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("l_fmpq_afe....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
{
fmpq_t x;
acb_t s, t, u;
dirichlet_group_t G;
dirichlet_char_t chi;
ulong q, k;
slong prec;
fmpq_init(x);
acb_init(s);
acb_init(t);
acb_init(u);
q = 1 + n_randint(state, 50);
prec = 2 + n_randint(state, 400);
k = n_randint(state, n_euler_phi(q));
dirichlet_group_init(G, q);
dirichlet_char_init(chi, G);
dirichlet_char_index(chi, G, k);
if (dirichlet_char_is_primitive(G, chi))
{
fmpq_randtest(x, state, 2 + n_randint(state, 8));
acb_set_fmpq(s, x, prec);
if (n_randint(state, 2))
acb_dirichlet_l_hurwitz(t, s, NULL, G, chi, prec);
else
acb_dirichlet_l(t, s, G, chi, prec);
acb_dirichlet_l_fmpq_afe(u, x, G, chi, prec);
if (!acb_overlaps(t, u))
{
flint_printf("FAIL: overlap\n\n");
flint_printf("iter = %wd q = %wu k = %wu prec = %wd\n\n", iter, q, k, prec);
flint_printf("x = "); fmpq_print(x); flint_printf("\n\n");
flint_printf("s = "); acb_printn(s, 100, 0); flint_printf("\n\n");
flint_printf("t = "); acb_printn(t, 100, 0); flint_printf("\n\n");
flint_printf("u = "); acb_printn(u, 100, 0); flint_printf("\n\n");
acb_sub(t, t, u, prec);
flint_printf("t - u = "); acb_printd(t, 30); flint_printf("\n\n");
flint_abort();
}
}
dirichlet_char_clear(chi);
dirichlet_group_clear(G);
fmpq_clear(x);
acb_clear(s);
acb_clear(t);
acb_clear(u);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

6
arb_hypgeom.h
View file

@ -172,6 +172,12 @@ void arb_hypgeom_central_bin_ui(arb_t res, ulong n, slong prec);
void arb_hypgeom_dilog(arb_t res, const arb_t z, slong prec);
slong _arb_hypgeom_gamma_upper_fmpq_inf_choose_N(mag_t err, const fmpq_t a, const arb_t z, const mag_t abs_tol);
void _arb_hypgeom_gamma_upper_fmpq_inf_bsplit(arb_t res, const fmpq_t a, const arb_t z, slong N, slong prec);
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);
#ifdef __cplusplus
}
#endif

353
arb_hypgeom/gamma_upper_fmpq.c Normal file
View file

@ -0,0 +1,353 @@
/*
Copyright (C) 2021 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 "arb_hypgeom.h"
static void
mag_div_fmpq(mag_t res, const mag_t x, const fmpq_t a)
{
mag_t t;
mag_init(t);
mag_set_fmpz_lower(t, fmpq_numref(a));
mag_div(res, x, t);
mag_set_fmpz(t, fmpq_denref(a));
mag_mul(res, res, t);
mag_clear(t);
}
static void
mag_pow_fmpq_fast(mag_t res, const mag_t x, const fmpq_t e)
{
fmpz_t b;
fmpz_init(b);
if (mag_cmp_2exp_si(x, 0) >= 0)
{
fmpz_cdiv_q(b, fmpq_numref(e), fmpq_denref(e));
mag_pow_fmpz(res, x, b);
}
else
{
fmpz_fdiv_q(b, fmpq_numref(e), fmpq_denref(e));
mag_pow_fmpz(res, x, b);
}
fmpz_clear(b);
}
slong
_arb_hypgeom_gamma_upper_fmpq_inf_choose_N(mag_t err, const fmpq_t a, const arb_t z, const mag_t abs_tol)
{
slong N, aa, ab;
fmpz_t az1, az2;
fmpq_t a1;
mag_t t, u;
fmpz_init(az1);
fmpz_init(az2);
fmpq_init(a1);
mag_init(t);
mag_init(u);
fmpz_fdiv_q(az1, fmpq_numref(a), fmpq_denref(a));
fmpz_cdiv_q(az2, fmpq_numref(a), fmpq_denref(a));
if (!fmpz_fits_si(az1) || !fmpz_fits_si(az2))
{
mag_inf(err);
N = -1;
}
else
{
aa = fmpz_get_si(az1);
ab = fmpz_get_si(az2);
/* prefactor z^(a-1) * exp(-z) */
arb_get_mag_lower(t, z);
mag_expinv(t, t);
fmpq_sub_ui(a1, a, 1);
arb_get_mag(u, z);
mag_pow_fmpq_fast(u, u, a1);
mag_mul(err, t, u);
arb_get_mag_lower(t, z);
mag_inv(t, t);
for (N = 1; ; N++)
{
mag_mul_ui(u, err, FLINT_MAX(FLINT_ABS(aa - N), FLINT_ABS(ab - N)));
mag_mul(u, u, t);
if (N >= ab - 1 && mag_cmp(u, abs_tol) < 0)
{
mag_swap(err, u);
break;
}
/* Stop if terms are increasing, unless a is a positive integer in
which case the series will terminate eventually. */
if (mag_cmp(u, err) > 0 && !(aa == ab && aa >= 1))
{
mag_inf(err);
N = -1;
break;
}
mag_swap(err, u);
}
}
fmpz_clear(az1);
fmpz_clear(az2);
mag_clear(t);
mag_clear(u);
fmpq_clear(a1);
return N;
}
static void
upper_bsplit(arb_t M, arb_t S, arb_t Q, const fmpz_t ap, const fmpz_t aq, 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;
fmpz_init_set(t, ap);
fmpz_submul_ui(t, aq, na + 1);
fmpz_neg(t, t);
arb_set_fmpz(M, t);
arb_mul_fmpz(S, z, aq, prec);
arb_neg(S, S);
arb_set(Q, S);
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);
upper_bsplit(M, S, Q, ap, aq, z, na, m, 1, prec);
upper_bsplit(M2, S2, Q2, ap, aq, z, m, nb, cont, prec);
/* todo: squaring opt; power table */
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_fmpq_inf_bsplit(arb_t res, const fmpq_t a, const arb_t z, slong N, slong prec)
{
arb_t M, S, Q;
fmpq_t a1;
arb_init(M);
arb_init(S);
arb_init(Q);
fmpq_init(a1);
N = FLINT_MAX(N, 0);
upper_bsplit(M, S, Q, fmpq_numref(a), fmpq_denref(a), z, 0, N, 0, prec);
arb_div(S, S, Q, prec);
fmpq_sub_ui(a1, a, 1);
arb_pow_fmpq(M, z, a1, prec);
arb_mul(S, S, M, prec);
arb_neg(M, z);
arb_exp(M, M, prec);
arb_mul(res, S, M, prec);
arb_clear(M);
arb_clear(S);
arb_clear(Q);
fmpq_clear(a1);
}
/* select N and bound error for z**a * exp(-z) / a * sum(z**n / rf(a+1, n)) */
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)
{
slong N, aa, ab, c;
fmpz_t az1, az2;
mag_t t, u;
fmpz_init(az1);
fmpz_init(az2);
mag_init(t);
mag_init(u);
fmpz_fdiv_q(az1, fmpq_numref(a), fmpq_denref(a));
fmpz_cdiv_q(az2, fmpq_numref(a), fmpq_denref(a));
if (!fmpz_fits_si(az1) || !fmpz_fits_si(az2))
{
mag_inf(err);
N = -1;
}
else
{
aa = fmpz_get_si(az1);
ab = fmpz_get_si(az2);
/* prefactor z^a * exp(-z) / a */
arb_get_mag_lower(t, z);
mag_expinv(t, t);
arb_get_mag(u, z);
mag_pow_fmpq_fast(u, u, a);
mag_mul(err, t, u);
mag_div_fmpq(err, err, a);
arb_get_mag(t, z);
for (N = 1; ; N++)
{
c = FLINT_MIN(FLINT_ABS(aa + N), FLINT_ABS(ab + N));
if (c == 0)
{
fmpq_t q;
fmpq_init(q);
fmpq_add_ui(q, a, N);
mag_div_fmpq(err, err, q);
fmpq_clear(q);
}
else
{
mag_div_ui(err, err, c);
}
mag_mul(err, err, t);
/* todo: condition can be relaxed */
/* todo: faster check (compare t) */
if ((aa + N) > 0 && mag_cmp(err, abs_tol) < 0)
{
mag_div_ui(u, t, aa + N);
mag_geom_series(u, u, 0);
mag_mul(u, err, u);
if (mag_cmp(u, abs_tol) < 0)
{
mag_swap(err, u);
break;
}
}
}
}
fmpz_clear(az1);
fmpz_clear(az2);
mag_clear(t);
mag_clear(u);
return N;
}
/* todo: squaring opt; power table */
static void
lower_bsplit(arb_t M, arb_t S, arb_t Q, const fmpz_t ap, const fmpz_t aq, 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;
fmpz_init_set(t, ap);
fmpz_addmul_ui(t, aq, na + 1);
arb_set_fmpz(S, t);
arb_set(Q, S);
arb_mul_fmpz(M, z, aq, prec);
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);
lower_bsplit(M, S, Q, ap, aq, z, na, m, 1, prec);
lower_bsplit(M2, S2, Q2, ap, aq, 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_lower_fmpq_0_bsplit(arb_t res, const fmpq_t a, 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);
lower_bsplit(M, S, Q, fmpq_numref(a), fmpq_denref(a), z, 0, N, 0, prec);
arb_div(S, S, Q, prec);
arb_pow_fmpq(M, z, a, prec);
arb_mul(S, S, M, prec);
arb_neg(M, z);
arb_exp(M, M, prec);
arb_mul(S, S, M, prec);
arb_div_fmpz(S, S, fmpq_numref(a), prec);
arb_mul_fmpz(res, S, fmpq_denref(a), prec);
arb_clear(M);
arb_clear(S);
arb_clear(Q);
}

316
arb_hypgeom/gamma_upper_fmpq_step_bsplit.c Normal file
View file

@ -0,0 +1,316 @@
/*
Copyright (C) 2021 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 "arb_mat.h"
#include "arb_hypgeom.h"
static void
taylor_M(mag_t M, const arb_t a, const arb_t z, const mag_t x, slong Rexp)
{
arb_t t, u;
arb_init(t);
arb_init(u);
arb_sub_ui(u, a, 1, 53);
arb_sgn(t, u);
arb_mul_2exp_si(t, t, Rexp);
arb_add(t, z, t, 53);
arb_pow(t, t, u, 53);
arb_one(u);
arb_mul_2exp_si(u, u, Rexp);
arb_sub(u, u, z, 53);
arb_exp(u, u, 53);
arb_mul(t, t, u, 53);
arb_get_mag(M, t);
arb_clear(t);
arb_clear(u);
}
/* choose N such that M * C^N / (1 - C) <= tol */
/* todo: fix */
static slong
mag_geom_choose_N(const mag_t M, const mag_t C, const mag_t tol)
{
mag_t t, u;
slong N;
/* N = log(M / ((1 - C) tol)) / log(1/C) */
mag_init(t);
mag_init(u);
mag_one(t);
mag_sub_lower(t, t, C);
mag_mul_lower(t, t, tol);
mag_div(t, M, t);
mag_log(t, t);
mag_inv_lower(u, C);
mag_log_lower(u, u);
mag_div(t, t, u);
N = ceil(mag_get_d(t));
N = FLINT_MAX(N, 1);
mag_clear(t);
mag_clear(u);
return N;
}
static void
taylor_bound(mag_t err, const arb_t a, const arb_t z, const mag_t x, slong Rexp, slong N)
{
mag_t C, M;
mag_init(C);
mag_init(M);
/* C = x / R */
mag_mul_2exp_si(C, x, -Rexp);
/* M R C^n / (1 - C) / N */
mag_geom_series(err, C, N);
if (!mag_is_inf(err))
{
taylor_M(M, a, z, x, Rexp);
mag_mul(err, err, M);
mag_mul_2exp_si(err, err, Rexp);
mag_div_ui(err, err, N);
}
mag_clear(C);
mag_clear(M);
}
static slong
taylor_N(const arb_t a, const arb_t z, const mag_t x, slong Rexp, const mag_t abs_tol)
{
mag_t C, M;
slong N;
mag_init(C);
mag_init(M);
/* C = x / R */
mag_mul_2exp_si(C, x, -Rexp);
if (mag_cmp_2exp_si(C, 0) < 0)
{
taylor_M(M, a, z, x, Rexp);
mag_mul_2exp_si(M, M, Rexp);
N = mag_geom_choose_N(M, C, abs_tol);
}
else
{
N = WORD_MAX;
}
mag_clear(C);
mag_clear(M);
return N;
}
static void
arb_hypgeom_gamma_upper_taylor_choose(slong * res_N, mag_t err, const arb_t a, const arb_t z, const mag_t x, const mag_t abs_tol)
{
slong N, New;
mag_t zlow;
slong Rexp;
mag_init(zlow);
arb_get_mag_lower(zlow, z);
Rexp = 0;
while (mag_cmp_2exp_si(zlow, Rexp + 1) < 0)
Rexp--;
N = taylor_N(a, z, x, Rexp, abs_tol);
while (N > 1 && mag_cmp_2exp_si(x, Rexp - 1) < 0)
{
New = taylor_N(a, z, x, Rexp - 1, abs_tol);
if (New <= N)
{
Rexp = Rexp - 1;
N = New;
}
else
{
break;
}
}
if (Rexp == 0)
{
while (N > 1 && mag_cmp_2exp_si(zlow, Rexp + 1) > 0)
{
New = taylor_N(a, z, x, Rexp + 1, abs_tol);
if (New <= N)
{
Rexp = Rexp + 1;
N = New;
}
else
{
break;
}
}
}
*res_N = N;
taylor_bound(err, a, z, x, Rexp, N);
if (mag_cmp(err, abs_tol) > 0)
{
printf("err = "); mag_printd(err, 10); printf("\n");
printf("abs_tol = "); mag_printd(abs_tol, 10); printf("\n");
printf("a = "); arb_printd(a, 10); printf("\n");
printf("z = "); arb_printd(z, 10); printf("\n");
printf("x = "); mag_printd(x, 10); printf("\n");
flint_abort();
}
mag_clear(zlow);
}
static void
gamma_upper_taylor_bsplit(arb_mat_t M, arb_t Q,
const fmpz_t ap, const fmpz_t aq, const arb_t z0, const arb_t x, const arb_t x2, slong a, slong b, int cont, slong prec)
{
if (b - a == 0)
{
arb_mat_one(M);
}
else if (b - a == 1)
{
slong n;
fmpz_t t;
n = a;
fmpz_init(t);
/* Q = -z0*(n+1)*(n+2)*aq */
fmpz_mul2_uiui(t, aq, n + 1, n + 2);
arb_mul_fmpz(Q, z0, t, prec);
arb_neg(Q, Q);
/* x Q */
arb_mul(arb_mat_entry(M, 0, 1), Q, x, prec);
/* aq n x */
fmpz_mul_ui(t, aq, n);
arb_mul_fmpz(arb_mat_entry(M, 1, 0), x, t, prec);
/* x*(-ap + aq*(n + z0 + 1))*(n + 1) */
arb_add_ui(arb_mat_entry(M, 1, 1), z0, n + 1, prec);
arb_mul_fmpz(arb_mat_entry(M, 1, 1), arb_mat_entry(M, 1, 1), aq, prec);
arb_sub_fmpz(arb_mat_entry(M, 1, 1), arb_mat_entry(M, 1, 1), ap, prec);
arb_mul_ui(arb_mat_entry(M, 1, 1), arb_mat_entry(M, 1, 1), n + 1, prec);
arb_mul(arb_mat_entry(M, 1, 1), arb_mat_entry(M, 1, 1), x, prec);
arb_set(arb_mat_entry(M, 2, 0), Q);
arb_set(arb_mat_entry(M, 2, 2), Q);
fmpz_clear(t);
}
else
{
arb_mat_t M1, M2;
arb_t Q2;
slong m;
arb_mat_init(M1, 3, 3);
arb_mat_init(M2, 3, 3);
arb_init(Q2);
m = a + (b - a) / 2;
gamma_upper_taylor_bsplit(M1, Q, ap, aq, z0, x, x2, a, m, 1, prec);
gamma_upper_taylor_bsplit(M2, Q2, ap, aq, z0, x, x2, m, b, cont, prec);
if (cont)
{
arb_mat_mul(M, M2, M1, prec);
}
else
{
arb_mat_transpose(M1, M1);
arb_dot(arb_mat_entry(M, 2, 0), NULL, 0, arb_mat_entry(M1, 0, 0), 1, arb_mat_entry(M2, 2, 0), 1, 3, prec);
arb_dot(arb_mat_entry(M, 2, 1), NULL, 0, arb_mat_entry(M1, 1, 0), 1, arb_mat_entry(M2, 2, 0), 1, 3, prec);
arb_dot(arb_mat_entry(M, 2, 2), NULL, 0, arb_mat_entry(M1, 2, 0), 1, arb_mat_entry(M2, 2, 0), 1, 3, prec);
}
arb_mul(Q, Q2, Q, prec);
arb_mat_clear(M1);
arb_mat_clear(M2);
arb_clear(Q2);
}
}
/*
Given Gz0 = Gamma(a, z0) and expmz0 = exp(-z0), compute Gz1 = Gamma(a, z1)
*/
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)
{
arb_t x, Q, a_real;
arb_mat_t M;
mag_t xmag, err;
slong N;
fmpq_t a1;
mag_init(xmag);
mag_init(err);
arb_init(x);
arb_init(Q);
arb_init(a_real);
fmpq_init(a1);
arb_mat_init(M, 3, 3);
arb_sub(x, z1, z0, prec);
arb_get_mag(xmag, x);
arb_set_fmpq(a_real, a, 53);
arb_hypgeom_gamma_upper_taylor_choose(&N, err, a_real, z0, xmag, abs_tol);
gamma_upper_taylor_bsplit(M, Q, fmpq_numref(a), fmpq_denref(a), z0, x, NULL, 0, N, 0, prec);
arb_mul(arb_mat_entry(M, 2, 0), arb_mat_entry(M, 2, 0), Gz0, prec);
fmpq_sub_ui(a1, a, 1);
arb_pow_fmpq(arb_mat_entry(M, 0, 0), z0, a1, prec);
arb_mul(arb_mat_entry(M, 0, 0), arb_mat_entry(M, 0, 0), expmz0, prec);
arb_submul(arb_mat_entry(M, 2, 0), arb_mat_entry(M, 2, 1), arb_mat_entry(M, 0, 0), prec);
arb_div(Gz1, arb_mat_entry(M, 2, 0), Q, prec);
arb_add_error_mag(Gz1, err);
mag_clear(xmag);
mag_clear(err);
arb_clear(x);
arb_clear(Q);
arb_clear(a_real);
fmpq_clear(a1);
arb_mat_clear(M);
}

148
arb_hypgeom/test/t-gamma_upper_fmpq.c Normal file
View file

@ -0,0 +1,148 @@
/*
Copyright (C) 2021 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 "arb_hypgeom.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("gamma_upper_fmpq....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
{
fmpq_t a;
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);
fmpq_init(a);
arb_init(z);
arb_init(r1);
arb_init(r2);
prec1 = 2 + n_randint(state, 100);
prec2 = 2 + n_randint(state, 100);
do {
fmpq_randtest(a, state, 5);
} while (fmpz_is_one(fmpq_denref(a)) && fmpz_sgn(fmpq_numref(a)) <= 0);
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_fmpq_inf_choose_N(err1, a, z, tol1);
N2 = _arb_hypgeom_gamma_upper_fmpq_inf_choose_N(err2, a, z, tol2);
_arb_hypgeom_gamma_upper_fmpq_inf_bsplit(r1, a, z, N1, prec1);
arb_add_error_mag(r1, err1);
_arb_hypgeom_gamma_upper_fmpq_inf_bsplit(r2, a, z, N2, prec2);
arb_add_error_mag(r2, err2);
if (!arb_overlaps(r1, r2))
{
flint_printf("FAIL: overlap\n\n");
flint_printf("a = "); fmpq_print(a); 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();
}
fmpq_clear(a);
arb_clear(z);
arb_clear(r1);
arb_clear(r2);
mag_clear(tol1);
mag_clear(tol2);
mag_clear(err1);
mag_clear(err2);
}
for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
{
fmpq_t a;
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);
fmpq_init(a);
arb_init(z);
arb_init(r1);
arb_init(r2);
prec1 = 2 + n_randint(state, 100);
prec2 = 2 + n_randint(state, 100);
do {
fmpq_randtest(a, state, 5);
} while (fmpz_is_one(fmpq_denref(a)) && fmpz_sgn(fmpq_numref(a)) <= 0);
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_lower_fmpq_0_choose_N(err1, a, z, tol1);
N2 = _arb_hypgeom_gamma_lower_fmpq_0_choose_N(err2, a, z, tol2);
_arb_hypgeom_gamma_lower_fmpq_0_bsplit(r1, a, z, N1, prec1);
arb_add_error_mag(r1, err1);
_arb_hypgeom_gamma_lower_fmpq_0_bsplit(r2, a, z, N2, prec2);
arb_add_error_mag(r2, err2);
if (!arb_overlaps(r1, r2))
{
flint_printf("FAIL: overlap\n\n");
flint_printf("a = "); fmpq_print(a); 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();
}
fmpq_clear(a);
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");
return EXIT_SUCCESS;
}

7
doc/source/acb_dirichlet.rst
View file

@ -541,6 +541,13 @@ Dirichlet L-functions
Computes `L(s,\chi)` using a default choice of algorithm.
.. function:: void acb_dirichlet_l_fmpq(acb_t res, const fmpq_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
void acb_dirichlet_l_fmpq_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
Computes `L(s,\chi)` where *s* is a rational number.
The *afe* version uses the approximate functional equation;
the default version chooses an algorithm automatically.
.. function:: void acb_dirichlet_l_vec_hurwitz(acb_ptr res, const acb_t s, const acb_dirichlet_hurwitz_precomp_t precomp, const dirichlet_group_t G, slong prec)
Compute all values `L(s,\chi)` for `\chi` mod `q`, using the

36
doc/source/arb_hypgeom.rst
View file

@ -285,6 +285,42 @@ Incomplete gamma and beta functions
The underscore method requires positive lengths and does not support
aliasing.
Internal evaluation functions
................................................................................
.. function:: slong _arb_hypgeom_gamma_upper_fmpq_inf_choose_N(mag_t err, const fmpq_t a, const arb_t z, const mag_t abs_tol)
Returns number of terms *N* and sets *err* to the truncation error for evaluating
`\Gamma(a,z)` using the asymptotic series at infinity, targeting an absolute
tolerance of *abs_tol*. The error may be set to *err* if the tolerance
cannot be achieved. Assumes that *z* is positive.
.. function:: void _arb_hypgeom_gamma_upper_fmpq_inf_bsplit(arb_t res, const fmpq_t a, const arb_t z, slong N, slong prec)
Sets *res* to the approximation of `\Gamma(a,z)` obtained by truncating
the asymptotic series at infinity before term *N*.
The truncation error bound has to be added separately.
.. function:: 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)
Returns number of terms *N* and sets *err* to the truncation error for evaluating
`\gamma(a,z)` using the Taylor series at zero, targeting an absolute
tolerance of *abs_tol*. Assumes that *z* is positive.
.. function:: void _arb_hypgeom_gamma_lower_fmpq_0_bsplit(arb_t res, const fmpq_t a, const arb_t z, slong N, slong prec)
Sets *res* to the approximation of `\gamma(a,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)`,
computes `\Gamma(a,z_1)` using the Taylor series at `z_0` evaluated
using binary splitting,
targeting an absolute error of *abs_tol*.
Assumes that `z_0` and `z_1` are positive.
Exponential and trigonometric integrals
-------------------------------------------------------------------------------