diff --git a/doc/source/credits.rst b/doc/source/credits.rst index a822c5d1..fad19a74 100644 --- a/doc/source/credits.rst +++ b/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 diff --git a/doc/source/fmprb.rst b/doc/source/fmprb.rst index d7a8c0ba..e35b08a3 100644 --- a/doc/source/fmprb.rst +++ b/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 ............................................................................... diff --git a/doc/source/fmprb_poly.rst b/doc/source/fmprb_poly.rst index 4a346777..345e40f8 100644 --- a/doc/source/fmprb_poly.rst +++ b/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. diff --git a/fmprb.h b/fmprb.h index 88be9822..ded705cc 100644 --- a/fmprb.h +++ b/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); diff --git a/fmprb/const_khinchin.c b/fmprb/const_khinchin.c new file mode 100644 index 00000000..09040ec9 --- /dev/null +++ b/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 +#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) + diff --git a/fmprb/test/t-const_khinchin.c b/fmprb/test/t-const_khinchin.c new file mode 100644 index 00000000..32f5b535 --- /dev/null +++ b/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; +} +