summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/add1.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/add1.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/add1.c104
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: