summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/sub1.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/sub1.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/sub1.c126
1 files changed, 87 insertions, 39 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/sub1.c b/Build/source/libs/mpfr/mpfr-src/src/sub1.c
index 4584e67d682..6ffdcd358a3 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/sub1.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/sub1.c
@@ -1,6 +1,6 @@
/* mpfr_sub1 -- internal function to perform a "real" subtraction
-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.
@@ -28,8 +28,6 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
a positive value otherwise.
*/
-/* TODO: check the code in case exp_b == MPFR_EXP_MAX. */
-
int
mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
@@ -55,8 +53,10 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
(void) MPFR_GET_PREC (c);
sign = mpfr_cmp2 (b, c, &cancel);
+
if (MPFR_UNLIKELY(sign == 0))
{
+ MPFR_LOG_MSG (("sign=0\n", 0));
if (rnd_mode == MPFR_RNDD)
MPFR_SET_NEG (a);
else
@@ -65,6 +65,10 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
MPFR_RET (0);
}
+ /* sign != 0, so that cancel has a valid value. */
+ MPFR_LOG_MSG (("sign=%d cancel=%Pd\n", sign, cancel));
+ MPFR_ASSERTD (cancel >= 0 && cancel <= MPFR_PREC_MAX);
+
/*
* If subtraction: sign(a) = sign * sign(b)
* If addition: sign(a) = sign of the larger argument in absolute value.
@@ -91,9 +95,29 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
if (MPFR_UNLIKELY (MPFR_IS_UBF (b) || MPFR_IS_UBF (c)))
{
- exp_b = MPFR_IS_UBF (b) ?
- mpfr_ubf_zexp2exp (MPFR_ZEXP (b)) : MPFR_GET_EXP (b);
+ exp_b = MPFR_UBF_GET_EXP (b);
+ /* Early underflow detection. Rare, but a test is needed anyway
+ since in the "MAX (aq, bq) + 2 <= diff_exp" branch, the exponent
+ may decrease and MPFR_EXP_MIN would yield an integer overflow. */
+ if (MPFR_UNLIKELY (exp_b < __gmpfr_emin - 1))
+ {
+ if (rnd_mode == MPFR_RNDN)
+ rnd_mode = MPFR_RNDZ;
+ return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
+ }
diff_exp = mpfr_ubf_diff_exp (b, c);
+ MPFR_LOG_MSG (("UBF: exp_b=%" MPFR_EXP_FSPEC "d%s "
+ "diff_exp=%" MPFR_EXP_FSPEC "d%s\n",
+ (mpfr_eexp_t) exp_b,
+ exp_b == MPFR_EXP_MAX ? "=MPFR_EXP_MAX" : "",
+ (mpfr_eexp_t) diff_exp,
+ diff_exp == MPFR_EXP_MAX ? "=MPFR_EXP_MAX" : ""));
+ /* If diff_exp == MPFR_EXP_MAX, the actual value can be larger,
+ but anyway, since mpfr_exp_t >= mp_size_t, this will be the
+ case c small below, and the exact value does not matter. */
+ /* mpfr_set4 below used with MPFR_RNDF does not support UBF. */
+ if (rnd_mode == MPFR_RNDF)
+ rnd_mode = MPFR_RNDN;
}
else
{
@@ -123,27 +147,24 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
if (rnd_mode == MPFR_RNDF)
return mpfr_set4 (a, b, MPFR_RNDZ, MPFR_SIGN (a));
- MPFR_EXP (a) = exp_b; /* may be up to MPFR_EXP_MAX */
+ exp_a = exp_b; /* may be any out-of-range value due to UBF */
MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), bq,
rnd_mode, MPFR_SIGN (a),
- if (MPFR_EXP (a) != MPFR_EXP_MAX)
- ++ MPFR_EXP (a));
+ if (exp_a != MPFR_EXP_MAX)
+ exp_a ++);
MPFR_LOG_MSG (("inexact=%d\n", inexact));
- if (inexact == 0)
- {
+ if (inexact == 0 &&
/* a = b, but the exact value of b - c is a bit below. Then,
except for directed rounding similar to toward zero and
before overflow checking: a is the correctly rounded value
and since |b| - |c| < |a|, the ternary value value is given
by the sign of a. */
- if (! MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a)))
- {
- inexact = MPFR_INT_SIGN (a);
- goto check_overflow;
- }
- }
- else /* inexact != 0 */
+ ! MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a)))
{
+ MPFR_LOG_MSG (("c small, case 1\n", 0));
+ inexact = MPFR_INT_SIGN (a);
+ }
+ else if (inexact != 0 &&
/* A.AAAAAAAAAAAAAA
= B.BBBBBBBBBBBBBBB
- C.CCCCCCCCCCCCC */
@@ -163,27 +184,41 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
1.BBBBBBBBBBBBB1+1 if inexact == -EVEN_INEX (x == 1)
which means we get a wrong rounded result if x == 1,
i.e. inexact == MPFR_EVEN_INEX (for positive numbers). */
- if (MPFR_LIKELY (inexact != MPFR_EVEN_INEX * MPFR_INT_SIGN (a)))
- goto check_overflow;
+ MPFR_LIKELY (inexact != MPFR_EVEN_INEX * MPFR_INT_SIGN (a)))
+ {
+ MPFR_LOG_MSG (("c small, case 2\n", 0));
+ /* nothing to do */
}
- /* We need to take the value preceding |a|. We can't use
- mpfr_nexttozero due to a possible out-of-range exponent.
- But this will allow us to have more specific code. */
- MPFR_LOG_MSG (("correcting the value of a\n", 0));
- sh = (mpfr_prec_t) an * GMP_NUMB_BITS - aq;
- mpn_sub_1 (ap, ap, an, MPFR_LIMB_ONE << sh);
- if (MPFR_UNLIKELY (MPFR_LIMB_MSB (ap[an-1]) == 0))
+ else
{
- MPFR_EXP (a) --;
- /* The following is valid whether an = 1 or an > 1. */
- ap[an-1] |= MPFR_LIMB_HIGHBIT;
+ /* We need to take the value preceding |a|. We can't use
+ mpfr_nexttozero due to a possible out-of-range exponent.
+ But this will allow us to have more specific code. */
+ MPFR_LOG_MSG (("c small, case 3: correcting the value of a\n", 0));
+ sh = (mpfr_prec_t) an * GMP_NUMB_BITS - aq;
+ mpn_sub_1 (ap, ap, an, MPFR_LIMB_ONE << sh);
+ if (MPFR_UNLIKELY (MPFR_LIMB_MSB (ap[an-1]) == 0))
+ {
+ exp_a --;
+ /* The following is valid whether an = 1 or an > 1. */
+ ap[an-1] |= MPFR_LIMB_HIGHBIT;
+ }
+ inexact = - MPFR_INT_SIGN (a);
}
- inexact = - MPFR_INT_SIGN (a);
- check_overflow:
- if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax))
+ /* The underflow case is possible only with UBF. The overflow case
+ is also possible with normal FP due to rounding. */
+ if (MPFR_UNLIKELY (exp_a > __gmpfr_emax))
return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
- else
- MPFR_RET (inexact);
+ if (MPFR_UNLIKELY (exp_a < __gmpfr_emin))
+ {
+ if (rnd_mode == MPFR_RNDN &&
+ (exp_a < __gmpfr_emin - 1 ||
+ (inexact * MPFR_INT_SIGN (a) >= 0 && mpfr_powerof2_raw (a))))
+ rnd_mode = MPFR_RNDZ;
+ return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
+ }
+ MPFR_SET_EXP (a, exp_a);
+ MPFR_RET (inexact);
}
/* reserve a space to store b aligned with the result, i.e. shifted by
@@ -656,6 +691,15 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
if (MPFR_LIKELY(cancel))
{
cancel -= add_exp; /* OK: add_exp is an int equal to 0 or 1 */
+ MPFR_ASSERTD (cancel >= 0);
+ /* Detect an underflow case to avoid a possible integer overflow
+ with UBF in the computation of exp_a. */
+ if (MPFR_UNLIKELY (exp_b < __gmpfr_emin - 1))
+ {
+ if (rnd_mode == MPFR_RNDN)
+ rnd_mode = MPFR_RNDZ;
+ return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
+ }
exp_a = exp_b - cancel;
/* The following assertion corresponds to a limitation of the MPFR
implementation. It may fail with a 32-bit ABI and huge precisions,
@@ -671,14 +715,17 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
rnd_mode = MPFR_RNDZ;
return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
}
- if (MPFR_UNLIKELY (exp_a > __gmpfr_emax))
+ /* We cannot have an overflow here, except for UBFs. Indeed:
+ exp_a = exp_b - cancel + add_exp <= emax - 1 + 1 <= emax.
+ For UBFs, we can have exp_b > emax. */
+ if (exp_a > __gmpfr_emax)
{
+ MPFR_ASSERTD (exp_b > __gmpfr_emax); /* since exp_b >= exp_a */
return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
}
}
else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */
{
-
/* in case cancel = 0, add_exp can still be 1, in case b is just
below a power of two, c is very small, prec(a) < prec(b),
and rnd=away or nearest */
@@ -687,12 +734,13 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
a subtraction below to avoid a potential integer overflow in
the case exp_b == MPFR_EXP_MAX. */
if (MPFR_UNLIKELY (exp_b > __gmpfr_emax - add_exp))
- {
- return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
- }
+ return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
exp_a = exp_b + add_exp;
+ /* Warning: an underflow can happen for UBFs, for example when
+ mpfr_add is called from mpfr_fmma or mpfr_fmms. */
if (MPFR_UNLIKELY (exp_a < __gmpfr_emin))
goto underflow;
+ MPFR_ASSERTD (exp_a >= __gmpfr_emin);
}
MPFR_SET_EXP (a, exp_a);
/* check that result is msb-normalized */