initial code for computing Khinchin's constant

doc/source/credits.rst

@@ -28,6 +28,8 @@ The following software has been helpful in the development of Arb.
 Bibliography
 -------------------------------------------------------------------------------
 
+.. [BBC1997] D.H. Bailey, J. M. Borwein and R. E. Crandall, "On the Khintchine constant", Mathematics of Computation 66 (1997) 417-431
+
 .. [BBC2000] J.Borwein, D. M. Bradley and R. E. Crandall, "Computational strategies for the Riemann zeta function", Journal of Computational and Applied Mathematics 121 (2000) 247-296
 
 .. [BZ1992]_ J.Borwein and I.Zucker, "Fast evaluation of the gamma function for small rational fractions using complete elliptic integrals of the first kind", IMA Journal of Numerical Analysis 12 (1992) 519-526

doc/source/fmprb.rst

@@ -550,6 +550,35 @@ Constants
     recently proved that this bound is correct, but without determining
     an explicit big-O factor [Bre2010]_.
 
+.. function:: void fmprb_const_khinchin(fmprb_t res, long prec)
+
+    Sets *res* to Khinchin's constant `K_0`, computed as
+
+    .. math ::
+
+        \log K_0 = \frac{1}{\log 2} \left[
+        \sum_{k=2}^{N-1} \log \left(\frac{k-1}{k} \right) \log \left(\frac{k+1}{k} \right)
+        + \sum_{n=1}^\infty 
+        \frac {\zeta (2n,N)}{n} \sum_{k=1}^{2n-1} \frac{(-1)^{k+1}}{k}
+        \right]
+
+    where `N \ge 2` is a free parameter that can be used for tuning [BBC1997]_.
+    If the infinite series is truncated after `n = M`, the remainder
+    is smaller in absolute value than
+
+    .. math ::
+
+        \sum_{n=M+1}^{\infty} \zeta(2n, N) = 
+        \sum_{n=M+1}^{\infty} \sum_{k=0}^{\infty} (k+N)^{-2n} \le
+        \sum_{n=M+1}^{\infty} \left( N^{-2n} + \int_0^{\infty} (t+N)^{-2n} dt \right)
+
+        = \sum_{n=M+1}^{\infty} \frac{1}{N^{2n}} \left(1 + \frac{N}{2n-1}\right)
+        \le \sum_{n=M+1}^{\infty} \frac{N+1}{N^{2n}} = \frac{1}{N^{2M} (N-1)}
+        \le \frac{1}{N^{2M}}.
+
+    Thus, for an error of at most `2^{-p}` in the series,
+    it is sufficient to choose `M \ge p / (2 \log_2 N)`.
+
 Riemann zeta function
 ...............................................................................
 

doc/source/fmprb_poly.rst

@@ -66,6 +66,8 @@ Conversions
 
 .. function:: void fmprb_poly_set_fmpq_poly(fmprb_poly_t poly, const fmpq_poly_t src, long prec)
 
+.. function:: void fmprb_poly_set_si(fmprb_poly_t poly, long src)
+
     Sets *poly* to *src*, rounding the coefficients to *prec* bits.
 
 

fmprb.h

