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