diff options
author | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2017-02-07 18:42:53 +0000 |
---|---|---|
committer | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2017-02-07 18:42:53 +0000 |
commit | dcb41575f892c942fa3004fff746468a23311c0c (patch) | |
tree | a1f16b5fea6942c0510a81caf74df40f7e95d408 /src | |
parent | 3a8ae11cb261aba272814572acf3cb0230ad9872 (diff) | |
download | mpfr-dcb41575f892c942fa3004fff746468a23311c0c.tar.gz |
[src/add1sp.c] added special code for p=GMP_NUMB_BITS
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@11271 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'src')
-rw-r--r-- | src/add1sp.c | 110 |
1 files changed, 110 insertions, 0 deletions
diff --git a/src/add1sp.c b/src/add1sp.c index ad64d379a..936f11b41 100644 --- a/src/add1sp.c +++ b/src/add1sp.c @@ -226,6 +226,111 @@ mpfr_add1sp1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, } } +/* same as mpfr_add1sp, but for p = GMP_NUMB_BITS */ +static int +mpfr_add1sp1n (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); + mp_limb_t rb; /* round bit */ + mp_limb_t sb; /* sticky bit */ + mp_limb_t a0; /* to store the result */ + 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]; + rb = a0 & MPFR_LIMB_ONE; + ap[0] = MPFR_LIMB_HIGHBIT | (a0 >> 1); + bx ++; + sb = 0; /* since b + c fits on p+1 bits, the sticky bit is zero */ + } + else if (bx > cx) + { + BGreater1: + d = (mpfr_uexp_t) bx - cx; + if (d < GMP_NUMB_BITS) /* 1 <= d < GMP_NUMB_BITS */ + { + a0 = bp[0] + (cp[0] >> d); + sb = cp[0] << (GMP_NUMB_BITS - d); /* bits from cp[-1] after shift */ + if (a0 < bp[0]) /* carry */ + { + ap[0] = MPFR_LIMB_HIGHBIT | (a0 >> 1); + rb = a0 & 1; + bx ++; + } + else /* no carry */ + { + ap[0] = a0; + rb = sb & MPFR_LIMB_HIGHBIT; + sb &= ~MPFR_LIMB_HIGHBIT; + } + } + else /* d >= GMP_NUMB_BITS */ + { + sb = d != GMP_NUMB_BITS || cp[0] != MPFR_LIMB_HIGHBIT; + ap[0] = bp[0]; + rb = d == GMP_NUMB_BITS; + } + } + else /* 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; + goto BGreater1; + } + + /* Note: we could keep the output significand in a0 for the rounding, + and only store it in ap[0] at the very end, but this seems slower + on average (but better for the worst case). */ + + /* now perform rounding */ + if (MPFR_UNLIKELY(bx > __gmpfr_emax)) + return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); + + MPFR_SET_EXP (a, bx); + if (rb == 0 && sb == 0) + MPFR_RET(0); + else if (rnd_mode == MPFR_RNDN) + { + /* the condition below should be rb == 0 || (rb != 0 && ...), but this + is equivalent to rb == 0 || (...) */ + 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; + /* no need to have MPFR_LIKELY here, since we are in a rare branch */ + if (bx + 1 <= __gmpfr_emax) + MPFR_SET_EXP (a, bx + 1); + else /* overflow */ + return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); + } + MPFR_RET(MPFR_SIGN(a)); + } +} + /* same as mpfr_add1sp, but for GMP_NUMB_BITS < p < 2*GMP_NUMB_BITS */ static int mpfr_add1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, @@ -540,6 +645,11 @@ mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) if (GMP_NUMB_BITS < p && p < 2 * GMP_NUMB_BITS) return mpfr_add1sp2 (a, b, c, rnd_mode, p); + /* we put this test after the mpfr_add1sp2() call, to avoid slowing down + the more frequent case GMP_NUMB_BITS < p < 2 * GMP_NUMB_BITS */ + if (p == GMP_NUMB_BITS) + return mpfr_add1sp1n (a, b, c, rnd_mode); + if (2 * GMP_NUMB_BITS < p && p < 3 * GMP_NUMB_BITS) return mpfr_add1sp3 (a, b, c, rnd_mode, p); |