1 /* mpfr_const_log2 -- compute natural logarithm of 2 2 3 Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramel 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 http://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 #define MPFR_NEED_LONGLONG_H 24 #include "mpfr-impl.h" 25 26 /* Declare the cache */ 27 #ifndef MPFR_USE_LOGGING 28 MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_log2, mpfr_const_log2_internal); 29 #else 30 MPFR_DECL_INIT_CACHE(__gmpfr_normal_log2, mpfr_const_log2_internal); 31 MPFR_DECL_INIT_CACHE(__gmpfr_logging_log2, mpfr_const_log2_internal); 32 mpfr_cache_ptr MPFR_THREAD_ATTR __gmpfr_cache_const_log2 = __gmpfr_normal_log2; 33 #endif 34 35 /* Set User interface */ 36 #undef mpfr_const_log2 37 int 38 mpfr_const_log2 (mpfr_ptr x, mpfr_rnd_t rnd_mode) { 39 return mpfr_cache (x, __gmpfr_cache_const_log2, rnd_mode); 40 } 41 42 /* Auxiliary function: Compute the terms from n1 to n2 (excluded) 43 3/4*sum((-1)^n*n!^2/2^n/(2*n+1)!, n = n1..n2-1). 44 Numerator is T[0], denominator is Q[0], 45 Compute P[0] only when need_P is non-zero. 46 Need 1+ceil(log(n2-n1)/log(2)) cells in T[],P[],Q[]. 47 */ 48 static void 49 S (mpz_t *T, mpz_t *P, mpz_t *Q, unsigned long n1, unsigned long n2, int need_P) 50 { 51 if (n2 == n1 + 1) 52 { 53 if (n1 == 0) 54 mpz_set_ui (P[0], 3); 55 else 56 { 57 mpz_set_ui (P[0], n1); 58 mpz_neg (P[0], P[0]); 59 } 60 if (n1 <= (ULONG_MAX / 4 - 1) / 2) 61 mpz_set_ui (Q[0], 4 * (2 * n1 + 1)); 62 else /* to avoid overflow in 4 * (2 * n1 + 1) */ 63 { 64 mpz_set_ui (Q[0], n1); 65 mpz_mul_2exp (Q[0], Q[0], 1); 66 mpz_add_ui (Q[0], Q[0], 1); 67 mpz_mul_2exp (Q[0], Q[0], 2); 68 } 69 mpz_set (T[0], P[0]); 70 } 71 else 72 { 73 unsigned long m = (n1 / 2) + (n2 / 2) + (n1 & 1UL & n2); 74 unsigned long v, w; 75 76 S (T, P, Q, n1, m, 1); 77 S (T + 1, P + 1, Q + 1, m, n2, need_P); 78 mpz_mul (T[0], T[0], Q[1]); 79 mpz_mul (T[1], T[1], P[0]); 80 mpz_add (T[0], T[0], T[1]); 81 if (need_P) 82 mpz_mul (P[0], P[0], P[1]); 83 mpz_mul (Q[0], Q[0], Q[1]); 84 85 /* remove common trailing zeroes if any */ 86 v = mpz_scan1 (T[0], 0); 87 if (v > 0) 88 { 89 w = mpz_scan1 (Q[0], 0); 90 if (w < v) 91 v = w; 92 if (need_P) 93 { 94 w = mpz_scan1 (P[0], 0); 95 if (w < v) 96 v = w; 97 } 98 /* now v = min(val(T), val(Q), val(P)) */ 99 if (v > 0) 100 { 101 mpz_fdiv_q_2exp (T[0], T[0], v); 102 mpz_fdiv_q_2exp (Q[0], Q[0], v); 103 if (need_P) 104 mpz_fdiv_q_2exp (P[0], P[0], v); 105 } 106 } 107 } 108 } 109 110 /* Don't need to save / restore exponent range: the cache does it */ 111 int 112 mpfr_const_log2_internal (mpfr_ptr x, mpfr_rnd_t rnd_mode) 113 { 114 unsigned long n = MPFR_PREC (x); 115 mpfr_prec_t w; /* working precision */ 116 unsigned long N; 117 mpz_t *T, *P, *Q; 118 mpfr_t t, q; 119 int inexact; 120 int ok = 1; /* ensures that the 1st try will give correct rounding */ 121 unsigned long lgN, i; 122 MPFR_ZIV_DECL (loop); 123 124 MPFR_LOG_FUNC ( 125 ("rnd_mode=%d", rnd_mode), 126 ("x[%Pu]=%.*Rg inex=%d", mpfr_get_prec(x), mpfr_log_prec, x, inexact)); 127 128 mpfr_init2 (t, MPFR_PREC_MIN); 129 mpfr_init2 (q, MPFR_PREC_MIN); 130 131 if (n < 1253) 132 w = n + 10; /* ensures correct rounding for the four rounding modes, 133 together with N = w / 3 + 1 (see below). */ 134 else if (n < 2571) 135 w = n + 11; /* idem */ 136 else if (n < 3983) 137 w = n + 12; 138 else if (n < 4854) 139 w = n + 13; 140 else if (n < 26248) 141 w = n + 14; 142 else 143 { 144 w = n + 15; 145 ok = 0; 146 } 147 148 MPFR_ZIV_INIT (loop, w); 149 for (;;) 150 { 151 N = w / 3 + 1; /* Warning: do not change that (even increasing N!) 152 without checking correct rounding in the above 153 ranges for n. */ 154 155 /* the following are needed for error analysis (see algorithms.tex) */ 156 MPFR_ASSERTD(w >= 3 && N >= 2); 157 158 lgN = MPFR_INT_CEIL_LOG2 (N) + 1; 159 T = (mpz_t *) (*__gmp_allocate_func) (3 * lgN * sizeof (mpz_t)); 160 P = T + lgN; 161 Q = T + 2*lgN; 162 for (i = 0; i < lgN; i++) 163 { 164 mpz_init (T[i]); 165 mpz_init (P[i]); 166 mpz_init (Q[i]); 167 } 168 169 S (T, P, Q, 0, N, 0); 170 171 mpfr_set_prec (t, w); 172 mpfr_set_prec (q, w); 173 174 mpfr_set_z (t, T[0], MPFR_RNDN); 175 mpfr_set_z (q, Q[0], MPFR_RNDN); 176 mpfr_div (t, t, q, MPFR_RNDN); 177 178 for (i = 0; i < lgN; i++) 179 { 180 mpz_clear (T[i]); 181 mpz_clear (P[i]); 182 mpz_clear (Q[i]); 183 } 184 (*__gmp_free_func) (T, 3 * lgN * sizeof (mpz_t)); 185 186 if (MPFR_LIKELY (ok != 0 187 || mpfr_can_round (t, w - 2, MPFR_RNDN, rnd_mode, n))) 188 break; 189 190 MPFR_ZIV_NEXT (loop, w); 191 } 192 MPFR_ZIV_FREE (loop); 193 194 inexact = mpfr_set (x, t, rnd_mode); 195 196 mpfr_clear (t); 197 mpfr_clear (q); 198 199 return inexact; 200 } 201