diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/add1.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/add1.c | 104 |
1 files changed, 73 insertions, 31 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/add1.c b/Build/source/libs/mpfr/mpfr-src/src/add1.c index 7ed4fb1e7c7..c4b053dd032 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/add1.c +++ b/Build/source/libs/mpfr/mpfr-src/src/add1.c @@ -22,27 +22,39 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #include "mpfr-impl.h" -/* compute sign(b) * (|b| + |c|), assuming b and c have same sign, - and are not NaN, Inf, nor zero. */ -int +/* compute sign(b) * (|b| + |c|), assuming that b and c + are not NaN, Inf, nor zero. Assumes EXP(b) >= EXP(c). +*/ +MPFR_HOT_FUNCTION_ATTR int mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) { mp_limb_t *ap, *bp, *cp; mpfr_prec_t aq, bq, cq, aq2; mp_size_t an, bn, cn; - mpfr_exp_t difw, exp; + mpfr_exp_t difw, exp, diff_exp; int sh, rb, fb, inex; - mpfr_uexp_t diff_exp; MPFR_TMP_DECL(marker); - MPFR_ASSERTD(MPFR_IS_PURE_FP(b)); - MPFR_ASSERTD(MPFR_IS_PURE_FP(c)); + MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); + MPFR_ASSERTD (MPFR_IS_PURE_UBF (c)); + MPFR_ASSERTD (! MPFR_UBF_EXP_LESS_P (b, c)); + + if (MPFR_UNLIKELY (MPFR_IS_UBF (b))) + { + exp = mpfr_ubf_zexp2exp (MPFR_ZEXP (b)); + if (exp > __gmpfr_emax) + return mpfr_overflow (a, rnd_mode, MPFR_SIGN (b));; + } + else + exp = MPFR_GET_EXP (b); + + MPFR_ASSERTD (exp <= __gmpfr_emax); MPFR_TMP_MARK(marker); - aq = MPFR_PREC(a); - bq = MPFR_PREC(b); - cq = MPFR_PREC(c); + aq = MPFR_GET_PREC (a); + bq = MPFR_GET_PREC (b); + cq = MPFR_GET_PREC (c); an = MPFR_PREC2LIMBS (aq); /* number of limbs of a */ aq2 = (mpfr_prec_t) an * GMP_NUMB_BITS; @@ -62,19 +74,24 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) if (ap == cp) { cp = bp; } } - else if (MPFR_UNLIKELY(ap == cp)) + else if (ap == cp) { cp = MPFR_TMP_LIMBS_ALLOC (cn); MPN_COPY(cp, ap, cn); } - exp = MPFR_GET_EXP (b); MPFR_SET_SAME_SIGN(a, b); MPFR_UPDATE2_RND_MODE(rnd_mode, MPFR_SIGN(b)); - /* now rnd_mode is either MPFR_RNDN, MPFR_RNDZ or MPFR_RNDA */ - /* Note: exponents can be negative, but the unsigned subtraction is - a modular subtraction, so that one gets the correct result. */ - diff_exp = (mpfr_uexp_t) exp - MPFR_GET_EXP(c); + /* now rnd_mode is either MPFR_RNDN, MPFR_RNDZ, MPFR_RNDA or MPFR_RNDF. */ + if (MPFR_UNLIKELY (MPFR_IS_UBF (c))) + { + MPFR_STAT_STATIC_ASSERT (MPFR_EXP_MAX > MPFR_PREC_MAX); + diff_exp = mpfr_ubf_diff_exp (b, c); + } + else + diff_exp = exp - MPFR_GET_EXP (c); + + MPFR_ASSERTD (diff_exp >= 0); /* * 1. Compute the significant part A', the non-significant bits of A @@ -153,7 +170,7 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) MPN_COPY(a2p - difn, cp + (cn - difn), difn); /* add b to a */ - cc = MPFR_UNLIKELY(an > bn) + cc = an > bn ? mpn_add_n(ap + (an - bn), ap + (an - bn), bp, bn) : mpn_add_n(ap, ap, bp + (bn - an), an); @@ -185,7 +202,7 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) } /* cc */ } /* aq2 > diff_exp */ - /* non-significant bits of a */ + /* zero the non-significant bits of a */ if (MPFR_LIKELY(rb < 0 && sh)) { mp_limb_t mask, bb; @@ -207,7 +224,11 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) } } - /* determine rounding and sticky bits (and possible carry) */ + /* Determine rounding and sticky bits (and possible carry). + In faithful rounding, we may stop two bits after ulp(a): + the approximation is regarded as the number formed by a, + the rounding bit rb and an additional bit fb; and the + corresponding error is < 1/2 ulp of the unrounded result. */ difw = (mpfr_exp_t) an - (mpfr_exp_t) (diff_exp / GMP_NUMB_BITS); /* difw is the number of limbs from b (regarded as having an infinite @@ -234,12 +255,13 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) MPFR_ASSERTD(fb != 0); if (fb > 0) { - if (bb != MP_LIMB_T_MAX) - { - fb = 1; /* c hasn't been taken into account - ==> sticky bit != 0 */ - goto rounding; - } + /* Note: Here, we can round to nearest, but the loop may still + be necessary to determine whether there is a carry from c, + which will have an effect on the ternary value. However, in + faithful rounding, we do not have to determine the ternary + value, so that we can end the loop here. */ + if (bb != MPFR_LIMB_MAX || rnd_mode == MPFR_RNDF) + goto rounding; } else /* fb not initialized yet */ { @@ -249,7 +271,7 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) bb |= MPFR_LIMB_HIGHBIT; } fb = 1; - if (bb != MP_LIMB_T_MAX) + if (bb != MPFR_LIMB_MAX) goto rounding; } @@ -318,10 +340,16 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) } fb = bb != 0; - if (fb && bb != MP_LIMB_T_MAX) + if (fb && bb != MPFR_LIMB_MAX) goto rounding; } /* fb < 0 */ + /* At least two bits after ulp(a) have been read, which is + sufficient for faithful rounding, as we do not need to + determine on which side of a breakpoint the result is. */ + if (rnd_mode == MPFR_RNDF) + goto rounding; + while (bk > 0) { mp_limb_t bb, cc; @@ -368,7 +396,7 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) fb = 1; goto rounding; } - if (fb && bb != MP_LIMB_T_MAX) + if (fb && bb != MPFR_LIMB_MAX) goto rounding; } /* while */ @@ -406,7 +434,7 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) } fb = bb != 0; } /* fb < 0 */ - if (fb) + if (fb || rnd_mode == MPFR_RNDF) goto rounding; while (bk) { @@ -458,6 +486,11 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) rb = cc >> (GMP_NUMB_BITS - 1); cc &= ~MPFR_LIMB_HIGHBIT; } + if (cc == 0 && rnd_mode == MPFR_RNDF) + { + fb = 0; + goto rounding; + } while (cc == 0) { if (ck == 0) @@ -473,8 +506,8 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) } /* fb != 1 */ rounding: - /* rnd_mode should be one of MPFR_RNDN, MPFR_RNDZ or MPFR_RNDA */ - if (MPFR_LIKELY(rnd_mode == MPFR_RNDN)) + /* rnd_mode should be one of MPFR_RNDN, MPFR_RNDF, MPFR_RNDZ or MPFR_RNDA */ + if (MPFR_LIKELY(rnd_mode == MPFR_RNDN || rnd_mode == MPFR_RNDF)) { if (fb == 0) { @@ -530,6 +563,15 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) } set_exponent: + if (MPFR_UNLIKELY (exp < __gmpfr_emin)) /* possible if b and c are UBF's */ + { + if (rnd_mode == MPFR_RNDN && + (exp < __gmpfr_emin - 1 || + (inex >= 0 && mpfr_powerof2_raw (a)))) + rnd_mode = MPFR_RNDZ; + inex = mpfr_underflow (a, rnd_mode, MPFR_SIGN(a)); + goto end_of_add; + } MPFR_SET_EXP (a, exp); end_of_add: |