1/* mpfr_mpn_exp -- auxiliary function for mpfr_get_str and mpfr_set_str 2 3Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao projects, INRIA. 5Contributed by Alain Delplanque and Paul Zimmermann. 6 7This file is part of the GNU MPFR Library. 8 9The GNU MPFR Library is free software; you can redistribute it and/or modify 10it under the terms of the GNU Lesser General Public License as published by 11the Free Software Foundation; either version 3 of the License, or (at your 12option) any later version. 13 14The GNU MPFR Library is distributed in the hope that it will be useful, but 15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17License for more details. 18 19You should have received a copy of the GNU Lesser General Public License 20along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 21http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2251 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 23 24 25#define MPFR_NEED_LONGLONG_H 26#include "mpfr-impl.h" 27 28/* this function computes an approximation of b^e in {a, n}, with exponent 29 stored in exp_r. The computed value is rounded toward zero (truncated). 30 It returns an integer f such that the final error is bounded by 2^f ulps, 31 that is: 32 a*2^exp_r <= b^e <= 2^exp_r (a + 2^f), 33 where a represents {a, n}, i.e. the integer 34 a[0] + a[1]*B + ... + a[n-1]*B^(n-1) where B=2^GMP_NUMB_BITS 35 36 Return -1 is the result is exact. 37 Return -2 if an overflow occurred in the computation of exp_r. 38*/ 39 40long 41mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n) 42{ 43 mp_limb_t *c, B; 44 mpfr_exp_t f, h; 45 int i; 46 unsigned long t; /* number of bits in e */ 47 unsigned long bits; 48 size_t n1; 49 unsigned int error; /* (number - 1) of loop a^2b inexact */ 50 /* error == t means no error */ 51 int err_s_a2 = 0; 52 int err_s_ab = 0; /* number of error when shift A^2, AB */ 53 MPFR_TMP_DECL(marker); 54 55 MPFR_ASSERTN(e > 0); 56 MPFR_ASSERTN((2 <= b) && (b <= 62)); 57 58 MPFR_TMP_MARK(marker); 59 60 /* initialization of a, b, f, h */ 61 62 /* normalize the base */ 63 B = (mp_limb_t) b; 64 count_leading_zeros (h, B); 65 66 bits = GMP_NUMB_BITS - h; 67 68 B = B << h; 69 h = - h; 70 71 /* allocate space for A and set it to B */ 72 c = (mp_limb_t*) MPFR_TMP_ALLOC(2 * n * BYTES_PER_MP_LIMB); 73 a [n - 1] = B; 74 MPN_ZERO (a, n - 1); 75 /* initial exponent for A: invariant is A = {a, n} * 2^f */ 76 f = h - (n - 1) * GMP_NUMB_BITS; 77 78 /* determine number of bits in e */ 79 count_leading_zeros (t, (mp_limb_t) e); 80 81 t = GMP_NUMB_BITS - t; /* number of bits of exponent e */ 82 83 error = t; /* error <= GMP_NUMB_BITS */ 84 85 MPN_ZERO (c, 2 * n); 86 87 for (i = t - 2; i >= 0; i--) 88 { 89 90 /* determine precision needed */ 91 bits = n * GMP_NUMB_BITS - mpn_scan1 (a, 0); 92 n1 = (n * GMP_NUMB_BITS - bits) / GMP_NUMB_BITS; 93 94 /* square of A : {c+2n1, 2(n-n1)} = {a+n1, n-n1}^2 */ 95 mpn_sqr_n (c + 2 * n1, a + n1, n - n1); 96 97 /* set {c+n, 2n1-n} to 0 : {c, n} = {a, n}^2*K^n */ 98 99 /* check overflow on f */ 100 if (MPFR_UNLIKELY(f < MPFR_EXP_MIN/2 || f > MPFR_EXP_MAX/2)) 101 { 102 overflow: 103 MPFR_TMP_FREE(marker); 104 return -2; 105 } 106 /* FIXME: Could f = 2*f + n * GMP_NUMB_BITS be used? */ 107 f = 2*f; 108 MPFR_SADD_OVERFLOW (f, f, n * GMP_NUMB_BITS, 109 mpfr_exp_t, mpfr_uexp_t, 110 MPFR_EXP_MIN, MPFR_EXP_MAX, 111 goto overflow, goto overflow); 112 if ((c[2*n - 1] & MPFR_LIMB_HIGHBIT) == 0) 113 { 114 /* shift A by one bit to the left */ 115 mpn_lshift (a, c + n, n, 1); 116 a[0] |= mpn_lshift (c + n - 1, c + n - 1, 1, 1); 117 f --; 118 if (error != t) 119 err_s_a2 ++; 120 } 121 else 122 MPN_COPY (a, c + n, n); 123 124 if ((error == t) && (2 * n1 <= n) && 125 (mpn_scan1 (c + 2 * n1, 0) < (n - 2 * n1) * GMP_NUMB_BITS)) 126 error = i; 127 128 if (e & ((mpfr_exp_t) 1 << i)) 129 { 130 /* multiply A by B */ 131 c[2 * n - 1] = mpn_mul_1 (c + n - 1, a, n, B); 132 f += h + GMP_NUMB_BITS; 133 if ((c[2 * n - 1] & MPFR_LIMB_HIGHBIT) == 0) 134 { /* shift A by one bit to the left */ 135 mpn_lshift (a, c + n, n, 1); 136 a[0] |= mpn_lshift (c + n - 1, c + n - 1, 1, 1); 137 f --; 138 } 139 else 140 { 141 MPN_COPY (a, c + n, n); 142 if (error != t) 143 err_s_ab ++; 144 } 145 if ((error == t) && (c[n - 1] != 0)) 146 error = i; 147 } 148 } 149 150 MPFR_TMP_FREE(marker); 151 152 *exp_r = f; 153 154 if (error == t) 155 return -1; /* result is exact */ 156 else /* error <= t-2 <= GMP_NUMB_BITS-2 157 err_s_ab, err_s_a2 <= t-1 */ 158 { 159 /* if there are p loops after the first inexact result, with 160 j shifts in a^2 and l shifts in a*b, then the final error is 161 at most 2^(p+ceil((j+1)/2)+l+1)*ulp(res). 162 This is bounded by 2^(5/2*t-1/2) where t is the number of bits of e. 163 */ 164 error = error + err_s_ab + err_s_a2 / 2 + 3; /* <= 5t/2-1/2 */ 165#if 0 166 if ((error - 1) >= ((n * GMP_NUMB_BITS - 1) / 2)) 167 error = n * GMP_NUMB_BITS; /* result is completely wrong: 168 this is very unlikely since error is 169 at most 5/2*log_2(e), and 170 n * GMP_NUMB_BITS is at least 171 3*log_2(e) */ 172#endif 173 return error; 174 } 175} 176