1 /* mpfr_exp2 -- power of 2 function 2^y
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 /* TODO: mpfr_get_exp_t is called 3 times, with 3 different directed
27 rounding modes. One could reduce it to only one call thanks to the
28 inexact flag, but is it worth? */
29
30 /* Convert x to an mpfr_eexp_t integer, with saturation at the minimum
31 and maximum values. Flags are unchanged. */
32 static mpfr_eexp_t
round_to_eexp_t(mpfr_srcptr x,mpfr_rnd_t rnd_mode)33 round_to_eexp_t (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
34 {
35 mpfr_flags_t flags = __gmpfr_flags;
36 mpfr_eexp_t e;
37
38 e = mpfr_get_exp_t (x, rnd_mode);
39 __gmpfr_flags = flags;
40 return e;
41 }
42
43 /* The computation of y = 2^z is done by *
44 * y = exp(z*log(2)). The result is exact iff z is an integer. */
45
46 int
mpfr_exp2(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)47 mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
48 {
49 int inexact;
50 mpfr_eexp_t xint; /* note: will fit in mpfr_exp_t */
51 mpfr_t xfrac;
52 MPFR_SAVE_EXPO_DECL (expo);
53
54 MPFR_LOG_FUNC
55 (("x[%Pd]=%.*Rg rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode),
56 ("y[%Pd]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y,
57 inexact));
58
59 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
60 {
61 if (MPFR_IS_NAN (x))
62 {
63 MPFR_SET_NAN (y);
64 MPFR_RET_NAN;
65 }
66 else if (MPFR_IS_INF (x))
67 {
68 if (MPFR_IS_POS (x))
69 MPFR_SET_INF (y);
70 else
71 MPFR_SET_ZERO (y);
72 MPFR_SET_POS (y);
73 MPFR_RET (0);
74 }
75 else /* 2^0 = 1 */
76 {
77 MPFR_ASSERTD (MPFR_IS_ZERO(x));
78 return mpfr_set_ui (y, 1, rnd_mode);
79 }
80 }
81
82 /* Since the smallest representable non-zero float is 1/2 * 2^emin,
83 if x <= emin - 2, the result is either 1/2 * 2^emin or 0.
84 Warning, for emin - 2 < x < emin - 1, we cannot conclude, since 2^x
85 might round to 2^(emin - 1) for rounding away or to nearest, and there
86 might be no underflow, since we consider underflow "after rounding". */
87
88 if (MPFR_UNLIKELY (round_to_eexp_t (x, MPFR_RNDU) <= __gmpfr_emin - 2))
89 return mpfr_underflow (y, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode, 1);
90
91 if (MPFR_UNLIKELY (round_to_eexp_t (x, MPFR_RNDD) >= __gmpfr_emax))
92 return mpfr_overflow (y, rnd_mode, 1);
93
94 /* We now know that emin - 2 < x < emax. Note that an underflow or
95 overflow is still possible (we have eliminated only easy cases). */
96
97 MPFR_SAVE_EXPO_MARK (expo);
98
99 /* 2^x = 1 + x*log(2) + O(x^2) for x near zero, and for |x| <= 1 we have
100 |2^x - 1| <= x < 2^EXP(x). If x > 0 we must round away from 0 (dir=1);
101 if x < 0 we must round toward 0 (dir=0). */
102 MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, __gmpfr_one, - MPFR_GET_EXP (x), 0,
103 MPFR_IS_POS (x), rnd_mode, expo, {});
104
105 xint = mpfr_get_exp_t (x, MPFR_RNDZ);
106 MPFR_ASSERTD (__gmpfr_emin - 2 < xint && xint < __gmpfr_emax);
107
108 mpfr_init2 (xfrac, MPFR_PREC (x));
109 MPFR_DBGRES (inexact = mpfr_frac (xfrac, x, MPFR_RNDN));
110 MPFR_ASSERTD (inexact == 0);
111
112 if (MPFR_IS_ZERO (xfrac))
113 {
114 /* Here, emin - 1 <= x <= emax - 1, so that an underflow or overflow
115 will not be possible. */
116 mpfr_set_ui (y, 1, MPFR_RNDN);
117 inexact = 0;
118 }
119 else
120 {
121 /* Declaration of the intermediary variable */
122 mpfr_t t;
123
124 /* Declaration of the size variable */
125 mpfr_prec_t Ny = MPFR_PREC(y); /* target precision */
126 mpfr_prec_t Nt; /* working precision */
127 mpfr_exp_t err; /* error */
128 MPFR_ZIV_DECL (loop);
129
130 /* compute the precision of intermediary variable */
131 /* the optimal number of bits : see algorithms.tex */
132 Nt = Ny + 5 + MPFR_INT_CEIL_LOG2 (Ny);
133
134 /* initialize of intermediary variable */
135 mpfr_init2 (t, Nt);
136
137 /* First computation */
138 MPFR_ZIV_INIT (loop, Nt);
139 for (;;)
140 {
141 /* compute exp(x*ln(2))*/
142 mpfr_const_log2 (t, MPFR_RNDU); /* ln(2) */
143 mpfr_mul (t, xfrac, t, MPFR_RNDU); /* xfrac * ln(2) */
144 err = Nt - (MPFR_GET_EXP (t) + 2); /* Estimate of the error */
145 mpfr_exp (t, t, MPFR_RNDN); /* exp(xfrac * ln(2)) */
146
147 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
148 break;
149
150 /* Actualisation of the precision */
151 MPFR_ZIV_NEXT (loop, Nt);
152 mpfr_set_prec (t, Nt);
153 }
154 MPFR_ZIV_FREE (loop);
155
156 inexact = mpfr_set (y, t, rnd_mode);
157
158 mpfr_clear (t);
159 }
160
161 mpfr_clear (xfrac);
162
163 if (MPFR_UNLIKELY (rnd_mode == MPFR_RNDN && xint == __gmpfr_emin - 1 &&
164 MPFR_GET_EXP (y) == 0 && mpfr_powerof2_raw (y)))
165 {
166 /* y was rounded down to 1/2 and the rounded value with an unbounded
167 exponent range would be 2^(emin-2), i.e. the midpoint between 0
168 and the smallest positive FP number. This is a double rounding
169 problem: we should not round to 0, but to (1/2) * 2^emin. */
170 MPFR_SET_EXP (y, __gmpfr_emin);
171 inexact = 1;
172 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_UNDERFLOW);
173 }
174 else
175 {
176 /* The following is OK due to early overflow/underflow checking.
177 the exponent may be slightly out-of-range, but this will be
178 handled by mpfr_check_range. */
179 MPFR_EXP (y) += xint;
180 }
181
182 MPFR_SAVE_EXPO_FREE (expo);
183 return mpfr_check_range (y, inexact, rnd_mode);
184 }
185