1/* mpfr_pow_ui-- compute the power of a floating-point 2 by a machine integer 3 4Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 5Contributed by the AriC and Caramel projects, INRIA. 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#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 */ 28int 29mpfr_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