1 /* mpfr_pow_ui -- compute the power of a floating-point by a machine integer 2 3 Copyright 1999-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 24 #include "mpfr-impl.h" 25 26 #ifndef POW_U 27 #define POW_U mpfr_pow_ui 28 #define MPZ_SET_U mpz_set_ui 29 #define UTYPE unsigned long int 30 #define FSPEC "l" 31 #endif 32 33 /* sets y to x^n, and return 0 if exact, non-zero otherwise */ 34 int 35 POW_U (mpfr_ptr y, mpfr_srcptr x, UTYPE n, mpfr_rnd_t rnd) 36 { 37 UTYPE m; 38 mpfr_t res; 39 mpfr_prec_t prec, err, nlen; 40 int inexact; 41 mpfr_rnd_t rnd1; 42 MPFR_SAVE_EXPO_DECL (expo); 43 MPFR_ZIV_DECL (loop); 44 MPFR_BLOCK_DECL (flags); 45 46 MPFR_LOG_FUNC 47 (("x[%Pd]=%.*Rg n=%" FSPEC "u rnd=%d", 48 mpfr_get_prec (x), mpfr_log_prec, x, n, rnd), 49 ("y[%Pd]=%.*Rg inexact=%d", 50 mpfr_get_prec (y), mpfr_log_prec, y, inexact)); 51 52 /* x^0 = 1 for any x, even a NaN */ 53 if (MPFR_UNLIKELY (n == 0)) 54 return mpfr_set_ui (y, 1, rnd); 55 56 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 57 { 58 if (MPFR_IS_NAN (x)) 59 { 60 MPFR_SET_NAN (y); 61 MPFR_RET_NAN; 62 } 63 else if (MPFR_IS_INF (x)) 64 { 65 /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */ 66 if (MPFR_IS_NEG (x) && (n & 1) == 1) 67 MPFR_SET_NEG (y); 68 else 69 MPFR_SET_POS (y); 70 MPFR_SET_INF (y); 71 MPFR_RET (0); 72 } 73 else /* x is zero */ 74 { 75 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 76 /* 0^n = 0 for any n */ 77 MPFR_SET_ZERO (y); 78 if (MPFR_IS_POS (x) || (n & 1) == 0) 79 MPFR_SET_POS (y); 80 else 81 MPFR_SET_NEG (y); 82 MPFR_RET (0); 83 } 84 } 85 else if (MPFR_UNLIKELY (n <= 2)) 86 { 87 if (n < 2) 88 /* x^1 = x */ 89 return mpfr_set (y, x, rnd); 90 else 91 /* x^2 = sqr(x) */ 92 return mpfr_sqr (y, x, rnd); 93 } 94 95 /* Augment exponent range */ 96 MPFR_SAVE_EXPO_MARK (expo); 97 98 for (m = n, nlen = 0; m != 0; nlen++, m >>= 1) 99 ; 100 /* 2^(nlen-1) <= n < 2^nlen */ 101 102 /* set up initial precision */ 103 prec = MPFR_PREC (y) + 3 + GMP_NUMB_BITS 104 + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)); 105 if (prec <= nlen) 106 prec = nlen + 1; 107 mpfr_init2 (res, prec); 108 109 rnd1 = MPFR_IS_POS (x) ? MPFR_RNDU : MPFR_RNDD; /* away */ 110 111 MPFR_ZIV_INIT (loop, prec); 112 for (;;) 113 { 114 int i; 115 116 MPFR_ASSERTD (prec > nlen); 117 err = prec - 1 - nlen; 118 /* First step: compute square from x */ 119 MPFR_BLOCK (flags, 120 inexact = mpfr_sqr (res, x, MPFR_RNDU); 121 MPFR_ASSERTD (nlen >= 2 && nlen <= INT_MAX); 122 i = nlen; 123 if (n & ((UTYPE) 1 << (i-2))) 124 inexact |= mpfr_mul (res, res, x, rnd1); 125 for (i -= 3; i >= 0 && !MPFR_BLOCK_EXCEP; i--) 126 { 127 inexact |= mpfr_sqr (res, res, MPFR_RNDU); 128 if (n & ((UTYPE) 1 << i)) 129 inexact |= mpfr_mul (res, res, x, rnd1); 130 }); 131 /* let r(n) be the number of roundings: we have r(2)=1, r(3)=2, 132 and r(2n)=2r(n)+1, r(2n+1)=2r(n)+2, thus r(n)=n-1. 133 Using Higham's method, to each rounding corresponds a factor 134 (1-theta) with 0 <= theta <= 2^(1-p), thus at the end the 135 absolute error is bounded by (n-1)*2^(1-p)*res <= 2*(n-1)*ulp(res) 136 since 2^(-p)*x <= ulp(x). Since n < 2^i, this gives a maximal 137 error of 2^(1+i)*ulp(res). 138 */ 139 if (MPFR_LIKELY (inexact == 0 140 || MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags) 141 || MPFR_CAN_ROUND (res, err, MPFR_PREC (y), rnd))) 142 break; 143 /* Actualisation of the precision */ 144 MPFR_ZIV_NEXT (loop, prec); 145 mpfr_set_prec (res, prec); 146 } 147 MPFR_ZIV_FREE (loop); 148 149 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags))) 150 { 151 mpz_t z; 152 153 /* Internal overflow or underflow. However, the approximation error has 154 * not been taken into account. So, let's solve this problem by using 155 * mpfr_pow_z, which can handle it. This case could be improved in the 156 * future, without having to use mpfr_pow_z. 157 */ 158 MPFR_LOG_MSG (("Internal overflow or underflow," 159 " let's use mpfr_pow_z.\n", 0)); 160 mpfr_clear (res); 161 MPFR_SAVE_EXPO_FREE (expo); 162 mpz_init (z); 163 MPZ_SET_U (z, n); 164 inexact = mpfr_pow_z (y, x, z, rnd); 165 mpz_clear (z); 166 return inexact; 167 } 168 169 inexact = mpfr_set (y, res, rnd); 170 mpfr_clear (res); 171 172 MPFR_SAVE_EXPO_FREE (expo); 173 return mpfr_check_range (y, inexact, rnd); 174 } 175