1 /* mpfr_pow_si -- power function x^y with y a signed int
3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao projects, INRIA.
6 This file is part of the GNU MPFR Library.
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 2.1 of the License, or (at your
11 option) any later version.
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.
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.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
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
32 mpfr_pow_si (mpfr_ptr y
, mpfr_srcptr x
, long int n
, mp_rnd_t rnd
)
34 MPFR_LOG_FUNC (("x[%#R]=%R n=%ld rnd=%d", x
, x
, n
, rnd
),
38 return mpfr_pow_ui (y
, x
, n
, rnd
);
41 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
48 else if (MPFR_IS_INF (x
))
51 if (MPFR_IS_POS (x
) || ((unsigned) n
& 1) == 0)
59 MPFR_ASSERTD (MPFR_IS_ZERO (x
));
61 if (MPFR_IS_POS (x
) || ((unsigned) n
& 1) == 0)
70 /* detect exact powers: x^(-n) is exact iff x is a power of 2 */
71 if (mpfr_cmp_si_2exp (x
, MPFR_SIGN(x
), MPFR_EXP(x
) - 1) == 0)
73 mp_exp_t expx
= MPFR_EXP (x
) - 1, expy
;
75 /* Warning: n * expx may overflow!
76 * Some systems (apparently alpha-freebsd) abort with
77 * LONG_MIN / 1, and LONG_MIN / -1 is undefined.
78 * Proof of the overflow checking. The expressions below are
79 * assumed to be on the rational numbers, but the word "overflow"
80 * still has its own meaning in the C context. / still denotes
81 * the integer (truncated) division, and // denotes the exact
83 * - First, (__gmpfr_emin - 1) / n and (__gmpfr_emax - 1) / n
84 * cannot overflow due to the constraints on the exponents of
86 * - If n = -1, then n * expx = - expx, which is representable
87 * because of the constraints on the exponents of MPFR numbers.
88 * - If expx = 0, then n * expx = 0, which is representable.
89 * - If n < -1 and expx > 0:
90 * + If expx > (__gmpfr_emin - 1) / n, then
91 * expx >= (__gmpfr_emin - 1) / n + 1
92 * > (__gmpfr_emin - 1) // n,
94 * n * expx < __gmpfr_emin - 1,
96 * n * expx <= __gmpfr_emin - 2.
97 * This corresponds to an underflow, with a null result in
98 * the rounding-to-nearest mode.
99 * + If expx <= (__gmpfr_emin - 1) / n, then n * expx cannot
100 * overflow since 0 < expx <= (__gmpfr_emin - 1) / n and
101 * 0 > n * expx >= n * ((__gmpfr_emin - 1) / n)
102 * >= __gmpfr_emin - 1.
103 * - If n < -1 and expx < 0:
104 * + If expx < (__gmpfr_emax - 1) / n, then
105 * expx <= (__gmpfr_emax - 1) / n - 1
106 * < (__gmpfr_emax - 1) // n,
108 * n * expx > __gmpfr_emax - 1,
110 * n * expx >= __gmpfr_emax.
111 * This corresponds to an overflow (2^(n * expx) has an
112 * exponent > __gmpfr_emax).
113 * + If expx >= (__gmpfr_emax - 1) / n, then n * expx cannot
114 * overflow since 0 > expx >= (__gmpfr_emax - 1) / n and
115 * 0 < n * expx <= n * ((__gmpfr_emax - 1) / n)
116 * <= __gmpfr_emax - 1.
117 * Note: one could use expx bounds based on MPFR_EXP_MIN and
118 * MPFR_EXP_MAX instead of __gmpfr_emin and __gmpfr_emax. The
119 * current bounds do not lead to noticeably slower code and
120 * allow us to avoid a bug in Sun's compiler for Solaris/x86
121 * (when optimizations are enabled).
124 n
!= -1 && expx
> 0 && expx
> (__gmpfr_emin
- 1) / n
?
125 MPFR_EMIN_MIN
- 2 /* Underflow */ :
126 n
!= -1 && expx
< 0 && expx
< (__gmpfr_emax
- 1) / n
?
127 MPFR_EMAX_MAX
/* Overflow */ : n
* expx
;
128 return mpfr_set_si_2exp (y
, n
% 2 ? MPFR_INT_SIGN (x
) : 1,
134 /* Declaration of the intermediary variable */
136 /* Declaration of the size variable */
137 mp_prec_t Ny
; /* target precision */
138 mp_prec_t Nt
; /* working precision */
143 MPFR_SAVE_EXPO_DECL (expo
);
144 MPFR_ZIV_DECL (loop
);
146 abs_n
= - (unsigned long) n
;
147 count_leading_zeros (size_n
, (mp_limb_t
) abs_n
);
148 size_n
= BITS_PER_MP_LIMB
- size_n
;
150 /* initial working precision */
152 Nt
= Ny
+ size_n
+ 3 + MPFR_INT_CEIL_LOG2 (Ny
);
154 MPFR_SAVE_EXPO_MARK (expo
);
156 /* initialise of intermediary variable */
159 /* We will compute rnd(rnd1(1/x) ^ |n|), where rnd1 is the rounding
160 toward sign(x), to avoid spurious overflow or underflow, as in
162 rnd1
= MPFR_EXP (x
) < 1 ? GMP_RNDZ
:
163 (MPFR_SIGN (x
) > 0 ? GMP_RNDU
: GMP_RNDD
);
165 MPFR_ZIV_INIT (loop
, Nt
);
168 MPFR_BLOCK_DECL (flags
);
170 /* compute (1/x)^|n| */
171 MPFR_BLOCK (flags
, mpfr_ui_div (t
, 1, x
, rnd1
));
172 MPFR_ASSERTD (! MPFR_UNDERFLOW (flags
));
173 /* t = (1/x)*(1+theta) where |theta| <= 2^(-Nt) */
174 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags
)))
176 MPFR_BLOCK (flags
, mpfr_pow_ui (t
, t
, abs_n
, rnd
));
177 /* t = (1/x)^|n|*(1+theta')^(|n|+1) where |theta'| <= 2^(-Nt).
178 If (|n|+1)*2^(-Nt) <= 1/2, which is satisfied as soon as
179 Nt >= bits(n)+2, then we can use Lemma \ref{lemma_graillat}
180 from algorithms.tex, which yields x^n*(1+theta) with
181 |theta| <= 2(|n|+1)*2^(-Nt), thus the error is bounded by
182 2(|n|+1) ulps <= 2^(bits(n)+2) ulps. */
183 if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags
)))
186 MPFR_ZIV_FREE (loop
);
188 MPFR_SAVE_EXPO_FREE (expo
);
189 MPFR_LOG_MSG (("overflow\n", 0));
190 return mpfr_overflow (y
, rnd
, abs_n
& 1 ?
191 MPFR_SIGN (x
) : MPFR_SIGN_POS
);
193 if (MPFR_UNLIKELY (MPFR_UNDERFLOW (flags
)))
195 MPFR_ZIV_FREE (loop
);
197 MPFR_LOG_MSG (("underflow\n", 0));
202 /* We cannot decide now whether the result should be
203 rounded toward zero or away from zero. So, like
204 in mpfr_pow_pos_z, let's use the general case of
205 mpfr_pow in precision 2. */
206 MPFR_ASSERTD (mpfr_cmp_si_2exp (x
, MPFR_SIGN (x
),
207 MPFR_EXP (x
) - 1) != 0);
209 mpfr_init2 (nn
, sizeof (long) * CHAR_BIT
);
210 inexact
= mpfr_set_si (nn
, n
, GMP_RNDN
);
211 MPFR_ASSERTN (inexact
== 0);
212 inexact
= mpfr_pow_general (y2
, x
, nn
, rnd
, 1,
213 (mpfr_save_expo_t
*) NULL
);
215 mpfr_set (y
, y2
, GMP_RNDN
);
217 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, MPFR_FLAGS_UNDERFLOW
);
222 MPFR_SAVE_EXPO_FREE (expo
);
223 return mpfr_underflow (y
, rnd
, abs_n
& 1 ?
224 MPFR_SIGN (x
) : MPFR_SIGN_POS
);
227 /* error estimate -- see pow function in algorithms.ps */
228 if (MPFR_LIKELY (MPFR_CAN_ROUND (t
, Nt
- size_n
- 2, Ny
, rnd
)))
231 /* actualisation of the precision */
232 MPFR_ZIV_NEXT (loop
, Nt
);
233 mpfr_set_prec (t
, Nt
);
235 MPFR_ZIV_FREE (loop
);
237 inexact
= mpfr_set (y
, t
, rnd
);
241 MPFR_SAVE_EXPO_FREE (expo
);
242 return mpfr_check_range (y
, inexact
, rnd
);