include/flint/bernoulli.h

0001 /*
0002     Copyright (C) 2012 Fredrik Johansson
0003
0004     This file is part of FLINT.
0005
0006     FLINT is free software: you can redistribute it and/or modify it under
0007     the terms of the GNU Lesser General Public License (LGPL) as published
0008     by the Free Software Foundation; either version 3 of the License, or
0009     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
0010 */
0011
0012 #ifndef BERNOULLI_H
0013 #define BERNOULLI_H
0014
0015 #include <math.h>
0016 #include "arith.h"
0017 #include "arb_types.h"
0018
0019 #ifdef __cplusplus
0020 extern "C" {
0021 #endif
0022
0023 extern slong FLINT_TLS_PREFIX bernoulli_cache_num;
0024
0025 extern FLINT_TLS_PREFIX fmpq * bernoulli_cache;
0026
0027 void bernoulli_cache_compute(slong n);
0028
0029 /*
0030 Crude bound for the bits in d(n) = denom(B_n).
0031 By von Staudt-Clausen, d(n) = prod_{p-1 | n} p
0032     <= prod_{k | n} 2k
0033     <= n^{sigma_0(n)}.
0034
0035 We get a more accurate estimate taking the square root of this.
0036 Further, at least for sufficiently large n,
0037 sigma_0(n) < exp(1.066 log(n) / log(log(n))).
0038 */
0039 FLINT_FORCE_INLINE slong bernoulli_denom_size(slong n)
0040 {
0041     return 0.5 * 1.4427 * log(n) * pow(n, 1.066 / log(log(n)));
0042 }
0043
0044 FLINT_FORCE_INLINE slong bernoulli_zeta_terms(ulong s, slong prec)
0045 {
0046     slong N;
0047     N = pow(2.0, (prec + 1.0) / (s - 1.0));
0048     N += ((N % 2) == 0);
0049     return N;
0050 }
0051
0052 FLINT_FORCE_INLINE slong bernoulli_power_prec(slong i, ulong s1, slong wp)
0053 {
0054     slong p = wp - s1 * log(i) * 1.44269504088896341;
0055     return FLINT_MAX(p, 10);
0056 }
0057
0058 /* we should technically add O(log(n)) guard bits, but this is unnecessary
0059    in practice since the denominator estimate is quite a bit larger
0060    than the true denominators
0061  */
0062 FLINT_FORCE_INLINE slong bernoulli_global_prec(ulong nmax)
0063 {
0064     return arith_bernoulli_number_size(nmax) + bernoulli_denom_size(nmax);
0065 }
0066
0067
0068 /* avoid potential numerical problems for very small n */
0069 #define BERNOULLI_REV_MIN 32
0070
0071 typedef struct
0072 {
0073     slong alloc;
0074     slong prec;
0075     slong max_power;
0076     fmpz * powers;
0077     fmpz_t pow_error;
0078     arb_t prefactor;
0079     arb_t two_pi_squared;
0080     ulong n;
0081 }
0082 bernoulli_rev_struct;
0083
0084 typedef bernoulli_rev_struct bernoulli_rev_t[1];
0085
0086 void bernoulli_rev_init(bernoulli_rev_t iter, ulong nmax);
0087
0088 void bernoulli_rev_next(fmpz_t numer, fmpz_t denom, bernoulli_rev_t iter);
0089
0090 void bernoulli_rev_clear(bernoulli_rev_t iter);
0091
0092 void bernoulli_fmpq_vec_no_cache(fmpq * res, ulong a, slong num);
0093
0094 #define BERNOULLI_ENSURE_CACHED(n) \
0095   do { \
0096     slong __n = (n); \
0097     if (__n >= bernoulli_cache_num) \
0098         bernoulli_cache_compute(__n + 1); \
0099   } while (0); \
0100
0101 slong bernoulli_bound_2exp_si(ulong n);
0102
0103 ulong bernoulli_mod_p_harvey(ulong k, ulong p);
0104
0105 void _bernoulli_fmpq_ui_multi_mod(fmpz_t num, fmpz_t den, ulong n, double alpha);
0106 void _bernoulli_fmpq_ui_zeta(fmpz_t num, fmpz_t den, ulong n);
0107 void _bernoulli_fmpq_ui(fmpz_t num, fmpz_t den, ulong n);
0108 void bernoulli_fmpq_ui(fmpq_t b, ulong n);
0109
0110 #ifdef __cplusplus
0111 }
0112 #endif
0113
0114 #endif