1 /* mpfr_pow_ui-- compute the power of a floating-point 2 by a machine integer 3 4 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 5 Contributed by the AriC and Caramel projects, INRIA. 6 7 This file is part of the GNU MPFR Library. 8 9 The GNU MPFR Library is free software; you can redistribute it and/or modify 10 it under the terms of the GNU Lesser General Public License as published by 11 the Free Software Foundation; either version 3 of the License, or (at your 12 option) any later version. 13 14 The GNU MPFR Library is distributed in the hope that it will be useful, but 15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17 License for more details. 18 19 You should have received a copy of the GNU Lesser General Public License 20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 23 24 #define MPFR_NEED_LONGLONG_H 25 #include "mpfr-impl.h" 26 27 /* sets y to x^n, and return 0 if exact, non-zero otherwise */ 28 int 29 mpfr_pow_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int n, mpfr_rnd_t rnd) 30 { 31 unsigned long m; 32 mpfr_t res; 33 mpfr_prec_t prec, err; 34 int inexact; 35 mpfr_rnd_t rnd1; 36 MPFR_SAVE_EXPO_DECL (expo); 37 MPFR_ZIV_DECL (loop); 38 MPFR_BLOCK_DECL (flags); 39 40 MPFR_LOG_FUNC 41 (("x[%Pu]=%.*Rg n=%lu rnd=%d", 42 mpfr_get_prec (x), mpfr_log_prec, x, n, rnd), 43 ("y[%Pu]=%.*Rg inexact=%d", 44 mpfr_get_prec (y), mpfr_log_prec, y, inexact)); 45 46 /* x^0 = 1 for any x, even a NaN */ 47 if (MPFR_UNLIKELY (n == 0)) 48 return mpfr_set_ui (y, 1, rnd); 49 50 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 51 { 52 if (MPFR_IS_NAN (x)) 53 { 54 MPFR_SET_NAN (y); 55 MPFR_RET_NAN; 56 } 57 else if (MPFR_IS_INF (x)) 58 { 59 /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */ 60 if (MPFR_IS_NEG (x) && (n & 1) == 1) 61 MPFR_SET_NEG (y); 62 else 63 MPFR_SET_POS (y); 64 MPFR_SET_INF (y); 65 MPFR_RET (0); 66 } 67 else /* x is zero */ 68 { 69 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 70 /* 0^n = 0 for any n */ 71 MPFR_SET_ZERO (y); 72 if (MPFR_IS_POS (x) || (n & 1) == 0) 73 MPFR_SET_POS (y); 74 else 75 MPFR_SET_NEG (y); 76 MPFR_RET (0); 77 } 78 } 79 else if (MPFR_UNLIKELY (n <= 2)) 80 { 81 if (n < 2) 82 /* x^1 = x */ 83 return mpfr_set (y, x, rnd); 84 else 85 /* x^2 = sqr(x) */ 86 return mpfr_sqr (y, x, rnd); 87 } 88 89 /* Augment exponent range */ 90 MPFR_SAVE_EXPO_MARK (expo); 91 92 /* setup initial precision */ 93 prec = MPFR_PREC (y) + 3 + GMP_NUMB_BITS 94 + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)); 95 mpfr_init2 (res, prec); 96 97 rnd1 = MPFR_IS_POS (x) ? MPFR_RNDU : MPFR_RNDD; /* away */ 98 99 MPFR_ZIV_INIT (loop, prec); 100 for (;;) 101 { 102 int i; 103 104 for (m = n, i = 0; m; i++, m >>= 1) 105 ; 106 /* now 2^(i-1) <= n < 2^i */ 107 MPFR_ASSERTD (prec > (mpfr_prec_t) i); 108 err = prec - 1 - (mpfr_prec_t) i; 109 /* First step: compute square from x */ 110 MPFR_BLOCK (flags, 111 inexact = mpfr_mul (res, x, x, MPFR_RNDU); 112 MPFR_ASSERTD (i >= 2); 113 if (n & (1UL << (i-2))) 114 inexact |= mpfr_mul (res, res, x, rnd1); 115 for (i -= 3; i >= 0 && !MPFR_BLOCK_EXCEP; i--) 116 { 117 inexact |= mpfr_mul (res, res, res, MPFR_RNDU); 118 if (n & (1UL << i)) 119 inexact |= mpfr_mul (res, res, x, rnd1); 120 }); 121 /* let r(n) be the number of roundings: we have r(2)=1, r(3)=2, 122 and r(2n)=2r(n)+1, r(2n+1)=2r(n)+2, thus r(n)=n-1. 123 Using Higham's method, to each rounding corresponds a factor 124 (1-theta) with 0 <= theta <= 2^(1-p), thus at the end the 125 absolute error is bounded by (n-1)*2^(1-p)*res <= 2*(n-1)*ulp(res) 126 since 2^(-p)*x <= ulp(x). Since n < 2^i, this gives a maximal 127 error of 2^(1+i)*ulp(res). 128 */ 129 if (MPFR_LIKELY (inexact == 0 130 || MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags) 131 || MPFR_CAN_ROUND (res, err, MPFR_PREC (y), rnd))) 132 break; 133 /* Actualisation of the precision */ 134 MPFR_ZIV_NEXT (loop, prec); 135 mpfr_set_prec (res, prec); 136 } 137 MPFR_ZIV_FREE (loop); 138 139 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags))) 140 { 141 mpz_t z; 142 143 /* Internal overflow or underflow. However the approximation error has 144 * not been taken into account. So, let's solve this problem by using 145 * mpfr_pow_z, which can handle it. This case could be improved in the 146 * future, without having to use mpfr_pow_z. 147 */ 148 MPFR_LOG_MSG (("Internal overflow or underflow," 149 " let's use mpfr_pow_z.\n", 0)); 150 mpfr_clear (res); 151 MPFR_SAVE_EXPO_FREE (expo); 152 mpz_init (z); 153 mpz_set_ui (z, n); 154 inexact = mpfr_pow_z (y, x, z, rnd); 155 mpz_clear (z); 156 return inexact; 157 } 158 159 inexact = mpfr_set (y, res, rnd); 160 mpfr_clear (res); 161 162 MPFR_SAVE_EXPO_FREE (expo); 163 return mpfr_check_range (y, inexact, rnd); 164 } 165