pow_ui.c revision 1.1.1.1
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 Free Software Foundation, Inc. 5Contributed by the Arenaire and Cacao 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 (("x[%#R]=%R n=%lu rnd=%d", x, x, n, rnd), 41 ("y[%#R]=%R inexact=%d", y, y, inexact)); 42 43 /* x^0 = 1 for any x, even a NaN */ 44 if (MPFR_UNLIKELY (n == 0)) 45 return mpfr_set_ui (y, 1, rnd); 46 47 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 48 { 49 if (MPFR_IS_NAN (x)) 50 { 51 MPFR_SET_NAN (y); 52 MPFR_RET_NAN; 53 } 54 else if (MPFR_IS_INF (x)) 55 { 56 /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */ 57 if (MPFR_IS_NEG (x) && (n & 1) == 1) 58 MPFR_SET_NEG (y); 59 else 60 MPFR_SET_POS (y); 61 MPFR_SET_INF (y); 62 MPFR_RET (0); 63 } 64 else /* x is zero */ 65 { 66 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 67 /* 0^n = 0 for any n */ 68 MPFR_SET_ZERO (y); 69 if (MPFR_IS_POS (x) || (n & 1) == 0) 70 MPFR_SET_POS (y); 71 else 72 MPFR_SET_NEG (y); 73 MPFR_RET (0); 74 } 75 } 76 else if (MPFR_UNLIKELY (n <= 2)) 77 { 78 if (n < 2) 79 /* x^1 = x */ 80 return mpfr_set (y, x, rnd); 81 else 82 /* x^2 = sqr(x) */ 83 return mpfr_sqr (y, x, rnd); 84 } 85 86 /* Augment exponent range */ 87 MPFR_SAVE_EXPO_MARK (expo); 88 89 /* setup initial precision */ 90 prec = MPFR_PREC (y) + 3 + GMP_NUMB_BITS 91 + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)); 92 mpfr_init2 (res, prec); 93 94 rnd1 = MPFR_IS_POS (x) ? MPFR_RNDU : MPFR_RNDD; /* away */ 95 96 MPFR_ZIV_INIT (loop, prec); 97 for (;;) 98 { 99 int i; 100 101 for (m = n, i = 0; m; i++, m >>= 1) 102 ; 103 /* now 2^(i-1) <= n < 2^i */ 104 MPFR_ASSERTD (prec > (mpfr_prec_t) i); 105 err = prec - 1 - (mpfr_prec_t) i; 106 /* First step: compute square from x */ 107 MPFR_BLOCK (flags, 108 inexact = mpfr_mul (res, x, x, MPFR_RNDU); 109 MPFR_ASSERTD (i >= 2); 110 if (n & (1UL << (i-2))) 111 inexact |= mpfr_mul (res, res, x, rnd1); 112 for (i -= 3; i >= 0 && !MPFR_BLOCK_EXCEP; i--) 113 { 114 inexact |= mpfr_mul (res, res, res, MPFR_RNDU); 115 if (n & (1UL << i)) 116 inexact |= mpfr_mul (res, res, x, rnd1); 117 }); 118 /* let r(n) be the number of roundings: we have r(2)=1, r(3)=2, 119 and r(2n)=2r(n)+1, r(2n+1)=2r(n)+2, thus r(n)=n-1. 120 Using Higham's method, to each rounding corresponds a factor 121 (1-theta) with 0 <= theta <= 2^(1-p), thus at the end the 122 absolute error is bounded by (n-1)*2^(1-p)*res <= 2*(n-1)*ulp(res) 123 since 2^(-p)*x <= ulp(x). Since n < 2^i, this gives a maximal 124 error of 2^(1+i)*ulp(res). 125 */ 126 if (MPFR_LIKELY (inexact == 0 127 || MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags) 128 || MPFR_CAN_ROUND (res, err, MPFR_PREC (y), rnd))) 129 break; 130 /* Actualisation of the precision */ 131 MPFR_ZIV_NEXT (loop, prec); 132 mpfr_set_prec (res, prec); 133 } 134 MPFR_ZIV_FREE (loop); 135 136 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags))) 137 { 138 mpz_t z; 139 140 /* Internal overflow or underflow. However the approximation error has 141 * not been taken into account. So, let's solve this problem by using 142 * mpfr_pow_z, which can handle it. This case could be improved in the 143 * future, without having to use mpfr_pow_z. 144 */ 145 MPFR_LOG_MSG (("Internal overflow or underflow," 146 " let's use mpfr_pow_z.\n", 0)); 147 mpfr_clear (res); 148 MPFR_SAVE_EXPO_FREE (expo); 149 mpz_init (z); 150 mpz_set_ui (z, n); 151 inexact = mpfr_pow_z (y, x, z, rnd); 152 mpz_clear (z); 153 return inexact; 154 } 155 156 inexact = mpfr_set (y, res, rnd); 157 mpfr_clear (res); 158 159 MPFR_SAVE_EXPO_FREE (expo); 160 return mpfr_check_range (y, inexact, rnd); 161} 162