1 /* mpfr_pow_si -- power function x^y with y a signed int
2
3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel 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 http://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 y = pow_si(x,n) is done by
27 * y = pow_ui(x,n) if n >= 0
28 * y = 1 / pow_ui(x,-n) if n < 0
29 */
30
31 int
mpfr_pow_si(mpfr_ptr y,mpfr_srcptr x,long int n,mpfr_rnd_t rnd)32 mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd)
33 {
34 MPFR_LOG_FUNC
35 (("x[%Pu]=%.*Rg n=%ld rnd=%d",
36 mpfr_get_prec (x), mpfr_log_prec, x, n, rnd),
37 ("y[%Pu]=%.*Rg", mpfr_get_prec (y), mpfr_log_prec, y));
38
39 if (n >= 0)
40 return mpfr_pow_ui (y, x, n, rnd);
41 else
42 {
43 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
44 {
45 if (MPFR_IS_NAN (x))
46 {
47 MPFR_SET_NAN (y);
48 MPFR_RET_NAN;
49 }
50 else
51 {
52 int positive = MPFR_IS_POS (x) || ((unsigned long) n & 1) == 0;
53 if (MPFR_IS_INF (x))
54 MPFR_SET_ZERO (y);
55 else /* x is zero */
56 {
57 MPFR_ASSERTD (MPFR_IS_ZERO (x));
58 MPFR_SET_INF (y);
59 mpfr_set_divby0 ();
60 }
61 if (positive)
62 MPFR_SET_POS (y);
63 else
64 MPFR_SET_NEG (y);
65 MPFR_RET (0);
66 }
67 }
68
69 /* detect exact powers: x^(-n) is exact iff x is a power of 2 */
70 if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0)
71 {
72 mpfr_exp_t expx = MPFR_EXP (x) - 1, expy;
73 MPFR_ASSERTD (n < 0);
74 /* Warning: n * expx may overflow!
75 *
76 * Some systems (apparently alpha-freebsd) abort with
77 * LONG_MIN / 1, and LONG_MIN / -1 is undefined.
78 * http://www.freebsd.org/cgi/query-pr.cgi?pr=72024
79 *
80 * Proof of the overflow checking. The expressions below are
81 * assumed to be on the rational numbers, but the word "overflow"
82 * still has its own meaning in the C context. / still denotes
83 * the integer (truncated) division, and // denotes the exact
84 * division.
85 * - First, (__gmpfr_emin - 1) / n and (__gmpfr_emax - 1) / n
86 * cannot overflow due to the constraints on the exponents of
87 * MPFR numbers.
88 * - If n = -1, then n * expx = - expx, which is representable
89 * because of the constraints on the exponents of MPFR numbers.
90 * - If expx = 0, then n * expx = 0, which is representable.
91 * - If n < -1 and expx > 0:
92 * + If expx > (__gmpfr_emin - 1) / n, then
93 * expx >= (__gmpfr_emin - 1) / n + 1
94 * > (__gmpfr_emin - 1) // n,
95 * and
96 * n * expx < __gmpfr_emin - 1,
97 * i.e.
98 * n * expx <= __gmpfr_emin - 2.
99 * This corresponds to an underflow, with a null result in
100 * the rounding-to-nearest mode.
101 * + If expx <= (__gmpfr_emin - 1) / n, then n * expx cannot
102 * overflow since 0 < expx <= (__gmpfr_emin - 1) / n and
103 * 0 > n * expx >= n * ((__gmpfr_emin - 1) / n)
104 * >= __gmpfr_emin - 1.
105 * - If n < -1 and expx < 0:
106 * + If expx < (__gmpfr_emax - 1) / n, then
107 * expx <= (__gmpfr_emax - 1) / n - 1
108 * < (__gmpfr_emax - 1) // n,
109 * and
110 * n * expx > __gmpfr_emax - 1,
111 * i.e.
112 * n * expx >= __gmpfr_emax.
113 * This corresponds to an overflow (2^(n * expx) has an
114 * exponent > __gmpfr_emax).
115 * + If expx >= (__gmpfr_emax - 1) / n, then n * expx cannot
116 * overflow since 0 > expx >= (__gmpfr_emax - 1) / n and
117 * 0 < n * expx <= n * ((__gmpfr_emax - 1) / n)
118 * <= __gmpfr_emax - 1.
119 * Note: one could use expx bounds based on MPFR_EXP_MIN and
120 * MPFR_EXP_MAX instead of __gmpfr_emin and __gmpfr_emax. The
121 * current bounds do not lead to noticeably slower code and
122 * allow us to avoid a bug in Sun's compiler for Solaris/x86
123 * (when optimizations are enabled); known affected versions:
124 * cc: Sun C 5.8 2005/10/13
125 * cc: Sun C 5.8 Patch 121016-02 2006/03/31
126 * cc: Sun C 5.8 Patch 121016-04 2006/10/18
127 */
128 expy =
129 n != -1 && expx > 0 && expx > (__gmpfr_emin - 1) / n ?
130 MPFR_EMIN_MIN - 2 /* Underflow */ :
131 n != -1 && expx < 0 && expx < (__gmpfr_emax - 1) / n ?
132 MPFR_EMAX_MAX /* Overflow */ : n * expx;
133 return mpfr_set_si_2exp (y, n % 2 ? MPFR_INT_SIGN (x) : 1,
134 expy, rnd);
135 }
136
137 /* General case */
138 {
139 /* Declaration of the intermediary variable */
140 mpfr_t t;
141 /* Declaration of the size variable */
142 mpfr_prec_t Ny; /* target precision */
143 mpfr_prec_t Nt; /* working precision */
144 mpfr_rnd_t rnd1;
145 int size_n;
146 int inexact;
147 unsigned long abs_n;
148 MPFR_SAVE_EXPO_DECL (expo);
149 MPFR_ZIV_DECL (loop);
150
151 abs_n = - (unsigned long) n;
152 count_leading_zeros (size_n, (mp_limb_t) abs_n);
153 size_n = GMP_NUMB_BITS - size_n;
154
155 /* initial working precision */
156 Ny = MPFR_PREC (y);
157 Nt = Ny + size_n + 3 + MPFR_INT_CEIL_LOG2 (Ny);
158
159 MPFR_SAVE_EXPO_MARK (expo);
160
161 /* initialise of intermediary variable */
162 mpfr_init2 (t, Nt);
163
164 /* We will compute rnd(rnd1(1/x) ^ |n|), where rnd1 is the rounding
165 toward sign(x), to avoid spurious overflow or underflow, as in
166 mpfr_pow_z. */
167 rnd1 = MPFR_EXP (x) < 1 ? MPFR_RNDZ :
168 (MPFR_SIGN (x) > 0 ? MPFR_RNDU : MPFR_RNDD);
169
170 MPFR_ZIV_INIT (loop, Nt);
171 for (;;)
172 {
173 MPFR_BLOCK_DECL (flags);
174
175 /* compute (1/x)^|n| */
176 MPFR_BLOCK (flags, mpfr_ui_div (t, 1, x, rnd1));
177 MPFR_ASSERTD (! MPFR_UNDERFLOW (flags));
178 /* t = (1/x)*(1+theta) where |theta| <= 2^(-Nt) */
179 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
180 goto overflow;
181 MPFR_BLOCK (flags, mpfr_pow_ui (t, t, abs_n, rnd));
182 /* t = (1/x)^|n|*(1+theta')^(|n|+1) where |theta'| <= 2^(-Nt).
183 If (|n|+1)*2^(-Nt) <= 1/2, which is satisfied as soon as
184 Nt >= bits(n)+2, then we can use Lemma \ref{lemma_graillat}
185 from algorithms.tex, which yields x^n*(1+theta) with
186 |theta| <= 2(|n|+1)*2^(-Nt), thus the error is bounded by
187 2(|n|+1) ulps <= 2^(bits(n)+2) ulps. */
188 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
189 {
190 overflow:
191 MPFR_ZIV_FREE (loop);
192 mpfr_clear (t);
193 MPFR_SAVE_EXPO_FREE (expo);
194 MPFR_LOG_MSG (("overflow\n", 0));
195 return mpfr_overflow (y, rnd, abs_n & 1 ?
196 MPFR_SIGN (x) : MPFR_SIGN_POS);
197 }
198 if (MPFR_UNLIKELY (MPFR_UNDERFLOW (flags)))
199 {
200 MPFR_ZIV_FREE (loop);
201 mpfr_clear (t);
202 MPFR_LOG_MSG (("underflow\n", 0));
203 if (rnd == MPFR_RNDN)
204 {
205 mpfr_t y2, nn;
206
207 /* We cannot decide now whether the result should be
208 rounded toward zero or away from zero. So, like
209 in mpfr_pow_pos_z, let's use the general case of
210 mpfr_pow in precision 2. */
211 MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_SIGN (x),
212 MPFR_EXP (x) - 1) != 0);
213 mpfr_init2 (y2, 2);
214 mpfr_init2 (nn, sizeof (long) * CHAR_BIT);
215 inexact = mpfr_set_si (nn, n, MPFR_RNDN);
216 MPFR_ASSERTN (inexact == 0);
217 inexact = mpfr_pow_general (y2, x, nn, rnd, 1,
218 (mpfr_save_expo_t *) NULL);
219 mpfr_clear (nn);
220 mpfr_set (y, y2, MPFR_RNDN);
221 mpfr_clear (y2);
222 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_UNDERFLOW);
223 goto end;
224 }
225 else
226 {
227 MPFR_SAVE_EXPO_FREE (expo);
228 return mpfr_underflow (y, rnd, abs_n & 1 ?
229 MPFR_SIGN (x) : MPFR_SIGN_POS);
230 }
231 }
232 /* error estimate -- see pow function in algorithms.ps */
233 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - size_n - 2, Ny, rnd)))
234 break;
235
236 /* actualisation of the precision */
237 MPFR_ZIV_NEXT (loop, Nt);
238 mpfr_set_prec (t, Nt);
239 }
240 MPFR_ZIV_FREE (loop);
241
242 inexact = mpfr_set (y, t, rnd);
243 mpfr_clear (t);
244
245 end:
246 MPFR_SAVE_EXPO_FREE (expo);
247 return mpfr_check_range (y, inexact, rnd);
248 }
249 }
250 }
251