1/* mpfr_const_euler -- Euler's constant 2 3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#define MPFR_NEED_LONGLONG_H 24#include "mpfr-impl.h" 25 26/* Declare the cache */ 27MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_euler, mpfr_const_euler_internal); 28 29/* Set User Interface */ 30#undef mpfr_const_euler 31int 32mpfr_const_euler (mpfr_ptr x, mpfr_rnd_t rnd_mode) { 33 return mpfr_cache (x, __gmpfr_cache_const_euler, rnd_mode); 34} 35 36 37static void mpfr_const_euler_S2 (mpfr_ptr, unsigned long); 38static void mpfr_const_euler_R (mpfr_ptr, unsigned long); 39 40int 41mpfr_const_euler_internal (mpfr_t x, mpfr_rnd_t rnd) 42{ 43 mpfr_prec_t prec = MPFR_PREC(x), m, log2m; 44 mpfr_t y, z; 45 unsigned long n; 46 int inexact; 47 MPFR_ZIV_DECL (loop); 48 49 log2m = MPFR_INT_CEIL_LOG2 (prec); 50 m = prec + 2 * log2m + 23; 51 52 mpfr_init2 (y, m); 53 mpfr_init2 (z, m); 54 55 MPFR_ZIV_INIT (loop, m); 56 for (;;) 57 { 58 mpfr_exp_t exp_S, err; 59 /* since prec >= 1, we have m >= 24 here, which ensures n >= 9 below */ 60 n = 1 + (unsigned long) ((double) m * LOG2 / 2.0); 61 MPFR_ASSERTD (n >= 9); 62 mpfr_const_euler_S2 (y, n); /* error <= 3 ulps */ 63 exp_S = MPFR_EXP(y); 64 mpfr_set_ui (z, n, MPFR_RNDN); 65 mpfr_log (z, z, MPFR_RNDD); /* error <= 1 ulp */ 66 mpfr_sub (y, y, z, MPFR_RNDN); /* S'(n) - log(n) */ 67 /* the error is less than 1/2 + 3*2^(exp_S-EXP(y)) + 2^(EXP(z)-EXP(y)) 68 <= 1/2 + 2^(exp_S+2-EXP(y)) + 2^(EXP(z)-EXP(y)) 69 <= 1/2 + 2^(1+MAX(exp_S+2,EXP(z))-EXP(y)) */ 70 err = 1 + MAX(exp_S + 2, MPFR_EXP(z)) - MPFR_EXP(y); 71 err = (err >= -1) ? err + 1 : 0; /* error <= 2^err ulp(y) */ 72 exp_S = MPFR_EXP(y); 73 mpfr_const_euler_R (z, n); /* err <= ulp(1/2) = 2^(-m) */ 74 mpfr_sub (y, y, z, MPFR_RNDN); 75 /* err <= 1/2 ulp(y) + 2^(-m) + 2^(err + exp_S - EXP(y)) ulp(y). 76 Since the result is between 0.5 and 1, ulp(y) = 2^(-m). 77 So we get 3/2*ulp(y) + 2^(err + exp_S - EXP(y)) ulp(y). 78 3/2 + 2^e <= 2^(e+1) for e>=1, and <= 2^2 otherwise */ 79 err = err + exp_S - MPFR_EXP(y); 80 err = (err >= 1) ? err + 1 : 2; 81 if (MPFR_LIKELY (MPFR_CAN_ROUND (y, m - err, prec, rnd))) 82 break; 83 MPFR_ZIV_NEXT (loop, m); 84 mpfr_set_prec (y, m); 85 mpfr_set_prec (z, m); 86 } 87 MPFR_ZIV_FREE (loop); 88 89 inexact = mpfr_set (x, y, rnd); 90 91 mpfr_clear (y); 92 mpfr_clear (z); 93 94 return inexact; /* always inexact */ 95} 96 97static void 98mpfr_const_euler_S2_aux (mpz_t P, mpz_t Q, mpz_t T, unsigned long n, 99 unsigned long a, unsigned long b, int need_P) 100{ 101 if (a + 1 == b) 102 { 103 mpz_set_ui (P, n); 104 if (a > 1) 105 mpz_mul_si (P, P, 1 - (long) a); 106 mpz_set (T, P); 107 mpz_set_ui (Q, a); 108 mpz_mul_ui (Q, Q, a); 109 } 110 else 111 { 112 unsigned long c = (a + b) / 2; 113 mpz_t P2, Q2, T2; 114 mpfr_const_euler_S2_aux (P, Q, T, n, a, c, 1); 115 mpz_init (P2); 116 mpz_init (Q2); 117 mpz_init (T2); 118 mpfr_const_euler_S2_aux (P2, Q2, T2, n, c, b, 1); 119 mpz_mul (T, T, Q2); 120 mpz_mul (T2, T2, P); 121 mpz_add (T, T, T2); 122 if (need_P) 123 mpz_mul (P, P, P2); 124 mpz_mul (Q, Q, Q2); 125 mpz_clear (P2); 126 mpz_clear (Q2); 127 mpz_clear (T2); 128 /* divide by 2 if possible */ 129 { 130 unsigned long v2; 131 v2 = mpz_scan1 (P, 0); 132 c = mpz_scan1 (Q, 0); 133 if (c < v2) 134 v2 = c; 135 c = mpz_scan1 (T, 0); 136 if (c < v2) 137 v2 = c; 138 if (v2) 139 { 140 mpz_tdiv_q_2exp (P, P, v2); 141 mpz_tdiv_q_2exp (Q, Q, v2); 142 mpz_tdiv_q_2exp (T, T, v2); 143 } 144 } 145 } 146} 147 148/* computes S(n) = sum(n^k*(-1)^(k-1)/k!/k, k=1..ceil(4.319136566 * n)) 149 using binary splitting. 150 We have S(n) = sum(f(k), k=1..N) with N=ceil(4.319136566 * n) 151 and f(k) = n^k*(-1)*(k-1)/k!/k, 152 thus f(k)/f(k-1) = -n*(k-1)/k^2 153*/ 154static void 155mpfr_const_euler_S2 (mpfr_t x, unsigned long n) 156{ 157 mpz_t P, Q, T; 158 unsigned long N = (unsigned long) (ALPHA * (double) n + 1.0); 159 mpz_init (P); 160 mpz_init (Q); 161 mpz_init (T); 162 mpfr_const_euler_S2_aux (P, Q, T, n, 1, N + 1, 0); 163 mpfr_set_z (x, T, MPFR_RNDN); 164 mpfr_div_z (x, x, Q, MPFR_RNDN); 165 mpz_clear (P); 166 mpz_clear (Q); 167 mpz_clear (T); 168} 169 170/* computes R(n) = exp(-n)/n * sum(k!/(-n)^k, k=0..n-2) 171 with error at most 4*ulp(x). Assumes n>=2. 172 Since x <= exp(-n)/n <= 1/8, then 4*ulp(x) <= ulp(1). 173*/ 174static void 175mpfr_const_euler_R (mpfr_t x, unsigned long n) 176{ 177 unsigned long k, m; 178 mpz_t a, s; 179 mpfr_t y; 180 181 MPFR_ASSERTN (n >= 2); /* ensures sum(k!/(-n)^k, k=0..n-2) >= 2/3 */ 182 183 /* as we multiply the sum by exp(-n), we need only PREC(x) - n/LOG2 bits */ 184 m = MPFR_PREC(x) - (unsigned long) ((double) n / LOG2); 185 186 mpz_init_set_ui (a, 1); 187 mpz_mul_2exp (a, a, m); 188 mpz_init_set (s, a); 189 190 for (k = 1; k <= n; k++) 191 { 192 mpz_mul_ui (a, a, k); 193 mpz_fdiv_q_ui (a, a, n); 194 /* the error e(k) on a is e(k) <= 1 + k/n*e(k-1) with e(0)=0, 195 i.e. e(k) <= k */ 196 if (k % 2) 197 mpz_sub (s, s, a); 198 else 199 mpz_add (s, s, a); 200 } 201 /* the error on s is at most 1+2+...+n = n*(n+1)/2 */ 202 mpz_fdiv_q_ui (s, s, n); /* err <= 1 + (n+1)/2 */ 203 MPFR_ASSERTN (MPFR_PREC(x) >= mpz_sizeinbase(s, 2)); 204 mpfr_set_z (x, s, MPFR_RNDD); /* exact */ 205 mpfr_div_2ui (x, x, m, MPFR_RNDD); 206 /* now x = 1/n * sum(k!/(-n)^k, k=0..n-2) <= 1/n */ 207 /* err(x) <= (n+1)/2^m <= (n+1)*exp(n)/2^PREC(x) */ 208 209 mpfr_init2 (y, m); 210 mpfr_set_si (y, -(long)n, MPFR_RNDD); /* assumed exact */ 211 mpfr_exp (y, y, MPFR_RNDD); /* err <= ulp(y) <= exp(-n)*2^(1-m) */ 212 mpfr_mul (x, x, y, MPFR_RNDD); 213 /* err <= ulp(x) + (n + 1 + 2/n) / 2^prec(x) 214 <= ulp(x) + (n + 1 + 2/n) ulp(x)/x since x*2^(-prec(x)) < ulp(x) 215 <= ulp(x) + (n + 1 + 2/n) 3/(2n) ulp(x) since x >= 2/3*n for n >= 2 216 <= 4 * ulp(x) for n >= 2 */ 217 mpfr_clear (y); 218 219 mpz_clear (a); 220 mpz_clear (s); 221} 222