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