1 /* mpfr_exp10m1 -- Compute 10^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 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. */ 33 static int 34 mpfr_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 60 int 61 mpfr_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[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), 73 ("y[%Pu]=%.*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