summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/pow.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/pow.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/pow.c151
1 files changed, 74 insertions, 77 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/pow.c b/Build/source/libs/mpfr/mpfr-src/src/pow.c
index 1c4e82118fe..d6ebc910663 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/pow.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/pow.c
@@ -1,6 +1,6 @@
/* mpfr_pow -- power function x^y
-Copyright 2001-2019 Free Software Foundation, Inc.
+Copyright 2001-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
@@ -23,6 +23,10 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
#define MPFR_NEED_LONGLONG_H
#include "mpfr-impl.h"
+#ifndef MPFR_POW_EXP_THRESHOLD
+# define MPFR_POW_EXP_THRESHOLD (MAX (sizeof(mpfr_exp_t) * CHAR_BIT, 256))
+#endif
+
/* return non zero iff x^y is exact.
Assumes x and y are ordinary numbers,
y is not an integer, x is not a power of 2 and x is positive
@@ -34,8 +38,7 @@ mpfr_pow_is_exact (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
mpfr_rnd_t rnd_mode, int *inexact)
{
mpz_t a, c;
- mpfr_exp_t d, b;
- unsigned long i;
+ mpfr_exp_t d, b, i;
int res;
MPFR_ASSERTD (!MPFR_IS_SINGULAR (y));
@@ -48,7 +51,9 @@ mpfr_pow_is_exact (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
if (MPFR_IS_NEG (y))
return 0; /* x is not a power of two => x^-y is not exact */
- /* compute d such that y = c*2^d with c odd integer */
+ /* Compute d such that y = c*2^d with c odd integer.
+ Since c comes from a regular MPFR number, due to the constraints on the
+ exponent and the precision, there can be no integer overflow below. */
mpz_init (c);
d = mpfr_get_z_2exp (c, y);
i = mpz_scan1 (c, 0);
@@ -58,7 +63,9 @@ mpfr_pow_is_exact (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
/* Since y is not an integer, d is necessarily < 0 */
MPFR_ASSERTD (d < 0);
- /* Compute a,b such that x=a*2^b */
+ /* Compute a,b such that x=a*2^b.
+ Since a comes from a regular MPFR number, due to the constrainst on the
+ exponent and the precision, there can be no integer overflow below. */
mpz_init (a);
b = mpfr_get_z_2exp (a, x);
i = mpz_scan1 (a, 0);
@@ -109,7 +116,8 @@ mpfr_pow_is_exact (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
/* Assumes that the exponent range has already been extended and if y is
an integer, then the result is not exact in unbounded exponent range.
- If x < 0, assumes y is an integer.
+ If y_is_integer is non-zero, y is an integer (always when x < 0).
+ expo is the saved exponent range and flags (at the call to mpfr_pow).
*/
int
mpfr_pow_general (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
@@ -126,7 +134,6 @@ mpfr_pow_general (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
mpfr_exp_t err; /* error */
MPFR_ZIV_DECL (ziv_loop);
-
MPFR_LOG_FUNC
(("x[%Pu]=%.*Rg y[%Pu]=%.*Rg rnd=%d",
mpfr_get_prec (x), mpfr_log_prec, x,
@@ -321,30 +328,13 @@ mpfr_pow_general (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
*/
if (rnd_mode == MPFR_RNDN && inexact < 0 && lk < 0 &&
MPFR_GET_EXP (z) == __gmpfr_emin - 1 - lk && mpfr_powerof2_raw (z))
- {
- /* Rounding to nearest, real result > z * 2^k = 2^(emin - 2),
- * underflow case:
- * (a) if the precision of z is > 1, we will obtain the correct
- * result and exceptions by replacing z by nextabove(z).
- * (b) if the precision of z is 1, we first copy z to zcopy of
- * precision 2 bits and perform nextabove(zcopy).
- */
- if (MPFR_PREC(z) >= 2)
- mpfr_nextabove (z);
- else
- {
- mpfr_t zcopy;
- mpfr_init2 (zcopy, MPFR_PREC(z) + 1);
- mpfr_set (zcopy, z, MPFR_RNDZ);
- mpfr_nextabove (zcopy);
- inex2 = mpfr_mul_2si (z, zcopy, lk, rnd_mode);
- mpfr_clear (zcopy);
- goto under_over;
- }
- }
+ /* Rounding to nearest, real result > z * 2^k = 2^(emin - 2),
+ * underflow case: we will obtain the correct result and exceptions
+ * by replacing z by nextabove(z).
+ */
+ mpfr_nextabove (z);
MPFR_CLEAR_FLAGS ();
inex2 = mpfr_mul_2si (z, z, lk, rnd_mode);
- under_over:
if (inex2) /* underflow or overflow */
{
inexact = inex2;
@@ -537,25 +527,19 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
if (cmp_x_1 * MPFR_SIGN (y) > 0)
{
- mpfr_t t;
+ mpfr_t t, x_abs;
int negative, overflow;
+ /* FIXME: since we round y*log2|x| toward zero, we could also do early
+ underflow detection */
MPFR_SAVE_EXPO_MARK (expo);
- mpfr_init2 (t, 53);
- /* we want a lower bound on y*log2|x|:
- (i) if x > 0, it suffices to round log2(x) toward zero, and
- to round y*o(log2(x)) toward zero too;
- (ii) if x < 0, we first compute t = o(-x), with rounding toward 1,
- and then follow as in case (1). */
- if (MPFR_IS_POS (x))
- mpfr_log2 (t, x, MPFR_RNDZ);
- else
- {
- mpfr_neg (t, x, (cmp_x_1 > 0) ? MPFR_RNDZ : MPFR_RNDU);
- mpfr_log2 (t, t, MPFR_RNDZ);
- }
+ mpfr_init2 (t, sizeof (mpfr_exp_t) * CHAR_BIT);
+ /* we want a lower bound on y*log2|x| */
+ MPFR_TMP_INIT_ABS (x_abs, x);
+ mpfr_log2 (t, x_abs, MPFR_RNDZ);
mpfr_mul (t, t, y, MPFR_RNDZ);
- overflow = mpfr_cmp_si (t, __gmpfr_emax) > 0;
+ overflow = mpfr_cmp_si (t, __gmpfr_emax) >= 0;
+ /* if t >= emax, then |z| >= 2^t >= 2^emax and we have overflow */
mpfr_clear (t);
MPFR_SAVE_EXPO_FREE (expo);
if (overflow)
@@ -617,8 +601,9 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
MBs may be necessary). Note that in such a case, either x = +/-2^b
(this case is handled below) or x^y cannot be represented exactly in
any precision supported by MPFR (the general case uses this property).
- */
- if (y_is_integer && (MPFR_GET_EXP (y) <= 256))
+ Note: the threshold of 256 should not be decreased too much, see the
+ comments about (-2^b)^y just below. */
+ if (y_is_integer && MPFR_GET_EXP (y) <= MPFR_POW_EXP_THRESHOLD)
{
mpz_t zi;
@@ -633,37 +618,49 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
/* Special case (+/-2^b)^Y which could be exact. If x is negative, then
necessarily y is a large integer. */
if (mpfr_powerof2_raw (x))
- {
- mpfr_exp_t b = MPFR_GET_EXP (x) - 1;
- mpfr_t tmp;
- int sgnx = MPFR_SIGN (x);
-
- MPFR_ASSERTN (b >= LONG_MIN && b <= LONG_MAX); /* FIXME... */
-
- MPFR_LOG_MSG (("special case (+/-2^b)^Y\n", 0));
- /* now x = +/-2^b, so x^y = (+/-1)^y*2^(b*y) is exact whenever b*y is
- an integer */
- MPFR_SAVE_EXPO_MARK (expo);
- mpfr_init2 (tmp, MPFR_PREC (y) + sizeof (long) * CHAR_BIT);
- inexact = mpfr_mul_si (tmp, y, b, MPFR_RNDN); /* exact */
- MPFR_ASSERTN (inexact == 0);
- /* Note: as the exponent range has been extended, an overflow is not
- possible (due to basic overflow and underflow checking above, as
- the result is ~ 2^tmp), and an underflow is not possible either
- because b is an integer (thus either 0 or >= 1). */
- MPFR_CLEAR_FLAGS ();
- inexact = mpfr_exp2 (z, tmp, rnd_mode);
- mpfr_clear (tmp);
- if (sgnx < 0 && mpfr_odd_p (y))
- {
- mpfr_neg (z, z, rnd_mode);
- inexact = -inexact;
- }
- /* Without the following, the overflows3 test in tpow.c fails. */
- MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
- MPFR_SAVE_EXPO_FREE (expo);
- return mpfr_check_range (z, inexact, rnd_mode);
- }
+ {
+ mpfr_exp_t b = MPFR_GET_EXP (x) - 1;
+ mpfr_t tmp;
+
+ MPFR_STAT_STATIC_ASSERT (MPFR_POW_EXP_THRESHOLD >=
+ sizeof(mpfr_exp_t) * CHAR_BIT);
+
+ /* For x < 0, we have EXP(y) > MPFR_POW_EXP_THRESHOLD, thus
+ EXP(y) > bitsize of mpfr_exp_t (1). Therefore, since |x| <> 1:
+ (a) either |x| >= 2, and we have an overflow due to (1), but
+ this was detected by the early overflow detection above,
+ i.e. this case is not possible;
+ (b) either |x| <= 1/2, and we have underflow. */
+ if (MPFR_SIGN (x) < 0)
+ {
+ MPFR_ASSERTD (MPFR_EXP (x) <= 0);
+ MPFR_ASSERTD (MPFR_EXP (y) > sizeof(mpfr_exp_t) * CHAR_BIT);
+ return mpfr_underflow (z,
+ rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode,
+ MPFR_IS_NEG (x) && mpfr_odd_p (y) ? -1 : 1);
+ }
+
+ MPFR_ASSERTN (b >= LONG_MIN && b <= LONG_MAX); /* FIXME... */
+
+ MPFR_LOG_MSG (("special case (+/-2^b)^Y\n", 0));
+ /* now x = +/-2^b, so x^y = (+/-1)^y*2^(b*y) is exact whenever b*y is
+ an integer */
+ MPFR_SAVE_EXPO_MARK (expo);
+ mpfr_init2 (tmp, MPFR_PREC (y) + sizeof (long) * CHAR_BIT);
+ inexact = mpfr_mul_si (tmp, y, b, MPFR_RNDN); /* exact */
+ MPFR_ASSERTN (inexact == 0);
+ /* Note: as the exponent range has been extended, an overflow is not
+ possible (due to basic overflow and underflow checking above, as
+ the result is ~ 2^tmp), and an underflow is not possible either
+ because b is an integer (thus either 0 or >= 1). */
+ MPFR_CLEAR_FLAGS ();
+ inexact = mpfr_exp2 (z, tmp, rnd_mode);
+ mpfr_clear (tmp);
+ /* Without the following, the overflows3 test in tpow.c fails. */
+ MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
+ MPFR_SAVE_EXPO_FREE (expo);
+ return mpfr_check_range (z, inexact, rnd_mode);
+ }
MPFR_SAVE_EXPO_MARK (expo);