1/* mpfr_ui_pow_ui -- compute the power between two machine integers 2 3Copyright 1999-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include "mpfr-impl.h" 24 25int 26mpfr_ui_pow_ui (mpfr_ptr x, unsigned long int k, unsigned long int n, 27 mpfr_rnd_t rnd) 28{ 29 mpfr_exp_t err; 30 unsigned long m; 31 mpfr_t res; 32 mpfr_prec_t prec; 33 int size_n; 34 int inexact; 35 MPFR_ZIV_DECL (loop); 36 MPFR_SAVE_EXPO_DECL (expo); 37 38 MPFR_LOG_FUNC 39 (("k=%lu n=%lu rnd=%d", k, n, rnd), 40 ("y[%Pd]=%.*Rg inexact=%d", 41 mpfr_get_prec (x), mpfr_log_prec, x, inexact)); 42 43 if (MPFR_UNLIKELY (n <= 1)) 44 { 45 if (n == 1) 46 return mpfr_set_ui (x, k, rnd); /* k^1 = k */ 47 else 48 return mpfr_set_ui (x, 1, rnd); /* k^0 = 1 for any k */ 49 } 50 else if (MPFR_UNLIKELY (k <= 1)) 51 { 52 if (k == 1) 53 return mpfr_set_ui (x, 1, rnd); /* 1^n = 1 for any n > 0 */ 54 else 55 return mpfr_set_ui (x, 0, rnd); /* 0^n = 0 for any n > 0 */ 56 } 57 58 for (size_n = 0, m = n; m != 0; size_n++, m >>= 1) 59 ; 60 61 MPFR_SAVE_EXPO_MARK (expo); 62 prec = MPFR_PREC (x) + 3 + size_n; 63 mpfr_init2 (res, prec); 64 65 MPFR_ZIV_INIT (loop, prec); 66 for (;;) 67 { 68 int i = size_n; 69 unsigned int inex_res; 70 71 inex_res = mpfr_set_ui (res, k, MPFR_RNDU); 72 err = 1; 73 /* now 2^(i-1) <= n < 2^i: i=1+floor(log2(n)) */ 74 for (i -= 2; i >= 0; i--) 75 { 76 inex_res |= mpfr_sqr (res, res, MPFR_RNDU); 77 err++; 78 if (n & (1UL << i)) 79 inex_res |= mpfr_mul_ui (res, res, k, MPFR_RNDU); 80 } 81 82 if (MPFR_UNLIKELY (MPFR_IS_INF (res))) 83 { 84 mpfr_t kf; 85 mpz_t z; 86 int size_k; 87 MPFR_BLOCK_DECL (flags); 88 89 /* Let's handle the overflow by calling mpfr_pow_z. 90 Alternatively, we could call mpfr_pow_ui; this would 91 need a bit shorter code below, but mpfr_pow_ui handles 92 the overflow by calling mpfr_pow_z, so that calling 93 mpfr_pow_z directly should be a bit more efficient. */ 94 95 MPFR_ZIV_FREE (loop); 96 mpfr_clear (res); 97 for (size_k = 0, m = k; m != 0; size_k++, m >>= 1) 98 ; 99 mpfr_init2 (kf, size_k); 100 inexact = mpfr_set_ui (kf, k, MPFR_RNDN); 101 MPFR_ASSERTD (inexact == 0); 102 mpz_init (z); 103 mpz_set_ui (z, n); 104 MPFR_BLOCK (flags, inexact = mpfr_pow_z (x, kf, z, rnd);); 105 mpz_clear (z); 106 mpfr_clear (kf); 107 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, flags); 108 goto end; 109 } 110 111 /* since the loop is executed floor(log2(n)) times, 112 we have err = 1+floor(log2(n)). 113 Since prec >= MPFR_PREC(x) + 4 + floor(log2(n)), prec > err */ 114 err = prec - err; 115 116 MPFR_LOG_VAR (res); 117 if (MPFR_LIKELY (!inex_res 118 || MPFR_CAN_ROUND (res, err, MPFR_PREC (x), rnd))) 119 break; 120 121 /* Actualisation of the precision */ 122 MPFR_ZIV_NEXT (loop, prec); 123 mpfr_set_prec (res, prec); 124 } 125 MPFR_ZIV_FREE (loop); 126 127 inexact = mpfr_set (x, res, rnd); 128 129 mpfr_clear (res); 130 131 end: 132 MPFR_SAVE_EXPO_FREE (expo); 133 return mpfr_check_range (x, inexact, rnd); 134} 135