summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/bernoulli.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/bernoulli.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/bernoulli.c265
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;
}