1 /* bernoulli -- internal function to compute Bernoulli numbers.
2
3 Copyright 2005-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #include "mpfr-impl.h"
24
25 /* assume p >= 5 and is odd */
26 static int
is_prime(unsigned long p)27 is_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. */
43 static void
mpfr_bernoulli_internal(mpz_t * b,unsigned long n)44 mpfr_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
204 static MPFR_THREAD_ATTR mpz_t *bernoulli_table = NULL;
205 static MPFR_THREAD_ATTR unsigned long bernoulli_size = 0;
206 static MPFR_THREAD_ATTR unsigned long bernoulli_alloc = 0;
207
208 mpz_srcptr
mpfr_bernoulli_cache(unsigned long n)209 mpfr_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
239 void
mpfr_bernoulli_freecache(void)240 mpfr_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