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[%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 }