use binary splitting to compute consecutive logarithms faster at high…

… precision
flintlib · Aug 3, 2013 · fea196d · fea196d
1 parent 2689439
commit fea196d
Show file tree

Hide file tree

Showing 4 changed files with 169 additions and 7 deletions.
diff --git a/zeta.h b/zeta.h
@@ -50,6 +50,8 @@ void zeta_series_em_bound(fmpr_t bound, const fmpcb_t s, const fmpcb_t a, long N
 void zeta_series_em_vec_bound(fmprb_ptr vec, const fmpcb_t s, const fmpcb_t a, ulong N, ulong M, long d, long wp);
 void zeta_series(fmpcb_ptr z, const fmpcb_t s, const fmpcb_t a, int deflate, long d, long prec);
 
+void zeta_log_ui_from_prev(fmprb_t s, ulong k, fmprb_t log_prev, ulong prev, long prec);
+
 void zeta_powsum_series_naive(fmpcb_ptr z, const fmpcb_t s, const fmpcb_t a, long n, long len, long prec);
 void zeta_powsum_one_series_sieved(fmpcb_ptr z, const fmpcb_t s, long n, long len, long prec);
 

diff --git a/zeta/log_ui_from_prev.c b/zeta/log_ui_from_prev.c
@@ -0,0 +1,72 @@
+/*=============================================================================
+
+    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) 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "zeta.h"
+#include "hypgeom.h"
+
+static void
+atanh_bsplit(fmprb_t s, ulong p, ulong q, long prec)
+{
+    fmprb_t t;
+    hypgeom_t series;
+    hypgeom_init(series);
+    fmprb_init(t);
+
+    fmpz_poly_set_ui(series->A, 1);
+    fmpz_poly_set_str(series->B, "2  1 2");
+    fmpz_poly_set_ui(series->P, p);
+    fmpz_mul(series->P->coeffs, series->P->coeffs, series->P->coeffs);
+    fmpz_poly_set_ui(series->Q, q);
+    fmpz_mul(series->Q->coeffs, series->Q->coeffs, series->Q->coeffs);
+
+    fmprb_hypgeom_infsum(s, t, series, prec, prec);
+    fmprb_mul_ui(s, s, p, prec);
+    fmprb_mul_ui(t, t, q, prec);
+    fmprb_div(s, s, t, prec);
+
+    fmprb_clear(t);
+    hypgeom_clear(series);
+}
+
+void
+zeta_log_ui_from_prev(fmprb_t s, ulong k, fmprb_t log_prev, ulong prev, long prec)
+{
+    if (k > 200 && prec > 5000)
+    {
+        fmprb_t t;
+        fmprb_init(t);
+
+        atanh_bsplit(t, k - prev, k + prev, prec);
+        fmprb_mul_2exp_si(t, t, 1);
+        fmprb_add(s, log_prev, t, prec);
+
+        fmprb_clear(t);
+    }
+    else
+    {
+        fmprb_log_ui(s, k, prec);
+    }
+}
+
diff --git a/zeta/powsum_one_series_sieved.c b/zeta/powsum_one_series_sieved.c
@@ -30,11 +30,10 @@
 void _fmpcb_poly_fmpcb_invpow_cpx(fmpcb_ptr res,
     const fmpcb_t N, const fmpcb_t c, long trunc, long prec);
 
-
 #define POWER(_k) (powers + (((_k)-1)/2) * (len))
 #define DIVISOR(_k) (divisors[((_k)-1)/2])
 
-#define COMPUTE_POWER(t, k) \
+#define COMPUTE_POWER(t, k, kprev) \
   do { \
     if (integer) \
     { \
@@ -44,13 +43,15 @@ void _fmpcb_poly_fmpcb_invpow_cpx(fmpcb_ptr res,
         fmprb_zero(fmpcb_imagref(t)); \
         if (len != 1) \
         { \
-            fmprb_log_ui(logk, k, prec); \
+            zeta_log_ui_from_prev(logk, k, logk, kprev, prec); \
+            kprev = k; \
             fmprb_neg(logk, logk); \
         } \
     } \
     else \
     { \
-        fmprb_log_ui(logk, k, prec); \
+        zeta_log_ui_from_prev(logk, k, logk, kprev, prec); \
+        kprev = k; \
         fmprb_neg(logk, logk); \
         fmprb_mul(w, logk, fmpcb_imagref(s), prec); \
         fmprb_sin_cos(fmpcb_imagref(t), fmpcb_realref(t), w, prec); \
@@ -71,14 +72,15 @@ void _fmpcb_poly_fmpcb_invpow_cpx(fmpcb_ptr res,
         fmpcb_mul_fmprb(t + i, t + i - 1, logk, prec); \
         fmpcb_div_ui(t + i, t + i, i, prec); \
     } \
+    fmprb_neg(logk, logk); \
   } while (0); \
 
 void
 zeta_powsum_one_series_sieved(fmpcb_ptr z, const fmpcb_t s, long n, long len, long prec)
 {
     long * divisors;
     long powers_alloc;
-    long i, j, k, power_of_two, horner_point;
+    long i, j, k, kprev, power_of_two, horner_point;
     int critical_line, integer;
 
     fmpcb_ptr powers;
@@ -113,14 +115,15 @@ zeta_powsum_one_series_sieved(fmpcb_ptr z, const fmpcb_t s, long n, long len, lo
 
     _fmpcb_vec_zero(z, len);
 
-    COMPUTE_POWER(x, 2);
+    kprev = 0;
+    COMPUTE_POWER(x, 2, kprev);
 
     for (k = 1; k <= n; k += 2)
     {
         /* t = k^(-s) */
         if (DIVISOR(k) == 0)
         {
-            COMPUTE_POWER(t, k);
+            COMPUTE_POWER(t, k, kprev);
         }
         else
         {

diff --git a/zeta/test/t-log_ui_from_prev.c b/zeta/test/t-log_ui_from_prev.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 "zeta.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("log_ui_from_prev....");
+    fflush(stdout);
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000; iter++)
+    {
+        fmprb_t z1, z2, z3;
+        ulong n1, n2;
+        long prec, accuracy;
+
+        prec = 2 + n_randint(state, 20000);
+        n1 = n_randint(state, 100000);
+        n2 = n1 + 1 + n_randint(state, 10);
+
+        fmprb_init(z1);
+        fmprb_init(z2);
+        fmprb_init(z3);
+
+        fmprb_log_ui(z1, n1, prec);
+        fmprb_log_ui(z2, n2, prec);
+
+        zeta_log_ui_from_prev(z3, n2, z1, n1, prec);
+
+        if (!fmprb_overlaps(z2, z3))
+        {
+            printf("FAIL: overlap\n\n");
+            printf("prec = %ld, n1 = %lu, n2 = %lu\n\n", prec, n1, n2);
+            printf("z1 = "); fmprb_printd(z1, prec / 3.33); printf("\n\n");
+            printf("z2 = "); fmprb_printd(z2, prec / 3.33); printf("\n\n");
+            printf("z3 = "); fmprb_printd(z3, prec / 3.33); printf("\n\n");
+            abort();
+        }
+
+        accuracy = fmprb_rel_accuracy_bits(z3);
+
+        if (accuracy < prec - 4)
+        {
+            printf("FAIL: accuracy = %ld, prec = %ld\n\n", accuracy, prec);
+            printf("n1 = %lu, n2 = %lu\n\n", n1, n2);
+            printf("z3 = "); fmprb_printd(z3, prec / 3.33); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(z1);
+        fmprb_clear(z2);
+        fmprb_clear(z3);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}