diff --git a/bernoulli.h b/bernoulli.h index 166ced66..df04ac0c 100644 --- a/bernoulli.h +++ b/bernoulli.h @@ -33,6 +33,7 @@ #include "fmpq.h" #include "arith.h" #include "fmprb.h" +#include "arb.h" #ifdef __cplusplus extern "C" { @@ -93,8 +94,8 @@ typedef struct long max_power; fmpz * powers; fmpz_t pow_error; - fmprb_t prefactor; - fmprb_t two_pi_squared; + arb_t prefactor; + arb_t two_pi_squared; ulong n; } bernoulli_rev_struct; diff --git a/bernoulli/rev_clear.c b/bernoulli/rev_clear.c index 7b52c476..d7dd838b 100644 --- a/bernoulli/rev_clear.c +++ b/bernoulli/rev_clear.c @@ -32,8 +32,8 @@ bernoulli_rev_clear(bernoulli_rev_t iter) { _fmpz_vec_clear(iter->powers, iter->alloc); fmpz_clear(iter->pow_error); - fmprb_clear(iter->prefactor); - fmprb_clear(iter->two_pi_squared); + arb_clear(iter->prefactor); + arb_clear(iter->two_pi_squared); } } diff --git a/bernoulli/rev_init.c b/bernoulli/rev_init.c index eca123d6..a354133a 100644 --- a/bernoulli/rev_init.c +++ b/bernoulli/rev_init.c @@ -30,7 +30,8 @@ bernoulli_rev_init(bernoulli_rev_t iter, ulong nmax) { long j; fmpz_t t; - fmprb_t x; + arb_t x; + arf_t u; int round1, round2; long wp; @@ -47,39 +48,42 @@ bernoulli_rev_init(bernoulli_rev_t iter, ulong nmax) iter->alloc = iter->max_power + 1; iter->powers = _fmpz_vec_init(iter->alloc); fmpz_init(iter->pow_error); - fmprb_init(iter->prefactor); - fmprb_init(iter->two_pi_squared); + arb_init(iter->prefactor); + arb_init(iter->two_pi_squared); - fmprb_init(x); + arb_init(x); fmpz_init(t); + arf_init(u); /* precompute powers */ for (j = 3; j <= iter->max_power; j += 2) { - fmprb_ui_pow_ui(x, j, nmax, power_prec(j, nmax, wp)); - fmprb_inv(x, x, power_prec(j, nmax, wp)); - round1 = fmpr_get_fmpz_fixed_si(t, fmprb_midref(x), -wp); + arb_ui_pow_ui(x, j, nmax, power_prec(j, nmax, wp)); + arb_inv(x, x, power_prec(j, nmax, wp)); + round1 = arf_get_fmpz_fixed_si(t, arb_midref(x), -wp); fmpz_set(iter->powers + j, t); /* error: the radius, plus two roundings */ - round2 = fmpr_get_fmpz_fixed_si(t, fmprb_radref(x), -wp); + arf_set_mag(u, arb_radref(x)); + round2 = arf_get_fmpz_fixed_si(t, u, -wp); fmpz_add_ui(t, t, (round1 != 0) + (round2 != 0)); if (fmpz_cmp(iter->pow_error, t) < 0) fmpz_set(iter->pow_error, t); } /* precompute (2pi)^2 and 2*(n!)/(2pi)^n */ - fmprb_fac_ui(iter->prefactor, nmax, wp); - fmprb_mul_2exp_si(iter->prefactor, iter->prefactor, 1); + arb_fac_ui(iter->prefactor, nmax, wp); + arb_mul_2exp_si(iter->prefactor, iter->prefactor, 1); - fmprb_const_pi(x, wp); - fmprb_mul_2exp_si(x, x, 1); - fmprb_mul(iter->two_pi_squared, x, x, wp); + arb_const_pi(x, wp); + arb_mul_2exp_si(x, x, 1); + arb_mul(iter->two_pi_squared, x, x, wp); - fmprb_pow_ui(x, iter->two_pi_squared, nmax / 2, wp); - fmprb_div(iter->prefactor, iter->prefactor, x, wp); + arb_pow_ui(x, iter->two_pi_squared, nmax / 2, wp); + arb_div(iter->prefactor, iter->prefactor, x, wp); fmpz_clear(t); - fmprb_clear(x); + arb_clear(x); + arf_clear(u); } diff --git a/bernoulli/rev_next.c b/bernoulli/rev_next.c index 22aede32..4f2d1042 100644 --- a/bernoulli/rev_next.c +++ b/bernoulli/rev_next.c @@ -25,14 +25,24 @@ #include "bernoulli.h" +static __inline__ void +mag_ui_div(mag_t z, ulong c, const mag_t x) +{ + mag_t t; + mag_init(t); + mag_set_ui(t, c); + mag_div(z, t, x); + mag_clear(t); +} + void bernoulli_rev_next(fmpz_t numer, fmpz_t denom, bernoulli_rev_t iter) { ulong n; long j, wp; fmpz_t sum; - fmpr_t err; - fmprb_t z, h; + mag_t err; + arb_t z, h; n = iter->n; wp = iter->prec; @@ -46,42 +56,42 @@ bernoulli_rev_next(fmpz_t numer, fmpz_t denom, bernoulli_rev_t iter) } fmpz_init(sum); - fmpr_init(err); - fmprb_init(z); - fmprb_init(h); + mag_init(err); + arb_init(z); + arb_init(h); /* add all odd powers */ fmpz_zero(sum); for (j = iter->max_power; j >= 3; j -= 2) fmpz_add(sum, sum, iter->powers + j); - fmprb_set_fmpz(z, sum); + arb_set_fmpz(z, sum); /* bound numerical error from the powers */ fmpz_mul_ui(sum, iter->pow_error, iter->max_power / 2); - fmpr_set_fmpz(err, sum); - fmprb_add_error_fmpr(z, err); - fmprb_mul_2exp_si(z, z, -wp); + mag_set_fmpz(err, sum); + mag_add(arb_radref(z), arb_radref(z), err); + arb_mul_2exp_si(z, z, -wp); - fmprb_add_ui(z, z, 1, wp); + arb_add_ui(z, z, 1, wp); /* add truncation error: sum_{k > N} 1/k^n <= 1/N^(i-1) */ - fmpr_set_ui(err, iter->max_power); - fmpr_pow_sloppy_ui(err, err, n - 1, FMPRB_RAD_PREC, FMPR_RND_DOWN); - fmpr_ui_div(err, 1, err, FMPRB_RAD_PREC, FMPR_RND_UP); - fmprb_add_error_fmpr(z, err); + mag_set_ui_lower(err, iter->max_power); + mag_pow_ui_lower(err, err, n - 1); + mag_ui_div(err, 1, err); + mag_add(arb_radref(z), arb_radref(z), err); /* convert zeta to Bernoulli number */ - fmprb_div_2expm1_ui(h, z, n, wp); - fmprb_add(z, z, h, wp); - fmprb_mul(z, z, iter->prefactor, wp); + arb_div_2expm1_ui(h, z, n, wp); + arb_add(z, z, h, wp); + arb_mul(z, z, iter->prefactor, wp); arith_bernoulli_number_denom(denom, n); - fmprb_mul_fmpz(z, z, denom, wp); + arb_mul_fmpz(z, z, denom, wp); if (n % 4 == 0) - fmprb_neg(z, z); + arb_neg(z, z); - /* printf("%ld: ", n); fmprb_printd(z, 5); printf("\n"); */ + /* printf("%ld: ", n); arb_printd(z, 5); printf("\n"); */ - if (!fmprb_get_unique_fmpz(numer, z)) + if (!arb_get_unique_fmpz(numer, z)) { printf("warning: insufficient precision for B_%ld\n", n); _bernoulli_fmpq_ui(numer, denom, n); @@ -90,9 +100,9 @@ bernoulli_rev_next(fmpz_t numer, fmpz_t denom, bernoulli_rev_t iter) /* update prefactor */ if (n > 0) { - fmprb_mul(iter->prefactor, iter->prefactor, iter->two_pi_squared, wp); - fmprb_div_ui(iter->prefactor, iter->prefactor, n, wp); - fmprb_div_ui(iter->prefactor, iter->prefactor, n - 1, wp); + arb_mul(iter->prefactor, iter->prefactor, iter->two_pi_squared, wp); + arb_div_ui(iter->prefactor, iter->prefactor, n, wp); + arb_div_ui(iter->prefactor, iter->prefactor, n - 1, wp); } /* update powers */ @@ -126,7 +136,7 @@ bernoulli_rev_next(fmpz_t numer, fmpz_t denom, bernoulli_rev_t iter) fmpz_add_ui(iter->pow_error, iter->pow_error, 1); /* speed improvement (could be skipped with better multiplication) */ - fmprb_set_round(iter->two_pi_squared, iter->two_pi_squared, new_prec); + arb_set_round(iter->two_pi_squared, iter->two_pi_squared, new_prec); iter->max_power = new_max; iter->prec = new_prec; @@ -136,8 +146,8 @@ bernoulli_rev_next(fmpz_t numer, fmpz_t denom, bernoulli_rev_t iter) iter->n -= 2; fmpz_clear(sum); - fmpr_clear(err); - fmprb_clear(z); - fmprb_clear(h); + mag_clear(err); + arb_clear(z); + arb_clear(h); }