summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/const_euler.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/const_euler.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/const_euler.c341
1 files changed, 198 insertions, 143 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/const_euler.c b/Build/source/libs/mpfr/mpfr-src/src/const_euler.c
index 7bd879d4b81..44af5b82375 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/const_euler.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/const_euler.c
@@ -1,7 +1,7 @@
/* mpfr_const_euler -- Euler's constant
Copyright 2001-2017 Free Software Foundation, Inc.
-Contributed by the AriC and Caramba projects, INRIA.
+Contributed by Fredrik Johansson.
This file is part of the GNU MPFR Library.
@@ -20,6 +20,9 @@ along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
+/* The approximation error bound uses Theorem 1 and Remark 2 in
+ http://arxiv.org/pdf/1312.0039v1.pdf */
+
#define MPFR_NEED_LONGLONG_H
#include "mpfr-impl.h"
@@ -34,188 +37,240 @@ mpfr_const_euler (mpfr_ptr x, mpfr_rnd_t rnd_mode) {
}
-static void mpfr_const_euler_S2 (mpfr_ptr, unsigned long);
-static void mpfr_const_euler_R (mpfr_ptr, unsigned long);
-
-int
-mpfr_const_euler_internal (mpfr_t x, mpfr_rnd_t rnd)
+typedef struct
{
- mpfr_prec_t prec = MPFR_PREC(x), m, log2m;
- mpfr_t y, z;
- unsigned long n;
- int inexact;
- MPFR_ZIV_DECL (loop);
+ mpz_t P;
+ mpz_t Q;
+ mpz_t T;
+ mpz_t C;
+ mpz_t D;
+ mpz_t V;
+} mpfr_const_euler_bs_struct;
+
+typedef mpfr_const_euler_bs_struct mpfr_const_euler_bs_t[1];
- log2m = MPFR_INT_CEIL_LOG2 (prec);
- m = prec + 2 * log2m + 23;
+static void
+mpfr_const_euler_bs_init (mpfr_const_euler_bs_t s)
+{
+ mpz_init (s->P);
+ mpz_init (s->Q);
+ mpz_init (s->T);
+ mpz_init (s->C);
+ mpz_init (s->D);
+ mpz_init (s->V);
+}
- mpfr_init2 (y, m);
- mpfr_init2 (z, m);
+static void
+mpfr_const_euler_bs_clear (mpfr_const_euler_bs_t s)
+{
+ mpz_clear (s->P);
+ mpz_clear (s->Q);
+ mpz_clear (s->T);
+ mpz_clear (s->C);
+ mpz_clear (s->D);
+ mpz_clear (s->V);
+}
- MPFR_ZIV_INIT (loop, m);
- for (;;)
+static void
+mpfr_const_euler_bs_1 (mpfr_const_euler_bs_t s,
+ unsigned long n1, unsigned long n2, unsigned long N,
+ int cont)
+{
+ if (n2 - n1 == 1)
{
- mpfr_exp_t exp_S, err;
- /* since prec >= 1, we have m >= 24 here, which ensures n >= 9 below */
- n = 1 + (unsigned long) ((double) m * LOG2 / 2.0);
- MPFR_ASSERTD (n >= 9);
- mpfr_const_euler_S2 (y, n); /* error <= 3 ulps */
- exp_S = MPFR_EXP(y);
- mpfr_set_ui (z, n, MPFR_RNDN);
- mpfr_log (z, z, MPFR_RNDD); /* error <= 1 ulp */
- mpfr_sub (y, y, z, MPFR_RNDN); /* S'(n) - log(n) */
- /* the error is less than 1/2 + 3*2^(exp_S-EXP(y)) + 2^(EXP(z)-EXP(y))
- <= 1/2 + 2^(exp_S+2-EXP(y)) + 2^(EXP(z)-EXP(y))
- <= 1/2 + 2^(1+MAX(exp_S+2,EXP(z))-EXP(y)) */
- err = 1 + MAX(exp_S + 2, MPFR_EXP(z)) - MPFR_EXP(y);
- err = (err >= -1) ? err + 1 : 0; /* error <= 2^err ulp(y) */
- exp_S = MPFR_EXP(y);
- mpfr_const_euler_R (z, n); /* err <= ulp(1/2) = 2^(-m) */
- mpfr_sub (y, y, z, MPFR_RNDN);
- /* err <= 1/2 ulp(y) + 2^(-m) + 2^(err + exp_S - EXP(y)) ulp(y).
- Since the result is between 0.5 and 1, ulp(y) = 2^(-m).
- So we get 3/2*ulp(y) + 2^(err + exp_S - EXP(y)) ulp(y).
- 3/2 + 2^e <= 2^(e+1) for e>=1, and <= 2^2 otherwise */
- err = err + exp_S - MPFR_EXP(y);
- err = (err >= 1) ? err + 1 : 2;
- if (MPFR_LIKELY (MPFR_CAN_ROUND (y, m - err, prec, rnd)))
- break;
- MPFR_ZIV_NEXT (loop, m);
- mpfr_set_prec (y, m);
- mpfr_set_prec (z, m);
+ mpz_set_ui (s->P, N);
+ mpz_mul (s->P, s->P, s->P);
+ mpz_set_ui (s->Q, n1 + 1);
+ mpz_mul (s->Q, s->Q, s->Q);
+ mpz_set_ui (s->C, 1);
+ mpz_set_ui (s->D, n1 + 1);
+ mpz_set (s->T, s->P);
+ mpz_set (s->V, s->P);
}
- MPFR_ZIV_FREE (loop);
+ else
+ {
+ mpfr_const_euler_bs_t L, R;
+ mpz_t t, u, v;
+ unsigned long m = (n1 + n2) / 2;
- inexact = mpfr_set (x, y, rnd);
+ mpfr_const_euler_bs_init (L);
+ mpfr_const_euler_bs_init (R);
+ mpfr_const_euler_bs_1 (L, n1, m, N, 1);
+ mpfr_const_euler_bs_1 (R, m, n2, N, 1);
- mpfr_clear (y);
- mpfr_clear (z);
+ mpz_init (t);
+ mpz_init (u);
+ mpz_init (v);
- return inexact; /* always inexact */
+ if (cont)
+ mpz_mul (s->P, L->P, R->P);
+
+ mpz_mul (s->Q, L->Q, R->Q);
+ mpz_mul (s->D, L->D, R->D);
+
+ /* T = LP RT + RQ LT*/
+ mpz_mul (t, L->P, R->T);
+ mpz_mul (v, R->Q, L->T);
+ mpz_add (s->T, t, v);
+
+ /* C = LC RD + RC LD */
+ if (cont)
+ {
+ mpz_mul (s->C, L->C, R->D);
+ mpz_addmul (s->C, R->C, L->D);
+ }
+
+ /* V = RD (RQ LV + LC LP RT) + LD LP RV */
+ mpz_mul (u, L->P, R->V);
+ mpz_mul (u, u, L->D);
+ mpz_mul (v, R->Q, L->V);
+ mpz_addmul (v, t, L->C);
+ mpz_mul (v, v, R->D);
+ mpz_add (s->V, u, v);
+
+ mpfr_const_euler_bs_clear (L);
+ mpfr_const_euler_bs_clear (R);
+ mpz_clear (t);
+ mpz_clear (u);
+ mpz_clear (v);
+ }
}
static void
-mpfr_const_euler_S2_aux (mpz_t P, mpz_t Q, mpz_t T, unsigned long n,
- unsigned long a, unsigned long b, int need_P)
+mpfr_const_euler_bs_2 (mpz_t P, mpz_t Q, mpz_t T,
+ unsigned long n1, unsigned long n2, unsigned long N,
+ int cont)
{
- if (a + 1 == b)
+ if (n2 - n1 == 1)
{
- mpz_set_ui (P, n);
- if (a > 1)
- mpz_mul_si (P, P, 1 - (long) a);
+ if (n1 == 0)
+ {
+ mpz_set_ui (P, 1);
+ mpz_set_ui (Q, 4 * N);
+ }
+ else
+ {
+ mpz_set_ui (P, 2 * n1 - 1);
+ mpz_pow_ui (P, P, 3);
+ mpz_set_ui (Q, 32 * n1);
+ mpz_mul_ui (Q, Q, N);
+ mpz_mul_ui (Q, Q, N);
+ }
mpz_set (T, P);
- mpz_set_ui (Q, a);
- mpz_mul_ui (Q, Q, a);
}
else
{
- unsigned long c = (a + b) / 2;
mpz_t P2, Q2, T2;
- mpfr_const_euler_S2_aux (P, Q, T, n, a, c, 1);
+ unsigned long m = (n1 + n2) / 2;
+
mpz_init (P2);
mpz_init (Q2);
mpz_init (T2);
- mpfr_const_euler_S2_aux (P2, Q2, T2, n, c, b, 1);
+ mpfr_const_euler_bs_2 (P, Q, T, n1, m, N, 1);
+ mpfr_const_euler_bs_2 (P2, Q2, T2, m, n2, N, 1);
mpz_mul (T, T, Q2);
mpz_mul (T2, T2, P);
mpz_add (T, T, T2);
- if (need_P)
+ if (cont)
mpz_mul (P, P, P2);
mpz_mul (Q, Q, Q2);
mpz_clear (P2);
mpz_clear (Q2);
mpz_clear (T2);
- /* divide by 2 if possible */
- {
- unsigned long v2;
- v2 = mpz_scan1 (P, 0);
- c = mpz_scan1 (Q, 0);
- if (c < v2)
- v2 = c;
- c = mpz_scan1 (T, 0);
- if (c < v2)
- v2 = c;
- if (v2)
- {
- mpz_tdiv_q_2exp (P, P, v2);
- mpz_tdiv_q_2exp (Q, Q, v2);
- mpz_tdiv_q_2exp (T, T, v2);
- }
- }
}
}
-/* computes S(n) = sum(n^k*(-1)^(k-1)/k!/k, k=1..ceil(4.319136566 * n))
- using binary splitting.
- We have S(n) = sum(f(k), k=1..N) with N=ceil(4.319136566 * n)
- and f(k) = n^k*(-1)*(k-1)/k!/k,
- thus f(k)/f(k-1) = -n*(k-1)/k^2
-*/
-static void
-mpfr_const_euler_S2 (mpfr_t x, unsigned long n)
-{
- mpz_t P, Q, T;
- unsigned long N = (unsigned long) (ALPHA * (double) n + 1.0);
- mpz_init (P);
- mpz_init (Q);
- mpz_init (T);
- mpfr_const_euler_S2_aux (P, Q, T, n, 1, N + 1, 0);
- mpfr_set_z (x, T, MPFR_RNDN);
- mpfr_div_z (x, x, Q, MPFR_RNDN);
- mpz_clear (P);
- mpz_clear (Q);
- mpz_clear (T);
-}
-
-/* computes R(n) = exp(-n)/n * sum(k!/(-n)^k, k=0..n-2)
- with error at most 4*ulp(x). Assumes n>=2.
- Since x <= exp(-n)/n <= 1/8, then 4*ulp(x) <= ulp(1).
-*/
-static void
-mpfr_const_euler_R (mpfr_t x, unsigned long n)
+int
+mpfr_const_euler_internal (mpfr_t x, mpfr_rnd_t rnd)
{
- unsigned long k, m;
- mpz_t a, s;
+ mpfr_const_euler_bs_t sum;
+ mpz_t t, u, v;
+ unsigned long n, N;
+ mpfr_prec_t prec, wp, magn;
mpfr_t y;
+ int inexact;
+ MPFR_ZIV_DECL (loop);
- MPFR_ASSERTN (n >= 2); /* ensures sum(k!/(-n)^k, k=0..n-2) >= 2/3 */
-
- /* as we multiply the sum by exp(-n), we need only PREC(x) - n/LOG2 bits */
- m = MPFR_PREC(x) - (unsigned long) ((double) n / LOG2);
+ prec = mpfr_get_prec (x);
+ wp = prec + 24;
- mpz_init_set_ui (a, 1);
- mpz_mul_2exp (a, a, m);
- mpz_init_set (s, a);
+ mpfr_init2 (y, wp);
+ mpfr_const_euler_bs_init (sum);
+ mpz_init (t);
+ mpz_init (u);
+ mpz_init (v);
- for (k = 1; k <= n; k++)
+ MPFR_ZIV_INIT (loop, wp);
+ for (;;)
{
- mpz_mul_ui (a, a, k);
- mpz_fdiv_q_ui (a, a, n);
- /* the error e(k) on a is e(k) <= 1 + k/n*e(k-1) with e(0)=0,
- i.e. e(k) <= k */
- if (k % 2)
- mpz_sub (s, s, a);
- else
- mpz_add (s, s, a);
+ /* The approximation error is bounded by 24 exp(-8n) when
+ n > 1, which is smaller than 2^-wp if
+ n > (wp + log_2(24)) * (log(2)/8).
+ Note log2(24) < 5 and log(2)/8 < 866434 / 10000000. */
+ mpz_set_ui (t, wp + 5);
+ mpz_mul_ui (t, t, 866434);
+ mpz_cdiv_q_ui (t, t, 10000000);
+ n = mpz_get_ui (t);
+
+ /* It is sufficient to take N >= alpha*n + 1
+ where alpha = 3/LambertW(3/e) = 4.970625759544... */
+ mpz_set_ui (t, n);
+ mpz_mul_ui (t, t, 4970626);
+ mpz_cdiv_q_ui (t, t, 1000000);
+ mpz_add_ui (t, t, 1);
+ N = mpz_get_ui (t);
+
+ /* V / ((T + Q) * D) = S / I
+ where S = sum_{k=0}^{N-1} H_k n^(2k) / (k!)^2,
+ I = sum_{k=0}^{N-1} n^(2k) / (k!)^2 */
+ mpfr_const_euler_bs_1 (sum, 0, N, n, 0);
+ mpz_add (sum->T, sum->T, sum->Q);
+ mpz_mul (t, sum->T, sum->D);
+ mpz_mul_2exp (u, sum->V, wp);
+ mpz_tdiv_q (v, u, t);
+ /* v * 2^-wp = S/I with error < 1 */
+
+ /* C / (D * V) = U where
+ U = (1/(4n)) sum_{k=0}^{2n-1} [(2k)!]^3 / ((k!)^4 8^(2k) (2n)^(2k)) */
+ mpfr_const_euler_bs_2 (sum->C, sum->D, sum->V, 0, 2*n, n, 0);
+ mpz_mul (t, sum->Q, sum->Q);
+ mpz_mul (t, t, sum->V);
+ mpz_mul (u, sum->T, sum->T);
+ mpz_mul (u, u, sum->D);
+ mpz_mul_2exp (t, t, wp);
+ mpz_tdiv_q (t, t, u);
+ /* t * 2^-wp = U/I^2 with error < 1 */
+
+ /* gamma = S/I - U/I^2 - log(n) with error at most 2^-wp */
+ mpz_sub (v, v, t);
+ /* v * 2^-wp now equals gamma + log(n) with error at most 3*2^-wp */
+
+ /* log(n) < 2^ceil(log2(n)) */
+ magn = MPFR_INT_CEIL_LOG2(n);
+ mpfr_set_prec (y, wp + magn);
+ mpfr_set_ui (y, n, MPFR_RNDZ); /* exact */
+ mpfr_log (y, y, MPFR_RNDZ); /* error < 2^-wp */
+
+ mpfr_mul_2exp (y, y, wp, MPFR_RNDZ);
+ mpfr_z_sub (y, v, y, MPFR_RNDZ);
+ mpfr_div_2exp (y, y, wp, MPFR_RNDZ);
+ /* rounding error from the last subtraction < 2^-wp */
+ /* so y = gamma with error < 5*2^-wp */
+
+ if (MPFR_LIKELY (MPFR_CAN_ROUND (y, wp - 3, prec, rnd)))
+ break;
+
+ MPFR_ZIV_NEXT (loop, wp);
}
- /* the error on s is at most 1+2+...+n = n*(n+1)/2 */
- mpz_fdiv_q_ui (s, s, n); /* err <= 1 + (n+1)/2 */
- MPFR_ASSERTN (MPFR_PREC(x) >= mpz_sizeinbase(s, 2));
- mpfr_set_z (x, s, MPFR_RNDD); /* exact */
- mpfr_div_2ui (x, x, m, MPFR_RNDD);
- /* now x = 1/n * sum(k!/(-n)^k, k=0..n-2) <= 1/n */
- /* err(x) <= (n+1)/2^m <= (n+1)*exp(n)/2^PREC(x) */
-
- mpfr_init2 (y, m);
- mpfr_set_si (y, -(long)n, MPFR_RNDD); /* assumed exact */
- mpfr_exp (y, y, MPFR_RNDD); /* err <= ulp(y) <= exp(-n)*2^(1-m) */
- mpfr_mul (x, x, y, MPFR_RNDD);
- /* err <= ulp(x) + (n + 1 + 2/n) / 2^prec(x)
- <= ulp(x) + (n + 1 + 2/n) ulp(x)/x since x*2^(-prec(x)) < ulp(x)
- <= ulp(x) + (n + 1 + 2/n) 3/(2n) ulp(x) since x >= 2/3*n for n >= 2
- <= 4 * ulp(x) for n >= 2 */
+
+ MPFR_ZIV_FREE (loop);
+ inexact = mpfr_set (x, y, rnd);
+
mpfr_clear (y);
+ mpz_clear (t);
+ mpz_clear (u);
+ mpz_clear (v);
+ mpfr_const_euler_bs_clear (sum);
- mpz_clear (a);
- mpz_clear (s);
+ return inexact; /* always inexact */
}