diff options
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.c | 69 |
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' */ |