1 /* mpfr_exp2m1 -- Compute 2^x-1 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 /* The computation of exp2m1 is done by expm1(x) = 2^x-1 */ 27 28 /* In case x is small in absolute value, 2^x - 1 ~ x*log(2). 29 If this is enough to deduce correct rounding, put in the auxiliary variable 30 t the approximation that will be rounded to get y, and return non-zero. 31 If we put 0 in t, it means underflow. 32 Otherwise return 0. */ 33 static int 34 mpfr_exp2m1_small (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode, mpfr_ptr t) 35 { 36 mpfr_prec_t prec; 37 mpfr_exp_t e; 38 MPFR_BLOCK_DECL (flags); 39 40 /* for |x| < 0.125, we have |2^x-1-x*log(2)| < x^2/4 */ 41 if (MPFR_EXP(x) > -3) 42 return 0; 43 /* now EXP(x) <= -3, thus x < 0.125 */ 44 prec = MPFR_PREC(t); 45 mpfr_const_log2 (t, MPFR_RNDN); 46 /* t = log(2)*(1 + theta) with |theta| <= 2^(-prec) */ 47 MPFR_BLOCK (flags, mpfr_mul (t, t, x, MPFR_RNDN)); 48 /* If an underflow occurs in log(2)*x, then return underflow. */ 49 if (MPFR_UNDERFLOW (flags)) 50 { 51 MPFR_SET_ZERO (t); 52 return 1; 53 } 54 /* t = x*log(2)*(1 + theta)^2 with |theta| <= 2^(-prec) */ 55 /* |t - x*log(2)| <= ((1 + theta)^2 - 1) * |t| <= 3*2^(-prec)*|t| */ 56 /* |t - x*log(2)| < 3*2^(EXP(t)-prec) */ 57 e = 2 * MPFR_GET_EXP (x) - 2 + prec - MPFR_GET_EXP(t); 58 /* |x^2/4| < 2^e*2^(EXP(t)-prec) thus 59 |t - exp2m1(x)| < (3+2^e)*2^(EXP(t)-prec) */ 60 e = (e <= 1) ? 2 + (e == 1) : e + 1; 61 /* now |t - exp2m1(x)| < 2^e*2^(EXP(t)-prec) */ 62 return MPFR_CAN_ROUND (t, prec - e, MPFR_PREC(y), rnd_mode); 63 } 64 65 int 66 mpfr_exp2m1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 67 { 68 int inexact, nloop; 69 mpfr_t t; 70 mpfr_prec_t Ny = MPFR_PREC(y); /* target precision */ 71 mpfr_prec_t Nt; /* working precision */ 72 mpfr_exp_t err, exp_te; /* error */ 73 MPFR_ZIV_DECL (loop); 74 MPFR_SAVE_EXPO_DECL (expo); 75 76 MPFR_LOG_FUNC 77 (("x[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), 78 ("y[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, 79 inexact)); 80 81 if (MPFR_IS_SINGULAR (x)) 82 return mpfr_expm1 (y, x, rnd_mode); /* singular cases are identical */ 83 84 MPFR_ASSERTN(!MPFR_IS_ZERO(x)); 85 86 MPFR_SAVE_EXPO_MARK (expo); 87 88 /* Check case where result is -1 or nextabove(-1) because x is a huge 89 negative number. */ 90 if (MPFR_IS_NEG(x) && mpfr_cmpabs_ui (x, MPFR_PREC(y) + 1) > 0) 91 { 92 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_INEXACT); 93 /* 1/2*ulp(-1) = 2^(-PREC(y)) thus 2^x < 1/4*ulp(-1): 94 result is -1 for RNDA,RNDD,RNDN, and nextabove(-1) for RNDZ,RNDU */ 95 mpfr_set_si (y, -1, MPFR_RNDZ); 96 if (!MPFR_IS_LIKE_RNDZ(rnd_mode,1)) 97 inexact = -1; 98 else 99 { 100 mpfr_nextabove (y); 101 inexact = 1; 102 } 103 goto end; 104 } 105 106 /* compute the precision of intermediary variable */ 107 /* the optimal number of bits : see algorithms.tex */ 108 Nt = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 6; 109 110 mpfr_init2 (t, Nt); 111 112 MPFR_ZIV_INIT (loop, Nt); 113 for (nloop = 0;; nloop++) 114 { 115 int inex1; 116 117 MPFR_BLOCK_DECL (flags); 118 119 /* 2^x may overflow and underflow */ 120 MPFR_BLOCK (flags, inex1 = mpfr_exp2 (t, x, MPFR_RNDN)); 121 122 if (MPFR_OVERFLOW (flags)) /* overflow case */ 123 { 124 inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN_POS); 125 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW); 126 goto clear; 127 } 128 129 /* integer case */ 130 if (inex1 == 0) 131 { 132 inexact = mpfr_sub_ui (y, t, 1, rnd_mode); 133 goto clear; 134 } 135 136 /* To get an underflow in 2^x, we need 2^x < 0.5*2^MPFR_EMIN_MIN 137 thus x < MPFR_EMIN_MIN-1. But in that case (huge negative x) 138 was already detected before Ziv's loop. */ 139 MPFR_ASSERTD(!MPFR_UNDERFLOW (flags)); 140 141 MPFR_ASSERTN(!MPFR_IS_ZERO(t)); 142 exp_te = MPFR_GET_EXP (t); 143 mpfr_sub_ui (t, t, 1, MPFR_RNDN); /* 2^x-1 */ 144 145 /* error estimate */ 146 /* err = __gmpfr_ceil_log2(1+pow(2,MPFR_EXP(te)-MPFR_EXP(t))) */ 147 if (!MPFR_IS_ZERO(t)) 148 { 149 err = MAX (exp_te - MPFR_GET_EXP (t), 0) + 1; 150 /* if inex1=0, this means that t=o(2^x) is exact, thus the correct 151 rounding is simply o(t-1) */ 152 if (inex1 == 0 || 153 MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - err, Ny, rnd_mode))) 154 break; 155 } 156 157 /* check small case: we need to do it at each step of Ziv's loop, 158 since the multiplication x*log(2) might not enable correct 159 rounding at the first loop */ 160 if (mpfr_exp2m1_small (y, x, rnd_mode, t)) 161 { 162 if (MPFR_IS_ZERO(t)) /* underflow */ 163 { 164 mpfr_clear (t); 165 MPFR_SAVE_EXPO_FREE (expo); 166 return mpfr_underflow (y, (rnd_mode == MPFR_RNDN) ? MPFR_RNDZ 167 : rnd_mode, 1); 168 } 169 break; 170 } 171 172 /* increase the precision */ 173 MPFR_ZIV_NEXT (loop, Nt); 174 mpfr_set_prec (t, Nt); 175 } 176 177 inexact = mpfr_set (y, t, rnd_mode); 178 clear: 179 MPFR_ZIV_FREE (loop); 180 mpfr_clear (t); 181 182 end: 183 MPFR_SAVE_EXPO_FREE (expo); 184 return mpfr_check_range (y, inexact, rnd_mode); 185 } 186
Indexes created Fri Jun 05 00:26:10 UTC 2026