diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/bernoulli.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/bernoulli.c | 265 |
1 files changed, 220 insertions, 45 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/bernoulli.c b/Build/source/libs/mpfr/mpfr-src/src/bernoulli.c index 26ecf25bd7c..399d31c729d 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/bernoulli.c +++ b/Build/source/libs/mpfr/mpfr-src/src/bernoulli.c @@ -22,59 +22,234 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #include "mpfr-impl.h" -/* assuming b[0]...b[2(n-1)] are computed, computes and stores B[2n]*(2n+1)! - - t/(exp(t)-1) = sum(B[j]*t^j/j!, j=0..infinity) - thus t = (exp(t)-1) * sum(B[j]*t^j/j!, n=0..infinity). - Taking the coefficient of degree n+1 > 1, we get: - 0 = sum(1/(n+1-k)!*B[k]/k!, k=0..n) - which gives: - B[n] = -sum(binomial(n+1,k)*B[k], k=0..n-1)/(n+1). - - Let C[n] = B[n]*(n+1)!. - Then C[n] = -sum(binomial(n+1,k)*C[k]*n!/(k+1)!, k=0..n-1), - which proves that the C[n] are integers. -*/ -mpz_t* +/* assume p >= 5 and is odd */ +static int +is_prime (unsigned long p) +{ + unsigned long q; + + MPFR_ASSERTD (p >= 5 && (p & 1) != 0); + for (q = 3; q * q <= p; q += 2) + if ((p % q) == 0) + return 0; + return 1; +} + +/* Computes and stores B[2n]*(2n+1)! in b[n] + using Von Staudt–Clausen theorem, which says that the denominator of B[n] + divides the product of all primes p such that p-1 divides n. + Since B[n] = zeta(n) * 2*n!/(2pi)^n, we compute an approximation of + d * zeta(n) * 2*n!/(2pi)^n and round it to the nearest integer. */ +static void mpfr_bernoulli_internal (mpz_t *b, unsigned long n) { + unsigned long p, err, zn; + mpz_t s, t, u, den; + mpz_ptr num; + mpfr_t y, z; + int ok; + /* Prec[n/2] is minimal precision so that result is correct for B[n] */ + mpfr_prec_t prec; + mpfr_prec_t Prec[] = {0, 5, 5, 6, 6, 9, 16, 10, 19, 23, 25, 27, 35, 31, + 42, 51, 51, 50, 73, 60, 76, 79, 83, 87, 101, 97, + 108, 113, 119, 125, 149, 133, 146}; + + mpz_init (b[n]); + if (n == 0) { - b = (mpz_t *) (*__gmp_allocate_func) (sizeof (mpz_t)); - mpz_init_set_ui (b[0], 1); + mpz_set_ui (b[0], 1); + return; + } + + /* compute denominator */ + num = b[n]; + n = 2 * n; + mpz_init_set_ui (den, 6); + for (p = 5; p <= n+1; p += 2) + { + if ((n % (p-1)) == 0 && is_prime (p)) + mpz_mul_ui (den, den, p); + } + if (n <= 64) + prec = Prec[n >> 1]; + else + { + /* evaluate the needed precision: zeta(n)*2*den*n!/(2*pi)^n <= + 3.3*den*(n/e/2/pi)^n*sqrt(2*pi*n) */ + prec = __gmpfr_ceil_log2 (7.0 * (double) n); /* bound 2*pi by 7 */ + prec = (prec + 1) >> 1; /* sqrt(2*pi*n) <= 2^prec */ + mpfr_init2 (z, 53); + mpfr_set_ui_2exp (z, 251469612, -32, MPFR_RNDU); /* 1/e/2/pi <= z */ + mpfr_mul_ui (z, z, n, MPFR_RNDU); + mpfr_log2 (z, z, MPFR_RNDU); + mpfr_mul_ui (z, z, n, MPFR_RNDU); + p = mpfr_get_ui (z, MPFR_RNDU); /* (n/e/2/pi)^n <= 2^p */ + mpfr_clear (z); + /* the +14 term ensures no rounding failure up to n=10000 */ + prec += p + mpz_sizeinbase (den, 2) + 14; + } + + try_again: + mpz_init (s); + mpz_init (t); + mpz_init (u); + mpz_set_ui (u, 1); + mpz_mul_2exp (u, u, prec); /* u = 2^prec */ + mpz_ui_pow_ui (t, 3, n); + mpz_fdiv_q (s, u, t); /* multiply all terms by 2^prec */ + /* we compute a lower bound of the series, thus the final result cannot + be too large */ + for (p = 4; mpz_cmp_ui (t, 0) > 0; p++) + { + mpz_ui_pow_ui (t, p, n); + mpz_fdiv_q (t, u, t); + /* 2^prec/p^n-1 < t <= 2^prec/p^n */ + mpz_add (s, s, t); } + /* sum(2^prec/q^n-1, q=3..p) < t <= sum(2^prec/q^n, q=3..p) + thus the error on the truncated series is at most p-2. + The neglected part of the series is R = sum(1/x^n, x=p+1..infinity) + with int(1/x^n, x=p+1..infinity) <= R <= int(1/x^n, x=p..infinity) + thus 1/(n-1)/(p+1)^(n-1) <= R <= 1/(n-1)/p^(n-1). The difference between + the lower and upper bound is bounded by p^(-n), which is bounded by + 2^(-prec) since t=0 in the above loop */ + mpz_ui_pow_ui (t, p, n - 1); + mpz_mul_ui (t, t, n - 1); + mpz_cdiv_q (t, u, t); + mpz_add (s, s, t); + /* now 2^prec * (zeta(n)-1-1/2^n) - p < s <= 2^prec * (zeta(n)-1-1/2^n) */ + /* add 1 which is 2^prec */ + mpz_add (s, s, u); + /* add 1/2^n which is 2^(prec-n) */ + mpz_cdiv_q_2exp (u, u, n); + mpz_add (s, s, u); + /* now 2^prec * zeta(n) - p < s <= 2^prec * zeta(n) */ + /* multiply by n! */ + mpz_fac_ui (t, n); + mpz_mul (s, s, t); + /* multiply by 2*den */ + mpz_mul (s, s, den); + mpz_mul_2exp (s, s, 1); + /* now convert to mpfr */ + mpfr_init2 (z, prec); + mpfr_set_z (z, s, MPFR_RNDZ); + /* now (2^prec * zeta(n) - p) * 2*den*n! - ulp(z) < z <= + 2^prec * zeta(n) * 2*den*n!. + Since z <= 2^prec * zeta(n) * 2*den*n!, + ulp(z) <= 2*zeta(n) * 2*den*n!, thus + (2^prec * zeta(n)-(p+1)) * 2*den*n! < z <= 2^prec * zeta(n) * 2*den*n! */ + mpfr_div_2exp (z, z, prec, MPFR_RNDZ); + /* now (zeta(n) - (p+1)/2^prec) * 2*den*n! < z <= zeta(n) * 2*den*n! */ + /* divide by (2pi)^n */ + mpfr_init2 (y, prec); + mpfr_const_pi (y, MPFR_RNDU); + /* pi <= y <= pi * (1 + 2^(1-prec)) */ + mpfr_mul_2exp (y, y, 1, MPFR_RNDU); + /* 2pi <= y <= 2pi * (1 + 2^(1-prec)) */ + mpfr_pow_ui (y, y, n, MPFR_RNDU); + /* (2pi)^n <= y <= (2pi)^n * (1 + 2^(1-prec))^(n+1) */ + mpfr_div (z, z, y, MPFR_RNDZ); + /* now (zeta(n) - (p+1)/2^prec) * 2*den*n! / (2pi)^n / (1+2^(1-prec))^(n+1) + <= z <= zeta(n) * 2*den*n! / (2pi)^n, and since zeta(n) >= 1: + den * B[n] * (1 - (p+1)/2^prec) / (1+2^(1-prec))^(n+1) + <= z <= den * B[n] + Since 1 / (1+2^(1-prec))^(n+1) >= (1 - 2^(1-prec))^(n+1) >= + 1 - (n+1) * 2^(1-prec): + den * B[n] / (2pi)^n * (1 - (p+1)/2^prec) * (1-(n+1)*2^(1-prec)) + <= z <= den * B[n] thus + den * B[n] * (1 - (2n+p+3)/2^prec) <= z <= den * B[n] */ + + /* the error is bounded by 2^(EXP(z)-prec) * (2n+p+3) */ + for (err = 0, p = 2 * n + p + 3; p > 1; err++, p = (p + 1) >> 1); + zn = MPFR_LIMB_SIZE(z) * GMP_NUMB_BITS; /* total number of bits of z */ + if (err >= prec) + ok = 0; else { - mpz_t t; - unsigned long k; - - b = (mpz_t *) (*__gmp_reallocate_func) - (b, n * sizeof (mpz_t), (n + 1) * sizeof (mpz_t)); - mpz_init (b[n]); - /* b[n] = -sum(binomial(2n+1,2k)*C[k]*(2n)!/(2k+1)!, k=0..n-1) */ - mpz_init_set_ui (t, 2 * n + 1); - mpz_mul_ui (t, t, 2 * n - 1); - mpz_mul_ui (t, t, 2 * n); - mpz_mul_ui (t, t, n); - mpz_fdiv_q_ui (t, t, 3); /* exact: t=binomial(2*n+1,2*k)*(2*n)!/(2*k+1)! - for k=n-1 */ - mpz_mul (b[n], t, b[n-1]); - for (k = n - 1; k-- > 0;) + err = prec - err; + /* now the absolute error is bounded by 2^(EXP(z) - err): + den * B[n] - 2^(EXP(z) - err) <= z <= den * B[n] + thus if subtracting 2^(EXP(z) - err) does not change the rounding + (up) we are ok */ + err = mpn_scan1 (MPFR_MANT(z), zn - err); + /* weight of this 1 bit is 2^(EXP(z) - zn + err) */ + ok = MPFR_EXP(z) < zn - err; + } + mpfr_get_z (num, z, MPFR_RNDU); + if ((n & 2) == 0) + mpz_neg (num, num); + + /* multiply by (n+1)! */ + mpz_mul_ui (t, t, n + 1); + mpz_divexact (t, t, den); /* t was still n! */ + mpz_mul (num, num, t); + mpz_set_ui (den, 1); + + mpfr_clear (y); + mpfr_clear (z); + mpz_clear (s); + mpz_clear (t); + mpz_clear (u); + + if (!ok) + { + prec += prec / 10; + goto try_again; + } + + mpz_clear (den); +} + +static MPFR_THREAD_ATTR mpz_t *bernoulli_table = NULL; +static MPFR_THREAD_ATTR unsigned long bernoulli_size = 0; +static MPFR_THREAD_ATTR unsigned long bernoulli_alloc = 0; + +mpz_srcptr +mpfr_bernoulli_cache (unsigned long n) +{ + unsigned long i; + + if (n >= bernoulli_size) + { + if (bernoulli_alloc == 0) + { + bernoulli_alloc = MAX(16, n + n/4); + bernoulli_table = (mpz_t *) + mpfr_allocate_func (bernoulli_alloc * sizeof (mpz_t)); + bernoulli_size = 0; + } + else if (n >= bernoulli_alloc) + { + bernoulli_table = (mpz_t *) mpfr_reallocate_func + (bernoulli_table, bernoulli_alloc * sizeof (mpz_t), + (n + n/4) * sizeof (mpz_t)); + bernoulli_alloc = n + n/4; + } + MPFR_ASSERTD (bernoulli_alloc > n); + MPFR_ASSERTD (bernoulli_size >= 0); + for (i = bernoulli_size; i <= n; i++) + mpfr_bernoulli_internal (bernoulli_table, i); + bernoulli_size = n+1; + } + MPFR_ASSERTD (bernoulli_size > n); + return bernoulli_table[n]; +} + +void +mpfr_bernoulli_freecache (void) +{ + unsigned long i; + + if (bernoulli_table != NULL) + { + for (i = 0; i < bernoulli_size; i++) { - mpz_mul_ui (t, t, 2 * k + 1); - mpz_mul_ui (t, t, 2 * k + 2); - mpz_mul_ui (t, t, 2 * k + 2); - mpz_mul_ui (t, t, 2 * k + 3); - mpz_fdiv_q_ui (t, t, 2 * (n - k) + 1); - mpz_fdiv_q_ui (t, t, 2 * (n - k)); - mpz_addmul (b[n], t, b[k]); + mpz_clear (bernoulli_table[i]); } - /* take into account C[1] */ - mpz_mul_ui (t, t, 2 * n + 1); - mpz_fdiv_q_2exp (t, t, 1); - mpz_sub (b[n], b[n], t); - mpz_neg (b[n], b[n]); - mpz_clear (t); + mpfr_free_func (bernoulli_table, bernoulli_alloc * sizeof (mpz_t)); + bernoulli_table = NULL; + bernoulli_alloc = 0; + bernoulli_size = 0; } - return b; } |