summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/round_prec.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/round_prec.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/round_prec.c69
1 files changed, 38 insertions, 31 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/round_prec.c b/Build/source/libs/mpfr/mpfr-src/src/round_prec.c
index 436adaff718..f03fbac7be7 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/round_prec.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/round_prec.c
@@ -1,7 +1,7 @@
/* mpfr_round_raw_generic, mpfr_round_raw2, mpfr_round_raw, mpfr_prec_round,
mpfr_can_round, mpfr_can_round_raw -- various rounding functions
-Copyright 1999-2019 Free Software Foundation, Inc.
+Copyright 1999-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
@@ -28,6 +28,7 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
#define use_inexp 1
#include "round_raw_generic.c"
+/* mpfr_round_raw_2 is called from mpfr_round_raw2 */
#define mpfr_round_raw_generic mpfr_round_raw_2
#define flag 1
#define use_inexp 0
@@ -194,20 +195,42 @@ mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err,
MPFR_ASSERT_SIGN(neg);
neg = MPFR_IS_NEG_SIGN(neg);
-
- /* Transform RNDD and RNDU to Zero / Away */
MPFR_ASSERTD (neg == 0 || neg == 1);
- /* transform RNDF to RNDN */
+
+ /* For rnd1 and rnd2, transform RNDF / RNDD / RNDU to RNDN / RNDZ / RNDA
+ (with a special case for rnd1 directed rounding, rnd2 = RNDF). */
+
if (rnd1 == MPFR_RNDF)
- rnd1 = MPFR_RNDN;
- if (rnd1 != MPFR_RNDN)
+ rnd1 = MPFR_RNDN; /* transform RNDF to RNDN */
+ else if (rnd1 != MPFR_RNDN)
rnd1 = MPFR_IS_LIKE_RNDZ(rnd1, neg) ? MPFR_RNDZ : MPFR_RNDA;
+
+ MPFR_ASSERTD (rnd1 == MPFR_RNDN ||
+ rnd1 == MPFR_RNDZ ||
+ rnd1 == MPFR_RNDA);
+
if (rnd2 == MPFR_RNDF)
- rnd2 = (rnd1 == MPFR_RNDN) ? MPFR_RNDN :
- MPFR_IS_LIKE_RNDZ(rnd2, neg) ? MPFR_RNDA : MPFR_RNDZ;
- if (rnd2 != MPFR_RNDN)
+ {
+ if (rnd1 == MPFR_RNDN)
+ rnd2 = MPFR_RNDN;
+ else
+ {
+ rnd2 = MPFR_IS_LIKE_RNDZ(rnd1, neg) ? MPFR_RNDA : MPFR_RNDZ;
+ /* Warning: in this case (rnd1 directed rounding, rnd2 = RNDF),
+ the specification of mpfr_can_round says that we should
+ return non-zero (i.e., we can round) when {bp, bn} is
+ exactly representable in precision prec. */
+ if (mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec) == 0)
+ return 1;
+ }
+ }
+ else if (rnd2 != MPFR_RNDN)
rnd2 = MPFR_IS_LIKE_RNDZ(rnd2, neg) ? MPFR_RNDZ : MPFR_RNDA;
+ MPFR_ASSERTD (rnd2 == MPFR_RNDN ||
+ rnd2 == MPFR_RNDZ ||
+ rnd2 == MPFR_RNDA);
+
/* For err < prec (+1 for rnd1=RNDN), we can never round correctly, since
the error is at least 2*ulp(b) >= ulp(round(b)).
However for err = prec (+1 for rnd1=RNDN), we can round correctly in some
@@ -324,24 +347,8 @@ mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err,
}
return 1;
}
- else if (rnd1 == rnd2)
- {
- if (rnd1 == MPFR_RNDN && prec < (mpfr_prec_t) bn * GMP_NUMB_BITS)
- {
- /* then rnd2 = RNDN, and for prec = bn * GMP_NUMB_BITS we cannot
- have b the middle of two representable numbers */
- k1 = MPFR_PREC2LIMBS (prec + 1);
- MPFR_UNSIGNED_MINUS_MODULO(s1, prec + 1);
- if (((bp[bn - k1] >> s1) & 1) &&
- mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec + 1) == 0)
- /* b is representable in precision prec+1 and ends with a 1 */
- return 0;
- else
- return 1;
- }
- else
- return 1;
- }
+ else if (rnd1 == rnd2) /* cases RNDZ RNDZ or RNDA RNDA: ok */
+ return 1;
else
return mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec) != 0;
}
@@ -376,7 +383,7 @@ mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err,
k++; /* since we work with k+1 everywhere */
tmp = MPFR_TMP_LIMBS_ALLOC (tn);
if (bn > k)
- MPN_COPY (tmp, bp, bn - k);
+ MPN_COPY (tmp, bp, bn - k); /* copy low bn-k limbs of b into tmp */
MPFR_ASSERTD (k > 0);
@@ -384,7 +391,7 @@ mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err,
{
case MPFR_RNDZ:
/* rnd1 = Round to Zero */
- cc = (bp[bn - 1] >> s1) & 1;
+ cc = (bp[bn - 1] >> s1) & 1; /* cc is the least significant bit of b */
/* mpfr_round_raw2 returns 1 if one should add 1 at ulp(b,prec),
and 0 otherwise */
cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec2);
@@ -394,7 +401,7 @@ mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err,
cy = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
/* propagate carry up to most significant limb */
for (tn = 0; tn + 1 < k1 && cy != 0; tn ++)
- cy = ~bp[bn + tn] == 0;
+ cy = bp[bn + tn] == MPFR_LIMB_MAX;
if (cy == 0 && err == prec)
{
res = 0;
@@ -424,7 +431,7 @@ mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err,
cy = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
/* propagate carry up to most significant limb */
for (tn = 0; tn + 1 < k1 && cy != 0; tn ++)
- cy = ~bp[bn + tn] == 0;
+ cy = bp[bn + tn] == MPFR_LIMB_MAX;
cc = (tmp[bn - 1] >> s1) & 1; /* gives 0 when cc=1 */
cc ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2);
/* cc is the new value of bit s1 in bp[bn-1]+eps after rounding 'rnd2' */