diff --git a/doc/source/fmprb.rst b/doc/source/fmprb.rst index 4dbea47d..4d32d92e 100644 --- a/doc/source/fmprb.rst +++ b/doc/source/fmprb.rst @@ -642,6 +642,13 @@ Constants \log 2 = \frac{3}{4} \sum_{k=0}^{\infty} \frac{(-1)^k (k!)^2}{2^k (2k+1)!} +.. function:: void fmprb_const_log10(fmprb_t s, long prec) + + Sets *x* to `\log 10`. The value is cached for repeated use. + Uses the generic hypergeometric series code to evaluate the + Machin-like formula + `\log 10 = 46 \operatorname{atanh}(1/31) + 34 \operatorname{atanh}(1/49) + 20 \operatorname{atanh}(1/161)`. + .. function:: void fmprb_const_e(fmprb_t s, long prec) Sets *x* to Euler's number `e = \sum_{n=0}^{\infty} 1/n!`, evaluated diff --git a/fmprb.h b/fmprb.h index 3fe35a2e..1c024cf0 100644 --- a/fmprb.h +++ b/fmprb.h @@ -307,6 +307,7 @@ void fmprb_const_pi(fmprb_t x, long prec); void fmprb_const_sqrt_pi(fmprb_t t, long prec); void fmprb_const_log_sqrt2pi(fmprb_t t, long prec); void fmprb_const_log2(fmprb_t s, long prec); +void fmprb_const_log10(fmprb_t s, long prec); void fmprb_const_catalan(fmprb_t s, long prec); void fmprb_const_e(fmprb_t s, long prec); diff --git a/fmprb/const_log10.c b/fmprb/const_log10.c new file mode 100644 index 00000000..69b3ff77 --- /dev/null +++ b/fmprb/const_log10.c @@ -0,0 +1,70 @@ +/*============================================================================= + + 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 "fmprb.h" +#include "hypgeom.h" + +static void +atanh_bsplit(fmprb_t s, ulong c, long a, long prec) +{ + fmprb_t t; + hypgeom_t series; + hypgeom_init(series); + fmprb_init(t); + + fmpz_poly_set_ui(series->A, 1); + fmpz_poly_set_coeff_ui(series->B, 0, 1); + fmpz_poly_set_coeff_ui(series->B, 1, 2); + fmpz_poly_set_ui(series->P, 1); + fmpz_poly_set_ui(series->Q, c * c); + + fmprb_hypgeom_infsum(s, t, series, prec, prec); + fmprb_mul_si(s, s, a, prec); + fmprb_mul_ui(t, t, c, prec); + fmprb_div(s, s, t, prec); + + fmprb_clear(t); + hypgeom_clear(series); +} + +void +fmprb_const_log10_eval(fmprb_t s, long prec) +{ + fmprb_t t; + fmprb_init(t); + + prec += FLINT_CLOG2(prec); + + atanh_bsplit(s, 31, 46, prec); + atanh_bsplit(t, 49, 34, prec); + fmprb_add(s, s, t, prec); + atanh_bsplit(t, 161, 20, prec); + fmprb_add(s, s, t, prec); + + fmprb_clear(t); +} + +DEF_CACHED_CONSTANT(fmprb_const_log10, fmprb_const_log10_eval) + diff --git a/fmprb/test/t-const_log10.c b/fmprb/test/t-const_log10.c new file mode 100644 index 00000000..38ec4139 --- /dev/null +++ b/fmprb/test/t-const_log10.c @@ -0,0 +1,79 @@ +/*============================================================================= + + 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_log10...."); + fflush(stdout); + flint_randinit(state); + + for (iter = 0; iter < 250; iter++) + { + fmprb_t r; + mpfr_t s; + long accuracy, prec; + + prec = 2 + n_randint(state, 1 << n_randint(state, 16)); + + fmprb_init(r); + mpfr_init2(s, prec + 1000); + + fmprb_const_log10(r, prec); + mpfr_set_ui(s, 10, MPFR_RNDN); + mpfr_log(s, s, MPFR_RNDN); + + if (!fmprb_contains_mpfr(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); + mpfr_clear(s); + } + + flint_randclear(state); + _fmpz_cleanup(); + mpfr_free_cache(); + printf("PASS\n"); + return EXIT_SUCCESS; +}