diff options
author | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2016-12-08 10:32:38 +0000 |
---|---|---|
committer | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2016-12-08 10:32:38 +0000 |
commit | 2480f0ab60393450888351a4f9a9c9703045f1f0 (patch) | |
tree | f131d4bb827e2c2de6a388697da70613006b4478 | |
parent | 9e3649c443a284c8727f9bb0f191f70d25c3146f (diff) | |
download | mpfr-2480f0ab60393450888351a4f9a9c9703045f1f0.tar.gz |
[sub1sp.c] improved further
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@10991 280ebfd0-de03-0410-8827-d642c229c3f4
-rw-r--r-- | src/sub1sp.c | 112 |
1 files changed, 56 insertions, 56 deletions
diff --git a/src/sub1sp.c b/src/sub1sp.c index 0d4599b4b..9f11d5ca0 100644 --- a/src/sub1sp.c +++ b/src/sub1sp.c @@ -307,7 +307,7 @@ mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, mpfr_prec_t cnt, INITIALIZED(sh); mp_limb_t rb; /* round bit */ mp_limb_t sb; /* sticky bit */ - mp_limb_t mask, b0, b1; + mp_limb_t mask, a0, a1; mpfr_uexp_t d; MPFR_ASSERTD(GMP_NUMB_BITS < p && p < 2 * GMP_NUMB_BITS); @@ -315,11 +315,9 @@ mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, /* save b0, b1 in case a and b are the same variable */ if (bx == cx) /* subtraction is exact in this case */ { - b0 = bp[0]; - b1 = bp[1]; - ap[0] = b0 - cp[0]; - ap[1] = b1 - cp[1] - (ap[0] > b0); - if (ap[1] == 0 && ap[0] == 0) /* result is zero */ + a1 = bp[1] - cp[1] - (bp[0] < cp[0]); + a0 = bp[0] - cp[0]; + if (a1 == 0 && a0 == 0) /* result is zero */ { if (rnd_mode == MPFR_RNDD) MPFR_SET_NEG(a); @@ -328,30 +326,30 @@ mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, MPFR_SET_ZERO(a); return 0; /* same as MPFR_RET(0) but faster */ } - else if (ap[1] >= b1) /* borrow: |c| > |b| */ + else if (a1 >= bp[1]) /* borrow: |c| > |b| */ { MPFR_SET_OPPOSITE_SIGN (a, b); /* a = b-c mod 2^(2*GMP_NUMB_BITS) */ - ap[1] = -ap[1] - (ap[0] != 0); - ap[0] = -ap[0]; + a1 = -a1 - (a0 != 0); + a0 = -a0; } else /* bp[0] > cp[0] */ MPFR_SET_SAME_SIGN (a, b); - if (ap[1] == 0) + if (a1 == 0) { - ap[1] = ap[0]; - ap[0] = 0; + a1 = a0; + a0 = 0; bx -= GMP_NUMB_BITS; } - /* now ap[1] != 0 */ - MPFR_ASSERTD(ap[1] != 0); - count_leading_zeros (cnt, ap[1]); + /* now a1 != 0 */ + MPFR_ASSERTD(a1 != 0); + count_leading_zeros (cnt, a1); if (cnt) { - ap[1] = (ap[1] << cnt) | (ap[0] >> (GMP_NUMB_BITS - cnt)); - ap[0] = ap[0] << cnt; + a1 = (a1 << cnt) | (a0 >> (GMP_NUMB_BITS - cnt)); + a0 = a0 << cnt; bx -= cnt; } rb = sb = 0; @@ -361,86 +359,85 @@ mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, { MPFR_SET_SAME_SIGN (a, b); BGreater2: - b0 = bp[0]; - b1 = bp[1]; d = (mpfr_uexp_t) bx - cx; sh = 2 * 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 */ - ap[0] = b0 - ((cp[0] >> d) | (cp[1] << (GMP_NUMB_BITS - d))); - ap[1] = b1 - (cp[1] >> d) - (ap[0] > b0); + a0 = bp[0] - ((cp[0] >> d) | (cp[1] << (GMP_NUMB_BITS - d))); + a1 = bp[1] - (cp[1] >> d) - (a0 > bp[0]); if (sb) { - ap[1] -= (ap[0] == 0); - ap[0] --; - /* a = ap[1],ap[0] cannot become zero here, since: - a) if d >= 2, then ap[1] >= 2^(w-1) - (2^(w-2)-1) with - w = GMP_NUMB_BITS, thus ap[1] - 1 >= 2^(w-2), + 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(ap[1] > 0 || ap[0] > 0); + MPFR_ASSERTD(a1 > 0 || a0 > 0); sb = -sb; /* 2^GMP_NUMB_BITS - sb */ } - if (ap[1] == 0) + if (a1 == 0) { - ap[1] = ap[0]; - ap[0] = sb; - sb = 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 ap[1] != 0 */ - MPFR_ASSERTD(ap[1] != 0); - count_leading_zeros (cnt, ap[1]); + /* now a1 != 0 */ + MPFR_ASSERTD(a1 != 0); + count_leading_zeros (cnt, a1); if (cnt) { - ap[1] = (ap[1] << cnt) | (ap[0] >> (GMP_NUMB_BITS - cnt)); - ap[0] = (ap[0] << cnt) | (sb >> (GMP_NUMB_BITS - cnt)); + a1 = (a1 << cnt) | (a0 >> (GMP_NUMB_BITS - cnt)); + a0 = (a0 << cnt) | (sb >> (GMP_NUMB_BITS - cnt)); sb <<= cnt; bx -= cnt; } /* sh > 0 since p < 2*GMP_NUMB_BITS */ MPFR_ASSERTD(sh > 0); - rb = ap[0] & (MPFR_LIMB_ONE << (sh - 1)); - sb |= (ap[0] & mask) ^ rb; - ap[0] = ap[0] & ~mask; + rb = a0 & (MPFR_LIMB_ONE << (sh - 1)); + sb |= (a0 & mask) ^ rb; + a0 = 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 ap[0] below */ + significant bit of a0 below */ sb = (d == GMP_NUMB_BITS) ? cp[0] : (cp[1] << (2*GMP_NUMB_BITS - d)) | (cp[0] != 0); - ap[0] = b0 - (cp[1] >> (d - GMP_NUMB_BITS)) - (sb != 0); - ap[1] = b1 - (ap[0] >= b0); /* since cp[1] has its most significant bit + a0 = bp[0] - (cp[1] >> (d - GMP_NUMB_BITS)) - (sb != 0); + a1 = bp[1] - (a0 >= bp[0]); /* since cp[1] has its most significant bit set, and d-GMP_NUMB_BITS < GMP_NUMB_BITS, - cp[1] >> (d - GMP_NUMB_BITS) is always - non-zero, thus a borrow occurs iff - ap[0] >= b0 */ + cp[1] >> (d - GMP_NUMB_BITS) is non-zero, + thus a borrow occurs iff a0 >= b0 */ sb = -sb; /* since has its most significant bit set, we can have an exponent decrease of at most one */ - if (ap[1] < MPFR_LIMB_HIGHBIT) + if (a1 < MPFR_LIMB_HIGHBIT) { - ap[1] = (ap[1] << 1) | (ap[0] >> (GMP_NUMB_BITS - 1)); - ap[0] = (ap[0] << 1) | (sb >> (GMP_NUMB_BITS - 1)); + a1 = (a1 << 1) | (a0 >> (GMP_NUMB_BITS - 1)); + a0 = (a0 << 1) | (sb >> (GMP_NUMB_BITS - 1)); bx --; } - rb = ap[0] & (MPFR_LIMB_ONE << (sh - 1)); - sb |= (ap[0] & mask) ^ rb; - ap[0] = ap[0] & ~mask; + rb = a0 & (MPFR_LIMB_ONE << (sh - 1)); + sb |= (a0 & mask) ^ rb; + a0 = 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. */ - ap[0] = b0 - (MPFR_LIMB_ONE << sh); - ap[1] = b1 - (ap[0] > b0); - if (ap[1] < MPFR_LIMB_HIGHBIT) + a0 = bp[0] - (MPFR_LIMB_ONE << sh); + a1 = bp[1] - (a0 > bp[0]); + if (a1 < MPFR_LIMB_HIGHBIT) { /* necessarily we had b = 1000...000 */ - ap[0] = ~mask; - ap[1] = ~0; + a0 = ~mask; + a1 = ~0; bx --; } rb = sb = 1; @@ -456,6 +453,9 @@ mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, goto BGreater2; } + ap[1] = a1; + ap[0] = a0; + /* now perform rounding */ /* Warning: MPFR considers underflow *after* rounding with an unbounded |