bernoulli.c revision 1.1.1.5
1/* bernoulli -- internal function to compute Bernoulli numbers. 2 3Copyright 2005-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include "mpfr-impl.h" 24 25/* assume p >= 5 and is odd */ 26static int 27is_prime (unsigned long p) 28{ 29 unsigned long q; 30 31 MPFR_ASSERTD (p >= 5 && (p & 1) != 0); 32 for (q = 3; q * q <= p; q += 2) 33 if ((p % q) == 0) 34 return 0; 35 return 1; 36} 37 38/* Computes and stores B[2n]*(2n+1)! in b[n] 39 using Von Staudt���Clausen theorem, which says that the denominator of B[n] 40 divides the product of all primes p such that p-1 divides n. 41 Since B[n] = zeta(n) * 2*n!/(2pi)^n, we compute an approximation of 42 (2n+1)! * zeta(n) * 2*n!/(2pi)^n and round it to the nearest integer. */ 43static void 44mpfr_bernoulli_internal (mpz_t *b, unsigned long n) 45{ 46 unsigned long p, err, zn; 47 mpz_t s, t, u, den; 48 mpz_ptr num; 49 mpfr_t y, z; 50 int ok; 51 /* Prec[n/2] is minimal precision so that result is correct for B[n] */ 52 mpfr_prec_t prec; 53 mpfr_prec_t Prec[] = {0, 5, 5, 6, 6, 9, 16, 10, 19, 23, 25, 27, 35, 31, 54 42, 51, 51, 50, 73, 60, 76, 79, 83, 87, 101, 97, 55 108, 113, 119, 125, 149, 133, 146}; 56 57 mpz_init (b[n]); 58 59 if (n == 0) 60 { 61 mpz_set_ui (b[0], 1); 62 return; 63 } 64 65 /* compute denominator */ 66 num = b[n]; 67 n = 2 * n; 68 mpz_init_set_ui (den, 6); 69 for (p = 5; p <= n+1; p += 2) 70 { 71 if ((n % (p-1)) == 0 && is_prime (p)) 72 mpz_mul_ui (den, den, p); 73 } 74 if (n <= 64) 75 prec = Prec[n >> 1]; 76 else 77 { 78 /* evaluate the needed precision: zeta(n)*2*den*n!/(2*pi)^n <= 79 3.3*den*(n/e/2/pi)^n*sqrt(2*pi*n) */ 80 prec = __gmpfr_ceil_log2 (7.0 * (double) n); /* bound 2*pi by 7 */ 81 prec = (prec + 1) >> 1; /* sqrt(2*pi*n) <= 2^prec */ 82 mpfr_init2 (z, 53); 83 mpfr_set_ui_2exp (z, 251469612, -32, MPFR_RNDU); /* 1/e/2/pi <= z */ 84 mpfr_mul_ui (z, z, n, MPFR_RNDU); 85 mpfr_log2 (z, z, MPFR_RNDU); 86 mpfr_mul_ui (z, z, n, MPFR_RNDU); 87 p = mpfr_get_ui (z, MPFR_RNDU); /* (n/e/2/pi)^n <= 2^p */ 88 mpfr_clear (z); 89 MPFR_INC_PREC (prec, p + mpz_sizeinbase (den, 2)); 90 /* the +2 term ensures no rounding failure up to n=10000 */ 91 MPFR_INC_PREC (prec, __gmpfr_ceil_log2 (prec) + 2); 92 } 93 94 try_again: 95 mpz_init (s); 96 mpz_init (t); 97 mpz_init (u); 98 mpz_set_ui (u, 1); 99 mpz_mul_2exp (u, u, prec); /* u = 2^prec */ 100 mpz_ui_pow_ui (t, 3, n); 101 mpz_fdiv_q (s, u, t); /* multiply all terms by 2^prec */ 102 /* we compute a lower bound of the series, thus the final result cannot 103 be too large */ 104 for (p = 4; mpz_cmp_ui (t, 0) > 0; p++) 105 { 106 mpz_ui_pow_ui (t, p, n); 107 mpz_fdiv_q (t, u, t); 108 /* 2^prec/p^n-1 < t <= 2^prec/p^n */ 109 mpz_add (s, s, t); 110 } 111 /* sum(2^prec/q^n-1, q=3..p) < t <= sum(2^prec/q^n, q=3..p) 112 thus the error on the truncated series is at most p-2. 113 The neglected part of the series is R = sum(1/x^n, x=p+1..infinity) 114 with int(1/x^n, x=p+1..infinity) <= R <= int(1/x^n, x=p..infinity) 115 thus 1/(n-1)/(p+1)^(n-1) <= R <= 1/(n-1)/p^(n-1). The difference between 116 the lower and upper bound is bounded by p^(-n), which is bounded by 117 2^(-prec) since t=0 in the above loop */ 118 mpz_ui_pow_ui (t, p, n - 1); 119 mpz_mul_ui (t, t, n - 1); 120 mpz_cdiv_q (t, u, t); 121 mpz_add (s, s, t); 122 /* now 2^prec * (zeta(n)-1-1/2^n) - p < s <= 2^prec * (zeta(n)-1-1/2^n) */ 123 /* add 1 which is 2^prec */ 124 mpz_add (s, s, u); 125 /* add 1/2^n which is 2^(prec-n) */ 126 mpz_cdiv_q_2exp (u, u, n); 127 mpz_add (s, s, u); 128 /* now 2^prec * zeta(n) - p < s <= 2^prec * zeta(n) */ 129 /* multiply by n! */ 130 mpz_fac_ui (t, n); 131 mpz_mul (s, s, t); 132 /* multiply by 2*den */ 133 mpz_mul (s, s, den); 134 mpz_mul_2exp (s, s, 1); 135 /* now convert to mpfr */ 136 mpfr_init2 (z, prec); 137 mpfr_set_z (z, s, MPFR_RNDZ); 138 /* now (2^prec * zeta(n) - p) * 2*den*n! - ulp(z) < z <= 139 2^prec * zeta(n) * 2*den*n!. 140 Since z <= 2^prec * zeta(n) * 2*den*n!, 141 ulp(z) <= 2*zeta(n) * 2*den*n!, thus 142 (2^prec * zeta(n)-(p+1)) * 2*den*n! < z <= 2^prec * zeta(n) * 2*den*n! */ 143 mpfr_div_2ui (z, z, prec, MPFR_RNDZ); 144 /* now (zeta(n) - (p+1)/2^prec) * 2*den*n! < z <= zeta(n) * 2*den*n! */ 145 /* divide by (2pi)^n */ 146 mpfr_init2 (y, prec); 147 mpfr_const_pi (y, MPFR_RNDU); 148 /* pi <= y <= pi * (1 + 2^(1-prec)) */ 149 mpfr_mul_2ui (y, y, 1, MPFR_RNDU); 150 /* 2pi <= y <= 2pi * (1 + 2^(1-prec)) */ 151 mpfr_pow_ui (y, y, n, MPFR_RNDU); 152 /* (2pi)^n <= y <= (2pi)^n * (1 + 2^(1-prec))^(n+1) */ 153 mpfr_div (z, z, y, MPFR_RNDZ); 154 /* now (zeta(n) - (p+1)/2^prec) * 2*den*n! / (2pi)^n / (1+2^(1-prec))^(n+1) 155 <= z <= zeta(n) * 2*den*n! / (2pi)^n, and since zeta(n) >= 1: 156 den * B[n] * (1 - (p+1)/2^prec) / (1+2^(1-prec))^(n+1) 157 <= z <= den * B[n] 158 Since 1 / (1+2^(1-prec))^(n+1) >= (1 - 2^(1-prec))^(n+1) >= 159 1 - (n+1) * 2^(1-prec): 160 den * B[n] / (2pi)^n * (1 - (p+1)/2^prec) * (1-(n+1)*2^(1-prec)) 161 <= z <= den * B[n] thus 162 den * B[n] * (1 - (2n+p+3)/2^prec) <= z <= den * B[n] */ 163 164 /* the error is bounded by 2^(EXP(z)-prec) * (2n+p+3) */ 165 for (err = 0, p = 2 * n + p + 3; p > 1; err++, p = (p + 1) >> 1); 166 zn = MPFR_LIMB_SIZE(z) * GMP_NUMB_BITS; /* total number of bits of z */ 167 if (err >= prec) 168 ok = 0; 169 else 170 { 171 err = prec - err; 172 /* now the absolute error is bounded by 2^(EXP(z) - err): 173 den * B[n] - 2^(EXP(z) - err) <= z <= den * B[n] 174 thus if subtracting 2^(EXP(z) - err) does not change the rounding 175 (up) we are ok */ 176 err = mpn_scan1 (MPFR_MANT(z), zn - err); 177 /* weight of this 1 bit is 2^(EXP(z) - zn + err) */ 178 ok = MPFR_EXP(z) < zn - err; 179 } 180 mpfr_get_z (num, z, MPFR_RNDU); 181 if ((n & 2) == 0) 182 mpz_neg (num, num); 183 184 /* multiply by (n+1)! */ 185 mpz_mul_ui (t, t, n + 1); 186 mpz_divexact (t, t, den); /* t was still n! */ 187 mpz_mul (num, num, t); 188 189 mpfr_clear (y); 190 mpfr_clear (z); 191 mpz_clear (s); 192 mpz_clear (t); 193 mpz_clear (u); 194 195 if (!ok) 196 { 197 MPFR_INC_PREC (prec, prec / 10); 198 goto try_again; 199 } 200 201 mpz_clear (den); 202} 203 204static MPFR_THREAD_ATTR mpz_t *bernoulli_table = NULL; 205static MPFR_THREAD_ATTR unsigned long bernoulli_size = 0; 206static MPFR_THREAD_ATTR unsigned long bernoulli_alloc = 0; 207 208mpz_srcptr 209mpfr_bernoulli_cache (unsigned long n) 210{ 211 unsigned long i; 212 213 if (n >= bernoulli_size) 214 { 215 if (bernoulli_alloc == 0) 216 { 217 bernoulli_alloc = MAX(16, n + n/4); 218 bernoulli_table = (mpz_t *) 219 mpfr_allocate_func (bernoulli_alloc * sizeof (mpz_t)); 220 bernoulli_size = 0; 221 } 222 else if (n >= bernoulli_alloc) 223 { 224 bernoulli_table = (mpz_t *) mpfr_reallocate_func 225 (bernoulli_table, bernoulli_alloc * sizeof (mpz_t), 226 (n + n/4) * sizeof (mpz_t)); 227 bernoulli_alloc = n + n/4; 228 } 229 MPFR_ASSERTD (bernoulli_alloc > n); 230 MPFR_ASSERTD (bernoulli_size >= 0); 231 for (i = bernoulli_size; i <= n; i++) 232 mpfr_bernoulli_internal (bernoulli_table, i); 233 bernoulli_size = n+1; 234 } 235 MPFR_ASSERTD (bernoulli_size > n); 236 return bernoulli_table[n]; 237} 238 239void 240mpfr_bernoulli_freecache (void) 241{ 242 unsigned long i; 243 244 if (bernoulli_table != NULL) 245 { 246 for (i = 0; i < bernoulli_size; i++) 247 { 248 mpz_clear (bernoulli_table[i]); 249 } 250 mpfr_free_func (bernoulli_table, bernoulli_alloc * sizeof (mpz_t)); 251 bernoulli_table = NULL; 252 bernoulli_alloc = 0; 253 bernoulli_size = 0; 254 } 255} 256