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.c1091
1 files changed, 625 insertions, 466 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/sub1sp.c b/Build/source/libs/mpfr/mpfr-src/src/sub1sp.c
index 2fc19db1d33..3b71e33a6de 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/sub1sp.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/sub1sp.c
@@ -1,7 +1,7 @@
/* mpfr_sub1sp -- internal function to perform a "real" subtraction
All the op must have the same precision
-Copyright 2003-2019 Free Software Foundation, Inc.
+Copyright 2003-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
@@ -24,8 +24,25 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
#define MPFR_NEED_LONGLONG_H
#include "mpfr-impl.h"
-/* Check if we have to check the result of mpfr_sub1sp with mpfr_sub1 */
+/* define MPFR_FULLSUB to use alternate code in mpfr_sub1sp2 and mpfr_sub1sp2n
+ (see comments in mpfr_sub1sp2) */
+/* #define MPFR_FULLSUB */
+
#if MPFR_WANT_ASSERT >= 2
+/* Check the result of mpfr_sub1sp with mpfr_sub1.
+
+ Note: mpfr_sub1sp itself has two algorithms: one always valid and one
+ faster for small precisions (up to 3 limbs). The latter one is disabled
+ if MPFR_GENERIC_ABI is defined. When MPFR_WANT_ASSERT >= 2, it could be
+ interesting to compare the results of these different algorithms. For
+ the time being, this is currently done by running the same code on the
+ same data with and without MPFR_GENERIC_ABI defined, where we have the
+ following comparisons in small precisions:
+ mpfr_sub1sp slow <-> mpfr_sub1 when MPFR_GENERIC_ABI is defined;
+ mpfr_sub1sp fast <-> mpfr_sub1 when MPFR_GENERIC_ABI is not defined.
+ By transitivity, the absence of failures implies that the 3 results are
+ the same.
+*/
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)
@@ -58,6 +75,10 @@ int mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
inexact = mpfr_sub1sp_ref (a, b, c, rnd_mode);
flags = __gmpfr_flags;
+ /* Convert the ternary values to (-1,0,1). */
+ inexact2 = VSIGN (inexact2);
+ inexact = VSIGN (inexact);
+
if (! mpfr_equal_p (tmpa, a) || inexact != inexact2 || flags != flags2)
{
fprintf (stderr, "sub1 & sub1sp return different values for %s\n"
@@ -73,9 +94,10 @@ int mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
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);
+ fprintf (stderr, "sub1 : ternary = %2d, flags =", inexact2);
+ flags_fout (stderr, flags2);
+ fprintf (stderr, "sub1sp: ternary = %2d, flags =", inexact);
+ flags_fout (stderr, flags);
MPFR_ASSERTN (0);
}
mpfr_clears (tmpa, tmpb, tmpc, (mpfr_ptr) 0);
@@ -86,6 +108,20 @@ int mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
#if !defined(MPFR_GENERIC_ABI)
+/* the sub1sp1_extracted.c is not ready yet */
+
+#if 0 && defined(MPFR_WANT_PROVEN_CODE) && GMP_NUMB_BITS == 64 && \
+ UINT_MAX == 0xffffffff && MPFR_PREC_BITS == 64 && \
+ _MPFR_PREC_FORMAT == 3 && _MPFR_EXP_FORMAT == _MPFR_PREC_FORMAT
+
+/* The code assumes that mp_limb_t has 64 bits exactly, unsigned int
+ has 32 bits exactly, mpfr_prec_t and mpfr_exp_t are of type long,
+ which has 64 bits exactly. */
+
+#include "sub1sp1_extracted.c"
+
+#else
+
/* 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,
@@ -107,8 +143,7 @@ mpfr_sub1sp1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
if (bx == cx)
{
- a0 = bp[0] - cp[0];
- if (a0 == 0) /* result is zero */
+ if (MPFR_UNLIKELY(bp[0] == cp[0])) /* result is zero */
{
if (rnd_mode == MPFR_RNDD)
MPFR_SET_NEG(a);
@@ -117,13 +152,16 @@ mpfr_sub1sp1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
MPFR_SET_ZERO(a);
MPFR_RET (0);
}
- else if (a0 > bp[0]) /* borrow: |c| > |b| */
+ else if (cp[0] > bp[0]) /* borrow: |c| > |b| */
{
+ a0 = cp[0] - bp[0];
MPFR_SET_OPPOSITE_SIGN (a, b);
- a0 = -a0;
}
else /* bp[0] > cp[0] */
- MPFR_SET_SAME_SIGN (a, b);
+ {
+ a0 = bp[0] - cp[0];
+ MPFR_SET_SAME_SIGN (a, b);
+ }
/* now a0 != 0 */
MPFR_ASSERTD(a0 != 0);
@@ -133,10 +171,21 @@ mpfr_sub1sp1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
rb = sb = 0;
/* Note: sh is not initialized, but will not be used in this case. */
}
- else if (bx > cx)
+ else
{
- MPFR_SET_SAME_SIGN (a, b);
- BGreater1:
+ if (bx < cx) /* swap b and c */
+ {
+ 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);
+ }
+ else
+ {
+ MPFR_SET_SAME_SIGN (a, b);
+ }
+ MPFR_ASSERTD (bx > cx);
d = (mpfr_uexp_t) bx - cx;
sh = GMP_NUMB_BITS - p;
mask = MPFR_LIMB_MASK(sh);
@@ -202,15 +251,6 @@ mpfr_sub1sp1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
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 */
@@ -266,6 +306,8 @@ mpfr_sub1sp1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
}
}
+#endif /* MPFR_WANT_PROVEN_CODE */
+
/* 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)
@@ -312,10 +354,21 @@ mpfr_sub1sp1n (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
bx -= cnt;
rb = sb = 0;
}
- else if (bx > cx)
+ else
{
- MPFR_SET_SAME_SIGN (a, b);
- BGreater1:
+ if (bx < cx) /* swap b and c */
+ {
+ 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);
+ }
+ else
+ {
+ MPFR_SET_SAME_SIGN (a, b);
+ }
+ MPFR_ASSERTD (bx > cx);
d = (mpfr_uexp_t) bx - cx;
if (d < GMP_NUMB_BITS)
{
@@ -369,8 +422,8 @@ mpfr_sub1sp1n (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
bx --;
if (d == GMP_NUMB_BITS && cp[0] > MPFR_LIMB_HIGHBIT)
{ /* case (b) */
- rb = (-cp[0]) >= (MPFR_LIMB_HIGHBIT >> 1);
- sb = (-cp[0]) << 2;
+ rb = MPFR_LIMB(-cp[0]) >= (MPFR_LIMB_HIGHBIT >> 1);
+ sb = MPFR_LIMB(-cp[0]) << 2;
ap[0] = -(MPFR_LIMB_ONE << 1);
}
else /* cases (a), (c), (d) and (e) */
@@ -388,15 +441,6 @@ mpfr_sub1sp1n (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
}
}
}
- 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 */
@@ -519,12 +563,23 @@ mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
rb = sb = 0;
/* Note: sh is not initialized, but will not be used in this case. */
}
- else if (bx > cx)
+ else
{
mp_limb_t t;
- MPFR_SET_SAME_SIGN (a, b);
- BGreater2:
+ if (bx < cx) /* swap b and c */
+ {
+ 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);
+ }
+ else
+ {
+ MPFR_SET_SAME_SIGN (a, b);
+ }
+ MPFR_ASSERTD (bx > cx);
d = (mpfr_uexp_t) bx - cx;
sh = 2 * GMP_NUMB_BITS - p;
mask = MPFR_LIMB_MASK(sh);
@@ -652,15 +707,6 @@ mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
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 */
@@ -712,6 +758,287 @@ mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
}
}
+/* special code for p = 2*GMP_NUMB_BITS */
+static int
+mpfr_sub1sp2n (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, a1;
+ mpfr_uexp_t d;
+
+/* this function is inspired by mpfr_sub1sp2 (for the operations of the
+ 2-limb arrays) and by mpfr_sub1sp1n (for the different cases to handle) */
+
+ if (bx == cx) /* subtraction is exact in this case */
+ {
+ 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);
+ }
+ /* since B/2 <= bp[1], cp[1] < B with B=2^GMP_NUMB_BITS,
+ if no borrow we have 0 <= bp[1] - cp[1] - x < B/2
+ where x = (bp[0] < cp[0]) is 0 or 1, thus a1 < B/2 <= bp[1] */
+ else if (a1 >= bp[1]) /* borrow: |c| > |b| */
+ {
+ MPFR_SET_OPPOSITE_SIGN (a, b);
+ /* negate [a1,a0] */
+ a0 = -a0;
+ a1 = -a1 - (a0 != 0);
+ }
+ else /* no borrow */
+ MPFR_SET_SAME_SIGN (a, b);
+
+ /* now [a1,a0] is the absolute value of b - c,
+ maybe not normalized */
+ 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)
+ {
+ /* shift [a1,a0] left by cnt bit and store in result */
+ 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; /* the subtraction is exact */
+ }
+ else
+ {
+ mp_limb_t t;
+
+ if (bx < cx) /* swap b and c */
+ {
+ 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);
+ }
+ else
+ {
+ MPFR_SET_SAME_SIGN (a, b);
+ }
+ MPFR_ASSERTD (bx > cx);
+ d = (mpfr_uexp_t) bx - cx;
+ if (d < GMP_NUMB_BITS)
+ {
+ t = (cp[1] << (GMP_NUMB_BITS - d)) | (cp[0] >> d);
+ /* t is the part that should be subtracted to bp[0]:
+ | a1 | a0 |
+ | bp[1] | bp[0] |
+ | cp[1]>>d | t | sb | */
+#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 */
+ /* now negate sb and subtract borrow to a0 if sb <> 0 */
+ if (sb)
+ {
+ a1 -= (a0 == 0);
+ a0 --;
+ /* a = a1,a0 can only be zero when d=1, b = 0.1000...000*2^bx,
+ and c = 0.111...111*2^(bx-1). In that case (where we have
+ sb = MPFR_LIMB_HIGHBIT below), the subtraction is exact, the
+ result is b/2^(2*GMP_NUMB_BITS). This case is dealt below. */
+ sb = -sb;
+ }
+#else
+ 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
+ /* now the result is formed of [a1,a0,sb], which might not be
+ normalized */
+ if (a1 == MPFR_LIMB_ZERO)
+ {
+ /* this implies d=1 */
+ MPFR_ASSERTD(d == 1);
+ a1 = a0;
+ a0 = sb;
+ sb = MPFR_LIMB_ZERO;
+ bx -= GMP_NUMB_BITS;
+ }
+ if (a1 == MPFR_LIMB_ZERO) /* case a = a1,a0 = 0 mentioned above */
+ {
+ MPFR_ASSERTD(a0 == MPFR_LIMB_HIGHBIT); /* was sb above */
+ a1 = a0;
+ a0 = sb;
+ bx -= GMP_NUMB_BITS;
+ sb = MPFR_LIMB_ZERO;
+ }
+ else
+ {
+ count_leading_zeros (cnt, a1);
+ if (cnt)
+ {
+ /* shift [a1,a0,sb] left by cnt bits and adjust exponent */
+ a1 = (a1 << cnt) | (a0 >> (GMP_NUMB_BITS - cnt));
+ a0 = (a0 << cnt) | (sb >> (GMP_NUMB_BITS - cnt));
+ sb <<= cnt;
+ bx -= cnt;
+ }
+ }
+ rb = sb & MPFR_LIMB_HIGHBIT;
+ sb = sb & ~MPFR_LIMB_HIGHBIT;
+ ap[1] = a1;
+ ap[0] = a0;
+ }
+ else if (d < 2 * GMP_NUMB_BITS)
+ { /* GMP_NUMB_BITS <= d < 2*GMP_NUMB_BITS:
+ compute t, the part to be subtracted to bp[0],
+ and sb, the neglected part of c:
+ | a1 | a0 |
+ | bp[1] | bp[0] |
+ | t | sb | */
+ /* warning: we should not ignore the low bits from cp[0]
+ in case d > GMP_NUMB_BITS */
+ sb = (d == GMP_NUMB_BITS) ? cp[0]
+ : (cp[1] << (2*GMP_NUMB_BITS - d))
+ | (cp[0] >> (d - GMP_NUMB_BITS))
+ | ((cp[0] << (2*GMP_NUMB_BITS - d)) != 0);
+ t = (cp[1] >> (d - GMP_NUMB_BITS)) + (sb != 0);
+ /* Warning: t might overflow to 0 if d=GMP_NUMB_BITS, sb <> 0,
+ and cp[1] = 111...111 */
+ a0 = bp[0] - t;
+ a1 = bp[1] - (bp[0] < t) - (t == 0 && sb != 0);
+ sb = -sb;
+ /* now the result is [a1,a0,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)
+ {
+ /* shift [a1,a0] left by 1 bit */
+ a1 = (a1 << 1) | (a0 >> (GMP_NUMB_BITS - 1));
+ MPFR_ASSERTD(a1 >= MPFR_LIMB_HIGHBIT);
+ a0 = (a0 << 1) | (sb >> (GMP_NUMB_BITS - 1));
+ sb <<= 1;
+ bx --;
+ }
+ ap[1] = a1;
+ ap[0] = a0;
+ rb = sb & MPFR_LIMB_HIGHBIT;
+ sb = sb & ~MPFR_LIMB_HIGHBIT;
+ }
+ else
+ { /* d >= 2*GMP_NUMB_BITS:
+ | a1 | a0 |
+ | bp[1] | bp[0] |
+ | cp[1] | cp[0] | */
+ /* we mimic the case d >= GMP_NUMB_BITS of mpfr_sub1sp1n */
+ int tst = cp[1] == MPFR_LIMB_HIGHBIT && cp[0] == MPFR_LIMB_ZERO;
+ /* if d = 2 * GMP_NUMB_BITS and tst=1, c = 1/2*ulp(b) */
+ if (bp[1] > MPFR_LIMB_HIGHBIT || bp[0] > MPFR_LIMB_ZERO)
+ {
+ /* no borrow in b - ulp(b) */
+ rb = d > 2 * GMP_NUMB_BITS || tst;
+ sb = d > 2 * GMP_NUMB_BITS || !tst;
+ ap[1] = bp[1] - (bp[0] == MPFR_LIMB_ZERO);
+ ap[0] = bp[0] - MPFR_LIMB_ONE;
+ }
+ else
+ {
+ /* b = 1000...000, thus subtracting c yields an exponent shift */
+ bx --;
+ if (d == 2 * GMP_NUMB_BITS && !tst) /* c > 1/2*ulp(b) */
+ {
+ t = -cp[1] - (cp[0] > MPFR_LIMB_ZERO);
+ /* the rounding bit is the 2nd most significant bit of t
+ (where the most significant bit of t is necessarily 0),
+ and the sticky bit is formed by the remaining bits of t,
+ and those from -cp[0] */
+ rb = t >= (MPFR_LIMB_HIGHBIT >> 1);
+ sb = (t << 2) | cp[0];
+ ap[1] = MPFR_LIMB_MAX;
+ ap[0] = -(MPFR_LIMB_ONE << 1);
+ }
+ else /* c <= 1/2*ulp(b) */
+ {
+ rb = d > 2 * GMP_NUMB_BITS + 1
+ || (d == 2 * GMP_NUMB_BITS + 1 && tst);
+ sb = d > 2 * GMP_NUMB_BITS + 1
+ || (d == 2 * GMP_NUMB_BITS + 1 && !tst);
+ ap[1] = -MPFR_LIMB_ONE;
+ ap[0] = -MPFR_LIMB_ONE;
+ }
+ }
+ }
+ }
+
+ /* 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) == 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;
+ 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,
@@ -791,10 +1118,21 @@ mpfr_sub1sp3 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
rb = sb = 0;
/* Note: sh is not initialized, but will not be used in this case. */
}
- else if (bx > cx)
+ else
{
- MPFR_SET_SAME_SIGN (a, b);
- BGreater2:
+ if (bx < cx) /* swap b and c */
+ {
+ 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);
+ }
+ else
+ {
+ MPFR_SET_SAME_SIGN (a, b);
+ }
+ MPFR_ASSERTD (bx > cx);
d = (mpfr_uexp_t) bx - cx;
sh = 3 * GMP_NUMB_BITS - p;
mask = MPFR_LIMB_MASK(sh);
@@ -976,15 +1314,6 @@ mpfr_sub1sp3 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
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 */
@@ -1097,14 +1426,21 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
mpfr_exp_t bx, cx;
mpfr_uexp_t d;
mpfr_prec_t p, sh, cnt;
- mp_size_t n;
- mp_limb_t *ap, *bp, *cp;
+ mp_size_t n, k;
+ mp_limb_t *ap = MPFR_MANT(a);
+ mp_limb_t *bp = MPFR_MANT(b);
+ mp_limb_t *cp = MPFR_MANT(c);
mp_limb_t limb;
int inexact;
- mp_limb_t bcp,bcp1; /* Cp and C'p+1 */
- mp_limb_t bbcp = MPFR_LIMB_MAX, bbcp1 = MPFR_LIMB_MAX; /* Cp+1 and C'p+2,
+ mp_limb_t rb, sb; /* round and sticky bits. They are interpreted as
+ negative, i.e., if rb <> 0, then we should subtract 1
+ at the round bit position, and of sb <> 0, we should
+ subtract something below the round bit position. */
+ mp_limb_t rbb = MPFR_LIMB_MAX, sbb = MPFR_LIMB_MAX; /* rbb is the next bit
+ after the round bit, and sbb the corresponding sticky bit.
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. */
+ int pow2;
MPFR_TMP_DECL(marker);
@@ -1132,6 +1468,9 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
/* 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);
+
+ if (p == 2 * GMP_NUMB_BITS)
+ return mpfr_sub1sp2n (a, b, c, rnd_mode);
#endif
n = MPFR_PREC2LIMBS (p);
@@ -1141,14 +1480,15 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
MPFR_TMP_MARK(marker);
+ k = n - 1;
+
if (bx == cx)
{
- mp_size_t k = n - 1;
/* Check mantissa since exponents are equal */
- bp = MPFR_MANT(b);
- cp = MPFR_MANT(c);
while (k >= 0 && MPFR_UNLIKELY(bp[k] == cp[k]))
k--;
+ /* now k = - 1 if b == c, otherwise k is the largest integer < n such
+ that bp[k] <> cp[k] */
if (k < 0)
/* b == c ! */
{
@@ -1173,10 +1513,13 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
/* Swap b and c and set sign */
mpfr_srcptr t;
mpfr_exp_t tx;
+ mp_limb_t *tp;
+
+ tx = bx; bx = cx; cx = tx;
CGreater:
MPFR_SET_OPPOSITE_SIGN(a,b);
t = b; b = c; c = t;
- tx = bx; bx = cx; cx = tx;
+ tp = bp; bp = cp; cp = tp;
}
else
{
@@ -1190,25 +1533,26 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
d = (mpfr_uexp_t) bx - cx;
/* printf ("New with diff=%lu\n", (unsigned long) d); */
- if (d <= 1)
+ /* FIXME: The goto's below are too complex (long backward) and make
+ the code unreadable. */
+
+ if (d == 0)
{
- if (d == 0)
- {
- /* <-- b -->
- <-- c --> : exact sub */
- ap = MPFR_MANT(a);
- mpn_sub_n (ap, MPFR_MANT(b), MPFR_MANT(c), n);
+ /* <-- b -->
+ <-- c --> : exact sub */
+ mpn_sub_n (ap, bp, cp, n);
/* Normalize */
ExactNormalize:
limb = ap[n-1];
if (MPFR_LIKELY (limb != 0))
{
/* First limb is not zero. */
- count_leading_zeros(cnt, limb);
- /* cnt could be == 0 <= SubD1Lose */
+ count_leading_zeros (cnt, limb);
+ /* Warning: cnt can be 0 when we come from the case SubD1Lose
+ with goto ExactNormalize */
if (MPFR_LIKELY(cnt))
{
- mpn_lshift(ap, ap, n, cnt); /* Normalize number */
+ mpn_lshift (ap, ap, n, cnt); /* Normalize number */
bx -= cnt; /* Update final expo */
}
/* Last limb should be OK */
@@ -1217,35 +1561,33 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
}
else
{
- /* First limb is zero */
- mp_size_t k = n-1, len;
+ /* First limb is zero: this can only occur for n >= 2 */
+ mp_size_t len;
/* Find the first limb not equal to zero. It necessarily exists
- since |b| > |c|. */
- do
- {
- MPFR_ASSERTD( k > 0 );
- limb = ap[--k];
- }
- while (limb == 0);
+ since |b| > |c|. We know that bp[k] > cp[k] and all upper
+ limbs are equal. */
+ while (ap[k] == 0)
+ k--;
+ limb = ap[k];
+ /* ap[k] is the non-zero limb of largest index, thus we have
+ to consider the k+1 least significant limbs */
MPFR_ASSERTD(limb != 0);
count_leading_zeros(cnt, limb);
k++;
- len = n - k; /* Number of last limb */
- MPFR_ASSERTD(k >= 0);
+ len = n - k; /* Number of most significant zero limbs */
+ MPFR_ASSERTD(k > 0);
if (cnt)
mpn_lshift (ap + len, ap, k, cnt); /* Normalize the High Limb*/
else
- {
- /* Must use copyd since src and dst may overlap & dst>=src */
- mpn_copyd (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 */
+ bx -= cnt + len * GMP_NUMB_BITS; /* update exponent */
+ /* ap[len] should have its low bits zero: it is bp[0]-cp[0] */
MPFR_ASSERTD(!(ap[len] & MPFR_LIMB_MASK((unsigned int) (-p)
% GMP_NUMB_BITS)));
}
- /* Check expo underflow */
+ /* Check exponent underflow (no overflow can happen) */
if (MPFR_UNLIKELY(bx < __gmpfr_emin))
{
MPFR_TMP_FREE(marker);
@@ -1261,122 +1603,110 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
}
MPFR_SET_EXP (a, bx);
/* No rounding is necessary since the result is exact */
- MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
+ MPFR_ASSERTD(ap[n-1] & MPFR_LIMB_HIGHBIT);
MPFR_TMP_FREE(marker);
return 0;
- }
- else /* if (d == 1) */
+ }
+ else if (d == 1)
{
/* | <-- b -->
| <-- c --> */
mp_limb_t c0, mask;
- mp_size_t k;
MPFR_UNSIGNED_MINUS_MODULO(sh, p);
/* If we lose at least one bit, compute 2*b-c (Exact)
* else compute b-c/2 */
- bp = MPFR_MANT(b);
- 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.
+ /* Let W = 2^GMP_NUMB_BITS:
+ we have |b|-|c| >= limb*W^k - (2*W^k-1)/2 >= limb*W^k - W^k + 1/2
+ thus if limb > W/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
+ minus the last significant bit called c0 below (in fact c0 is that
+ bit shifted by sh bits), 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)
+ if (limb > MPFR_LIMB_HIGHBIT) /* case limb > W/2 */
{
+ mp_limb_t *tp;
/* The exponent cannot decrease: compute b-c/2 */
/* Shift c in the allocated temporary block */
SubD1NoLose:
c0 = cp[0] & (MPFR_LIMB_ONE << sh);
mask = ~MPFR_LIMB_MASK(sh);
- ap = MPFR_MANT(a);
- cp = MPFR_TMP_LIMBS_ALLOC (n);
+ tp = MPFR_TMP_LIMBS_ALLOC (n);
/* 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]);
+ mpn_rshift (tp, cp, n, 1);
+ tp[0] &= mask; /* Zero last bit of c if set */
+ mpn_sub_n (ap, bp, tp, n);
+ MPFR_SET_EXP(a, bx); /* No exponent overflow! */
+ MPFR_ASSERTD(ap[n-1] & MPFR_LIMB_HIGHBIT);
if (MPFR_LIKELY(c0 == 0))
{
/* Result is exact: no need of rounding! */
MPFR_TMP_FREE(marker);
return 0;
}
+ /* c0 is non-zero, thus we have to subtract 1/2*ulp(a),
+ however we know (see analysis above) that this cannot
+ make the exponent decrease */
MPFR_ASSERTD( !(ap[0] & ~mask) ); /* Check last bits */
/* No normalize is needed */
- /* Rounding is necessary since c0 = 1 */
- /* Cp =-1 and C'p+1=0 */
- bcp = 1; bcp1 = 0;
-
- 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 last bit of a. */
- if (MPFR_LIKELY( (ap[0] & (MPFR_LIMB_ONE << sh)) == 0) )
- goto truncate;
- else
- goto sub_one_ulp;
- }
- MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
- if (rnd_mode == MPFR_RNDZ)
- goto sub_one_ulp;
- else
- goto truncate;
+
+ /* Rounding is necessary since c0 is non-zero */
+ /* we have to subtract 1 at the round bit position,
+ and 0 for the lower bits */
+ rb = 1; rbb = sbb = 0;
}
else if (MPFR_LIKELY(limb < MPFR_LIMB_HIGHBIT))
{
+ mp_limb_t *tp;
/* |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)
+#if defined(WANT_GMP_INTERNALS) && defined(HAVE___GMPN_RSBLSH1_N)
/* {ap, n} = 2*{bp, n} - {cp, n} */
- __gmpn_rsblsh_n (ap, MPFR_MANT(c), MPFR_MANT(b), n, 1);
+ /* mpn_rsblsh1_n -- rp[] = (vp[] << 1) - up[] */
+ __gmpn_rsblsh1_n (ap, cp, bp, n);
#else
- bp = MPFR_TMP_LIMBS_ALLOC (n);
+ tp = MPFR_TMP_LIMBS_ALLOC (n);
/* Shift b in the allocated temporary block */
- mpn_lshift (bp, MPFR_MANT(b), n, 1);
- mpn_sub_n (ap, bp, cp, n);
+ mpn_lshift (tp, bp, n, 1);
+ mpn_sub_n (ap, tp, cp, n);
#endif
bx--;
+ MPFR_ASSERTD(k == n-1);
goto ExactNormalize;
}
- else
+ else /* limb = MPFR_LIMB_HIGHBIT */
{
/* Case: limb = 100000000000 */
- /* Check while b[k] == c'[k] (C' is C shifted by 1) */
- /* If b[k]<c'[k] => We lose at least one bit*/
- /* If b[k]>c'[k] => We don't lose any bit */
- /* If k==-1 => We don't lose any bit
+ /* Check while b[l] == c'[l] (C' is C shifted by 1) */
+ /* If b[l]<c'[l] => We lose at least one bit */
+ /* If b[l]>c'[l] => We don't lose any bit */
+ /* If l==-1 => We don't lose any bit
AND the result is 100000000000 0000000000 00000000000 */
- mp_limb_t carry;
+ mp_size_t l = n - 1;
+ mp_limb_t cl_shifted;
do
{
- carry = cp[k] << (GMP_NUMB_BITS - 1);
- if (--k < 0)
+ /* the first loop will compare b[n-2] and c'[n-2] */
+ cl_shifted = cp[l] << (GMP_NUMB_BITS - 1);
+ if (--l < 0)
break;
- carry += cp[k] >> 1;
+ cl_shifted += cp[l] >> 1;
}
- while (bp[k] == carry);
- if (MPFR_UNLIKELY(k < 0))
+ while (bp[l] == cl_shifted);
+ if (MPFR_UNLIKELY(l < 0))
{
- ap = MPFR_MANT (a);
- if (MPFR_UNLIKELY(carry))
+ if (MPFR_UNLIKELY(cl_shifted))
{
- /* 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. */
+ /* Since cl_shifted is what should be subtracted
+ from ap[-1], if non-zero then necessarily the
+ precision is a 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
@@ -1384,138 +1714,79 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
}
else
{
- /* No carry: result is a power of 2. */
+ /* cl_shifted=0: 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! */
+ MPFR_SET_EXP (a, bx); /* No exponent 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)
+ /* cl_shifted is the shifted value c'[l] */
+ else if (bp[l] > cl_shifted)
goto SubD1NoLose; /* |b|-|c| >= 1/2*W^n */
else
{
- MPFR_ASSERTD(bp[k] < carry);
+ /* we cannot have bp[l] = cl_shifted since the only way we
+ can exit the while loop above is when bp[l] <> cl_shifted
+ or l < 0, and the case l < 0 was already treated above. */
+ MPFR_ASSERTD(bp[l] < cl_shifted);
goto SubD1Lose; /* |b|-|c| <= 1/2*W^n-1 and is exact */
}
}
}
- }
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);
/* We can't set A before since we use cp for rounding... */
/* Perform rounding: check if a=b or a=b-ulp(b) */
if (MPFR_UNLIKELY(d == p))
{
- /* cp == -1 and c'p+1 = ? */
- bcp = 1;
- /* We need Cp+1 later for a very improbable case. */
- bbcp = (MPFR_MANT(c)[n-1] & (MPFR_LIMB_ONE<<(GMP_NUMB_BITS-2)));
- /* We need also C'p+1 for an even more unprobable case... */
- if (MPFR_LIKELY( bbcp ))
- bcp1 = 1;
- else
- {
- cp = MPFR_MANT(c);
- if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
- {
- mp_size_t k = n-1;
- do
- k--;
- while (k >= 0 && cp[k] == 0);
- bcp1 = (k >= 0);
- }
- else
- bcp1 = 1;
- }
- /* 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_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) */
- {
- 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)
- {
- MPN_COPY(ap, bp, n);
- goto truncate;
- }
- MPN_COPY(ap, bp, n);
- goto sub_one_ulp;
- }
- MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
- if (rnd_mode == MPFR_RNDZ)
- {
- MPN_COPY(ap, bp, n);
- goto sub_one_ulp;
- }
- else
- {
- MPN_COPY(ap, bp, n);
- goto truncate;
- }
+ /* since c is normalized, we need to subtract 1/2 ulp(b) */
+ rb = 1;
+ /* rbb is the bit of weight 1/4 ulp(b) in c. We assume a limb has
+ at least 2 bits. If the precision is 1, we read in the unused
+ bits, which should be zero, and this is what we want. */
+ rbb = cp[n-1] & (MPFR_LIMB_HIGHBIT >> 1);
+
+ /* We need also sbb */
+ sbb = cp[n-1] & MPFR_LIMB_MASK(GMP_NUMB_BITS - 2);
+ for (k = n-1; sbb == 0 && k > 0;)
+ sbb = cp[--k];
}
else
{
- /* Cp=0, Cp+1=-1 if d==p+1, C'p+1=-1 */
- bcp = 0; bbcp = (d==p+1); bcp1 = 1;
- /* 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))
+ rb = 0;
+ if (d == p + 1)
{
- cp = MPFR_MANT(c);
- if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
- {
- mp_size_t k = n-1;
- do
- k--;
- while (k >= 0 && cp[k] == 0);
- bbcp1 = (k >= 0);
- }
- else
- bbcp1 = 1;
- /* printf("(D>P) C'p+2=%d\n", bbcp1!=0); */
+ rbb = 1;
+ sbb = cp[n-1] & MPFR_LIMB_MASK(GMP_NUMB_BITS - 1);
+ for (k = n-1; sbb == 0 && k > 0;)
+ sbb = cp[--k];
+ }
+ else
+ {
+ rbb = 0;
+ sbb = 1; /* since C is non-zero */
}
- /* Copy mantissa B in A */
- MPN_COPY(ap, MPFR_MANT(b), n);
- /* Round */
- 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)
- goto sub_one_ulp;
- else /* rnd_mode = AWAY */
- goto truncate;
}
+ /* Copy mantissa B in A */
+ MPN_COPY(ap, bp, n);
}
else /* case 2 <= d < p */
{
mpfr_uexp_t dm;
mp_size_t m;
- mp_limb_t mask;
+ mp_limb_t mask, *tp;
MPFR_UNSIGNED_MINUS_MODULO(sh, p);
- cp = MPFR_TMP_LIMBS_ALLOC (n);
+ tp = MPFR_TMP_LIMBS_ALLOC (n);
/* Shift c in temporary allocated place */
dm = d % GMP_NUMB_BITS;
@@ -1524,276 +1795,164 @@ mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
/* dm = 0 and m > 0: Just copy */
MPFR_ASSERTD(m != 0);
- MPN_COPY(cp, MPFR_MANT(c)+m, n-m);
- MPN_ZERO(cp+n-m, m);
+ MPN_COPY(tp, cp + m, n - m);
+ MPN_ZERO(tp + n - m, m);
}
else if (MPFR_LIKELY(m == 0))
{
/* dm >=2 and m == 0: just shift */
MPFR_ASSERTD(dm >= 2);
- mpn_rshift(cp, MPFR_MANT(c), n, dm);
+ mpn_rshift (tp, cp, n, dm);
}
else
{
/* dm > 0 and m > 0: shift and zero */
- mpn_rshift(cp, MPFR_MANT(c)+m, n-m, dm);
- MPN_ZERO(cp+n-m, m);
+ mpn_rshift (tp, cp + m, n - m, dm);
+ MPN_ZERO (tp + n - m, m);
}
+ /* FIXME: Instead of doing "cp = tp;", keep using tp to avoid
+ confusion? Thus in the block below, we don't need
+ "mp_limb_t *cp = MPFR_MANT(c);". In short, cp should always
+ be MPFR_MANT(c) defined earlier, possibly after the swap. */
+ cp = tp;
/* 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 (cp[0] & MPFR_LIMB_MASK(sh-1))
- bcp1 = 1;
- else
- {
- /* We can't compute C'p+1 from C'. Compute it from C */
- /* 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 (x > p)
- /* We are already looked at all the bits of c, so C'p+1 = 0*/
- bcp1 = 0;
- else
- {
- 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);
- /* 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;
- else
- {
- /*kx += (sx==0);*/
- /*If sx==0, tp[kx] hasn't been checked*/
- do
- kx--;
- while (kx >= 0 && tp[kx] == 0);
- bcp1 = (kx >= 0);
- }
- }
- }
- }
- else
- {
- /* Compute Cp and C'p+1 from C with sh=0 */
- mp_limb_t *tp = MPFR_MANT(c);
- /* Start from bit x=p-d in mantissa C */
- mpfr_prec_t x = p-d;
- 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);
- bcp = (tp[kx] & (MPFR_LIMB_ONE<<sx));
- /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
- if (tp[kx] & MPFR_LIMB_MASK(sx))
- bcp1 = 1;
- else
- {
- /*kx += (sx==0);*/ /*If sx==0, tp[kx] hasn't been checked*/
- do
- kx--;
- while (kx >= 0 && tp[kx] == 0);
- bcp1 = (kx >= 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))
- {
- /* 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 */
- if (MPFR_LIKELY(bcp1 == 0))
- {
- bbcp = 0;
- bbcp1 = 0;
- }
- else /* bcp1 != 0 */
- {
- /* We can lose a bit:
- compute Cp+1 and C'p+2 from mantissa C */
- 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);
-
- MPFR_ASSERTD(p > d);
- /* 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))))
- bbcp1 = 1;
- else
- {
- /*kx += (sx==0);*/ /*If sx==0, tp[kx] hasn't been checked*/
- do
- kx--;
- while (kx >= 0 && tp[kx] == 0);
- bbcp1 = (kx >= 0);
- /* printf ("(Pre) Scan done for %ld\n", (long) kx); */
- }
- } /*End of Bcp1 != 0*/
- /* printf("(Pre) Cp+1=%d C'p+2=%d\n", bbcp!=0, bbcp1!=0); */
- } /* End of "can lose a bit" */
+ /* Compute rb=Cp and sb=C'p+1 */
+ {
+ /* Compute rb and rbb from C */
+ mp_limb_t *cp = MPFR_MANT(c);
+ /* The round bit is bit p-d in C, assuming the most significant bit
+ of C is bit 0 */
+ mpfr_prec_t x = p - d;
+ mp_size_t kx = n - 1 - (x / GMP_NUMB_BITS);
+ mpfr_prec_t sx = GMP_NUMB_BITS - 1 - (x % GMP_NUMB_BITS);
+ /* the round bit is in cp[kx], at position sx */
+ MPFR_ASSERTD(p >= d);
+ rb = cp[kx] & (MPFR_LIMB_ONE << sx);
+
+ /* Now compute rbb: since d >= 2 it always exists in C */
+ if (sx == 0) /* rbb is in the next limb */
+ {
+ kx --;
+ sx = GMP_NUMB_BITS - 1;
+ }
+ else
+ sx --; /* rb and rbb are in the same limb */
+ rbb = cp[kx] & (MPFR_LIMB_ONE << sx);
+
+ /* Now look at the remaining low bits of C to determine sbb */
+ sbb = cp[kx] & MPFR_LIMB_MASK(sx);
+ while (sbb == 0 && kx > 0)
+ sbb = cp[--kx];
+ }
+ /* printf("sh=%lu Cp=%d C'p+1=%d\n", sh, rb!=0, sb!=0); */
/* Clean shifted C' */
mask = ~MPFR_LIMB_MASK (sh);
cp[0] &= mask;
/* Subtract the mantissa c from b in a */
- ap = MPFR_MANT(a);
mpn_sub_n (ap, bp, cp, n);
/* mpfr_print_mant_binary("Sub= ", ap, p); */
- /* Normalize: we lose at max one bit*/
+ /* Normalize: we lose at most one bit */
if (MPFR_UNLIKELY(MPFR_LIMB_MSB(ap[n-1]) == 0))
{
/* High bit is not set and we have to fix it! */
/* Ap >= 010000xxx001 */
- mpn_lshift(ap, ap, n, 1);
+ mpn_lshift (ap, ap, n, 1);
/* Ap >= 100000xxx010 */
- if (MPFR_UNLIKELY(bcp != 0)) /* Check if Cp = -1 */
+ if (MPFR_UNLIKELY(rb != 0)) /* Check if Cp = -1 */
/* Since Cp == -1, we have to subtract one more */
{
- mpn_sub_1(ap, ap, n, MPFR_LIMB_ONE<<sh);
+ mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
MPFR_ASSERTD(MPFR_LIMB_MSB(ap[n-1]) != 0);
}
/* Ap >= 10000xxx001 */
/* Final exponent -1 since we have shifted the mantissa */
bx--;
- /* Update bcp and bcp1 */
- MPFR_ASSERTD(bbcp != MPFR_LIMB_MAX);
- MPFR_ASSERTD(bbcp1 != MPFR_LIMB_MAX);
- bcp = bbcp;
- bcp1 = bbcp1;
+ /* Update rb and sb */
+ rb = rbb;
+ rbb = sbb;
/* 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_RNDF))
- goto truncate;
- else if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
- {
- if (MPFR_LIKELY(bcp == 0))
- goto truncate;
- else if (bcp1 != 0 || (ap[0] & (MPFR_LIMB_ONE << sh)) != 0)
- goto sub_one_ulp;
- else
- goto truncate;
- }
+ rounding:
+ /* at this point 'a' contains b - high(c), normalized,
+ and we have to subtract rb * 1/2 ulp(a), rbb * 1/4 ulp(a),
+ and sbb * 1/8 ulp(a), interpreting rb/rbb/sbb as 1 if non-zero. */
- /* Update rounding mode */
- MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
- if (rnd_mode == MPFR_RNDZ && MPFR_LIKELY (bcp != 0 || bcp1 != 0))
- goto sub_one_ulp;
- goto truncate;
+ sb = rbb | sbb;
+
+ if (rb == 0 && sb == 0)
+ {
+ inexact = 0;
+ goto end_of_sub;
}
- MPFR_RET_NEVER_GO_HERE ();
-
- /* Sub one ulp to the result */
- sub_one_ulp:
- mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
- /* Result should be smaller than exact value: inexact=-1 */
- inexact = -1;
- /* Check normalization */
- if (MPFR_UNLIKELY(ap[n-1] < MPFR_LIMB_HIGHBIT))
+
+ pow2 = mpfr_powerof2_raw (a);
+ if (pow2 && rb != 0) /* subtract 1 ulp */
{
- /* ap was a power of 2, and we lose a bit */
- /* Now it is 0111111111111111111[00000 */
- /* The following 2 lines are equivalent to: mpn_lshift(ap, ap, n, 1); */
- ap[0] <<= 1;
+ mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
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 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_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;
- if (rnd_mode == MPFR_RNDN)
- inexact = 1;
- /* 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 (rnd_mode == MPFR_RNDZ && bcp1 == 0)
- {
- /* printf("(SubOneUlp) Exact result\n"); */
- inexact = 0;
- }
+ rb = rbb;
+ rbb = sbb;
+ sbb = 0;
+ /* Note: for p=1 this case can only happen with d=1, but then we will
+ have rb=sb=0 at the next round. */
+ goto rounding;
}
- goto end_of_sub;
-
- truncate:
- /* Check if the result is an exact power of 2: 100000000000
- in which cases, we could have to do sub_one_ulp due to some nasty reasons:
- If Result is a Power of 2:
- + If rnd = AWAY,
- | If Cp=-1 and C'p+1 = 0, SubOneUlp and the result is EXACT.
- If Cp=-1 and C'p+1 =-1, SubOneUlp and the result is above.
- Otherwise truncate
- + If rnd = NEAREST,
- If Cp= 0 and Cp+1 =-1 and C'p+2=-1, SubOneUlp and the result is above
- If cp=-1 and C'p+1 = 0, SubOneUlp and the result is exact.
- Otherwise truncate.
- X bit should always be set if SubOneUlp*/
- if (MPFR_UNLIKELY(ap[n-1] == MPFR_LIMB_HIGHBIT))
+ /* now if a is a power of two, necessary rb = 0,
+ which means the exact result is always in (pred(a), a),
+ and the bounds cannot be attained */
+
+ if (rnd_mode == MPFR_RNDF)
+ inexact = 0;
+ else if (rnd_mode == MPFR_RNDN)
{
- mp_size_t k = n-1;
- do
- k--;
- while (k >= 0 && ap[k] == 0);
- if (MPFR_UNLIKELY (k < 0))
+ if (pow2)
{
- /* It is a power of 2! */
- /* Compute Cp+1 if it isn't already compute (ie d==1) */
- /* 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))
+ MPFR_ASSERTD(rb == 0);
+ /* since we have at the end of the binade, we have in fact rb = rbb
+ and sb = sbb */
+ rb = rbb;
+ sb = sbb;
+ }
+ /* Warning: for p=1, the significand is always odd: the "even" rule
+ rounds to the value with largest magnitude, thus we have to check
+ that case separately */
+ if (rb == 0 ||
+ (rb != 0 && sb == 0 &&
+ ((ap[0] & (MPFR_LIMB_ONE << sh)) == 0 || p == 1)))
+ inexact = 1; /* round to a and return 1 */
+ else /* round to pred(a) and return -1 */
+ {
+ subtract:
+ mpn_sub_1 (ap, ap, n, MPFR_LIMB_ONE << sh);
+ if (pow2) /* deal with cancellation */
{
- /* printf("(Truncate) Do sub\n"); */
- mpn_sub_1 (ap, ap, n, 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;
- goto end_of_sub;
}
+ inexact = -1;
}
}
-
- /* Calcul of Inexact flag.*/
- inexact = (bcp != 0 || bcp1 != 0);
+ else /* directed rounding */
+ {
+ MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(a));
+ if (rnd_mode == MPFR_RNDZ)
+ goto subtract;
+ else
+ inexact = 1;
+ }
end_of_sub:
/* Update Exponent */