diff options
-rw-r--r-- | src/add1sp.c | 40 |
1 files changed, 22 insertions, 18 deletions
diff --git a/src/add1sp.c b/src/add1sp.c index 55affc356..176fe403a 100644 --- a/src/add1sp.c +++ b/src/add1sp.c @@ -88,7 +88,7 @@ int mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) mpfr_clears (tmpa, tmpb, tmpc, (mpfr_ptr) 0); return inexact; } -# define mpfr_add1sp mpfr_add1sp2 +# define mpfr_add1sp mpfr_add1sp_ref #endif /* MPFR_WANT_ASSERT >= 2 */ /* Debugging support */ @@ -230,13 +230,13 @@ mpfr_add1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, mp_limb_t mask; mpfr_uexp_t d; - MPFR_ASSERTD(p < GMP_NUMB_BITS); + MPFR_ASSERTD(GMP_NUMB_BITS < p && p < 2 * GMP_NUMB_BITS); if (bx == cx) { /* since bp[1], cp[1] >= MPFR_LIMB_HIGHBIT, a carry always occurs */ a0 = bp[0] + cp[0]; - a1 = bp[1] + cp[1] + (sb < bp[0]); + a1 = bp[1] + cp[1] + (a0 < bp[0]); ap[0] = (a0 >> 1) | (a1 << (GMP_NUMB_BITS - 1)); ap[1] = MPFR_LIMB_HIGHBIT | (a1 >> 1); bx ++; @@ -246,7 +246,7 @@ mpfr_add1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, } else if (bx > cx) { - BGreater1: + BGreater2: d = (mpfr_uexp_t) bx - cx; mask = MPFR_LIMB_MASK(sh); if (d < GMP_NUMB_BITS) /* 0 < d < GMP_NUMB_BITS */ @@ -254,30 +254,34 @@ mpfr_add1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, a0 = bp[0]; a1 = bp[1]; sb = cp[0] << (GMP_NUMB_BITS - d); /* bits from cp[-1] after shift */ - ap[0] = a0 + ((cp[1] << (GMP_NUMB_BITS - d)) | (cp[0] >> d)); + a0 += (cp[1] << (GMP_NUMB_BITS - d)) | (cp[0] >> d); ap[1] = a1 + (cp[1] >> d); - if (ap[0] < a0) /* carry in low word */ + if (a0 < bp[0]) /* carry in low word */ ap[1] ++; if (ap[1] < a1) /* carry in high word */ { - sb |= ap[0] & 1; + exponent_shift: + sb |= a0 & 1; /* shift a by 1 */ - ap[0] = (ap[1] << (GMP_NUMB_BITS - 1)) | (ap[0] >> 1); + a0 = (ap[1] << (GMP_NUMB_BITS - 1)) | (a0 >> 1); ap[1] = MPFR_LIMB_HIGHBIT | (ap[1] >> 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; + ap[0] = a0 & ~mask; } else if (d < 2*GMP_NUMB_BITS) /* GMP_NUMB_BITS <= d < 2*GMP_NUMB_BITS */ { - sb = cp[0] | (cp[1] << (2*GMP_NUMB_BITS-d)); - ap[1] = bp[1]; - ap[0] = bp[0] + (cp[1] >> (d - GMP_NUMB_BITS)); - rb = ap[0] & (MPFR_LIMB_ONE << (sh - 1)); - sb |= (ap[0] & mask) ^ rb; - ap[0] = ap[0] & ~mask; + sb = (d == GMP_NUMB_BITS) ? cp[0] + : cp[0] | (cp[1] << (2*GMP_NUMB_BITS-d)); + a0 = bp[0] + (cp[1] >> (d - GMP_NUMB_BITS)); + ap[1] = bp[1] + (a0 < bp[0]); + if (ap[1] == 0) + goto exponent_shift; + rb = a0 & (MPFR_LIMB_ONE << (sh - 1)); + sb |= (a0 & mask) ^ rb; + ap[0] = a0 & ~mask; } else /* d >= 2*GMP_NUMB_BITS */ { @@ -293,7 +297,7 @@ mpfr_add1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, mp_limb_t *tp; tx = bx; bx = cx; cx = tx; tp = bp; bp = cp; cp = tp; - goto BGreater1; + goto BGreater2; } /* now perform rounding */ |