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
mpfr_exp2m1_small(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode,mpfr_ptr t)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
mpfr_exp2m1(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)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