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

use the fast zeta multi-evaluation algorithm

Browse source
This commit is contained in:
Fredrik Johansson 2012-09-10 10:35:11 +02:00
parent eda22bb13f
commit da38e48fe5
4 changed files with 69 additions and 17 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
fmprb.h
View file

@ -230,6 +230,10 @@ void fmprb_zeta_ui_bernoulli(fmprb_t x, ulong n, long prec);
void fmprb_zeta_ui_vec_borwein(fmprb_struct * z, ulong start, long num, ulong step, long prec);
void fmprb_zeta_ui(fmprb_t x, ulong n, long prec);
void fmprb_zeta_ui_vec_even(fmprb_struct * x, ulong start, long num, long prec);
void fmprb_zeta_ui_vec_odd(fmprb_struct * x, ulong start, long num, long prec);
void fmprb_zeta_ui_vec(fmprb_struct * x, ulong start, long num, long prec);
static __inline__ void
fmprb_add_error_2exp_si(fmprb_t x, long err)
{

57
fmprb/zeta_ui.c
View file

@ -40,7 +40,7 @@ fmprb_zeta_ui(fmprb_t x, ulong n, long prec)
printf("exception: zeta_ui(1)\n");
abort();
}
/* asymptotic case */
/* fast detection of asymptotic case */
else if (n > 6 && n > 0.7 * prec)
{
fmprb_zeta_ui_euler_product(x, n, prec);
@ -70,3 +70,58 @@ fmprb_zeta_ui(fmprb_t x, ulong n, long prec)
fmprb_zeta_ui_vec_borwein(x, n, 1, 0, prec);
}
}
void
fmprb_zeta_ui_vec_even(fmprb_struct * x, ulong start, long num, long prec)
{
long i;
for (i = 0; i < num; i++)
fmprb_zeta_ui(x + i, start + 2 * i, prec);
}
void
fmprb_zeta_ui_vec_odd(fmprb_struct * x, ulong start, long num, long prec)
{
long i, num_borwein;
ulong end, cutoff;
end = start + num * 2;
cutoff = 40 + 0.3 * prec;
if (cutoff > start)
{
num_borwein = 1 + (cutoff - start) / 2;
num_borwein = FLINT_MIN(num_borwein, num);
}
else
num_borwein = 0;
fmprb_zeta_ui_vec_borwein(x, start, num_borwein, 2, prec);
for (i = num_borwein; i < num; i++)
fmprb_zeta_ui(x + i, start + 2 * i, prec);
}
void
fmprb_zeta_ui_vec(fmprb_struct * x, ulong start, long num, long prec)
{
long i, num_odd, num_even, start_odd;
fmprb_struct * tmp;
num_odd = num / 2 + (start & num & 1);
num_even = num - num_odd;
start_odd = start % 2;
fmprb_zeta_ui_vec_even(x, start + start_odd, num_even, prec);
fmprb_zeta_ui_vec_odd(x + num_even, start + !start_odd, num_even, prec);
/* interleave */
tmp = flint_malloc(sizeof(fmprb_struct) * num);
for (i = 0; i < num_even; i++) tmp[i] = x[i];
for (i = 0; i < num_odd; i++) tmp[num_even + i] = x[num_even + i];
for (i = 0; i < num_even; i++) x[start_odd + 2 * i] = tmp[i];
for (i = 0; i < num_odd; i++) x[!start_odd + 2 * i] = tmp[num_even + i];
flint_free(tmp);
}

3
fmprb/zeta_ui_vec_borwein.c
View file

@ -74,6 +74,9 @@ fmprb_zeta_ui_vec_borwein(fmprb_struct * z, ulong start, long num,
fmpz_t c, d, t, u;
fmpz * zeta;
if (num < 1)
return;
wp = prec + FLINT_BIT_COUNT(prec);
n = wp / 2.5431066063272239453 + 1;

22
fmprb_poly/log_gamma_series.c
View file

@ -36,22 +36,12 @@ fmprb_poly_log_gamma_series(fmprb_poly_t z, long n, long prec)
fmprb_poly_fit_length(z, n);
_fmprb_poly_set_length(z, n);
for (i = 0; i < n; i++)
{
if (i == 0)
{
fmprb_zero(z->coeffs + i);
}
else if (i == 1)
{
fmprb_const_euler_brent_mcmillan(z->coeffs + i, prec);
}
else
{
fmprb_zeta_ui(z->coeffs + i, i, prec);
fmprb_div_ui(z->coeffs + i, z->coeffs + i, i, prec);
}
}
if (n > 0) fmprb_zero(z->coeffs);
if (n > 1) fmprb_const_euler_brent_mcmillan(z->coeffs + 1, prec);
if (n > 2) fmprb_zeta_ui_vec(z->coeffs + 2, 2, n - 2, prec);
for (i = 2; i < n; i++)
fmprb_div_ui(z->coeffs + i, z->coeffs + i, i, prec);
_fmprb_poly_normalise(z);
}