1 /* mpfr_pow_si -- power function x^y with y a signed int
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 /* for MPFR_INT_CEIL_LOG2 */
24 #include "mpfr-impl.h"
25
26 /* The computation of y = pow_si/sj(x,n) is done by
27 * y = pow_ui/uj(x,n) if n >= 0
28 * y = 1 / pow_ui/uj(x,-n) if n < 0
29 */
30
31 #ifndef POW_S
32 #define POW_S mpfr_pow_si
33 #define POW_U mpfr_pow_ui
34 #define SET_S mpfr_set_si
35 #define SET_S_2EXP mpfr_set_si_2exp
36 #define NBITS_UTYPE mpfr_nbits_ulong
37 #define TYPE long int
38 #define UTYPE unsigned long
39 #define FSPEC "l"
40 #endif
41
42 int
POW_S(mpfr_ptr y,mpfr_srcptr x,TYPE n,mpfr_rnd_t rnd)43 POW_S (mpfr_ptr y, mpfr_srcptr x, TYPE n, mpfr_rnd_t rnd)
44 {
45 MPFR_LOG_FUNC
46 (("x[%Pd]=%.*Rg n=%" FSPEC "d rnd=%d",
47 mpfr_get_prec (x), mpfr_log_prec, x, n, rnd),
48 ("y[%Pd]=%.*Rg", mpfr_get_prec (y), mpfr_log_prec, y));
49
50 if (n >= 0)
51 return POW_U (y, x, n, rnd);
52 else
53 {
54 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
55 {
56 if (MPFR_IS_NAN (x))
57 {
58 MPFR_SET_NAN (y);
59 MPFR_RET_NAN;
60 }
61 else
62 {
63 int positive = MPFR_IS_POS (x) || ((UTYPE) n & 1) == 0;
64 if (MPFR_IS_INF (x))
65 MPFR_SET_ZERO (y);
66 else /* x is zero */
67 {
68 MPFR_ASSERTD (MPFR_IS_ZERO (x));
69 MPFR_SET_INF (y);
70 MPFR_SET_DIVBY0 ();
71 }
72 if (positive)
73 MPFR_SET_POS (y);
74 else
75 MPFR_SET_NEG (y);
76 MPFR_RET (0);
77 }
78 }
79
80 /* detect exact powers: x^(-n) is exact iff x is a power of 2 */
81 if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0)
82 {
83 mpfr_exp_t expx = MPFR_EXP (x) - 1, expy;
84 MPFR_ASSERTD (n < 0);
85 /* Warning: n * expx may overflow!
86 *
87 * Some systems (apparently alpha-freebsd) abort with
88 * LONG_MIN / 1, and LONG_MIN / -1 is undefined.
89 * https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=72024
90 *
91 * Proof of the overflow checking. The expressions below are
92 * assumed to be on the rational numbers, but the word "overflow"
93 * still has its own meaning in the C context. / still denotes
94 * the integer (truncated) division, and // denotes the exact
95 * division.
96 * - First, (__gmpfr_emin - 1) / n and (__gmpfr_emax - 1) / n
97 * cannot overflow due to the constraints on the exponents of
98 * MPFR numbers.
99 * - If n = -1, then n * expx = - expx, which is representable
100 * because of the constraints on the exponents of MPFR numbers.
101 * - If expx = 0, then n * expx = 0, which is representable.
102 * - If n < -1 and expx > 0:
103 * + If expx > (__gmpfr_emin - 1) / n, then
104 * expx >= (__gmpfr_emin - 1) / n + 1
105 * > (__gmpfr_emin - 1) // n,
106 * and
107 * n * expx < __gmpfr_emin - 1,
108 * i.e.
109 * n * expx <= __gmpfr_emin - 2.
110 * This corresponds to an underflow, with a null result in
111 * the rounding-to-nearest mode.
112 * + If expx <= (__gmpfr_emin - 1) / n, then n * expx cannot
113 * overflow since 0 < expx <= (__gmpfr_emin - 1) / n and
114 * 0 > n * expx >= n * ((__gmpfr_emin - 1) / n)
115 * >= __gmpfr_emin - 1.
116 * - If n < -1 and expx < 0:
117 * + If expx < (__gmpfr_emax - 1) / n, then
118 * expx <= (__gmpfr_emax - 1) / n - 1
119 * < (__gmpfr_emax - 1) // n,
120 * and
121 * n * expx > __gmpfr_emax - 1,
122 * i.e.
123 * n * expx >= __gmpfr_emax.
124 * This corresponds to an overflow (2^(n * expx) has an
125 * exponent > __gmpfr_emax).
126 * + If expx >= (__gmpfr_emax - 1) / n, then n * expx cannot
127 * overflow since 0 > expx >= (__gmpfr_emax - 1) / n and
128 * 0 < n * expx <= n * ((__gmpfr_emax - 1) / n)
129 * <= __gmpfr_emax - 1.
130 * Note: one could use expx bounds based on MPFR_EXP_MIN and
131 * MPFR_EXP_MAX instead of __gmpfr_emin and __gmpfr_emax. The
132 * current bounds do not lead to noticeably slower code and
133 * allow us to avoid a bug in Sun's compiler for Solaris/x86
134 * (when optimizations are enabled); known affected versions:
135 * cc: Sun C 5.8 2005/10/13
136 * cc: Sun C 5.8 Patch 121016-02 2006/03/31
137 * cc: Sun C 5.8 Patch 121016-04 2006/10/18
138 */
139 expy =
140 n != -1 && expx > 0 && expx > (__gmpfr_emin - 1) / n ?
141 MPFR_EMIN_MIN - 2 /* Underflow */ :
142 n != -1 && expx < 0 && expx < (__gmpfr_emax - 1) / n ?
143 MPFR_EMAX_MAX /* Overflow */ : n * expx;
144 return SET_S_2EXP (y, n % 2 ? MPFR_INT_SIGN (x) : 1, expy, rnd);
145 }
146
147 /* General case */
148 {
149 /* Declaration of the intermediary variable */
150 mpfr_t t;
151 /* Declaration of the size variable */
152 mpfr_prec_t Ny; /* target precision */
153 mpfr_prec_t Nt; /* working precision */
154 mpfr_rnd_t rnd1;
155 int size_n;
156 int inexact;
157 UTYPE abs_n;
158 MPFR_SAVE_EXPO_DECL (expo);
159 MPFR_ZIV_DECL (loop);
160
161 abs_n = - (UTYPE) n;
162 size_n = NBITS_UTYPE (abs_n);
163
164 /* initial working precision */
165 Ny = MPFR_PREC (y);
166 Nt = Ny + size_n + 3 + MPFR_INT_CEIL_LOG2 (Ny);
167
168 MPFR_SAVE_EXPO_MARK (expo);
169
170 /* initialize of intermediary variable */
171 mpfr_init2 (t, Nt);
172
173 /* We will compute rnd(rnd1(1/x) ^ |n|), where rnd1 is the rounding
174 toward sign(x), to avoid spurious overflow or underflow, as in
175 mpfr_pow_z. */
176 rnd1 = MPFR_EXP (x) < 1 ? MPFR_RNDZ :
177 (MPFR_IS_POS (x) ? MPFR_RNDU : MPFR_RNDD);
178
179 /* The following ensures that 1/x cannot underflow.
180 Since |x| < 2^emax, |1/x| > 2^(-emax) >= 2^emin. */
181 MPFR_STAT_STATIC_ASSERT (MPFR_EMIN_MIN + MPFR_EMAX_MAX <= 0);
182
183 MPFR_ZIV_INIT (loop, Nt);
184 for (;;)
185 {
186 MPFR_BLOCK_DECL (flags);
187
188 /* TODO: Compute POW_U before the division (instead of after)
189 in order to reduce the error in the intermediate result?
190 POW_U, whose condition number is |n|, which may be large,
191 would be called on an exact value. This may be important
192 in very small precisions.
193 In this case, if x^|n| underflows, then |x^n| > 2^emax
194 (real overflow, and we can return the result); and if
195 x^|n| overflows, then the result underflows or is very
196 close to the underflow threshold, so that we should use
197 mpfr_pow_general (as already done for MPFR_RNDN), which
198 can handle such a case.
199 So the advantage of computing POW_U before the division
200 is that the code would be slightly faster is the general
201 case, but it could be noticeably slower in very uncommon
202 cases (and only with the extended exponent range). */
203
204 /* compute (1/x)^|n| */
205 MPFR_BLOCK (flags, mpfr_ui_div (t, 1, x, rnd1));
206 MPFR_ASSERTD (! MPFR_UNDERFLOW (flags));
207 /* t = (1/x)*(1+theta) where |theta| <= 2^(-Nt) */
208 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
209 goto overflow;
210 MPFR_BLOCK (flags, POW_U (t, t, abs_n, rnd));
211 /* t = (1/x)^|n|*(1+theta')^(|n|+1) where |theta'| <= 2^(-Nt).
212 If (|n|+1)*2^(-Nt) <= 1/2, which is satisfied as soon as
213 Nt >= bits(n)+2, then we can use Lemma \ref{lemma_graillat}
214 from algorithms.tex, which yields x^n*(1+theta) with
215 |theta| <= 2(|n|+1)*2^(-Nt), thus the error is bounded by
216 2(|n|+1) ulps <= 2^(bits(n)+2) ulps. */
217 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
218 {
219 overflow:
220 MPFR_ZIV_FREE (loop);
221 mpfr_clear (t);
222 MPFR_SAVE_EXPO_FREE (expo);
223 MPFR_LOG_MSG (("overflow\n", 0));
224 return mpfr_overflow (y, rnd, abs_n & 1 ?
225 MPFR_SIGN (x) : MPFR_SIGN_POS);
226 }
227 if (MPFR_UNLIKELY (MPFR_UNDERFLOW (flags)))
228 {
229 MPFR_ZIV_FREE (loop);
230 mpfr_clear (t);
231 MPFR_LOG_MSG (("underflow\n", 0));
232 if (rnd == MPFR_RNDN)
233 {
234 mpfr_t y2, nn;
235
236 /* We cannot decide now whether the result should be
237 rounded toward zero or away from zero. So, like
238 in mpfr_pow_pos_z, let's use the general case of
239 mpfr_pow in precision 2. */
240 MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_SIGN (x),
241 MPFR_EXP (x) - 1) != 0);
242 mpfr_init2 (y2, 2);
243 mpfr_init2 (nn, sizeof (TYPE) * CHAR_BIT);
244 inexact = SET_S (nn, n, MPFR_RNDN);
245 MPFR_ASSERTN (inexact == 0);
246 inexact = mpfr_pow_general (y2, x, nn, rnd, 1,
247 (mpfr_save_expo_t *) NULL);
248 mpfr_clear (nn);
249 mpfr_set (y, y2, MPFR_RNDN);
250 mpfr_clear (y2);
251 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_UNDERFLOW);
252 goto end;
253 }
254 else
255 {
256 MPFR_SAVE_EXPO_FREE (expo);
257 return mpfr_underflow (y, rnd, abs_n & 1 ?
258 MPFR_SIGN (x) : MPFR_SIGN_POS);
259 }
260 }
261 /* error estimate -- see pow function in algorithms.ps */
262 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - size_n - 2, Ny, rnd)))
263 break;
264
265 /* actualization of the precision */
266 MPFR_ZIV_NEXT (loop, Nt);
267 mpfr_set_prec (t, Nt);
268 }
269 MPFR_ZIV_FREE (loop);
270
271 inexact = mpfr_set (y, t, rnd);
272 mpfr_clear (t);
273
274 end:
275 MPFR_SAVE_EXPO_FREE (expo);
276 return mpfr_check_range (y, inexact, rnd);
277 }
278 }
279 }
280