summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/sub1sp.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/sub1sp.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/sub1sp.c1418
1 files changed, 1211 insertions, 207 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/sub1sp.c b/Build/source/libs/mpfr/mpfr-src/src/sub1sp.c
index 67d597677bb..bc776ee1b55 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/sub1sp.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/sub1sp.c
@@ -1,4 +1,4 @@
-/* mpfr_sub1sp -- internal function to perform a "real" substraction
+/* mpfr_sub1sp -- internal function to perform a "real" subtraction
All the op must have the same precision
Copyright 2003-2017 Free Software Foundation, Inc.
@@ -25,15 +25,20 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
#include "mpfr-impl.h"
/* Check if we have to check the result of mpfr_sub1sp with mpfr_sub1 */
-#ifdef MPFR_WANT_ASSERT
-# if MPFR_WANT_ASSERT >= 2
+#if MPFR_WANT_ASSERT >= 2
-int mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode);
+int mpfr_sub1sp_ref (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t);
int mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
mpfr_t tmpa, tmpb, tmpc;
+ mpfr_flags_t old_flags, flags, flags2;
int inexb, inexc, inexact, inexact2;
+ if (rnd_mode == MPFR_RNDF)
+ return mpfr_sub1sp_ref (a, b, c, rnd_mode);
+
+ old_flags = __gmpfr_flags;
+
mpfr_init2 (tmpa, MPFR_PREC (a));
mpfr_init2 (tmpb, MPFR_PREC (b));
mpfr_init2 (tmpc, MPFR_PREC (c));
@@ -44,40 +49,993 @@ int mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
inexc = mpfr_set (tmpc, c, MPFR_RNDN);
MPFR_ASSERTN (inexc == 0);
+ MPFR_ASSERTN (__gmpfr_flags == old_flags);
+
inexact2 = mpfr_sub1 (tmpa, tmpb, tmpc, rnd_mode);
- inexact = mpfr_sub1sp2(a, b, c, rnd_mode);
+ flags2 = __gmpfr_flags;
- if (mpfr_cmp (tmpa, a) || inexact != inexact2)
+ __gmpfr_flags = old_flags;
+ inexact = mpfr_sub1sp_ref (a, b, c, rnd_mode);
+ flags = __gmpfr_flags;
+
+ if (! mpfr_equal_p (tmpa, a) || inexact != inexact2 || flags != flags2)
{
fprintf (stderr, "sub1 & sub1sp return different values for %s\n"
"Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nB = ",
- mpfr_print_rnd_mode (rnd_mode), (unsigned long) MPFR_PREC (a),
- (unsigned long) MPFR_PREC (b), (unsigned long) MPFR_PREC (c));
- mpfr_fprint_binary (stderr, tmpb);
- fprintf (stderr, "\nC = ");
- mpfr_fprint_binary (stderr, tmpc);
- fprintf (stderr, "\nSub1 : ");
- mpfr_fprint_binary (stderr, tmpa);
- fprintf (stderr, "\nSub1sp: ");
- mpfr_fprint_binary (stderr, a);
- fprintf (stderr, "\nInexact sp = %d | Inexact = %d\n",
- inexact, inexact2);
+ mpfr_print_rnd_mode (rnd_mode),
+ (unsigned long) MPFR_PREC (a),
+ (unsigned long) MPFR_PREC (b),
+ (unsigned long) MPFR_PREC (c));
+ mpfr_fdump (stderr, tmpb);
+ fprintf (stderr, "C = ");
+ mpfr_fdump (stderr, tmpc);
+ fprintf (stderr, "sub1 : ");
+ mpfr_fdump (stderr, tmpa);
+ fprintf (stderr, "sub1sp: ");
+ mpfr_fdump (stderr, a);
+ fprintf (stderr, "Inexact sp = %d | Inexact = %d\n"
+ "Flags sp = %u | Flags = %u\n",
+ inexact, inexact2, flags, flags2);
MPFR_ASSERTN (0);
}
mpfr_clears (tmpa, tmpb, tmpc, (mpfr_ptr) 0);
return inexact;
}
-# define mpfr_sub1sp mpfr_sub1sp2
-# endif
-#endif
+# define mpfr_sub1sp mpfr_sub1sp_ref
+#endif /* MPFR_WANT_ASSERT >= 2 */
+
+#if !defined(MPFR_GENERIC_ABI)
+
+/* special code for p < GMP_NUMB_BITS */
+static int
+mpfr_sub1sp1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
+ mpfr_prec_t p)
+{
+ mpfr_exp_t bx = MPFR_GET_EXP (b);
+ mpfr_exp_t cx = MPFR_GET_EXP (c);
+ mp_limb_t *ap = MPFR_MANT(a);
+ mp_limb_t *bp = MPFR_MANT(b);
+ mp_limb_t *cp = MPFR_MANT(c);
+ mpfr_prec_t cnt, INITIALIZED(sh);
+ mp_limb_t rb; /* round bit */
+ mp_limb_t sb; /* sticky bit */
+ mp_limb_t a0;
+ mp_limb_t mask;
+ mpfr_uexp_t d;
+
+ MPFR_ASSERTD(p < GMP_NUMB_BITS);
+
+ if (bx == cx)
+ {
+ a0 = bp[0] - cp[0];
+ if (a0 == 0) /* result is zero */
+ {
+ if (rnd_mode == MPFR_RNDD)
+ MPFR_SET_NEG(a);
+ else
+ MPFR_SET_POS(a);
+ MPFR_SET_ZERO(a);
+ MPFR_RET (0);
+ }
+ else if (a0 > bp[0]) /* borrow: |c| > |b| */
+ {
+ MPFR_SET_OPPOSITE_SIGN (a, b);
+ a0 = -a0;
+ }
+ else /* bp[0] > cp[0] */
+ MPFR_SET_SAME_SIGN (a, b);
+
+ /* now a0 != 0 */
+ MPFR_ASSERTD(a0 != 0);
+ count_leading_zeros (cnt, a0);
+ ap[0] = a0 << cnt;
+ bx -= cnt;
+ rb = sb = 0;
+ /* Note: sh is not initialized, but will not be used in this case. */
+ }
+ else if (bx > cx)
+ {
+ MPFR_SET_SAME_SIGN (a, b);
+ BGreater1:
+ d = (mpfr_uexp_t) bx - cx;
+ sh = GMP_NUMB_BITS - p;
+ mask = MPFR_LIMB_MASK(sh);
+ if (d < GMP_NUMB_BITS)
+ {
+ sb = - (cp[0] << (GMP_NUMB_BITS - d)); /* neglected part of -c */
+ /* Note that "a0 = bp[0] - (cp[0] >> d) - (sb != 0);" is faster
+ on some other machines and has no immediate dependencies for
+ the first subtraction. In the future, make sure that the code
+ is recognized as a *single* subtraction with borrow and/or use
+ a builtin when available (currently provided by Clang, but not
+ by GCC); create a new macro for that. See the TODO later.
+ Note also that with Clang, the constant 0 for the first
+ subtraction instead of a variable breaks the optimization:
+ https://llvm.org/bugs/show_bug.cgi?id=31754 */
+ a0 = bp[0] - (sb != 0) - (cp[0] >> d);
+ /* a0 cannot be zero here since:
+ a) if d >= 2, then a0 >= 2^(w-1) - (2^(w-2)-1) with
+ w = GMP_NUMB_BITS, thus a0 - 1 >= 2^(w-2),
+ b) if d = 1, then since p < GMP_NUMB_BITS we have sb=0.
+ */
+ MPFR_ASSERTD(a0 > 0);
+ count_leading_zeros (cnt, a0);
+ if (cnt)
+ a0 = (a0 << cnt) | (sb >> (GMP_NUMB_BITS - cnt));
+ sb <<= cnt;
+ bx -= cnt;
+ /* sh > 0 since p < GMP_NUMB_BITS */
+ MPFR_ASSERTD(sh > 0);
+ rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
+ sb |= (a0 & mask) ^ rb;
+ ap[0] = a0 & ~mask;
+ }
+ else /* d >= GMP_NUMB_BITS */
+ {
+ if (bp[0] > MPFR_LIMB_HIGHBIT)
+ {
+ /* We compute b - ulp(b), and the remainder ulp(b) - c satisfies:
+ 1/2 ulp(b) < ulp(b) - c < ulp(b), thus rb = sb = 1. */
+ ap[0] = bp[0] - (MPFR_LIMB_ONE << sh);
+ rb = 1;
+ }
+ else
+ {
+ /* Warning: since we have an exponent decrease, when
+ p = GMP_NUMB_BITS - 1 and d = GMP_NUMB_BITS, the round bit
+ corresponds to the upper bit of -c. In that case rb = 0 and
+ sb = 1, except when c0 = MPFR_LIMB_HIGHBIT where rb = 1 and
+ sb = 0. */
+ rb = sh > 1 || d > GMP_NUMB_BITS || cp[0] == MPFR_LIMB_HIGHBIT;
+ /* sb=1 below is incorrect when p = GMP_NUMB_BITS - 1,
+ d = GMP_NUMB_BITS and c0 = MPFR_LIMB_HIGHBIT, but in
+ that case the even rule wound round up too. */
+ ap[0] = ~mask;
+ bx --;
+ /* Warning: if d = GMP_NUMB_BITS and c0 = 1000...000, then
+ b0 - c0 = |0111...111|1000...000|, which after the shift
+ becomes |111...111|000...000| thus if p = GMP_NUMB_BITS-1
+ we have rb = 1 but sb = 0. However in this case the round
+ even rule will round up, which is what we get with sb = 1:
+ the final result will be correct, while sb is incorrect. */
+ }
+ sb = 1;
+ }
+ }
+ else /* cx > bx */
+ {
+ mpfr_exp_t tx;
+ mp_limb_t *tp;
+ tx = bx; bx = cx; cx = tx;
+ tp = bp; bp = cp; cp = tp;
+ MPFR_SET_OPPOSITE_SIGN (a, b);
+ goto BGreater1;
+ }
+
+ /* now perform rounding */
+
+ /* Warning: MPFR considers underflow *after* rounding with an unbounded
+ exponent range. However since b and c have same precision p, they are
+ multiples of 2^(emin-p), likewise for b-c. Thus if bx < emin, the
+ subtraction (with an unbounded exponent range) is exact, so that bx is
+ also the exponent after rounding with an unbounded exponent range. */
+ if (MPFR_UNLIKELY(bx < __gmpfr_emin))
+ {
+ /* For RNDN, mpfr_underflow always rounds away, thus for |a| <= 2^(emin-2)
+ we have to change to RNDZ. This corresponds to:
+ (a) either bx < emin - 1
+ (b) or bx = emin - 1 and ap[0] = 1000....000 (in this case necessarily
+ rb = sb = 0 since b-c is multiple of 2^(emin-p) */
+ if (rnd_mode == MPFR_RNDN &&
+ (bx < __gmpfr_emin - 1 || ap[0] == MPFR_LIMB_HIGHBIT))
+ {
+ MPFR_ASSERTD(rb == 0 && sb == 0);
+ rnd_mode = MPFR_RNDZ;
+ }
+ return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
+ }
+
+ MPFR_SET_EXP (a, bx);
+ if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF)
+ MPFR_RET (0);
+ else if (rnd_mode == MPFR_RNDN)
+ {
+ if (rb == 0 || (sb == 0 && (ap[0] & (MPFR_LIMB_ONE << sh)) == 0))
+ goto truncate;
+ else
+ goto add_one_ulp;
+ }
+ else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
+ {
+ truncate:
+ MPFR_RET(-MPFR_SIGN(a));
+ }
+ else /* round away from zero */
+ {
+ add_one_ulp:
+ ap[0] += MPFR_LIMB_ONE << sh;
+ if (MPFR_UNLIKELY(ap[0] == 0))
+ {
+ ap[0] = MPFR_LIMB_HIGHBIT;
+ /* Note: bx+1 cannot exceed __gmpfr_emax, since |a| <= |b|, thus
+ bx+1 is at most equal to the original exponent of b. */
+ MPFR_ASSERTD(bx + 1 <= __gmpfr_emax);
+ MPFR_SET_EXP (a, bx + 1);
+ }
+ MPFR_RET(MPFR_SIGN(a));
+ }
+}
+
+/* special code for p = GMP_NUMB_BITS */
+static int
+mpfr_sub1sp1n (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
+{
+ mpfr_exp_t bx = MPFR_GET_EXP (b);
+ mpfr_exp_t cx = MPFR_GET_EXP (c);
+ mp_limb_t *ap = MPFR_MANT(a);
+ mp_limb_t *bp = MPFR_MANT(b);
+ mp_limb_t *cp = MPFR_MANT(c);
+ mpfr_prec_t cnt;
+ mp_limb_t rb; /* round bit */
+ mp_limb_t sb; /* sticky bit */
+ mp_limb_t a0;
+ mpfr_uexp_t d;
+
+ MPFR_ASSERTD(MPFR_PREC(a) == GMP_NUMB_BITS);
+ MPFR_ASSERTD(MPFR_PREC(b) == GMP_NUMB_BITS);
+ MPFR_ASSERTD(MPFR_PREC(c) == GMP_NUMB_BITS);
+
+ if (bx == cx)
+ {
+ a0 = bp[0] - cp[0];
+ if (a0 == 0) /* result is zero */
+ {
+ if (rnd_mode == MPFR_RNDD)
+ MPFR_SET_NEG(a);
+ else
+ MPFR_SET_POS(a);
+ MPFR_SET_ZERO(a);
+ MPFR_RET (0);
+ }
+ else if (a0 > bp[0]) /* borrow: |c| > |b| */
+ {
+ MPFR_SET_OPPOSITE_SIGN (a, b);
+ a0 = -a0;
+ }
+ else /* bp[0] > cp[0] */
+ MPFR_SET_SAME_SIGN (a, b);
+
+ /* now a0 != 0 */
+ MPFR_ASSERTD(a0 != 0);
+ count_leading_zeros (cnt, a0);
+ ap[0] = a0 << cnt;
+ bx -= cnt;
+ rb = sb = 0;
+ }
+ else if (bx > cx)
+ {
+ MPFR_SET_SAME_SIGN (a, b);
+ BGreater1:
+ d = (mpfr_uexp_t) bx - cx;
+ if (d < GMP_NUMB_BITS)
+ {
+ sb = - (cp[0] << (GMP_NUMB_BITS - d)); /* neglected part of -c */
+ a0 = bp[0] - (sb != 0) - (cp[0] >> d);
+ /* a0 can only be zero when d=1, b0 = B/2, and c0 = B-1, where
+ B = 2^GMP_NUMB_BITS, thus b0 - c0/2 = 1/2 */
+ if (a0 == MPFR_LIMB_ZERO)
+ {
+ bx -= GMP_NUMB_BITS;
+ ap[0] = MPFR_LIMB_HIGHBIT;
+ rb = sb = 0;
+ }
+ else
+ {
+ count_leading_zeros (cnt, a0);
+ if (cnt)
+ a0 = (a0 << cnt) | (sb >> (GMP_NUMB_BITS - cnt));
+ sb <<= cnt;
+ bx -= cnt;
+ rb = sb & MPFR_LIMB_HIGHBIT;
+ sb &= ~MPFR_LIMB_HIGHBIT;
+ ap[0] = a0;
+ }
+ }
+ else /* d >= GMP_NUMB_BITS */
+ {
+ /* We compute b - ulp(b) */
+ if (bp[0] > MPFR_LIMB_HIGHBIT)
+ {
+ /* If d = GMP_NUMB_BITS, rb = 0 and sb = 1,
+ unless c0 = MPFR_LIMB_HIGHBIT in which case rb = 1 and sb = 0.
+ If d > GMP_NUMB_BITS, rb = sb = 1. */
+ rb = d > GMP_NUMB_BITS || cp[0] == MPFR_LIMB_HIGHBIT;
+ sb = d > GMP_NUMB_BITS || cp[0] != MPFR_LIMB_HIGHBIT;
+ ap[0] = bp[0] - MPFR_LIMB_ONE;
+ }
+ else
+ {
+ /* Warning: in this case a0 is shifted by one!
+ If d=GMP_NUMB_BITS:
+ (a) if c0 = MPFR_LIMB_HIGHBIT, a0 = 111...111, rb = sb = 0
+ (b) otherwise, a0 = 111...110, rb = -c0 >= 01000...000,
+ sb = (-c0) << 2
+ If d=GMP_NUMB_BITS+1: a0 = 111...111
+ (c) if c0 = MPFR_LIMB_HIGHBIT, rb = 1 and sb = 0
+ (d) otherwise rb = 0 and sb = 1
+ If d > GMP_NUMB_BITS+1:
+ (e) a0 = 111...111, rb = sb = 1
+ */
+ bx --;
+ if (d == GMP_NUMB_BITS && cp[0] > MPFR_LIMB_HIGHBIT)
+ { /* case (b) */
+ rb = (-cp[0]) >= (MPFR_LIMB_HIGHBIT >> 1);
+ sb = (-cp[0]) << 2;
+ ap[0] = -(MPFR_LIMB_ONE << 1);
+ }
+ else /* cases (a), (c), (d) and (e) */
+ {
+ ap[0] = -MPFR_LIMB_ONE;
+ /* rb=1 in case (e) and case (c) */
+ rb = d > GMP_NUMB_BITS + 1
+ || (d == GMP_NUMB_BITS + 1 && cp[0] == MPFR_LIMB_HIGHBIT);
+ /* sb = 1 in case (d) and (e) */
+ sb = d > GMP_NUMB_BITS + 1
+ || (d == GMP_NUMB_BITS + 1 && cp[0] > MPFR_LIMB_HIGHBIT);
+ }
+ }
+ }
+ }
+ else /* cx > bx */
+ {
+ mpfr_exp_t tx;
+ mp_limb_t *tp;
+ tx = bx; bx = cx; cx = tx;
+ tp = bp; bp = cp; cp = tp;
+ MPFR_SET_OPPOSITE_SIGN (a, b);
+ goto BGreater1;
+ }
+
+ /* now perform rounding */
+
+ /* Warning: MPFR considers underflow *after* rounding with an unbounded
+ exponent range. However since b and c have same precision p, they are
+ multiples of 2^(emin-p), likewise for b-c. Thus if bx < emin, the
+ subtraction (with an unbounded exponent range) is exact, so that bx is
+ also the exponent after rounding with an unbounded exponent range. */
+ if (MPFR_UNLIKELY(bx < __gmpfr_emin))
+ {
+ /* For RNDN, mpfr_underflow always rounds away, thus for |a| <= 2^(emin-2)
+ we have to change to RNDZ. This corresponds to:
+ (a) either bx < emin - 1
+ (b) or bx = emin - 1 and ap[0] = 1000....000 (in this case necessarily
+ rb = sb = 0 since b-c is multiple of 2^(emin-p) */
+ if (rnd_mode == MPFR_RNDN &&
+ (bx < __gmpfr_emin - 1 || ap[0] == MPFR_LIMB_HIGHBIT))
+ {
+ MPFR_ASSERTD(rb == 0 && sb == 0);
+ rnd_mode = MPFR_RNDZ;
+ }
+ return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
+ }
-/* Debugging support */
-#ifdef DEBUG
-# undef DEBUG
-# define DEBUG(x) (x)
+ MPFR_SET_EXP (a, bx);
+ if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF)
+ MPFR_RET (0);
+ else if (rnd_mode == MPFR_RNDN)
+ {
+ if (rb == 0 || (sb == 0 && (ap[0] & MPFR_LIMB_ONE) == 0))
+ goto truncate;
+ else
+ goto add_one_ulp;
+ }
+ else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
+ {
+ truncate:
+ MPFR_RET(-MPFR_SIGN(a));
+ }
+ else /* round away from zero */
+ {
+ add_one_ulp:
+ ap[0] += MPFR_LIMB_ONE;
+ if (MPFR_UNLIKELY(ap[0] == 0))
+ {
+ ap[0] = MPFR_LIMB_HIGHBIT;
+ /* Note: bx+1 cannot exceed __gmpfr_emax, since |a| <= |b|, thus
+ bx+1 is at most equal to the original exponent of b. */
+ MPFR_ASSERTD(bx + 1 <= __gmpfr_emax);
+ MPFR_SET_EXP (a, bx + 1);
+ }
+ MPFR_RET(MPFR_SIGN(a));
+ }
+}
+
+/* special code for GMP_NUMB_BITS < p < 2*GMP_NUMB_BITS */
+static int
+mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
+ mpfr_prec_t p)
+{
+ mpfr_exp_t bx = MPFR_GET_EXP (b);
+ mpfr_exp_t cx = MPFR_GET_EXP (c);
+ mp_limb_t *ap = MPFR_MANT(a);
+ mp_limb_t *bp = MPFR_MANT(b);
+ mp_limb_t *cp = MPFR_MANT(c);
+ mpfr_prec_t cnt, INITIALIZED(sh);
+ mp_limb_t rb; /* round bit */
+ mp_limb_t sb; /* sticky bit */
+ mp_limb_t mask, a0, a1;
+ mpfr_uexp_t d;
+
+ MPFR_ASSERTD(GMP_NUMB_BITS < p && p < 2 * GMP_NUMB_BITS);
+
+ if (bx == cx) /* subtraction is exact in this case */
+ {
+ /* first compute a0: if the compiler is smart enough, it will use the generated
+ borrow to get for free the term (bp[0] < cp[0]) */
+ a0 = bp[0] - cp[0];
+ a1 = bp[1] - cp[1] - (bp[0] < cp[0]);
+ if (a1 == 0 && a0 == 0) /* result is zero */
+ {
+ if (rnd_mode == MPFR_RNDD)
+ MPFR_SET_NEG(a);
+ else
+ MPFR_SET_POS(a);
+ MPFR_SET_ZERO(a);
+ MPFR_RET (0);
+ }
+ else if (a1 >= bp[1]) /* borrow: |c| > |b| */
+ {
+ MPFR_SET_OPPOSITE_SIGN (a, b);
+ /* a = b-c mod 2^(2*GMP_NUMB_BITS) */
+ a0 = -a0;
+ a1 = -a1 - (a0 != 0);
+ }
+ else /* bp[0] > cp[0] */
+ MPFR_SET_SAME_SIGN (a, b);
+
+ if (a1 == 0)
+ {
+ a1 = a0;
+ a0 = 0;
+ bx -= GMP_NUMB_BITS;
+ }
+
+ /* now a1 != 0 */
+ MPFR_ASSERTD(a1 != 0);
+ count_leading_zeros (cnt, a1);
+ if (cnt)
+ {
+ ap[1] = (a1 << cnt) | (a0 >> (GMP_NUMB_BITS - cnt));
+ ap[0] = a0 << cnt;
+ bx -= cnt;
+ }
+ else
+ {
+ ap[1] = a1;
+ ap[0] = a0;
+ }
+ rb = sb = 0;
+ /* Note: sh is not initialized, but will not be used in this case. */
+ }
+ else if (bx > cx)
+ {
+ mp_limb_t t;
+
+ MPFR_SET_SAME_SIGN (a, b);
+ BGreater2:
+ d = (mpfr_uexp_t) bx - cx;
+ sh = 2 * GMP_NUMB_BITS - p;
+ mask = MPFR_LIMB_MASK(sh);
+ if (d < GMP_NUMB_BITS)
+ {
+ t = (cp[1] << (GMP_NUMB_BITS - d)) | (cp[0] >> d);
+ /* TODO: Change the code to generate a full subtraction with borrow,
+ avoiding the test on sb and the corresponding correction. Note
+ that Clang has builtins:
+ http://clang.llvm.org/docs/LanguageExtensions.html#multiprecision-arithmetic-builtins
+ but the generated code may not be good:
+ https://llvm.org/bugs/show_bug.cgi?id=20748
+ With the current source code, Clang generates on x86_64:
+ 1. sub %rsi,%rbx for the first subtraction in a1;
+ 2. sub %rdi,%rax for the subtraction in a0;
+ 3. sbb $0x0,%rbx for the second subtraction in a1, i.e. just
+ subtracting the borrow out from (2).
+ So, Clang recognizes the borrow, but doesn't merge (1) and (3).
+ Bug: https://llvm.org/bugs/show_bug.cgi?id=25858
+ For GCC: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=79173
+ */
+#ifndef MPFR_FULLSUB
+ a0 = bp[0] - t;
+ a1 = bp[1] - (cp[1] >> d) - (bp[0] < t);
+ sb = cp[0] << (GMP_NUMB_BITS - d); /* neglected part of c */
+ if (sb)
+ {
+ a1 -= (a0 == 0);
+ a0 --;
+ /* a = a1,a0 cannot become zero here, since:
+ a) if d >= 2, then a1 >= 2^(w-1) - (2^(w-2)-1) with
+ w = GMP_NUMB_BITS, thus a1 - 1 >= 2^(w-2),
+ b) if d = 1, then since p < 2*GMP_NUMB_BITS we have sb=0. */
+ MPFR_ASSERTD(a1 > 0 || a0 > 0);
+ sb = -sb; /* 2^GMP_NUMB_BITS - sb */
+ }
#else
-# define DEBUG(x) /**/
+ sb = - (cp[0] << (GMP_NUMB_BITS - d));
+ a0 = bp[0] - t - (sb != 0);
+ a1 = bp[1] - (cp[1] >> d) - (bp[0] < t || (bp[0] == t && sb != 0));
#endif
+ if (a1 == 0)
+ {
+ /* this implies d=1, which in turn implies sb=0 */
+ MPFR_ASSERTD(sb == 0);
+ a1 = a0;
+ a0 = 0; /* should be a0 = sb */
+ /* since sb=0 already, no need to set it to 0 */
+ bx -= GMP_NUMB_BITS;
+ }
+ /* now a1 != 0 */
+ MPFR_ASSERTD(a1 != 0);
+ count_leading_zeros (cnt, a1);
+ if (cnt)
+ {
+ ap[1] = (a1 << cnt) | (a0 >> (GMP_NUMB_BITS - cnt));
+ a0 = (a0 << cnt) | (sb >> (GMP_NUMB_BITS - cnt));
+ sb <<= cnt;
+ bx -= cnt;
+ }
+ else
+ ap[1] = a1;
+ /* sh > 0 since p < 2*GMP_NUMB_BITS */
+ MPFR_ASSERTD(sh > 0);
+ rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
+ sb |= (a0 & mask) ^ rb;
+ ap[0] = a0 & ~mask;
+ }
+ else if (d < 2 * GMP_NUMB_BITS)
+ { /* GMP_NUMB_BITS <= d < 2*GMP_NUMB_BITS */
+ /* warning: the most significant bit of sb might become the least
+ significant bit of a0 below */
+ sb = (d == GMP_NUMB_BITS) ? cp[0]
+ : (cp[1] << (2*GMP_NUMB_BITS - d)) | (cp[0] != 0);
+ t = (cp[1] >> (d - GMP_NUMB_BITS)) + (sb != 0);
+ /* warning: t might overflow to 0 if d=GMP_NUMB_BITS and sb <> 0 */
+ a0 = bp[0] - t;
+ a1 = bp[1] - (bp[0] < t) - (t == 0 && sb != 0);
+ sb = -sb;
+ /* since bp[1] has its most significant bit set, we can have an
+ exponent decrease of at most one */
+ if (a1 < MPFR_LIMB_HIGHBIT)
+ {
+ ap[1] = (a1 << 1) | (a0 >> (GMP_NUMB_BITS - 1));
+ a0 = (a0 << 1) | (sb >> (GMP_NUMB_BITS - 1));
+ sb <<= 1;
+ bx --;
+ }
+ else
+ ap[1] = a1;
+ rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
+ sb |= (a0 & mask) ^ rb;
+ ap[0] = a0 & ~mask;
+ }
+ else /* d >= 2*GMP_NUMB_BITS */
+ {
+ /* We compute b - ulp(b), and the remainder ulp(b) - c satisfies:
+ 1/2 ulp(b) < ulp(b) - c < ulp(b), thus rb = sb = 1, unless we
+ had an exponent decrease. */
+ t = MPFR_LIMB_ONE << sh;
+ a0 = bp[0] - t;
+ a1 = bp[1] - (bp[0] < t);
+ if (a1 < MPFR_LIMB_HIGHBIT)
+ {
+ /* necessarily we had b = 1000...000 */
+ /* Warning: since we have an exponent decrease, when
+ p = 2*GMP_NUMB_BITS - 1 and d = 2*GMP_NUMB_BITS, the round bit
+ corresponds to the upper bit of -c. In that case rb = 0 and
+ sb = 1, except when c = 1000...000 where rb = 1 and sb = 0. */
+ rb = sh > 1 || d > 2 * GMP_NUMB_BITS
+ || (cp[1] == MPFR_LIMB_HIGHBIT && cp[0] == MPFR_LIMB_ZERO);
+ /* sb=1 below is incorrect when p = 2*GMP_NUMB_BITS - 1,
+ d = 2*GMP_NUMB_BITS and c = 1000...000, but in
+ that case the even rule wound round up too. */
+ ap[0] = ~mask;
+ ap[1] = MPFR_LIMB_MAX;
+ bx --;
+ }
+ else
+ {
+ ap[0] = a0;
+ ap[1] = a1;
+ rb = 1;
+ }
+ sb = 1;
+ }
+ }
+ else /* cx > bx */
+ {
+ mpfr_exp_t tx;
+ mp_limb_t *tp;
+ tx = bx; bx = cx; cx = tx;
+ tp = bp; bp = cp; cp = tp;
+ MPFR_SET_OPPOSITE_SIGN (a, b);
+ goto BGreater2;
+ }
+
+ /* now perform rounding */
+
+ /* Warning: MPFR considers underflow *after* rounding with an unbounded
+ exponent range. However since b and c have same precision p, they are
+ multiples of 2^(emin-p), likewise for b-c. Thus if bx < emin, the
+ subtraction (with an unbounded exponent range) is exact, so that bx is
+ also the exponent after rounding with an unbounded exponent range. */
+ if (MPFR_UNLIKELY(bx < __gmpfr_emin))
+ {
+ /* for RNDN, mpfr_underflow always rounds away, thus for |a|<=2^(emin-2)
+ we have to change to RNDZ */
+ if (rnd_mode == MPFR_RNDN &&
+ (bx < __gmpfr_emin - 1 ||
+ (ap[1] == MPFR_LIMB_HIGHBIT && ap[0] == 0)))
+ rnd_mode = MPFR_RNDZ;
+ return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
+ }
+
+ MPFR_SET_EXP (a, bx);
+ if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF)
+ MPFR_RET (0);
+ else if (rnd_mode == MPFR_RNDN)
+ {
+ if (rb == 0 || (sb == 0 && (ap[0] & (MPFR_LIMB_ONE << sh)) == 0))
+ goto truncate;
+ else
+ goto add_one_ulp;
+ }
+ else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
+ {
+ truncate:
+ MPFR_RET(-MPFR_SIGN(a));
+ }
+ else /* round away from zero */
+ {
+ add_one_ulp:
+ ap[0] += MPFR_LIMB_ONE << sh;
+ ap[1] += (ap[0] == 0);
+ if (MPFR_UNLIKELY(ap[1] == 0))
+ {
+ ap[1] = MPFR_LIMB_HIGHBIT;
+ /* Note: bx+1 cannot exceed __gmpfr_emax, since |a| <= |b|, thus
+ bx+1 is at most equal to the original exponent of b. */
+ MPFR_ASSERTD(bx + 1 <= __gmpfr_emax);
+ MPFR_SET_EXP (a, bx + 1);
+ }
+ MPFR_RET(MPFR_SIGN(a));
+ }
+}
+
+/* special code for 2*GMP_NUMB_BITS < p < 3*GMP_NUMB_BITS */
+static int
+mpfr_sub1sp3 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
+ mpfr_prec_t p)
+{
+ mpfr_exp_t bx = MPFR_GET_EXP (b);
+ mpfr_exp_t cx = MPFR_GET_EXP (c);
+ mp_limb_t *ap = MPFR_MANT(a);
+ mp_limb_t *bp = MPFR_MANT(b);
+ mp_limb_t *cp = MPFR_MANT(c);
+ mpfr_prec_t cnt, INITIALIZED(sh);
+ mp_limb_t rb; /* round bit */
+ mp_limb_t sb; /* sticky bit */
+ mp_limb_t mask, a0, a1, a2;
+ mpfr_uexp_t d;
+
+ MPFR_ASSERTD(2 * GMP_NUMB_BITS < p && p < 3 * GMP_NUMB_BITS);
+
+ if (bx == cx) /* subtraction is exact in this case */
+ {
+ a0 = bp[0] - cp[0];
+ a1 = bp[1] - cp[1] - (bp[0] < cp[0]);
+ /* a borrow is generated for a when either bp[1] < cp[1],
+ or bp[1] = cp[1] and bp[0] < cp[0] */
+ a2 = bp[2] - cp[2]
+ - (bp[1] < cp[1] || (bp[1] == cp[1] && bp[0] < cp[0]));
+ if (a2 == 0 && a1 == 0 && a0 == 0) /* result is zero */
+ {
+ if (rnd_mode == MPFR_RNDD)
+ MPFR_SET_NEG(a);
+ else
+ MPFR_SET_POS(a);
+ MPFR_SET_ZERO(a);
+ MPFR_RET (0);
+ }
+ else if (a2 >= bp[2]) /* borrow: |c| > |b| */
+ {
+ MPFR_SET_OPPOSITE_SIGN (a, b);
+ /* a = b-c mod 2^(3*GMP_NUMB_BITS) */
+ a0 = -a0;
+ a1 = -a1 - (a0 != 0);
+ a2 = -a2 - (a0 != 0 || a1 != 0);
+ }
+ else /* bp[0] > cp[0] */
+ MPFR_SET_SAME_SIGN (a, b);
+
+ if (a2 == 0)
+ {
+ a2 = a1;
+ a1 = a0;
+ a0 = 0;
+ bx -= GMP_NUMB_BITS;
+ if (a2 == 0)
+ {
+ a2 = a1;
+ a1 = 0;
+ bx -= GMP_NUMB_BITS;
+ }
+ }
+
+ /* now a2 != 0 */
+ MPFR_ASSERTD(a2 != 0);
+ count_leading_zeros (cnt, a2);
+ if (cnt)
+ {
+ ap[2] = (a2 << cnt) | (a1 >> (GMP_NUMB_BITS - cnt));
+ ap[1] = (a1 << cnt) | (a0 >> (GMP_NUMB_BITS - cnt));
+ ap[0] = a0 << cnt;
+ bx -= cnt;
+ }
+ else
+ {
+ ap[2] = a2;
+ ap[1] = a1;
+ ap[0] = a0;
+ }
+ rb = sb = 0;
+ /* Note: sh is not initialized, but will not be used in this case. */
+ }
+ else if (bx > cx)
+ {
+ MPFR_SET_SAME_SIGN (a, b);
+ BGreater2:
+ d = (mpfr_uexp_t) bx - cx;
+ sh = 3 * GMP_NUMB_BITS - p;
+ mask = MPFR_LIMB_MASK(sh);
+ if (d < GMP_NUMB_BITS)
+ {
+ mp_limb_t cy;
+ /* warning: we must have the most significant bit of sb correct
+ since it might become the round bit below */
+ sb = cp[0] << (GMP_NUMB_BITS - d); /* neglected part of c */
+ a0 = bp[0] - ((cp[1] << (GMP_NUMB_BITS - d)) | (cp[0] >> d));
+ a1 = bp[1] - ((cp[2] << (GMP_NUMB_BITS - d)) | (cp[1] >> d))
+ - (a0 > bp[0]);
+ cy = a1 > bp[1] || (a1 == bp[1] && a0 > bp[0]); /* borrow in a1 */
+ a2 = bp[2] - (cp[2] >> d) - cy;
+ /* if sb is non-zero, subtract 1 from a2, a1, a0 since we want a
+ non-negative neglected part */
+ if (sb)
+ {
+ a2 -= (a1 == 0 && a0 == 0);
+ a1 -= (a0 == 0);
+ a0 --;
+ /* a = a2,a1,a0 cannot become zero here, since:
+ a) if d >= 2, then a2 >= 2^(w-1) - (2^(w-2)-1) with
+ w = GMP_NUMB_BITS, thus a2 - 1 >= 2^(w-2),
+ b) if d = 1, then since p < 3*GMP_NUMB_BITS we have sb=0. */
+ MPFR_ASSERTD(a2 > 0 || a1 > 0 || a0 > 0);
+ sb = -sb; /* 2^GMP_NUMB_BITS - sb */
+ }
+ if (a2 == 0)
+ {
+ /* this implies d=1, which in turn implies sb=0 */
+ MPFR_ASSERTD(sb == 0);
+ a2 = a1;
+ a1 = a0;
+ a0 = 0; /* should be a0 = sb */
+ /* since sb=0 already, no need to set it to 0 */
+ bx -= GMP_NUMB_BITS;
+ if (a2 == 0)
+ {
+ a2 = a1;
+ a1 = 0; /* should be a1 = a0 */
+ bx -= GMP_NUMB_BITS;
+ }
+ }
+ /* now a1 != 0 */
+ MPFR_ASSERTD(a2 != 0);
+ count_leading_zeros (cnt, a2);
+ if (cnt)
+ {
+ ap[2] = (a2 << cnt) | (a1 >> (GMP_NUMB_BITS - cnt));
+ ap[1] = (a1 << cnt) | (a0 >> (GMP_NUMB_BITS - cnt));
+ a0 = (a0 << cnt) | (sb >> (GMP_NUMB_BITS - cnt));
+ sb <<= cnt;
+ bx -= cnt;
+ }
+ else
+ {
+ ap[2] = a2;
+ ap[1] = a1;
+ }
+ /* sh > 0 since p < 2*GMP_NUMB_BITS */
+ MPFR_ASSERTD(sh > 0);
+ rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
+ sb |= (a0 & mask) ^ rb;
+ ap[0] = a0 & ~mask;
+ }
+ else if (d < 2 * GMP_NUMB_BITS)
+ {
+ mp_limb_t c0shifted;
+ /* warning: we must have the most significant bit of sb correct
+ since it might become the round bit below */
+ sb = (d == GMP_NUMB_BITS) ? cp[0]
+ : (cp[1] << (2*GMP_NUMB_BITS - d)) | (cp[0] != 0);
+ c0shifted = (d == GMP_NUMB_BITS) ? cp[1]
+ : (cp[2] << (2*GMP_NUMB_BITS-d)) | (cp[1] >> (d - GMP_NUMB_BITS));
+ a0 = bp[0] - c0shifted;
+ /* TODO: add a non-regression test for cp[2] == MPFR_LIMB_MAX,
+ d == GMP_NUMB_BITS and a0 > bp[0]. */
+ a1 = bp[1] - (cp[2] >> (d - GMP_NUMB_BITS)) - (a0 > bp[0]);
+ a2 = bp[2] - (a1 > bp[1] || (a1 == bp[1] && a0 > bp[0]));
+ /* if sb is non-zero, subtract 1 from a2, a1, a0 since we want a
+ non-negative neglected part */
+ if (sb)
+ {
+ a2 -= (a1 == 0 && a0 == 0);
+ a1 -= (a0 == 0);
+ a0 --;
+ /* a = a2,a1,a0 cannot become zero here, since:
+ a) if d >= 2, then a2 >= 2^(w-1) - (2^(w-2)-1) with
+ w = GMP_NUMB_BITS, thus a2 - 1 >= 2^(w-2),
+ b) if d = 1, then since p < 3*GMP_NUMB_BITS we have sb=0. */
+ MPFR_ASSERTD(a2 > 0 || a1 > 0 || a0 > 0);
+ sb = -sb; /* 2^GMP_NUMB_BITS - sb */
+ }
+ /* since bp[2] has its most significant bit set, we can have an
+ exponent decrease of at most one */
+ if (a2 < MPFR_LIMB_HIGHBIT)
+ {
+ ap[2] = (a2 << 1) | (a1 >> (GMP_NUMB_BITS - 1));
+ ap[1] = (a1 << 1) | (a0 >> (GMP_NUMB_BITS - 1));
+ a0 = (a0 << 1) | (sb >> (GMP_NUMB_BITS - 1));
+ sb <<= 1;
+ bx --;
+ }
+ else
+ {
+ ap[2] = a2;
+ ap[1] = a1;
+ }
+ rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
+ sb |= (a0 & mask) ^ rb;
+ ap[0] = a0 & ~mask;
+ }
+ else if (d < 3 * GMP_NUMB_BITS) /* 2*GMP_NUMB_BITS<=d<3*GMP_NUMB_BITS */
+ {
+ MPFR_ASSERTD (2*GMP_NUMB_BITS <= d && d < 3*GMP_NUMB_BITS);
+ /* warning: we must have the most significant bit of sb correct
+ since it might become the round bit below */
+ if (d == 2 * GMP_NUMB_BITS)
+ sb = cp[1] | (cp[0] != 0);
+ else
+ sb = cp[2] << (3*GMP_NUMB_BITS - d) | (cp[1] != 0) | (cp[0] != 0);
+ sb = -sb;
+ /* TODO: add a non-regression test for cp[2] == MPFR_LIMB_MAX,
+ d == 2*GMP_NUMB_BITS and sb != 0. */
+ a0 = bp[0] - (cp[2] >> (d - 2*GMP_NUMB_BITS)) - (sb != 0);
+ a1 = bp[1] - (a0 > bp[0] || (a0 == bp[0] && sb != 0));
+ a2 = bp[2] - (a1 > bp[1]);
+ if (a2 < MPFR_LIMB_HIGHBIT)
+ {
+ ap[2] = (a2 << 1) | (a1 >> (GMP_NUMB_BITS - 1));
+ ap[1] = (a1 << 1) | (a0 >> (GMP_NUMB_BITS - 1));
+ a0 = (a0 << 1) | (sb >> (GMP_NUMB_BITS - 1));
+ sb <<= 1;
+ bx --;
+ }
+ else
+ {
+ ap[2] = a2;
+ ap[1] = a1;
+ }
+ rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
+ sb |= (a0 & mask) ^ rb;
+ ap[0] = a0 & ~mask;
+ }
+ else /* d >= 3*GMP_NUMB_BITS */
+ {
+ /* We compute b - ulp(b), and the remainder ulp(b) - c satisfies:
+ 1/2 ulp(b) < ulp(b) - c < ulp(b), thus rb = sb = 1. */
+ mp_limb_t t = MPFR_LIMB_ONE << sh;
+ a0 = bp[0] - t;
+ a1 = bp[1] - (bp[0] < t);
+ a2 = bp[2] - (a1 > bp[1]);
+ if (a2 < MPFR_LIMB_HIGHBIT)
+ {
+ /* necessarily we had b = 1000...000 */
+ /* Warning: since we have an exponent decrease, when
+ p = 3*GMP_NUMB_BITS - 1 and d = 3*GMP_NUMB_BITS, the round bit
+ corresponds to the upper bit of -c. In that case rb = 0 and
+ sb = 1, except when c = 1000...000 where rb = 1 and sb = 0. */
+ rb = sh > 1 || d > 3 * GMP_NUMB_BITS
+ || (cp[2] == MPFR_LIMB_HIGHBIT && cp[1] == MPFR_LIMB_ZERO &&
+ cp[0] == MPFR_LIMB_ZERO);
+ /* sb=1 below is incorrect when p = 2*GMP_NUMB_BITS - 1,
+ d = 2*GMP_NUMB_BITS and c = 1000...000, but in
+ that case the even rule wound round up too. */
+ ap[0] = ~mask;
+ ap[1] = MPFR_LIMB_MAX;
+ ap[2] = MPFR_LIMB_MAX;
+ bx --;
+ }
+ else
+ {
+ ap[0] = a0;
+ ap[1] = a1;
+ ap[2] = a2;
+ rb = 1;
+ }
+ sb = 1;
+ }
+ }
+ else /* cx > bx */
+ {
+ mpfr_exp_t tx;
+ mp_limb_t *tp;
+ tx = bx; bx = cx; cx = tx;
+ tp = bp; bp = cp; cp = tp;
+ MPFR_SET_OPPOSITE_SIGN (a, b);
+ goto BGreater2;
+ }
+
+ /* now perform rounding */
+
+ /* Warning: MPFR considers underflow *after* rounding with an unbounded
+ exponent range. However since b and c have same precision p, they are
+ multiples of 2^(emin-p), likewise for b-c. Thus if bx < emin, the
+ subtraction (with an unbounded exponent range) is exact, so that bx is
+ also the exponent after rounding with an unbounded exponent range. */
+ if (MPFR_UNLIKELY(bx < __gmpfr_emin))
+ {
+ /* for RNDN, mpfr_underflow always rounds away, thus for |a|<=2^(emin-2)
+ we have to change to RNDZ */
+ if (rnd_mode == MPFR_RNDN &&
+ (bx < __gmpfr_emin - 1 ||
+ (ap[2] == MPFR_LIMB_HIGHBIT && ap[1] == 0 && ap[0] == 0)))
+ rnd_mode = MPFR_RNDZ;
+ return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
+ }
+
+ MPFR_SET_EXP (a, bx);
+ if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF)
+ MPFR_RET (0);
+ else if (rnd_mode == MPFR_RNDN)
+ {
+ if (rb == 0 || (sb == 0 && (ap[0] & (MPFR_LIMB_ONE << sh)) == 0))
+ goto truncate;
+ else
+ goto add_one_ulp;
+ }
+ else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
+ {
+ truncate:
+ MPFR_RET(-MPFR_SIGN(a));
+ }
+ else /* round away from zero */
+ {
+ add_one_ulp:
+ ap[0] += MPFR_LIMB_ONE << sh;
+ ap[1] += (ap[0] == 0);
+ ap[2] += (ap[1] == 0 && ap[0] == 0);
+ if (MPFR_UNLIKELY(ap[2] == 0))
+ {
+ ap[2] = MPFR_LIMB_HIGHBIT;
+ /* Note: bx+1 cannot exceed __gmpfr_emax, since |a| <= |b|, thus
+ bx+1 is at most equal to the original exponent of b. */
+ MPFR_ASSERTD(bx + 1 <= __gmpfr_emax);
+ MPFR_SET_EXP (a, bx + 1);
+ }
+ MPFR_RET(MPFR_SIGN(a));
+ }
+}
+
+#endif /* !defined(MPFR_GENERIC_ABI) */
/* Rounding Sub */
@@ -131,10 +1089,10 @@ int mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
*
*/
-int
+MPFR_HOT_FUNCTION_ATTR int
mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
- mpfr_exp_t bx,cx;
+ mpfr_exp_t bx, cx;
mpfr_uexp_t d;
mpfr_prec_t p, sh, cnt;
mp_size_t n;
@@ -142,34 +1100,54 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
mp_limb_t limb;
int inexact;
mp_limb_t bcp,bcp1; /* Cp and C'p+1 */
- mp_limb_t bbcp = (mp_limb_t) -1, bbcp1 = (mp_limb_t) -1; /* Cp+1 and C'p+2,
+ mp_limb_t bbcp = MPFR_LIMB_MAX, bbcp1 = MPFR_LIMB_MAX; /* Cp+1 and C'p+2,
gcc claims that they might be used uninitialized. We fill them with invalid
values, which should produce a failure if so. See README.dev file. */
MPFR_TMP_DECL(marker);
- MPFR_TMP_MARK(marker);
-
MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c));
MPFR_ASSERTD(MPFR_IS_PURE_FP(b));
MPFR_ASSERTD(MPFR_IS_PURE_FP(c));
/* Read prec and num of limbs */
- p = MPFR_PREC (b);
- n = MPFR_PREC2LIMBS (p);
+ p = MPFR_GET_PREC (b);
+
+#if !defined(MPFR_GENERIC_ABI)
+ /* special case for p < GMP_NUMB_BITS */
+ if (p < GMP_NUMB_BITS)
+ return mpfr_sub1sp1 (a, b, c, rnd_mode, p);
+
+ /* special case for GMP_NUMB_BITS < p < 2*GMP_NUMB_BITS */
+ if (GMP_NUMB_BITS < p && p < 2 * GMP_NUMB_BITS)
+ return mpfr_sub1sp2 (a, b, c, rnd_mode, p);
+
+ /* special case for p = GMP_NUMB_BITS: we put it *after* mpfr_sub1sp2,
+ in order not to slow down mpfr_sub1sp2, which should be more frequent */
+ if (p == GMP_NUMB_BITS)
+ return mpfr_sub1sp1n (a, b, c, rnd_mode);
- /* Fast cmp of |b| and |c|*/
+ /* special case for 2*GMP_NUMB_BITS < p < 3*GMP_NUMB_BITS */
+ if (2 * GMP_NUMB_BITS < p && p < 3 * GMP_NUMB_BITS)
+ return mpfr_sub1sp3 (a, b, c, rnd_mode, p);
+#endif
+
+ n = MPFR_PREC2LIMBS (p);
+ /* Fast cmp of |b| and |c| */
bx = MPFR_GET_EXP (b);
cx = MPFR_GET_EXP (c);
- if (MPFR_UNLIKELY(bx == cx))
+
+ MPFR_TMP_MARK(marker);
+
+ if (bx == cx)
{
mp_size_t k = n - 1;
- /* Check mantissa since exponent are equals */
+ /* Check mantissa since exponents are equal */
bp = MPFR_MANT(b);
cp = MPFR_MANT(c);
- while (k>=0 && MPFR_UNLIKELY(bp[k] == cp[k]))
+ while (k >= 0 && MPFR_UNLIKELY(bp[k] == cp[k]))
k--;
- if (MPFR_UNLIKELY(k < 0))
+ if (k < 0)
/* b == c ! */
{
/* Return exact number 0 */
@@ -184,11 +1162,11 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
goto BGreater;
else
{
- MPFR_ASSERTD(bp[k]<cp[k]);
+ MPFR_ASSERTD(bp[k] < cp[k]);
goto CGreater;
}
}
- else if (MPFR_UNLIKELY(bx < cx))
+ else if (bx < cx)
{
/* Swap b and c and set sign */
mpfr_srcptr t;
@@ -200,19 +1178,19 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
}
else
{
- /* b > c */
+ /* |b| > |c| */
BGreater:
MPFR_SET_SAME_SIGN(a,b);
}
- /* Now b > c */
+ /* Now |b| > |c| */
MPFR_ASSERTD(bx >= cx);
d = (mpfr_uexp_t) bx - cx;
- DEBUG (printf ("New with diff=%lu\n", (unsigned long) d));
+ /* printf ("New with diff=%lu\n", (unsigned long) d); */
- if (MPFR_UNLIKELY(d <= 1))
+ if (d <= 1)
{
- if (MPFR_LIKELY(d < 1))
+ if (d == 0)
{
/* <-- b -->
<-- c --> : exact sub */
@@ -221,7 +1199,7 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
/* Normalize */
ExactNormalize:
limb = ap[n-1];
- if (MPFR_LIKELY(limb))
+ if (MPFR_LIKELY (limb != 0))
{
/* First limb is not zero. */
count_leading_zeros(cnt, limb);
@@ -231,7 +1209,7 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
mpn_lshift(ap, ap, n, cnt); /* Normalize number */
bx -= cnt; /* Update final expo */
}
- /* Last limb should be ok */
+ /* Last limb should be OK */
MPFR_ASSERTD(!(ap[0] & MPFR_LIMB_MASK((unsigned int) (-p)
% GMP_NUMB_BITS)));
}
@@ -239,8 +1217,8 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
/* First limb is zero */
mp_size_t k = n-1, len;
- /* Find the first limb not equal to zero.
- FIXME:It is assume it exists (since |b| > |c| and same prec)*/
+ /* Find the first limb not equal to zero. It necessarily exists
+ since |b| > |c|. */
do
{
MPFR_ASSERTD( k > 0 );
@@ -252,28 +1230,30 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
k++;
len = n - k; /* Number of last limb */
MPFR_ASSERTD(k >= 0);
- if (MPFR_LIKELY(cnt))
- mpn_lshift(ap+len, ap, k, cnt); /* Normalize the High Limb*/
+ if (cnt)
+ mpn_lshift (ap + len, ap, k, cnt); /* Normalize the High Limb*/
else
{
- /* Must use DECR since src and dest may overlap & dest>=src*/
- MPN_COPY_DECR(ap+len, ap, k);
+ /* Must use copyd since src and dst may overlap & dst>=src */
+ mpn_copyd (ap+len, ap, k);
}
MPN_ZERO(ap, len); /* Zeroing the last limbs */
bx -= cnt + len*GMP_NUMB_BITS; /* Update Expo */
- /* Last limb should be ok */
- MPFR_ASSERTD(!(ap[len]&MPFR_LIMB_MASK((unsigned int) (-p)
- % GMP_NUMB_BITS)));
+ /* Last limb should be OK */
+ MPFR_ASSERTD(!(ap[len] & MPFR_LIMB_MASK((unsigned int) (-p)
+ % GMP_NUMB_BITS)));
}
/* Check expo underflow */
if (MPFR_UNLIKELY(bx < __gmpfr_emin))
{
MPFR_TMP_FREE(marker);
- /* inexact=0 */
- DEBUG( printf("(D==0 Underflow)\n") );
+ /* since b and c have same sign, exponent and precision, the
+ subtraction is exact */
+ /* printf("(D==0 Underflow)\n"); */
+ /* for MPFR_RNDN, mpfr_underflow always rounds away from zero,
+ thus for |a| <= 2^(emin-2) we change to RNDZ. */
if (rnd_mode == MPFR_RNDN &&
- (bx < __gmpfr_emin - 1 ||
- (/*inexact >= 0 &&*/ mpfr_powerof2_raw (a))))
+ (bx < __gmpfr_emin - 1 || mpfr_powerof2_raw (a)))
rnd_mode = MPFR_RNDZ;
return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
}
@@ -296,41 +1276,48 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
cp = MPFR_MANT(c);
k = n-1;
limb = bp[k] - cp[k]/2;
+ /* we have |b|-|c| >= limb*W^k - (2*W^k-1)/2 >= limb*W^k - W^k + 1/2
+ thus if limb > W^k/2, |b|-|c| >= 1/2*W^n.
+ Moreover if trunc(|c|) represents the first p-1 bits of |c|,
+ minus the last significant bit called c0 below, then we have
+ |b|-trunc(|c|) >= 1/2*W^n+1, thus the two mpn_sub_n calls
+ below necessarily yield a > 1/2*W^n. */
if (limb > MPFR_LIMB_HIGHBIT)
{
- /* We can't lose precision: compute b-c/2 */
+ /* The exponent cannot decrease: compute b-c/2 */
/* Shift c in the allocated temporary block */
SubD1NoLose:
- c0 = cp[0] & (MPFR_LIMB_ONE<<sh);
+ c0 = cp[0] & (MPFR_LIMB_ONE << sh);
+ mask = ~MPFR_LIMB_MASK(sh);
+ ap = MPFR_MANT(a);
cp = MPFR_TMP_LIMBS_ALLOC (n);
- mpn_rshift(cp, MPFR_MANT(c), n, 1);
+ /* FIXME: it might be faster to have one function shifting c by 1
+ to the right and adding with b to a, which would read c once
+ only, and avoid a temporary allocation. */
+ mpn_rshift (cp, MPFR_MANT(c), n, 1);
+ cp[0] &= mask; /* Zero last bit of c if set */
+ mpn_sub_n (ap, bp, cp, n);
+ MPFR_SET_EXP(a, bx); /* No expo overflow! */
+ MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
if (MPFR_LIKELY(c0 == 0))
{
/* Result is exact: no need of rounding! */
- ap = MPFR_MANT(a);
- mpn_sub_n (ap, bp, cp, n);
- MPFR_SET_EXP(a, bx); /* No expo overflow! */
- /* No truncate or normalize is needed */
- MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
- /* No rounding is necessary since the result is exact */
MPFR_TMP_FREE(marker);
return 0;
}
- ap = MPFR_MANT(a);
- mask = ~MPFR_LIMB_MASK(sh);
- cp[0] &= mask; /* Delete last bit of c */
- mpn_sub_n (ap, bp, cp, n);
- MPFR_SET_EXP(a, bx); /* No expo overflow! */
MPFR_ASSERTD( !(ap[0] & ~mask) ); /* Check last bits */
/* No normalize is needed */
- MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
- /* Rounding is necessary since c0 = 1*/
+ /* Rounding is necessary since c0 = 1 */
/* Cp =-1 and C'p+1=0 */
bcp = 1; bcp1 = 0;
- if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
+
+ if (rnd_mode == MPFR_RNDF)
+ goto truncate; /* low(b) = 0 and low(c) is 0 or 1/2 ulp(b), thus
+ low(b) - low(c) = 0 or -1/2 ulp(b) */
+ else if (rnd_mode == MPFR_RNDN)
{
- /* Even Rule apply: Check Ap-1 */
- if (MPFR_LIKELY( (ap[0] & (MPFR_LIMB_ONE<<sh)) == 0) )
+ /* Even Rule apply: Check last bit of a. */
+ if (MPFR_LIKELY( (ap[0] & (MPFR_LIMB_ONE << sh)) == 0) )
goto truncate;
else
goto sub_one_ulp;
@@ -343,14 +1330,23 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
}
else if (MPFR_LIKELY(limb < MPFR_LIMB_HIGHBIT))
{
- /* We lose at least one bit of prec */
- /* Calcul of 2*b-c (Exact) */
- /* Shift b in the allocated temporary block */
+ /* |b| - |c| <= (W/2-1)*W^k + W^k-1 = 1/2*W^n - 1 */
+ /* The exponent decreases by one. */
SubD1Lose:
+ ap = MPFR_MANT(a);
+ /* Compute 2*b-c (Exact) */
+ /* Experiments with __gmpn_rsblsh_n show that it is not always
+ faster than mpn_lshift + mpn_sub_n, thus we don't enable it
+ for now (HAVE___GMPN_RSBLSH_N -> HAVE___GMPN_RSBLSH_Nxxx). */
+#if defined(WANT_GMP_INTERNALS) && defined(HAVE___GMPN_RSBLSH_Nxxx)
+ /* {ap, n} = 2*{bp, n} - {cp, n} */
+ __gmpn_rsblsh_n (ap, MPFR_MANT(c), MPFR_MANT(b), n, 1);
+#else
bp = MPFR_TMP_LIMBS_ALLOC (n);
+ /* Shift b in the allocated temporary block */
mpn_lshift (bp, MPFR_MANT(b), n, 1);
- ap = MPFR_MANT(a);
mpn_sub_n (ap, bp, cp, n);
+#endif
bx--;
goto ExactNormalize;
}
@@ -363,42 +1359,55 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
/* If k==-1 => We don't lose any bit
AND the result is 100000000000 0000000000 00000000000 */
mp_limb_t carry;
- do {
- carry = cp[k]&MPFR_LIMB_ONE;
- k--;
- } while (k>=0 &&
- bp[k]==(carry=cp[k]/2+(carry<<(GMP_NUMB_BITS-1))));
- if (MPFR_UNLIKELY(k<0))
+ do
+ {
+ carry = cp[k] << (GMP_NUMB_BITS - 1);
+ if (--k < 0)
+ break;
+ carry += cp[k] >> 1;
+ }
+ while (bp[k] == carry);
+ if (MPFR_UNLIKELY(k < 0))
{
- /*If carry then (sh==0 and Virtual c'[-1] > Virtual b[-1]) */
- if (MPFR_UNLIKELY(carry)) /* carry = cp[0]&MPFR_LIMB_ONE */
+ ap = MPFR_MANT (a);
+ if (MPFR_UNLIKELY(carry))
{
- /* FIXME: Can be faster? */
- MPFR_ASSERTD(sh == 0);
- goto SubD1Lose;
+ /* If carry then necessarily the precision is an exact
+ multiple of GMP_NUMB_BITS, and we lose one bit,
+ thus the (exact) result is a power of 2 minus 1. */
+ memset (ap, -1, n * MPFR_BYTES_PER_MP_LIMB);
+ MPFR_SET_EXP (a, bx - 1);
+ /* No underflow is possible since cx = bx-1 is a valid
+ exponent. */
}
- /* Result is a power of 2 */
- ap = MPFR_MANT (a);
- MPN_ZERO (ap, n);
- ap[n-1] = MPFR_LIMB_HIGHBIT;
- MPFR_SET_EXP (a, bx); /* No expo overflow! */
- /* No Normalize is needed*/
+ else
+ {
+ /* No carry: result is a power of 2. */
+ MPN_ZERO (ap, n - 1);
+ ap[n-1] = MPFR_LIMB_HIGHBIT;
+ MPFR_SET_EXP (a, bx); /* No expo overflow! */
+ }
+ /* No Normalize is needed */
/* No Rounding is needed */
MPFR_TMP_FREE (marker);
return 0;
}
/* carry = cp[k]/2+(cp[k-1]&1)<<(GMP_NUMB_BITS-1) = c'[k]*/
else if (bp[k] > carry)
- goto SubD1NoLose;
+ goto SubD1NoLose; /* |b|-|c| >= 1/2*W^n */
else
{
- MPFR_ASSERTD(bp[k]<carry);
- goto SubD1Lose;
+ MPFR_ASSERTD(bp[k] < carry);
+ goto SubD1Lose; /* |b|-|c| <= 1/2*W^n-1 and is exact */
}
}
}
}
- else if (MPFR_UNLIKELY(d >= p))
+ else if (MPFR_UNLIKELY(d >= p)) /* the difference of exponents is larger
+ than the precision of all operands, thus
+ the result is either b or b - 1 ulp,
+ with a possible exact result when
+ d = p, b = 2^e and c = 1/2 ulp(b) */
{
ap = MPFR_MANT(a);
MPFR_UNSIGNED_MINUS_MODULO(sh, p);
@@ -419,24 +1428,31 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
{
mp_size_t k = n-1;
- do {
+ do
k--;
- } while (k>=0 && cp[k]==0);
- bcp1 = (k>=0);
+ while (k >= 0 && cp[k] == 0);
+ bcp1 = (k >= 0);
}
else
bcp1 = 1;
}
- DEBUG( printf("(D=P) Cp=-1 Cp+1=%d C'p+1=%d \n", bbcp!=0, bcp1!=0) );
+ /* printf("(D=P) Cp=-1 Cp+1=%d C'p+1=%d \n", bbcp!=0, bcp1!=0); */
bp = MPFR_MANT (b);
- /* Even if src and dest overlap, it is ok using MPN_COPY */
- if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
+ /* Even if src and dest overlap, it is OK using MPN_COPY */
+ if (MPFR_LIKELY(rnd_mode == MPFR_RNDF))
+ /* then d = p, and subtracting one ulp of b is ok even in the
+ exact case b = 2^e and c = 1/2 ulp(b) */
{
- if (MPFR_UNLIKELY( bcp && bcp1==0 ))
+ MPN_COPY(ap, bp, n);
+ goto sub_one_ulp;
+ }
+ else if (rnd_mode == MPFR_RNDN)
+ {
+ if (MPFR_UNLIKELY (bcp != 0 && bcp1 == 0))
/* Cp=-1 and C'p+1=0: Even rule Apply! */
/* Check Ap-1 = Bp-1 */
- if ((bp[0] & (MPFR_LIMB_ONE<<sh)) == 0)
+ if ((bp[0] & (MPFR_LIMB_ONE << sh)) == 0)
{
MPN_COPY(ap, bp, n);
goto truncate;
@@ -460,7 +1476,7 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
/* Cp=0, Cp+1=-1 if d==p+1, C'p+1=-1 */
bcp = 0; bbcp = (d==p+1); bcp1 = 1;
- DEBUG( printf("(D>P) Cp=%d Cp+1=%d C'p+1=%d\n", bcp!=0,bbcp!=0,bcp1!=0) );
+ /* printf("(D>P) Cp=%d Cp+1=%d C'p+1=%d\n", bcp!=0,bbcp!=0,bcp1!=0); */
/* Need to compute C'p+2 if d==p+1 and if rnd_mode=NEAREST
(Because of a very improbable case) */
if (MPFR_UNLIKELY(d==p+1 && rnd_mode==MPFR_RNDN))
@@ -469,19 +1485,19 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
{
mp_size_t k = n-1;
- do {
+ do
k--;
- } while (k>=0 && cp[k]==0);
- bbcp1 = (k>=0);
+ while (k >= 0 && cp[k] == 0);
+ bbcp1 = (k >= 0);
}
else
bbcp1 = 1;
- DEBUG( printf("(D>P) C'p+2=%d\n", bbcp1!=0) );
+ /* printf("(D>P) C'p+2=%d\n", bbcp1!=0); */
}
/* Copy mantissa B in A */
MPN_COPY(ap, MPFR_MANT(b), n);
/* Round */
- if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
+ if (rnd_mode == MPFR_RNDF || rnd_mode == MPFR_RNDN)
goto truncate;
MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
if (rnd_mode == MPFR_RNDZ)
@@ -490,13 +1506,12 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
goto truncate;
}
}
- else
+ else /* case 2 <= d < p */
{
mpfr_uexp_t dm;
mp_size_t m;
mp_limb_t mask;
- /* General case: 2 <= d < p */
MPFR_UNSIGNED_MINUS_MODULO(sh, p);
cp = MPFR_TMP_LIMBS_ALLOC (n);
@@ -506,7 +1521,7 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
if (MPFR_UNLIKELY(dm == 0))
{
/* dm = 0 and m > 0: Just copy */
- MPFR_ASSERTD(m!=0);
+ MPFR_ASSERTD(m != 0);
MPN_COPY(cp, MPFR_MANT(c)+m, n-m);
MPN_ZERO(cp+n-m, m);
}
@@ -523,16 +1538,16 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
MPN_ZERO(cp+n-m, m);
}
- DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c), p) );
- DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) );
- DEBUG( mpfr_print_mant_binary("After ", cp, p) );
+ /* mpfr_print_mant_binary("Before", MPFR_MANT(c), p); */
+ /* mpfr_print_mant_binary("B= ", MPFR_MANT(b), p); */
+ /* mpfr_print_mant_binary("After ", cp, p); */
/* Compute bcp=Cp and bcp1=C'p+1 */
if (MPFR_LIKELY(sh))
{
/* Try to compute them from C' rather than C (FIXME: Faster?) */
bcp = (cp[0] & (MPFR_LIMB_ONE<<(sh-1))) ;
- if (MPFR_LIKELY( cp[0] & MPFR_LIMB_MASK(sh-1) ))
+ if (cp[0] & MPFR_LIMB_MASK(sh-1))
bcp1 = 1;
else
{
@@ -540,7 +1555,7 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
/* Start from bit x=p-d+sh in mantissa C
(+sh since we have already looked sh bits in C'!) */
mpfr_prec_t x = p-d+sh-1;
- if (MPFR_LIKELY(x>p))
+ if (x > p)
/* We are already looked at all the bits of c, so C'p+1 = 0*/
bcp1 = 0;
else
@@ -548,9 +1563,8 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
mp_limb_t *tp = MPFR_MANT(c);
mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
- DEBUG (printf ("(First) x=%lu Kx=%ld Sx=%lu\n",
- (unsigned long) x, (long) kx,
- (unsigned long) sx));
+ /* printf ("(First) x=%lu Kx=%ld Sx=%lu\n",
+ (unsigned long) x, (long) kx, (unsigned long) sx); */
/* Looks at the last bits of limb kx (if sx=0 does nothing)*/
if (tp[kx] & MPFR_LIMB_MASK(sx))
bcp1 = 1;
@@ -558,9 +1572,9 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
/*kx += (sx==0);*/
/*If sx==0, tp[kx] hasn't been checked*/
- do {
+ do
kx--;
- } while (kx>=0 && tp[kx]==0);
+ while (kx >= 0 && tp[kx] == 0);
bcp1 = (kx >= 0);
}
}
@@ -582,17 +1596,17 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
else
{
/*kx += (sx==0);*/ /*If sx==0, tp[kx] hasn't been checked*/
- do {
+ do
kx--;
- } while (kx>=0 && tp[kx]==0);
- bcp1 = (kx>=0);
+ while (kx >= 0 && tp[kx] == 0);
+ bcp1 = (kx >= 0);
}
}
- DEBUG( printf("sh=%lu Cp=%d C'p+1=%d\n", sh, bcp!=0, bcp1!=0) );
+ /* printf("sh=%lu Cp=%d C'p+1=%d\n", sh, bcp!=0, bcp1!=0); */
/* Check if we can lose a bit, and if so compute Cp+1 and C'p+2 */
bp = MPFR_MANT(b);
- if (MPFR_UNLIKELY((bp[n-1]-cp[n-1]) <= MPFR_LIMB_HIGHBIT))
+ if (MPFR_UNLIKELY (bp[n-1] - cp[n-1] <= MPFR_LIMB_HIGHBIT))
{
/* We can lose a bit so we precompute Cp+1 and C'p+2 */
/* Test for trivial case: since C'p+1=0, Cp+1=0 and C'p+2 =0 */
@@ -608,28 +1622,28 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
mp_limb_t *tp = MPFR_MANT(c);
/* Start from bit x=(p+1)-d in mantissa C */
mpfr_prec_t x = p+1-d;
- mp_size_t kx = n-1 - (x/GMP_NUMB_BITS);
- mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
+ mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
+ mpfr_prec_t sx = GMP_NUMB_BITS-1 - (x % GMP_NUMB_BITS);
+
MPFR_ASSERTD(p > d);
- DEBUG (printf ("(pre) x=%lu Kx=%ld Sx=%lu\n",
- (unsigned long) x, (long) kx,
- (unsigned long) sx));
+ /* printf ("(pre) x=%lu Kx=%ld Sx=%lu\n",
+ (unsigned long) x, (long) kx, (unsigned long) sx); */
bbcp = (tp[kx] & (MPFR_LIMB_ONE<<sx)) ;
/* Looks at the last bits of limb kx (If sx=0, does nothing)*/
/* If Cp+1=0, since C'p+1!=0, C'p+2=1 ! */
- if (MPFR_LIKELY(bbcp==0 || (tp[kx]&MPFR_LIMB_MASK(sx))))
+ if (MPFR_LIKELY (bbcp == 0 || (tp[kx] & MPFR_LIMB_MASK(sx))))
bbcp1 = 1;
else
{
/*kx += (sx==0);*/ /*If sx==0, tp[kx] hasn't been checked*/
- do {
+ do
kx--;
- } while (kx>=0 && tp[kx]==0);
- bbcp1 = (kx>=0);
- DEBUG (printf ("(Pre) Scan done for %ld\n", (long) kx));
+ while (kx >= 0 && tp[kx] == 0);
+ bbcp1 = (kx >= 0);
+ /* printf ("(Pre) Scan done for %ld\n", (long) kx); */
}
} /*End of Bcp1 != 0*/
- DEBUG( printf("(Pre) Cp+1=%d C'p+2=%d\n", bbcp!=0, bbcp1!=0) );
+ /* printf("(Pre) Cp+1=%d C'p+2=%d\n", bbcp!=0, bbcp1!=0); */
} /* End of "can lose a bit" */
/* Clean shifted C' */
@@ -639,7 +1653,7 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
/* Subtract the mantissa c from b in a */
ap = MPFR_MANT(a);
mpn_sub_n (ap, bp, cp, n);
- DEBUG( mpfr_print_mant_binary("Sub= ", ap, p) );
+ /* mpfr_print_mant_binary("Sub= ", ap, p); */
/* Normalize: we lose at max one bit*/
if (MPFR_UNLIKELY(MPFR_LIMB_MSB(ap[n-1]) == 0))
@@ -648,8 +1662,8 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
/* Ap >= 010000xxx001 */
mpn_lshift(ap, ap, n, 1);
/* Ap >= 100000xxx010 */
- if (MPFR_UNLIKELY(bcp!=0)) /* Check if Cp = -1 */
- /* Since Cp == -1, we have to substract one more */
+ if (MPFR_UNLIKELY(bcp != 0)) /* Check if Cp = -1 */
+ /* Since Cp == -1, we have to subtract one more */
{
mpn_sub_1(ap, ap, n, MPFR_LIMB_ONE<<sh);
MPFR_ASSERTD(MPFR_LIMB_MSB(ap[n-1]) != 0);
@@ -658,21 +1672,23 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
/* Final exponent -1 since we have shifted the mantissa */
bx--;
/* Update bcp and bcp1 */
- MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
- MPFR_ASSERTN(bbcp1 != (mp_limb_t) -1);
+ MPFR_ASSERTD(bbcp != MPFR_LIMB_MAX);
+ MPFR_ASSERTD(bbcp1 != MPFR_LIMB_MAX);
bcp = bbcp;
bcp1 = bbcp1;
- /* We dont't have anymore a valid Cp+1!
+ /* We don't have anymore a valid Cp+1!
But since Ap >= 100000xxx001, the final sub can't unnormalize!*/
}
MPFR_ASSERTD( !(ap[0] & ~mask) );
/* Rounding */
- if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
+ if (MPFR_LIKELY(rnd_mode == MPFR_RNDF))
+ goto truncate;
+ else if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
{
- if (MPFR_LIKELY(bcp==0))
+ if (MPFR_LIKELY(bcp == 0))
goto truncate;
- else if ((bcp1) || ((ap[0] & (MPFR_LIMB_ONE<<sh)) != 0))
+ else if (bcp1 != 0 || (ap[0] & (MPFR_LIMB_ONE << sh)) != 0)
goto sub_one_ulp;
else
goto truncate;
@@ -680,7 +1696,7 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
/* Update rounding mode */
MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
- if (rnd_mode == MPFR_RNDZ && (MPFR_LIKELY(bcp || bcp1)))
+ if (rnd_mode == MPFR_RNDZ && MPFR_LIKELY (bcp != 0 || bcp1 != 0))
goto sub_one_ulp;
goto truncate;
}
@@ -691,37 +1707,40 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
/* Result should be smaller than exact value: inexact=-1 */
inexact = -1;
- /* Check normalisation */
- if (MPFR_UNLIKELY(MPFR_LIMB_MSB(ap[n-1]) == 0))
+ /* Check normalization */
+ if (MPFR_UNLIKELY(ap[n-1] < MPFR_LIMB_HIGHBIT))
{
/* ap was a power of 2, and we lose a bit */
/* Now it is 0111111111111111111[00000 */
- mpn_lshift(ap, ap, n, 1);
+ /* The following 2 lines are equivalent to: mpn_lshift(ap, ap, n, 1); */
+ ap[0] <<= 1;
+ ap[n-1] |= MPFR_LIMB_HIGHBIT;
bx--;
/* And the lost bit x depends on Cp+1, and Cp */
- /* Compute Cp+1 if it isn't already compute (ie d==1) */
- /* FIXME: Is this case possible? */
- if (MPFR_UNLIKELY(d == 1))
- bbcp = 0;
- DEBUG( printf("(SubOneUlp)Cp=%d, Cp+1=%d C'p+1=%d\n", bcp!=0,bbcp!=0,bcp1!=0));
+ /* Compute Cp+1 if it isn't already computed (ie d==1) */
+ /* Note: we can't have d = 1 here, since the only "goto sub_one_ulp"
+ for d = 1 are in the "SubD1NoLose" case, and in that case
+ |b|-|c| >= 1/2*W^n, thus round(|b|-|c|) >= 1/2*W^n, and ap[n-1]
+ cannot go below MPFR_LIMB_HIGHBIT. */
+ /* printf("(SubOneUlp)Cp=%d, Cp+1=%d C'p+1=%d\n", bcp!=0,bbcp!=0,bcp1!=0); */
/* Compute the last bit (Since we have shifted the mantissa)
- we need one more bit!*/
- MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
- if ( (rnd_mode == MPFR_RNDZ && bcp==0)
- || (rnd_mode==MPFR_RNDN && bbcp==0)
- || (bcp && bcp1==0) ) /*Exact result*/
+ we need one more bit! */
+ MPFR_ASSERTD(bbcp != MPFR_LIMB_MAX);
+ if ((rnd_mode == MPFR_RNDZ && bcp == 0) ||
+ (rnd_mode == MPFR_RNDN && bbcp == 0) ||
+ (bcp != 0 && bcp1 == 0)) /* Exact result */
{
- ap[0] |= MPFR_LIMB_ONE<<sh;
+ ap[0] |= MPFR_LIMB_ONE << sh;
if (rnd_mode == MPFR_RNDN)
inexact = 1;
- DEBUG( printf("(SubOneUlp) Last bit set\n") );
+ /* printf("(SubOneUlp) Last bit set\n"); */
}
/* Result could be exact if C'p+1 = 0 and rnd == Zero
since we have had one more bit to the result */
/* Fixme: rnd_mode == MPFR_RNDZ needed ? */
- if (bcp1==0 && rnd_mode==MPFR_RNDZ)
+ if (rnd_mode == MPFR_RNDZ && bcp1 == 0)
{
- DEBUG( printf("(SubOneUlp) Exact result\n") );
+ /* printf("(SubOneUlp) Exact result\n"); */
inexact = 0;
}
}
@@ -744,46 +1763,43 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
if (MPFR_UNLIKELY(ap[n-1] == MPFR_LIMB_HIGHBIT))
{
mp_size_t k = n-1;
- do {
+ do
k--;
- } while (k>=0 && ap[k]==0);
- if (MPFR_UNLIKELY(k<0))
+ while (k >= 0 && ap[k] == 0);
+ if (MPFR_UNLIKELY (k < 0))
{
/* It is a power of 2! */
/* Compute Cp+1 if it isn't already compute (ie d==1) */
- /* FIXME: Is this case possible? */
- if (d == 1)
- bbcp=0;
- DEBUG( printf("(Truncate) Cp=%d, Cp+1=%d C'p+1=%d C'p+2=%d\n", \
- bcp!=0, bbcp!=0, bcp1!=0, bbcp1!=0) );
- MPFR_ASSERTN(bbcp != (mp_limb_t) -1);
- MPFR_ASSERTN((rnd_mode != MPFR_RNDN) || (bcp != 0) || (bbcp == 0) || (bbcp1 != (mp_limb_t) -1));
- if (((rnd_mode != MPFR_RNDZ) && bcp)
- ||
- ((rnd_mode == MPFR_RNDN) && (bcp == 0) && (bbcp) && (bbcp1)))
- {
- DEBUG( printf("(Truncate) Do sub\n") );
+ /* Note: if d=1, we have {a, n} > 1/2*W^n, thus we cannot have k < 0. */
+ /* printf("(Truncate) Cp=%d, Cp+1=%d C'p+1=%d C'p+2=%d\n",
+ bcp!=0, bbcp!=0, bcp1!=0, bbcp1!=0); */
+ MPFR_ASSERTD(bbcp != MPFR_LIMB_MAX);
+ MPFR_ASSERTD(rnd_mode != MPFR_RNDN || bcp != 0 ||
+ bbcp == 0 || bbcp1 != MPFR_LIMB_MAX);
+ if ((rnd_mode != MPFR_RNDZ && bcp != 0) ||
+ (rnd_mode == MPFR_RNDN && bcp == 0 && bbcp != 0 && bbcp1 != 0))
+ {
+ /* printf("(Truncate) Do sub\n"); */
mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
- mpn_lshift(ap, ap, n, 1);
- ap[0] |= MPFR_LIMB_ONE<<sh;
+ ap[n-1] |= MPFR_LIMB_HIGHBIT;
bx--;
/* FIXME: Explain why it works (or why not)... */
- inexact = (bcp1 == 0) ? 0 : (rnd_mode==MPFR_RNDN) ? -1 : 1;
+ inexact = (bcp1 == 0) ? 0 : (rnd_mode == MPFR_RNDN) ? -1 : 1;
goto end_of_sub;
}
}
}
/* Calcul of Inexact flag.*/
- inexact = MPFR_LIKELY(bcp || bcp1) ? 1 : 0;
+ inexact = (bcp != 0 || bcp1 != 0);
end_of_sub:
- /* Update Expo */
- /* FIXME: Is this test really useful?
+ /* Update Exponent */
+ /* bx >= emin. Proof:
If d==0 : Exact case. This is never called.
if 1 < d < p : bx=MPFR_EXP(b) or MPFR_EXP(b)-1 > MPFR_EXP(c) > emin
if d == 1 : bx=MPFR_EXP(b). If we could lose any bits, the exact
- normalisation is called.
+ normalization is called.
if d >= p : bx=MPFR_EXP(b) >= MPFR_EXP(c) + p > emin
After SubOneUlp, we could have one bit less.
if 1 < d < p : bx >= MPFR_EXP(b)-2 >= MPFR_EXP(c) > emin
@@ -791,18 +1807,6 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
if d >= p : bx >= MPFR_EXP(b)-1 > emin since p>=2.
*/
MPFR_ASSERTD( bx >= __gmpfr_emin);
- /*
- if (MPFR_UNLIKELY(bx < __gmpfr_emin))
- {
- DEBUG( printf("(Final Underflow)\n") );
- if (rnd_mode == MPFR_RNDN &&
- (bx < __gmpfr_emin - 1 ||
- (inexact >= 0 && mpfr_powerof2_raw (a))))
- rnd_mode = MPFR_RNDZ;
- MPFR_TMP_FREE(marker);
- return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
- }
- */
MPFR_SET_EXP (a, bx);
MPFR_TMP_FREE(marker);
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-12 05:08:30 +0000