use new types for bernoulli iteration code

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;

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);
     }
 }
 

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);
 }
 

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);
 }