1 /* mpfr_compound_si --- compound(x,n) = (1+x)^n
2
3 Copyright 2021-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 /* needed for MPFR_INT_CEIL_LOG2 */
24 #include "mpfr-impl.h"
25
26 /* assuming |(1+x)^n - 1| < 1/4*ulp(1), return correct rounding,
27 where s is the sign of n*log2(1+x) */
28 static int
mpfr_compound_near_one(mpfr_ptr y,int s,mpfr_rnd_t rnd_mode)29 mpfr_compound_near_one (mpfr_ptr y, int s, mpfr_rnd_t rnd_mode)
30 {
31 mpfr_set_ui (y, 1, rnd_mode); /* exact */
32 if (rnd_mode == MPFR_RNDN || rnd_mode == MPFR_RNDF
33 || (s > 0 && (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDD))
34 || (s < 0 && (rnd_mode == MPFR_RNDA || rnd_mode == MPFR_RNDU)))
35 {
36 /* round toward 1 */
37 return -s;
38 }
39 else if (s > 0) /* necessarily RNDA or RNDU */
40 {
41 /* round toward +Inf */
42 mpfr_nextabove (y);
43 return +1;
44 }
45 else /* necessarily s < 0 and RNDZ or RNDD */
46 {
47 /* round toward 0 */
48 mpfr_nextbelow (y);
49 return -1;
50 }
51 }
52
53 /* put in y the correctly rounded value of (1+x)^n */
54 int
mpfr_compound_si(mpfr_ptr y,mpfr_srcptr x,long n,mpfr_rnd_t rnd_mode)55 mpfr_compound_si (mpfr_ptr y, mpfr_srcptr x, long n, mpfr_rnd_t rnd_mode)
56 {
57 int inexact, compared, k, nloop;
58 mpfr_t t, u;
59 mpfr_prec_t py, prec, extra;
60 mpfr_rnd_t rnd1;
61 MPFR_ZIV_DECL (loop);
62 MPFR_SAVE_EXPO_DECL (expo);
63
64 MPFR_LOG_FUNC
65 (("x[%Pd]=%.*Rg n=%ld rnd=%d",
66 mpfr_get_prec(x), mpfr_log_prec, x, n, rnd_mode),
67 ("y[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inexact));
68
69 /* Special cases */
70 if (MPFR_IS_SINGULAR (x))
71 {
72 if (MPFR_IS_INF (x) && MPFR_IS_NEG (x))
73 {
74 /* compound(-Inf,n) is NaN */
75 MPFR_SET_NAN (y);
76 MPFR_RET_NAN;
77 }
78 else if (n == 0 || MPFR_IS_ZERO (x))
79 {
80 /* compound(x,0) = 1 for x >= -1 or NaN (the only special value
81 of x that is not concerned is -Inf, already handled);
82 compound(0,n) = 1 */
83 return mpfr_set_ui (y, 1, rnd_mode);
84 }
85 else if (MPFR_IS_NAN (x))
86 {
87 /* compound(NaN,n) is NaN, except for n = 0, already handled. */
88 MPFR_SET_NAN (y);
89 MPFR_RET_NAN;
90 }
91 else if (MPFR_IS_INF (x)) /* x = +Inf */
92 {
93 MPFR_ASSERTD (MPFR_IS_POS (x));
94 if (n < 0) /* (1+Inf)^n = +0 for n < 0 */
95 MPFR_SET_ZERO (y);
96 else /* n > 0: (1+Inf)^n = +Inf */
97 MPFR_SET_INF (y);
98 MPFR_SET_POS (y);
99 MPFR_RET (0); /* exact 0 or infinity */
100 }
101 }
102
103 /* (1+x)^n = NaN for x < -1 */
104 compared = mpfr_cmp_si (x, -1);
105 if (compared < 0)
106 {
107 MPFR_SET_NAN (y);
108 MPFR_RET_NAN;
109 }
110
111 /* compound(x,0) gives 1 for x >= 1 */
112 if (n == 0)
113 return mpfr_set_ui (y, 1, rnd_mode);
114
115 if (compared == 0)
116 {
117 if (n < 0)
118 {
119 /* compound(-1,n) = +Inf with divide-by-zero exception */
120 MPFR_SET_INF (y);
121 MPFR_SET_POS (y);
122 MPFR_SET_DIVBY0 ();
123 MPFR_RET (0);
124 }
125 else
126 {
127 /* compound(-1,n) = +0 */
128 MPFR_SET_ZERO (y);
129 MPFR_SET_POS (y);
130 MPFR_RET (0);
131 }
132 }
133
134 if (n == 1)
135 return mpfr_add_ui (y, x, 1, rnd_mode);
136
137 MPFR_SAVE_EXPO_MARK (expo);
138
139 py = MPFR_GET_PREC (y);
140 prec = py + MPFR_INT_CEIL_LOG2 (py) + 6;
141
142 mpfr_init2 (t, prec);
143 mpfr_init2 (u, prec);
144
145 k = MPFR_INT_CEIL_LOG2(SAFE_ABS (unsigned long, n)); /* thus |n| <= 2^k */
146
147 /* We compute u=log2p1(x) with prec+extra bits, since we lose some bits
148 in 2^u. */
149 extra = 0;
150 rnd1 = VSIGN (n) == MPFR_SIGN (x) ? MPFR_RNDD : MPFR_RNDU;
151
152 MPFR_ZIV_INIT (loop, prec);
153 for (nloop = 0; ; nloop++)
154 {
155 unsigned int inex;
156 mpfr_exp_t e, e2, ex;
157 mpfr_prec_t precu = MPFR_ADD_PREC (prec, extra);
158 mpfr_prec_t new_extra;
159 mpfr_rnd_t rnd2;
160
161 /* We compute (1+x)^n as 2^(n*log2p1(x)),
162 and we round toward 1, thus we round n*log2p1(x) toward 0,
163 thus for x*n > 0 we round log2p1(x) toward -Inf, and for x*n < 0
164 we round log2p1(x) toward +Inf. */
165 inex = mpfr_log2p1 (u, x, rnd1) != 0;
166 e = MPFR_GET_EXP (u);
167 /* |u - log2(1+x)| <= ulp(t) = 2^(e-precu) */
168 inex |= mpfr_mul_si (u, u, n, MPFR_RNDZ) != 0;
169 e2 = MPFR_GET_EXP (u);
170 /* |u - n*log2(1+x)| <= 2^(e2-precu) + |n|*2^(e-precu)
171 <= 2^(e2-precu) + 2^(e+k-precu) <= 2^(e+k+1-precu)
172 where |n| <= 2^k, and e2 is the new exponent of u. */
173 MPFR_ASSERTD (e2 <= e + k);
174 e += k + 1;
175 MPFR_ASSERTN (e2 <= MPFR_PREC_MAX);
176 new_extra = e2 > 0 ? e2 : 0;
177 /* |u - n*log2(1+x)| <= 2^(e-precu) */
178 /* detect overflow: since we rounded n*log2p1(x) toward 0,
179 if n*log2p1(x) >= __gmpfr_emax, we are sure there is overflow. */
180 if (mpfr_cmp_si (u, __gmpfr_emax) >= 0)
181 {
182 MPFR_ZIV_FREE (loop);
183 mpfr_clear (t);
184 mpfr_clear (u);
185 MPFR_SAVE_EXPO_FREE (expo);
186 return mpfr_overflow (y, rnd_mode, 1);
187 }
188 /* detect underflow: similarly, since we rounded n*log2p1(x) toward 0,
189 if n*log2p1(x) < __gmpfr_emin-1, we are sure there is underflow. */
190 if (mpfr_cmp_si (u, __gmpfr_emin - 1) < 0)
191 {
192 MPFR_ZIV_FREE (loop);
193 mpfr_clear (t);
194 mpfr_clear (u);
195 MPFR_SAVE_EXPO_FREE (expo);
196 return mpfr_underflow (y,
197 rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode, 1);
198 }
199 /* Detect cases where result is 1 or 1+ulp(1) or 1-1/2*ulp(1):
200 |2^u - 1| = |exp(u*log(2)) - 1| <= |u|*log(2) < |u| */
201 if (nloop == 0 && MPFR_GET_EXP(u) < - py)
202 {
203 /* since ulp(1) = 2^(1-py), we have |u| < 1/4*ulp(1) */
204 /* mpfr_compound_near_one must be called in the extended
205 exponent range, so that 1 is representable. */
206 inexact = mpfr_compound_near_one (y, MPFR_SIGN (u), rnd_mode);
207 goto end;
208 }
209 /* round 2^u toward 1 */
210 rnd2 = MPFR_IS_POS (u) ? MPFR_RNDD : MPFR_RNDU;
211 inex |= mpfr_exp2 (t, u, rnd2) != 0;
212 /* we had |u - n*log2(1+x)| < 2^(e-precu)
213 thus u = n*log2(1+x) + delta with |delta| < 2^(e-precu)
214 then 2^u = (1+x)^n * 2^delta with |delta| < 2^(e-precu).
215 For |delta| < 0.5, |2^delta - 1| <= |delta| thus
216 |t - (1+x)^n| <= ulp(t) + |t|*2^(e-precu)
217 < 2^(EXP(t)-prec) + 2^(EXP(t)+e-precu) */
218 e = (precu - prec >= e) ? 1 : e + 1 - (precu - prec);
219 /* now |t - (1+x)^n| < 2^(EXP(t)+e-prec) */
220
221 if (MPFR_LIKELY (!inex || MPFR_CAN_ROUND (t, prec - e, py, rnd_mode)))
222 break;
223
224 /* If t fits in the target precision (or with 1 more bit), then we can
225 round, assuming the working precision is large enough, but the above
226 MPFR_CAN_ROUND() will fail because we cannot determine the ternary
227 value. However since we rounded t toward 1, we can determine it.
228 Since the error in the approximation t is at most 2^e ulp(t),
229 this error should be less than 1/2 ulp(y), thus we should have
230 prec - py >= e + 1. */
231 if (mpfr_min_prec (t) <= py + 1 && prec - py >= e + 1)
232 {
233 /* we add/subtract one ulp to get the correct rounding */
234 if (rnd2 == MPFR_RNDD) /* t was rounded downwards */
235 mpfr_nextabove (t);
236 else
237 mpfr_nextbelow (t);
238 break;
239 }
240
241 /* Detect particular cases where Ziv's strategy may take too much
242 memory and be too long, i.e. when x^n fits in the target precision
243 (+ 1 additional bit for rounding to nearest) and the exact result
244 (1+x)^n is very close to x^n.
245 Necessarily, x is a large even integer and n > 0 (thus n > 1).
246 Since this does not depend on the working precision, we only
247 check this at the first iteration (nloop == 0).
248 Hence the first "if" below and the kx < ex test of the second "if"
249 (x is an even integer iff its least bit 1 has exponent >= 1).
250 The second test of the second "if" corresponds to another simple
251 condition that implies that x^n fits in the target precision.
252 Here are the details:
253 Let k be the minimum length of the significand of x, and x' the odd
254 (integer) significand of x. This means that 2^(k-1) <= x' < 2^k.
255 Thus 2^(n*(k-1)) <= (x')^n < 2^(k*n), and x^n has between n*(k-1)+1
256 and k*n bits. So x^n can fit into p bits only if p >= n*(k-1)+1,
257 i.e. n*(k-1) <= p-1.
258 Note that x >= 2^k, so that x^n >= 2^(k*n). Since raw overflow
259 has already been detected, k*n cannot overflow if computed with
260 the mpfr_exp_t type. Hence the second test of the second "if",
261 which cannot overflow. */
262 MPFR_ASSERTD (n < 0 || n > 1);
263 if (nloop == 0 && n > 1 && (ex = MPFR_GET_EXP (x)) >= 17)
264 {
265 mpfr_prec_t kx = mpfr_min_prec (x);
266 mpfr_prec_t p = py + (rnd_mode == MPFR_RNDN);
267
268 MPFR_LOG_MSG (("Check if x^n fits... n=%ld kx=%Pd p=%Pd\n",
269 n, kx, p));
270 if (kx < ex && n * (mpfr_exp_t) (kx - 1) <= p - 1)
271 {
272 mpfr_t v;
273
274 /* Check whether x^n really fits into p bits. */
275 mpfr_init2 (v, p);
276 inexact = mpfr_pow_ui (v, x, n, MPFR_RNDZ);
277 if (inexact == 0)
278 {
279 MPFR_LOG_MSG (("x^n fits into p bits\n", 0));
280 /* (x+1)^n = x^n * (1 + 1/x)^n
281 For directed rounding, we can round when (1 + 1/x)^n
282 < 1 + 2^-p, and then the result is x^n,
283 except for rounding up. Indeed, if (1 + 1/x)^n < 1 + 2^-p,
284 1 <= (x+1)^n < x^n * (1 + 2^-p) = x^n + x^n/2^p
285 < x^n + ulp(x^n).
286 For rounding to nearest, we can round when (1 + 1/x)^n
287 < 1 + 2^-p, and then the result is x^n when x^n fits
288 into p-1 bits, and nextabove(x^n) otherwise. */
289 mpfr_ui_div (t, 1, x, MPFR_RNDU);
290 mpfr_add_ui (t, t, 1, MPFR_RNDU);
291 mpfr_pow_ui (t, t, n, MPFR_RNDU);
292 mpfr_sub_ui (t, t, 1, MPFR_RNDU);
293 /* t cannot be zero */
294 if (MPFR_GET_EXP(t) < - py)
295 {
296 mpfr_set (y, v, MPFR_RNDZ);
297 if ((rnd_mode == MPFR_RNDN && mpfr_min_prec (v) == p)
298 || rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA)
299 {
300 /* round up */
301 mpfr_nextabove (y);
302 inexact = 1;
303 }
304 else
305 inexact = -1;
306 mpfr_clear (v);
307 goto end;
308 }
309 }
310 mpfr_clear (v);
311 }
312 }
313
314 /* Exact cases like compound(0.5,2) = 9/4 must be detected, since
315 except for 1+x power of 2, the log2p1 above will be inexact,
316 so that in the Ziv test, inexact != 0 and MPFR_CAN_ROUND will
317 fail (even for RNDN, as the ternary value cannot be determined),
318 yielding an infinite loop.
319 For an exact case in precision prec(y), 1+x will necessarily
320 be exact in precision prec(y), thus also in prec(t), where
321 prec(t) >= prec(y), and we can use mpfr_pow_si under this
322 condition (which will also evaluate some non-exact cases). */
323 if (mpfr_add_ui (t, x, 1, MPFR_RNDZ) == 0)
324 {
325 inexact = mpfr_pow_si (y, t, n, rnd_mode);
326 goto end;
327 }
328
329 MPFR_ZIV_NEXT (loop, prec);
330 mpfr_set_prec (t, prec);
331 extra = new_extra;
332 mpfr_set_prec (u, MPFR_ADD_PREC (prec, extra));
333 }
334
335 inexact = mpfr_set (y, t, rnd_mode);
336
337 end:
338 MPFR_ZIV_FREE (loop);
339 mpfr_clear (t);
340 mpfr_clear (u);
341
342 MPFR_SAVE_EXPO_FREE (expo);
343 return mpfr_check_range (y, inexact, rnd_mode);
344 }
345