code for gamma and zeta(n); some bugfixes

fmpr.h

@@ -131,6 +131,18 @@ static __inline__ void fmpr_nan(fmpr_t x) {
 
 
 
+static __inline__ int
+fmpr_is_one(const fmpr_t x)
+{
+    return fmpz_is_one(fmpr_manref(x)) && fmpz_is_zero(fmpr_expref(x));
+}
+
+static __inline__ void
+fmpr_one(fmpr_t x)
+{
+    fmpz_one(fmpr_manref(x));
+    fmpz_zero(fmpr_expref(x));
+}
 
 /* ------------------------------------------------------------------------ */
 
@@ -350,6 +362,7 @@ long fmpr_submul(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t
 
 long fmpr_sqrt(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd);
 
+long fmpr_log(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd);
 
 
 static __inline__ long fmpr_sqrt_ui(fmpr_t z, ulong x, long prec, fmpr_rnd_t rnd)

fmpr/log.c

@@ -0,0 +1,85 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpr.h"
+
+static mpfr_rnd_t fmpr_rnd_to_mpfr(fmpr_rnd_t rnd)
+{
+    if (rnd == FMPR_RND_NEAR) return MPFR_RNDN;
+    if (rnd == FMPR_RND_FLOOR) return MPFR_RNDD;
+    if (rnd == FMPR_RND_CEIL) return MPFR_RNDU;
+    if (rnd == FMPR_RND_DOWN) return MPFR_RNDZ;
+    if (rnd == FMPR_RND_UP) return MPFR_RNDA;
+    abort();
+}
+
+long
+fmpr_log(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
+{
+    if (fmpr_is_special(x))
+    {
+        if (fmpr_is_zero(x))
+            fmpr_neg_inf(y);
+        else if (fmpr_is_pos_inf(x))
+            fmpr_pos_inf(y);
+        else
+            fmpr_nan(y);
+
+        return FMPR_RESULT_EXACT;
+    }
+
+    if (fmpr_sgn(x) < 0)
+    {
+        fmpr_nan(y);
+        return FMPR_RESULT_EXACT;
+    }
+
+    if (fmpr_is_one(x))
+    {
+        fmpr_zero(y);
+        return FMPR_RESULT_EXACT;
+    }
+
+    {
+        mpfr_t t, u;
+        long r;
+
+        mpfr_init2(t, 1 + fmpz_bits(fmpr_manref(x)));
+        mpfr_init2(u, FLINT_MAX(2, prec));
+
+        fmpr_get_mpfr(t, x, MPFR_RNDD);
+
+        mpfr_log(u, t, fmpr_rnd_to_mpfr(rnd));
+
+        fmpr_set_mpfr(y, u);
+
+        r = prec - fmpz_bits(fmpr_manref(y));
+
+        mpfr_clear(t);
+        mpfr_clear(u);
+
+        return r;
+    }
+}

fmpr/sqrt.c

@@ -92,26 +92,3 @@ fmpr_sqrt(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
 
     return _fmpr_normalise(fmpr_manref(y), fmpr_expref(y), prec, rnd);
 }
-
-
-
-/*
-    if exp & 1:
-        exp -= 1
-        man <<= 1
-        bc += 1
-    elif man == 1:
-        return normalize1(sign, man, exp//2, bc, prec, rnd)
-    shift = max(4, 2*prec-bc+4)
-    shift += shift & 1
-    if rnd in 'fd':
-        man = isqrt(man<<shift)
-    else:
-        man, rem = sqrtrem(man<<shift)
-        # Perturb up
-        if rem:
-            man = (man<<1)+1
-            shift += 2
-    return from_man_exp(man, (exp-shift)//2, prec, rnd)
-*/
-

fmpr/test/t-log.c

@@ -0,0 +1,103 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpr.h"
+
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("log....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        long bits;
+        fmpr_t x, z, w;
+        mpfr_t X, Z;
+
+        bits = 2 + n_randint(state, 200);
+
+        fmpr_init(x);
+        fmpr_init(z);
+        fmpr_init(w);
+
+        mpfr_init2(X, bits + 100);
+        mpfr_init2(Z, bits);
+
+        fmpr_randtest_special(x, state, bits + n_randint(state, 100), 10);
+        fmpr_randtest_special(z, state, bits + n_randint(state, 100), 10);
+
+        fmpr_get_mpfr(X, x, MPFR_RNDN);
+
+        switch (n_randint(state, 4))
+        {
+            case 0:
+                mpfr_log(Z, X, MPFR_RNDZ);
+                fmpr_log(z, x, bits, FMPR_RND_DOWN);
+                break;
+            case 1:
+                mpfr_log(Z, X, MPFR_RNDA);
+                fmpr_log(z, x, bits, FMPR_RND_UP);
+                break;
+            case 2:
+                mpfr_log(Z, X, MPFR_RNDD);
+                fmpr_log(z, x, bits, FMPR_RND_FLOOR);
+                break;
+            case 3:
+                mpfr_log(Z, X, MPFR_RNDU);
+                fmpr_log(z, x, bits, FMPR_RND_CEIL);
+                break;
+        }
+
+        fmpr_set_mpfr(w, Z);
+
+        if (!fmpr_equal(z, w))
+        {
+            printf("FAIL\n\n");
+            printf("bits = %ld\n", bits);
+            printf("x = "); fmpr_print(x); printf("\n\n");
+            printf("z = "); fmpr_print(z); printf("\n\n");
+            printf("w = "); fmpr_print(w); printf("\n\n");
+            abort();
+        }
+
+        fmpr_clear(x);
+        fmpr_clear(z);
+        fmpr_clear(w);
+
+        mpfr_clear(X);
+        mpfr_clear(Z);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

fmprb.h

@@ -64,6 +64,13 @@ fmprb_is_exact(const fmprb_t x)
     return fmpr_is_zero(fmprb_radref(x));
 }
 
+static __inline__ void
+fmprb_zero(fmprb_t x)
+{
+    fmpr_zero(fmprb_midref(x));
+    fmpr_zero(fmprb_radref(x));
+}
+
 static __inline__ void
 fmprb_set(fmprb_t x, const fmprb_t y)
 {
@@ -150,7 +157,29 @@ void fmprb_submul_ui(fmprb_t z, const fmprb_t x, ulong y, long prec);
 void fmprb_submul_si(fmprb_t z, const fmprb_t x, long y, long prec);
 void fmprb_submul_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec);
 
+void fmprb_ui_pow_ui(fmprb_t y, ulong b, ulong e, long prec);
+void fmprb_si_pow_ui(fmprb_t y, long b, ulong e, long prec);
+
+void fmprb_log(fmprb_t z, const fmprb_t x, long prec);
+void fmprb_log_ui(fmprb_t z, ulong x, long prec);
+void fmprb_log_fmpz(fmprb_t z, const fmpz_t x, long prec);
+
 void fmprb_const_pi_chudnovsky(fmprb_t x, long prec);
+void fmprb_const_euler_brent_mcmillan(fmprb_t res, long prec);
+void fmprb_const_zeta3_bsplit(fmprb_t x, long prec);
+
+void fmprb_zeta_ui_bsplit(fmprb_t x, ulong s, long prec);
+
+static __inline__ void
+fmprb_add_error_2exp_si(fmprb_t x, long err)
+{
+    fmpr_t t;
+    fmpr_init(t);
+    fmpz_set_ui(fmpr_manref(t), 1);
+    fmpz_set_si(fmpr_expref(t), err);
+    fmpr_add(fmprb_radref(x), fmprb_radref(x), t, FMPRB_RAD_PREC, FMPR_RND_UP);
+    fmpr_clear(t);
+}
 
 static __inline__ void
 fmprb_printd(const fmprb_t x, long digits)

fmprb/add.c

@@ -30,7 +30,7 @@ fmprb_add(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
 {
     long r;
 
-    fmpr_add(fmprb_radref(z), fmprb_radref(z), fmprb_radref(z), FMPRB_RAD_PREC, FMPR_RND_UP);
+    fmpr_add(fmprb_radref(z), fmprb_radref(x), fmprb_radref(y), FMPRB_RAD_PREC, FMPR_RND_UP);
     r = fmpr_add(fmprb_midref(z), fmprb_midref(x), fmprb_midref(y), prec, FMPR_RND_DOWN);
 
     fmpr_add_error_result(fmprb_radref(z), fmprb_radref(z), fmprb_midref(z),

fmprb/const_euler_brent_mcmillan.c

@@ -0,0 +1,235 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb.h"
+
+typedef struct
+{
+    fmprb_t P;
+    fmprb_t Q;
+    fmprb_t T;
+    fmprb_t C;
+    fmprb_t D;
+    fmprb_t V;
+} euler_bsplit_struct;
+
+typedef euler_bsplit_struct euler_bsplit_t[1];
+
+static void euler_bsplit_init(euler_bsplit_t s)
+{
+    fmprb_init(s->P);
+    fmprb_init(s->Q);
+    fmprb_init(s->T);
+    fmprb_init(s->C);
+    fmprb_init(s->D);
+    fmprb_init(s->V);
+}
+
+static void euler_bsplit_clear(euler_bsplit_t s)
+{
+    fmprb_clear(s->P);
+    fmprb_clear(s->Q);
+    fmprb_clear(s->T);
+    fmprb_clear(s->C);
+    fmprb_clear(s->D);
+    fmprb_clear(s->V);
+}
+
+static void
+euler_bsplit_1_merge(euler_bsplit_t s, euler_bsplit_t L, euler_bsplit_t R,
+                    long wp, int cont)
+{
+    fmprb_t t, u, v;
+
+    fmprb_init(t);
+    fmprb_init(u);
+    fmprb_init(v);
+
+    if (cont)
+        fmprb_mul(s->P, L->P, R->P, wp);
+
+    fmprb_mul(s->Q, L->Q, R->Q, wp);
+    fmprb_mul(s->D, L->D, R->D, wp);
+
+    /* T = LP RT + RQ LT*/
+    fmprb_mul(t, L->P, R->T, wp);
+    fmprb_mul(v, R->Q, L->T, wp);
+    fmprb_add(s->T, t, v, wp);
+
+    /* C = LC RD + RC LD */
+    if (cont)
+    {
+        fmprb_mul(s->C, L->C, R->D, wp);
+        fmprb_addmul(s->C, R->C, L->D, wp);
+    }
+
+    /* V = RD (RQ LV + LC LP RT) + LD LP RV */
+    fmprb_mul(u, L->P, R->V, wp);
+    fmprb_mul(u, u, L->D, wp);
+    fmprb_mul(v, R->Q, L->V, wp);
+    fmprb_addmul(v, t, L->C, wp);
+    fmprb_mul(v, v, R->D, wp);
+    fmprb_add(s->V, u, v, wp);
+
+    fmprb_clear(t);
+    fmprb_clear(u);
+    fmprb_clear(v);
+}
+
+void
+euler_bsplit_1(euler_bsplit_t s, long n1, long n2, long N, long wp, int cont)
+{
+    if (n2 - n1 == 1)
+    {
+        fmprb_set_si(s->P, N); /* p = N^2 todo: shift optimization */
+        fmprb_mul(s->P, s->P, s->P, wp);
+        fmprb_set_si(s->Q, n1 + 1); /* q = (k + 1)^2 */
+        fmprb_mul(s->Q, s->Q, s->Q, wp);
+        fmprb_set_si(s->C, 1);
+        fmprb_set_si(s->D, n1 + 1);
+        fmprb_set(s->T, s->P);
+        fmprb_set(s->V, s->P);
+    }
+    else
+    {
+        euler_bsplit_t L, R;
+        long m = (n1 + n2) / 2;
+
+        euler_bsplit_init(L);
+        euler_bsplit_init(R);
+        euler_bsplit_1(L, n1, m, N, wp, 1);
+        euler_bsplit_1(R, m, n2, N, wp, 1);
+        euler_bsplit_1_merge(s, L, R, wp, cont);
+        euler_bsplit_clear(L);
+        euler_bsplit_clear(R);
+    }
+}
+
+void
+euler_bsplit_2(fmprb_t P, fmprb_t Q, fmprb_t T, long n1, long n2,
+                        long N, long wp, int cont)
+{
+    if (n2 - n1 == 1)
+    {
+        if (n1 == 0)
+        {
+            fmprb_set_si(P, 1);
+            fmprb_set_si(Q, 4 * N);
+            fmprb_set_si(T, 1);
+        }
+        else
+        {
+            fmprb_si_pow_ui(P, 1 - 2*n1, 3, wp);
+            fmprb_neg(P, P);
+
+            fmprb_set_si(Q, 32 * n1);
+            fmprb_mul_ui(Q, Q, N, wp);
+            fmprb_mul_ui(Q, Q, N, wp);
+        }
+
+        fmprb_set(T, P);
+    }
+    else
+    {
+        fmprb_t P2, Q2, T2;
+        long m = (n1 + n2) / 2;
+
+        fmprb_init(P2);
+        fmprb_init(Q2);
+        fmprb_init(T2);
+
+        euler_bsplit_2(P, Q, T, n1, m, N, wp, 1);
+        euler_bsplit_2(P2, Q2, T2, m, n2, N, wp, 1);
+
+        fmprb_mul(T, T, Q2, wp);
+        fmprb_mul(T2, T2, P, wp);
+        fmprb_add(T, T, T2, wp);
+
+        if (cont)
+            fmprb_mul(P, P, P2, wp);
+
+        fmprb_mul(Q, Q, Q2, wp);
+
+        fmprb_clear(P2);
+        fmprb_clear(Q2);
+        fmprb_clear(T2);
+    }
+}
+
+void
+fmprb_const_euler_brent_mcmillan(fmprb_t res, long prec)
+{
+    euler_bsplit_t sum;
+    fmprb_t t, u, v, P2, T2, Q2;
+    long bits, wp, n, nterms1, nterms2;
+
+    bits = prec + 20;
+    n = 0.08665 * bits + 1;
+
+    nterms1 = 4.9706258 * n + 1;
+    nterms2 = 2 * n + 1;
+    wp = bits + FLINT_BIT_COUNT(n);
+
+    euler_bsplit_init(sum);
+    fmprb_init(P2);
+    fmprb_init(T2);
+    fmprb_init(Q2);
+    fmprb_init(t);
+    fmprb_init(u);
+    fmprb_init(v);
+
+    /* Compute S0 = V / (Q * D), I0 = 1 + T / Q */
+    euler_bsplit_1(sum, 0, nterms1, n, wp, 0);
+
+    /* Compute K0 = T2 / Q2 */
+    euler_bsplit_2(P2, Q2, T2, 0, nterms2, n, wp, 0);
+
+    /* Compute (S0/I0 + K0/I0^2) = (Q2*(Q+T)*V - D*Q^2*T2)/(D*Q2*(Q+T)^2) */
+    fmprb_add(v, sum->Q, sum->T, wp);
+    fmprb_mul(t, v, Q2, wp);
+    fmprb_mul(u, sum->Q, sum->Q, wp);
+    fmprb_mul(u, u, T2, wp);
+    fmprb_mul(u, u, sum->D, wp);
+    fmprb_mul(sum->V, t, sum->V, wp);
+    fmprb_sub(sum->V, sum->V, u, wp);
+    fmprb_mul(u, sum->D, t, wp);
+    fmprb_mul(u, u, v, wp);
+    fmprb_div(t, sum->V, u, wp);
+
+    /* subtract log(n) */
+    fmprb_log_ui(u, n, wp);
+    fmprb_sub(res, t, u, wp);
+
+    /* TODO: add error term */
+
+    fmprb_clear(P2);
+    fmprb_clear(T2);
+    fmprb_clear(Q2);
+    fmprb_clear(t);
+    fmprb_clear(u);
+    fmprb_clear(v);
+
+    euler_bsplit_clear(sum);
+}

fmprb/const_zeta3_bsplit.c

@@ -0,0 +1,90 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb.h"
+
+static void
+zeta3_bsplit(fmprb_t P, fmprb_t Q, fmprb_t T, long a, long b, long wp, int cont)
+{
+    if (b - a == 1)
+    {
+        fmprb_ui_pow_ui(P, a ? a : 1, 5, wp);
+
+        fmprb_ui_pow_ui(Q, 2*a + 1, 5, wp);
+        fmprb_mul_ui(Q, Q, a ? 32 : 64, wp);
+
+        fmprb_set_ui(T, 205*a*a + 250*a + 77); /* xxx: fix overflow */
+        if (a % 2)
+            fmprb_neg(T, T);
+        fmprb_mul(T, T, P, wp);
+    }
+    else
+    {
+        long m;
+        fmprb_t P2, Q2, T2;
+
+        m = (a + b) / 2;
+
+        fmprb_init(P2);
+        fmprb_init(Q2);
+        fmprb_init(T2);
+
+        zeta3_bsplit(P, Q, T, a, m, wp, 1);
+        zeta3_bsplit(P2, Q2, T2, m, b, wp, 1);
+
+        fmprb_mul(T, T, Q2, wp);
+        fmprb_addmul(T, P, T2, wp);
+        fmprb_mul(Q, Q, Q2, wp);
+        if (cont)
+            fmprb_mul(P, P, P2, wp);
+
+        fmprb_clear(P2);
+        fmprb_clear(Q2);
+        fmprb_clear(T2);
+    }
+}
+
+void
+fmprb_const_zeta3_bsplit(fmprb_t x, long prec)
+{
+    long wp, N;
+    fmprb_t P, Q, T;
+
+    wp = prec + 10 + FLINT_BIT_COUNT(prec);
+    N = prec / 10 + 2;
+
+    fmprb_init(P);
+    fmprb_init(Q);
+    fmprb_init(T);
+
+    zeta3_bsplit(P, Q, T, 0, N, wp, 0);
+    fmprb_div(x, T, Q, prec);
+
+    fmprb_add_error_2exp_si(x, -10 * (N - 1));
+
+    fmprb_clear(P);
+    fmprb_clear(Q);
+    fmprb_clear(T);
+}

fmprb/log.c

@@ -0,0 +1,61 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb.h"
+
+void
+fmprb_log(fmprb_t z, const fmprb_t x, long prec)
+{
+    long r;
+
+    if (!fmprb_is_exact(x))
+        abort();
+
+    r = fmpr_log(fmprb_midref(z), fmprb_midref(x), prec, FMPR_RND_DOWN);
+
+    fmpr_set_error_result(fmprb_radref(z), fmprb_midref(z), r);
+
+    fmprb_adjust(z);
+}
+
+void
+fmprb_log_ui(fmprb_t z, ulong x, long prec)
+{
+    fmprb_t t;
+    fmprb_init(t);
+    fmprb_set_ui(t, x);
+    fmprb_log(z, t, prec);
+    fmprb_clear(t);
+}
+
+void
+fmprb_log_fmpz(fmprb_t z, const fmpz_t x, long prec)
+{
+    fmprb_t t;
+    fmprb_init(t);
+    fmprb_set_fmpz(t, x);
+    fmprb_log(z, t, prec);
+    fmprb_clear(t);
+}

fmprb/pow.c

@@ -0,0 +1,66 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb.h"
+
+void
+fmprb_ui_pow_ui(fmprb_t y, ulong b, ulong e, long prec)
+{
+    long i;
+
+    if (e <= 1)
+    {
+        fmprb_set_ui(y, e == 0 ? 1 : b);
+        return;
+    }
+
+    fmprb_set_ui(y, b);
+    for (i = FLINT_BIT_COUNT(e) - 2; i >= 0; i--)
+    {
+        fmprb_mul(y, y, y, prec);
+        if (e & (1UL<<i))
+            fmprb_mul_ui(y, y, b, prec);
+    }
+}
+
+void
+fmprb_si_pow_ui(fmprb_t y, long b, ulong e, long prec)
+{
+    long i;
+
+    if (e <= 1)
+    {
+        fmprb_set_si(y, e == 0 ? 1 : b);
+        return;
+    }
+
+    fmprb_set_si(y, b);
+    for (i = FLINT_BIT_COUNT(e) - 2; i >= 0; i--)
+    {
+        fmprb_mul(y, y, y, prec);
+        if (e & (1UL<<i))
+            fmprb_mul_si(y, y, b, prec);
+    }
+}

fmprb/sub.c

@@ -30,7 +30,8 @@ fmprb_sub(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
 {
     long r;
 
-    fmpr_add(fmprb_radref(z), fmprb_radref(z), fmprb_radref(z), FMPRB_RAD_PREC, FMPR_RND_UP);
+    fmpr_add(fmprb_radref(z), fmprb_radref(x), fmprb_radref(y), FMPRB_RAD_PREC, FMPR_RND_UP);
+
     r = fmpr_sub(fmprb_midref(z), fmprb_midref(x), fmprb_midref(y), prec, FMPR_RND_DOWN);
 
     fmpr_add_error_result(fmprb_radref(z), fmprb_radref(z), fmprb_midref(z),

fmprb/zeta_ui_bsplit.c

@@ -0,0 +1,235 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmprb.h"
+
+/* With parameter n, the error is bounded by 3/(3+sqrt(8))^n */
+#define ERROR_A 1.5849625007211561815 /* log2(3) */
+#define ERROR_B 2.5431066063272239453 /* log2(3+sqrt(8)) */
+
+typedef struct
+{
+    fmprb_t A;
+    fmprb_t B;
+    fmprb_t C;
+    fmprb_t D;
+    fmprb_t E;
+    fmprb_t Q1;
+    fmprb_t Q2;
+    fmprb_t Q3;
+}
+zeta_bsplit_state;
+
+typedef zeta_bsplit_state zeta_bsplit_t[1];
+
+static __inline__ void
+zeta_bsplit_init(zeta_bsplit_t S)
+{
+    fmprb_init(S->A);
+    fmprb_init(S->B);
+    fmprb_init(S->C);
+    fmprb_init(S->D);
+    fmprb_init(S->E);
+    fmprb_init(S->Q1);
+    fmprb_init(S->Q2);
+    fmprb_init(S->Q3);
+}
+
+static __inline__ void
+zeta_bsplit_clear(zeta_bsplit_t S)
+{
+    fmprb_clear(S->A);
+    fmprb_clear(S->B);
+    fmprb_clear(S->C);
+    fmprb_clear(S->D);
+    fmprb_clear(S->E);
+    fmprb_clear(S->Q1);
+    fmprb_clear(S->Q2);
+    fmprb_clear(S->Q3);
+}
+
+
+static __inline__ void
+zeta_coeff_k(zeta_bsplit_t S, long k, long n, long s)
+{
+    if (k + 1 < 0)
+    {
+        fmprb_set_si(S->D, 1);
+        fmprb_set_si(S->Q1, 1);
+    }
+    else if (k + 1 > n)
+    {
+        fmprb_zero(S->D);
+        fmprb_set_si(S->Q1, 1);
+    }
+    else
+    {
+        fmprb_set_si(S->D, 2 * (n + (k + 1) - 1));
+        fmprb_mul_si(S->D, S->D, n + 1 - (k + 1), FMPR_PREC_EXACT);
+        fmprb_set_si(S->Q1, k + 1);
+        fmprb_mul_si(S->Q1, S->Q1, 2*(k + 1) - 1, FMPR_PREC_EXACT);
+    }
+
+    if (k - 1 < 0)
+    {
+        fmprb_zero(S->E);
+        fmprb_set_si(S->Q2, 1);
+    }
+    else if (k - 1 >= n)
+    {
+        fmprb_set_si(S->E, 1);
+        fmprb_set_si(S->Q2, 1);
+    }
+    else
+    {
+        fmprb_set_si(S->E, ((k - 1) % 2) ? -1 : 1);
+        fmprb_set_si(S->Q2, k);
+        fmprb_ui_pow_ui(S->Q2, k, s, FMPR_PREC_EXACT);   /* XXX */
+    }
+
+    fmprb_mul(S->Q3, S->Q1, S->Q2, FMPR_PREC_EXACT);     /* XXX */
+    fmprb_mul(S->A, S->E, S->Q1, FMPR_PREC_EXACT);
+    fmprb_zero(S->B);
+    fmprb_set(S->C, S->Q1);
+}
+
+static void
+zeta_bsplit(zeta_bsplit_t L, long a, long b,
+    long n, long s, int cont, long bits)
+{
+    if (a + 1 == b)
+    {
+        zeta_coeff_k(L, a, n, s);
+    }
+    else
+    {
+        zeta_bsplit_t R;
+
+        long m = (a + b) / 2;
+
+        zeta_bsplit(L, m, b, n, s, 1, bits);
+
+        zeta_bsplit_init(R);
+        zeta_bsplit(R, a, m, n, s, 1, bits);
+
+        fmprb_mul(L->E, L->E, R->Q2, bits);
+        fmprb_addmul(L->E, R->E, L->Q2, bits);
+
+        fmprb_mul(L->B, L->B, R->D, bits);
+        fmprb_addmul(L->B, L->A, R->C, bits);
+
+        fmprb_mul(L->B, L->B, R->Q2, bits);
+        fmprb_addmul(L->B, R->B, L->Q3, bits);
+
+        if (cont)
+        {
+            fmprb_mul(L->A, L->A, R->Q3, bits);
+            fmprb_addmul(L->A, R->A, L->Q3, bits);
+        }
+
+        fmprb_mul(L->C, L->C, R->D, bits);
+        fmprb_addmul(L->C, R->C, L->Q1, bits);
+        fmprb_mul(L->Q2, L->Q2, R->Q2, bits);
+
+        if (cont)
+        {
+            fmprb_mul(L->D, L->D, R->D, bits);
+            fmprb_mul(L->Q1, L->Q1, R->Q1, bits);
+            fmprb_mul(L->Q3, L->Q3, R->Q3, bits);
+        }
+
+        zeta_bsplit_clear(R);
+    }
+}
+
+void
+fmprb_zeta_ui_bsplit(fmprb_t x, ulong s, long prec)
+{
+    zeta_bsplit_t sum;
+    long wp, n;
+    long i;
+
+    /* zeta(0) = -1/2 */
+    if (s == 0)
+    {
+        /* XXX */
+        fmpz_set_si(fmpr_manref(fmprb_midref(x)), -1);
+        fmpz_set_si(fmpr_expref(fmprb_midref(x)), -1);
+        fmpr_zero(fmprb_radref(x));
+        return;
+    }
+
+    if (s == 1)
+    {
+        printf("fmprb_zeta_ui_bsplit: zeta(1)");
+        abort();
+    }
+
+    /* for s > p, zeta(s) - 1 < 2^(-p) */
+    if (s != 2 && s > prec)
+    {
+        fmprb_set_ui(x, 1UL);
+
+        /* XXX: could make this even smaller when s is extremely large */
+        fmpz_one(fmpr_manref(fmprb_radref(x)));
+        fmpz_set_si(fmpr_expref(fmprb_radref(x)), -prec);
+
+        return;
+    }
+
+    n = prec / ERROR_B + 2;
+    wp = prec + 30;
+
+    zeta_bsplit_init(sum);
+    zeta_bsplit(sum, 0, n + 1, n, s, 0, wp);
+
+    fmprb_mul(sum->E, sum->E, sum->C, wp);
+    fmprb_sub(sum->E, sum->E, sum->B, wp);
+    fmprb_mul(sum->Q2, sum->Q2, sum->C, wp);
+    fmprb_div(sum->C, sum->E, sum->Q2, wp);
+
+    /* D = 1/(1 - 2^(1-s)) */
+    {
+        fmpz_t t;
+        fmpz_init(t);
+        for (i = wp; i >= 0; i -= (s - 1))
+            fmpz_setbit(t, i);
+
+        fmprb_set_fmpz(sum->D, t);
+        fmpz_sub_ui(fmpr_expref(fmprb_midref(sum->D)),
+            fmpr_expref(fmprb_midref(sum->D)), wp);
+
+        fmprb_add_error_2exp_si(sum->D, -wp);
+
+        fmpz_clear(t);
+    }
+
+    fmprb_mul(x, sum->C, sum->D, wp);
+
+    /* The truncation error is bounded by 3/(3+sqrt(8))^n */
+    fmprb_add_error_2exp_si(x, (long) (ERROR_A - ERROR_B * n + 1));
+
+    zeta_bsplit_clear(sum);
+}