w32tex
About: TeX Live provides a comprehensive TeX system including all the major TeX-related programs, macro packages, and fonts that are free software. Windows sources.
  Fossies Dox: w32tex-src.tar.xz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

bernoulli.c
Go to the documentation of this file.
1 /* bernoulli -- internal function to compute Bernoulli numbers.
2 
3 Copyright 2005-2020 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
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
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);
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  }
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 
205 static MPFR_THREAD_ATTR unsigned long bernoulli_size = 0;
206 static MPFR_THREAD_ATTR unsigned long bernoulli_alloc = 0;
207 
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 *)
220  bernoulli_size = 0;
221  }
222  else if (n >= bernoulli_alloc)
223  {
226  (n + n/4) * sizeof (mpz_t));
227  bernoulli_alloc = n + n/4;
228  }
231  for (i = bernoulli_size; i <= n; i++)
233  bernoulli_size = n+1;
234  }
236  return bernoulli_table[n];
237 }
238 
239 void
241 {
242  unsigned long i;
243 
244  if (bernoulli_table != NULL)
245  {
246  for (i = 0; i < bernoulli_size; i++)
247  {
249  }
252  bernoulli_alloc = 0;
253  bernoulli_size = 0;
254  }
255 }
q
Definition: afm2pl.c:2287
static unsigned long bernoulli_size
Definition: bernoulli.c:205
mpz_srcptr mpfr_bernoulli_cache(unsigned long n)
Definition: bernoulli.c:209
static unsigned long bernoulli_alloc
Definition: bernoulli.c:206
static void mpfr_bernoulli_internal(mpz_t *b, unsigned long n)
Definition: bernoulli.c:44
void mpfr_bernoulli_freecache(void)
Definition: bernoulli.c:240
static int is_prime(unsigned long p)
Definition: bernoulli.c:27
static mpz_t * bernoulli_table
Definition: bernoulli.c:204
#define n
Definition: t4ht.c:1290
#define b
Definition: jpegint.h:372
int mpfr_div_2ui(mpfr_ptr y, mpfr_srcptr x, unsigned long n, mpfr_rnd_t rnd_mode)
Definition: div_2ui.c:26
int z
Definition: dviconv.c:26
#define s
Definition: afcover.h:80
#define t
Definition: afcover.h:96
int mpfr_get_z(mpz_ptr z, mpfr_srcptr f, mpfr_rnd_t rnd)
Definition: get_z.c:27
#define mpz_init
Definition: gmp.h:851
#define mpz_neg
Definition: gmp.h:942
#define mpz_clear
Definition: gmp.h:679
#define mpz_add
Definition: gmp.h:628
#define mpz_cmp_ui(Z, UI)
Definition: gmp.h:2268
#define mpz_cdiv_q_2exp
Definition: gmp.h:655
#define mpz_set_ui
Definition: gmp.h:1035
#define mpz_fdiv_q
Definition: gmp.h:757
__mpz_struct mpz_t[1]
Definition: gmp.h:164
#define mpz_init_set_ui
Definition: gmp.h:872
#define mpn_scan1
Definition: gmp.h:1591
#define mpz_mul
Definition: gmp.h:930
#define mpz_mul_2exp
Definition: gmp.h:933
#define mpz_cdiv_q
Definition: gmp.h:652
#define GMP_NUMB_BITS
Definition: gmp.h:46
#define mpz_sizeinbase
Definition: gmp.h:1046
#define mpz_ui_pow_ui
Definition: gmp.h:1103
#define mpz_fac_ui
Definition: gmp.h:745
#define mpz_mul_ui
Definition: gmp.h:939
#define mpz_divexact
Definition: gmp.h:724
#define MAX(a, b)
Definition: jpegint.h:267
#define NULL
Definition: ftobjs.h:61
small capitals from c petite p scientific f u
Definition: afcover.h:88
small capitals from c petite p
Definition: afcover.h:72
small capitals from c petite p scientific i
Definition: afcover.h:80
kerning y
Definition: ttdriver.c:212
int den
Definition: dvi.c:19
int num
Definition: disdvi.c:621
int mpfr_log2(mpfr_ptr r, mpfr_srcptr a, mpfr_rnd_t rnd_mode)
Definition: log2.c:30
void mpfr_free_func(void *ptr, size_t size)
Definition: mpfr-gmp.c:330
void * mpfr_reallocate_func(void *ptr, size_t old_size, size_t new_size)
Definition: mpfr-gmp.c:319
void * mpfr_allocate_func(size_t alloc_size)
Definition: mpfr-gmp.c:308
long __gmpfr_ceil_log2(double)
Definition: uceil_log2.c:30
#define MPFR_INC_PREC(P, X)
Definition: mpfr-impl.h:2047
#define MPFR_LIMB_SIZE(x)
Definition: mpfr-impl.h:963
#define mpfr_const_pi(_d, _r)
Definition: mpfr-impl.h:1384
#define MPFR_ASSERTD(expr)
Definition: mpfr-impl.h:516
#define MPFR_EXP(x)
Definition: mpfr-impl.h:959
#define MPFR_MANT(x)
Definition: mpfr-impl.h:960
#define MPFR_THREAD_ATTR
Definition: mpfr-thread.h:44
void mpfr_clear(mpfr_ptr m)
Definition: clear.c:26
int mpfr_div(mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
Definition: div.c:748
unsigned long mpfr_get_ui(mpfr_srcptr f, mpfr_rnd_t rnd)
Definition: get_ui.c:26
void mpfr_init2(mpfr_ptr x, mpfr_prec_t p)
Definition: init2.c:26
int mpfr_pow_ui(mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t)
long mpfr_prec_t
Definition: mpfr.h:168
__mpfr_struct mpfr_t[1]
Definition: mpfr.h:250
int mpfr_mul_2ui(mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t)
int mpfr_set_ui_2exp(mpfr_ptr, unsigned long, mpfr_exp_t, mpfr_rnd_t)
Definition: set_ui_2exp.c:28
@ MPFR_RNDU
Definition: mpfr.h:105
@ MPFR_RNDZ
Definition: mpfr.h:104
int mpfr_set_z(mpfr_ptr, mpz_srcptr, mpfr_rnd_t)
Definition: set_z.c:27
int mpfr_mul_ui(mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t)
#define Prec
Definition: sqrtrem.c:168
@ err
Definition: mtxline.h:24
Definition: dvips.h:235