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