1 /* mpfr_exp2 -- power of 2 function 2^y 2 3 Copyright 2001-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 /* TODO: mpfr_get_exp_t is called 3 times, with 3 different directed 27 rounding modes. One could reduce it to only one call thanks to the 28 inexact flag, but is it worth? */ 29 30 /* Convert x to an mpfr_eexp_t integer, with saturation at the minimum 31 and maximum values. Flags are unchanged. */ 32 static mpfr_eexp_t 33 round_to_eexp_t (mpfr_srcptr x, mpfr_rnd_t rnd_mode) 34 { 35 mpfr_flags_t flags = __gmpfr_flags; 36 mpfr_eexp_t e; 37 38 e = mpfr_get_exp_t (x, rnd_mode); 39 __gmpfr_flags = flags; 40 return e; 41 } 42 43 /* The computation of y = 2^z is done by * 44 * y = exp(z*log(2)). The result is exact iff z is an integer. */ 45 46 int 47 mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 48 { 49 int inexact; 50 mpfr_eexp_t xint; /* note: will fit in mpfr_exp_t */ 51 mpfr_t xfrac; 52 MPFR_SAVE_EXPO_DECL (expo); 53 54 MPFR_LOG_FUNC 55 (("x[%Pd]=%.*Rg rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode), 56 ("y[%Pd]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y, 57 inexact)); 58 59 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 60 { 61 if (MPFR_IS_NAN (x)) 62 { 63 MPFR_SET_NAN (y); 64 MPFR_RET_NAN; 65 } 66 else if (MPFR_IS_INF (x)) 67 { 68 if (MPFR_IS_POS (x)) 69 MPFR_SET_INF (y); 70 else 71 MPFR_SET_ZERO (y); 72 MPFR_SET_POS (y); 73 MPFR_RET (0); 74 } 75 else /* 2^0 = 1 */ 76 { 77 MPFR_ASSERTD (MPFR_IS_ZERO(x)); 78 return mpfr_set_ui (y, 1, rnd_mode); 79 } 80 } 81 82 /* Since the smallest representable non-zero float is 1/2 * 2^emin, 83 if x <= emin - 2, the result is either 1/2 * 2^emin or 0. 84 Warning, for emin - 2 < x < emin - 1, we cannot conclude, since 2^x 85 might round to 2^(emin - 1) for rounding away or to nearest, and there 86 might be no underflow, since we consider underflow "after rounding". */ 87 88 if (MPFR_UNLIKELY (round_to_eexp_t (x, MPFR_RNDU) <= __gmpfr_emin - 2)) 89 return mpfr_underflow (y, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode, 1); 90 91 if (MPFR_UNLIKELY (round_to_eexp_t (x, MPFR_RNDD) >= __gmpfr_emax)) 92 return mpfr_overflow (y, rnd_mode, 1); 93 94 /* We now know that emin - 2 < x < emax. Note that an underflow or 95 overflow is still possible (we have eliminated only easy cases). */ 96 97 MPFR_SAVE_EXPO_MARK (expo); 98 99 /* 2^x = 1 + x*log(2) + O(x^2) for x near zero, and for |x| <= 1 we have 100 |2^x - 1| <= x < 2^EXP(x). If x > 0 we must round away from 0 (dir=1); 101 if x < 0 we must round toward 0 (dir=0). */ 102 MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, __gmpfr_one, - MPFR_GET_EXP (x), 0, 103 MPFR_IS_POS (x), rnd_mode, expo, {}); 104 105 xint = mpfr_get_exp_t (x, MPFR_RNDZ); 106 MPFR_ASSERTD (__gmpfr_emin - 2 < xint && xint < __gmpfr_emax); 107 108 mpfr_init2 (xfrac, MPFR_PREC (x)); 109 MPFR_DBGRES (inexact = mpfr_frac (xfrac, x, MPFR_RNDN)); 110 MPFR_ASSERTD (inexact == 0); 111 112 if (MPFR_IS_ZERO (xfrac)) 113 { 114 /* Here, emin - 1 <= x <= emax - 1, so that an underflow or overflow 115 will not be possible. */ 116 mpfr_set_ui (y, 1, MPFR_RNDN); 117 inexact = 0; 118 } 119 else 120 { 121 /* Declaration of the intermediary variable */ 122 mpfr_t t; 123 124 /* Declaration of the size variable */ 125 mpfr_prec_t Ny = MPFR_PREC(y); /* target precision */ 126 mpfr_prec_t Nt; /* working precision */ 127 mpfr_exp_t err; /* error */ 128 MPFR_ZIV_DECL (loop); 129 130 /* compute the precision of intermediary variable */ 131 /* the optimal number of bits : see algorithms.tex */ 132 Nt = Ny + 5 + MPFR_INT_CEIL_LOG2 (Ny); 133 134 /* initialize of intermediary variable */ 135 mpfr_init2 (t, Nt); 136 137 /* First computation */ 138 MPFR_ZIV_INIT (loop, Nt); 139 for (;;) 140 { 141 /* compute exp(x*ln(2))*/ 142 mpfr_const_log2 (t, MPFR_RNDU); /* ln(2) */ 143 mpfr_mul (t, xfrac, t, MPFR_RNDU); /* xfrac * ln(2) */ 144 err = Nt - (MPFR_GET_EXP (t) + 2); /* Estimate of the error */ 145 mpfr_exp (t, t, MPFR_RNDN); /* exp(xfrac * ln(2)) */ 146 147 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) 148 break; 149 150 /* Actualisation of the precision */ 151 MPFR_ZIV_NEXT (loop, Nt); 152 mpfr_set_prec (t, Nt); 153 } 154 MPFR_ZIV_FREE (loop); 155 156 inexact = mpfr_set (y, t, rnd_mode); 157 158 mpfr_clear (t); 159 } 160 161 mpfr_clear (xfrac); 162 163 if (MPFR_UNLIKELY (rnd_mode == MPFR_RNDN && xint == __gmpfr_emin - 1 && 164 MPFR_GET_EXP (y) == 0 && mpfr_powerof2_raw (y))) 165 { 166 /* y was rounded down to 1/2 and the rounded value with an unbounded 167 exponent range would be 2^(emin-2), i.e. the midpoint between 0 168 and the smallest positive FP number. This is a double rounding 169 problem: we should not round to 0, but to (1/2) * 2^emin. */ 170 MPFR_SET_EXP (y, __gmpfr_emin); 171 inexact = 1; 172 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_UNDERFLOW); 173 } 174 else 175 { 176 /* The following is OK due to early overflow/underflow checking. 177 the exponent may be slightly out-of-range, but this will be 178 handled by mpfr_check_range. */ 179 MPFR_EXP (y) += xint; 180 } 181 182 MPFR_SAVE_EXPO_FREE (expo); 183 return mpfr_check_range (y, inexact, rnd_mode); 184 } 185
Indexes created Tue Sep 01 00:24:37 UTC 2026