1/* mpfr_pow_si -- power function x^y with y a signed int 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/* 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 31int 32mpfr_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