@@ -302,6 +302,8 @@ void fmprb_const_log_sqrt2pi(fmprb_t t, 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_const_khinchin(fmprb_t K, long prec);
+
 void fmprb_zeta_ui_asymp(fmprb_t x, ulong s, long prec);
 void fmprb_zeta_ui_bsplit(fmprb_t x, ulong s, long prec);
 void fmprb_zeta_ui_euler_product(fmprb_t z, ulong s, long prec);

fmprb/const_khinchin.c

@@ -0,0 +1,124 @@
+/*=============================================================================
+
+    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 <math.h>
+#include "fmprb.h"
+
+void
+fmprb_const_khinchin_eval_param(fmprb_t s, ulong N, ulong M, long prec)
+{
+    fmprb_t t, u, h;
+    fmprb_struct * pows;
+    long k, n;
+
+    fmprb_init(t);
+    fmprb_init(u);
+    fmprb_init(h);
+
+    if (N < 2) abort();
+    if (M < 0) abort();
+
+    pows = _fmprb_vec_init(N - 2);
+
+    /* sum of logarithms */
+    fmprb_zero(s);
+    for (k = 2; k < N; k++)
+    {
+        fmprb_set_ui(t, k - 1);
+        fmprb_div_ui(t, t, k, prec);
+        fmprb_log(t, t, prec);
+
+        fmprb_set_ui(u, k + 1);
+        fmprb_div_ui(u, u, k, prec);
+        fmprb_log(u, u, prec);
+
+        fmprb_mul(t, t, u, prec);
+        fmprb_sub(s, s, t, prec);
+    }
+
+    /* alternating harmonic numbers */
+    fmprb_one(h);
+
+    /* powers */
+    for (k = 0; k < N - 2; k++)
+        fmprb_one(pows + k);
+
+    /* sum of zetas */
+    for (n = 1; n <= M; n++)
+    {
+        /* zeta(2n,N) / n */
+        fmprb_zeta_ui(t, 2 * n, prec);
+        fmprb_sub_ui(t, t, 1, prec);
+        for (k = 0; k < N - 2; k++)
+        {
+            fmprb_div_ui(pows + k, pows + k, (k + 2) * (k + 2), prec);
+            fmprb_sub(t, t, pows + k, prec);
+        }
+        fmprb_div_ui(t, t, n, prec);
+
+        fmprb_mul(t, t, h, prec);
+        fmprb_add(s, s, t, prec);
+
+        /* forward h by two */
+        fmprb_set_ui(u, 2 * n);
+        fmprb_mul_ui(u, u, 2 * n + 1, prec);
+        fmprb_ui_div(u, 1, u, prec);
+        fmprb_sub(h, h, u, prec);
+    }
+
+    /* error bound 1/N^(2M) */
+    fmprb_set_ui(t, N);
+    fmprb_pow_ui(t, t, 2 * M, FMPRB_RAD_PREC);
+    fmprb_ui_div(t, 1, t, FMPRB_RAD_PREC);
+    fmprb_add_error(s, t);
+
+    fmprb_log_ui(t, 2, prec);
+    fmprb_div(s, s, t, prec);
+    fmprb_exp(s, s, prec);
+
+    _fmprb_vec_clear(pows, N - 2);
+
+    fmprb_clear(t);
+    fmprb_clear(u);
+    fmprb_clear(h);
+}
+
+void
+fmprb_const_khinchin_eval(fmprb_t K, long prec)
+{
+    ulong N, M;
+
+    prec += 10 + 2 * FLINT_BIT_COUNT(prec);
+
+    /* heuristic */
+    N = pow(prec, 0.35);
+
+    M = (prec * 0.69314718055994530942) / (2 * log(N));
+
+    fmprb_const_khinchin_eval_param(K, N, M, prec);
+}
+
+DEF_CACHED_CONSTANT(fmprb_const_khinchin, fmprb_const_khinchin_eval)
+

fmprb/test/t-const_khinchin.c

@@ -0,0 +1,77 @@
+/*=============================================================================
+
+    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"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("const_khinchin....");
+    fflush(stdout);
+    flint_randinit(state);
+
+    for (iter = 0; iter < 250; iter++)
+    {
+        fmprb_t r, s;
+        long accuracy, prec;
+
+        prec = 2 + n_randint(state, 1 << n_randint(state, 14));
+
+        fmprb_init(r);
+        fmprb_init(s);
+
+        fmprb_const_khinchin(r, prec);
+        fmprb_const_khinchin(s, prec + 1000);
+
+        if (!fmprb_overlaps(r, s))
+        {
+            printf("FAIL: containment\n\n");
+            printf("prec = %ld\n", prec);
+            printf("r = "); fmprb_printd(r, prec / 3.33); printf("\n\n");
+            abort();
+        }
+
+        accuracy = fmprb_rel_accuracy_bits(r);
+
+        if (accuracy < prec - 4)
+        {
+            printf("FAIL: poor accuracy\n\n");
+            printf("prec = %ld\n", prec);
+            printf("r = "); fmprb_printd(r, prec / 3.33); printf("\n\n");
+            abort();
+        }
+
+        fmprb_clear(r);
+        fmprb_clear(s);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+