1 /* mpfr_pow_si -- power function x^y with y a signed int 2 3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4 Contributed by the Arenaire 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 /* The computation of y = pow_si(x,n) is done by 27 * y = pow_ui(x,n) if n >= 0 28 * y = 1 / pow_ui(x,-n) if n < 0 29 */ 30 31 int 32 mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd) 33 { 34 MPFR_LOG_FUNC 35 (("x[%Pu]=%.*Rg n=%ld rnd=%d", 36 mpfr_get_prec (x), mpfr_log_prec, x, n, rnd), 37 ("y[%Pu]=%.*Rg", mpfr_get_prec (y), mpfr_log_prec, y)); 38 39 if (n >= 0) 40 return mpfr_pow_ui (y, x, n, rnd); 41 else 42 { 43 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 44 { 45 if (MPFR_IS_NAN (x)) 46 { 47 MPFR_SET_NAN (y); 48 MPFR_RET_NAN; 49 } 50 else 51 { 52 int positive = MPFR_IS_POS (x) || ((unsigned long) n & 1) == 0; 53 if (MPFR_IS_INF (x)) 54 MPFR_SET_ZERO (y); 55 else /* x is zero */ 56 { 57 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 58 MPFR_SET_INF (y); 59 mpfr_set_divby0 (); 60 } 61 if (positive) 62 MPFR_SET_POS (y); 63 else 64 MPFR_SET_NEG (y); 65 MPFR_RET (0); 66 } 67 } 68 69 /* detect exact powers: x^(-n) is exact iff x is a power of 2 */ 70 if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0) 71 { 72 mpfr_exp_t expx = MPFR_EXP (x) - 1, expy; 73 MPFR_ASSERTD (n < 0); 74 /* Warning: n * expx may overflow! 75 * 76 * Some systems (apparently alpha-freebsd) abort with 77 * LONG_MIN / 1, and LONG_MIN / -1 is undefined. 78 * http://www.freebsd.org/cgi/query-pr.cgi?pr=72024 79 * 80 * Proof of the overflow checking. The expressions below are 81 * assumed to be on the rational numbers, but the word "overflow" 82 * still has its own meaning in the C context. / still denotes 83 * the integer (truncated) division, and // denotes the exact 84 * division. 85 * - First, (__gmpfr_emin - 1) / n and (__gmpfr_emax - 1) / n 86 * cannot overflow due to the constraints on the exponents of 87 * MPFR numbers. 88 * - If n = -1, then n * expx = - expx, which is representable 89 * because of the constraints on the exponents of MPFR numbers. 90 * - If expx = 0, then n * expx = 0, which is representable. 91 * - If n < -1 and expx > 0: 92 * + If expx > (__gmpfr_emin - 1) / n, then 93 * expx >= (__gmpfr_emin - 1) / n + 1 94 * > (__gmpfr_emin - 1) // n, 95 * and 96 * n * expx < __gmpfr_emin - 1, 97 * i.e. 98 * n * expx <= __gmpfr_emin - 2. 99 * This corresponds to an underflow, with a null result in 100 * the rounding-to-nearest mode. 101 * + If expx <= (__gmpfr_emin - 1) / n, then n * expx cannot 102 * overflow since 0 < expx <= (__gmpfr_emin - 1) / n and 103 * 0 > n * expx >= n * ((__gmpfr_emin - 1) / n) 104 * >= __gmpfr_emin - 1. 105 * - If n < -1 and expx < 0: 106 * + If expx < (__gmpfr_emax - 1) / n, then 107 * expx <= (__gmpfr_emax - 1) / n - 1 108 * < (__gmpfr_emax - 1) // n, 109 * and 110 * n * expx > __gmpfr_emax - 1, 111 * i.e. 112 * n * expx >= __gmpfr_emax. 113 * This corresponds to an overflow (2^(n * expx) has an 114 * exponent > __gmpfr_emax). 115 * + If expx >= (__gmpfr_emax - 1) / n, then n * expx cannot 116 * overflow since 0 > expx >= (__gmpfr_emax - 1) / n and 117 * 0 < n * expx <= n * ((__gmpfr_emax - 1) / n) 118 * <= __gmpfr_emax - 1. 119 * Note: one could use expx bounds based on MPFR_EXP_MIN and 120 * MPFR_EXP_MAX instead of __gmpfr_emin and __gmpfr_emax. The 121 * current bounds do not lead to noticeably slower code and 122 * allow us to avoid a bug in Sun's compiler for Solaris/x86 123 * (when optimizations are enabled); known affected versions: 124 * cc: Sun C 5.8 2005/10/13 125 * cc: Sun C 5.8 Patch 121016-02 2006/03/31 126 * cc: Sun C 5.8 Patch 121016-04 2006/10/18 127 */ 128 expy = 129 n != -1 && expx > 0 && expx > (__gmpfr_emin - 1) / n ? 130 MPFR_EMIN_MIN - 2 /* Underflow */ : 131 n != -1 && expx < 0 && expx < (__gmpfr_emax - 1) / n ? 132 MPFR_EMAX_MAX /* Overflow */ : n * expx; 133 return mpfr_set_si_2exp (y, n % 2 ? MPFR_INT_SIGN (x) : 1, 134 expy, rnd); 135 } 136 137 /* General case */ 138 { 139 /* Declaration of the intermediary variable */ 140 mpfr_t t; 141 /* Declaration of the size variable */ 142 mpfr_prec_t Ny; /* target precision */ 143 mpfr_prec_t Nt; /* working precision */ 144 mpfr_rnd_t rnd1; 145 int size_n; 146 int inexact; 147 unsigned long abs_n; 148 MPFR_SAVE_EXPO_DECL (expo); 149 MPFR_ZIV_DECL (loop); 150 151 abs_n = - (unsigned long) n; 152 count_leading_zeros (size_n, (mp_limb_t) abs_n); 153 size_n = GMP_NUMB_BITS - size_n; 154 155 /* initial working precision */ 156 Ny = MPFR_PREC (y); 157 Nt = Ny + size_n + 3 + MPFR_INT_CEIL_LOG2 (Ny); 158 159 MPFR_SAVE_EXPO_MARK (expo); 160 161 /* initialise of intermediary variable */ 162 mpfr_init2 (t, Nt); 163 164 /* We will compute rnd(rnd1(1/x) ^ |n|), where rnd1 is the rounding 165 toward sign(x), to avoid spurious overflow or underflow, as in 166 mpfr_pow_z. */ 167 rnd1 = MPFR_EXP (x) < 1 ? MPFR_RNDZ : 168 (MPFR_SIGN (x) > 0 ? MPFR_RNDU : MPFR_RNDD); 169 170 MPFR_ZIV_INIT (loop, Nt); 171 for (;;) 172 { 173 MPFR_BLOCK_DECL (flags); 174 175 /* compute (1/x)^|n| */ 176 MPFR_BLOCK (flags, mpfr_ui_div (t, 1, x, rnd1)); 177 MPFR_ASSERTD (! MPFR_UNDERFLOW (flags)); 178 /* t = (1/x)*(1+theta) where |theta| <= 2^(-Nt) */ 179 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags))) 180 goto overflow; 181 MPFR_BLOCK (flags, mpfr_pow_ui (t, t, abs_n, rnd)); 182 /* t = (1/x)^|n|*(1+theta')^(|n|+1) where |theta'| <= 2^(-Nt). 183 If (|n|+1)*2^(-Nt) <= 1/2, which is satisfied as soon as 184 Nt >= bits(n)+2, then we can use Lemma \ref{lemma_graillat} 185 from algorithms.tex, which yields x^n*(1+theta) with 186 |theta| <= 2(|n|+1)*2^(-Nt), thus the error is bounded by 187 2(|n|+1) ulps <= 2^(bits(n)+2) ulps. */ 188 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags))) 189 { 190 overflow: 191 MPFR_ZIV_FREE (loop); 192 mpfr_clear (t); 193 MPFR_SAVE_EXPO_FREE (expo); 194 MPFR_LOG_MSG (("overflow\n", 0)); 195 return mpfr_overflow (y, rnd, abs_n & 1 ? 196 MPFR_SIGN (x) : MPFR_SIGN_POS); 197 } 198 if (MPFR_UNLIKELY (MPFR_UNDERFLOW (flags))) 199 { 200 MPFR_ZIV_FREE (loop); 201 mpfr_clear (t); 202 MPFR_LOG_MSG (("underflow\n", 0)); 203 if (rnd == MPFR_RNDN) 204 { 205 mpfr_t y2, nn; 206 207 /* We cannot decide now whether the result should be 208 rounded toward zero or away from zero. So, like 209 in mpfr_pow_pos_z, let's use the general case of 210 mpfr_pow in precision 2. */ 211 MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_SIGN (x), 212 MPFR_EXP (x) - 1) != 0); 213 mpfr_init2 (y2, 2); 214 mpfr_init2 (nn, sizeof (long) * CHAR_BIT); 215 inexact = mpfr_set_si (nn, n, MPFR_RNDN); 216 MPFR_ASSERTN (inexact == 0); 217 inexact = mpfr_pow_general (y2, x, nn, rnd, 1, 218 (mpfr_save_expo_t *) NULL); 219 mpfr_clear (nn); 220 mpfr_set (y, y2, MPFR_RNDN); 221 mpfr_clear (y2); 222 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_UNDERFLOW); 223 goto end; 224 } 225 else 226 { 227 MPFR_SAVE_EXPO_FREE (expo); 228 return mpfr_underflow (y, rnd, abs_n & 1 ? 229 MPFR_SIGN (x) : MPFR_SIGN_POS); 230 } 231 } 232 /* error estimate -- see pow function in algorithms.ps */ 233 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - size_n - 2, Ny, rnd))) 234 break; 235 236 /* actualisation of the precision */ 237 MPFR_ZIV_NEXT (loop, Nt); 238 mpfr_set_prec (t, Nt); 239 } 240 MPFR_ZIV_FREE (loop); 241 242 inexact = mpfr_set (y, t, rnd); 243 mpfr_clear (t); 244 245 end: 246 MPFR_SAVE_EXPO_FREE (expo); 247 return mpfr_check_range (y, inexact, rnd); 248 } 249 } 250 } 251