--- /dev/null
+/* mpfr_add1 -- internal function to perform a "real" addition
+
+Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
+Contributed by the Arenaire and Cacao projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to
+the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
+MA 02110-1301, USA. */
+
+#include "mpfr-impl.h"
+
+/* compute sign(b) * (|b| + |c|) */
+int
+mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
+{
+ mp_limb_t *ap, *bp, *cp;
+ mp_prec_t aq, bq, cq, aq2;
+ mp_size_t an, bn, cn;
+ mp_exp_t difw, exp;
+ int sh, rb, fb, inex;
+ mp_exp_unsigned_t diff_exp;
+ MPFR_TMP_DECL(marker);
+
+ MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c));
+
+ MPFR_TMP_MARK(marker);
+
+ aq = MPFR_PREC(a);
+ bq = MPFR_PREC(b);
+ cq = MPFR_PREC(c);
+
+ an = (aq-1)/BITS_PER_MP_LIMB+1; /* number of limbs of a */
+ aq2 = (mp_prec_t) an * BITS_PER_MP_LIMB;
+ sh = aq2 - aq; /* non-significant bits in low limb */
+
+ bn = (bq-1)/BITS_PER_MP_LIMB+1; /* number of limbs of b */
+ cn = (cq-1)/BITS_PER_MP_LIMB+1; /* number of limbs of c */
+
+ ap = MPFR_MANT(a);
+ bp = MPFR_MANT(b);
+ cp = MPFR_MANT(c);
+
+ if (MPFR_UNLIKELY(ap == bp))
+ {
+ bp = (mp_ptr) MPFR_TMP_ALLOC (bn * BYTES_PER_MP_LIMB);
+ MPN_COPY (bp, ap, bn);
+ if (ap == cp)
+ { cp = bp; }
+ }
+ else if (MPFR_UNLIKELY(ap == cp))
+ {
+ cp = (mp_ptr) MPFR_TMP_ALLOC (cn * BYTES_PER_MP_LIMB);
+ MPN_COPY(cp, ap, cn);
+ }
+
+ exp = MPFR_GET_EXP (b);
+ MPFR_SET_SAME_SIGN(a, b);
+ diff_exp = (mp_exp_unsigned_t) exp - MPFR_GET_EXP(c);
+
+ /*
+ * 1. Compute the significant part A', the non-significant bits of A
+ * are taken into account.
+ *
+ * 2. Perform the rounding. At each iteration, we remember:
+ * _ r = rounding bit
+ * _ f = following bits (same value)
+ * where the result has the form: [number A]rfff...fff + a remaining
+ * value in the interval [0,2) ulp. We consider the most significant
+ * bits of the remaining value to update the result; a possible carry
+ * is immediately taken into account and A is updated accordingly. As
+ * soon as the bits f don't have the same value, A can be rounded.
+ * Variables:
+ * _ rb = rounding bit (0 or 1).
+ * _ fb = following bits (0 or 1), then sticky bit.
+ * If fb == 0, the only thing that can change is the sticky bit.
+ */
+
+ rb = fb = -1; /* means: not initialized */
+
+ if (MPFR_UNLIKELY(aq2 <= diff_exp))
+ { /* c does not overlap with a' */
+ if (MPFR_UNLIKELY(an > bn))
+ { /* a has more limbs than b */
+ /* copy b to the most significant limbs of a */
+ MPN_COPY(ap + (an - bn), bp, bn);
+ /* zero the least significant limbs of a */
+ MPN_ZERO(ap, an - bn);
+ }
+ else /* an <= bn */
+ {
+ /* copy the most significant limbs of b to a */
+ MPN_COPY(ap, bp + (bn - an), an);
+ }
+ }
+ else /* aq2 > diff_exp */
+ { /* c overlaps with a' */
+ mp_limb_t *a2p;
+ mp_limb_t cc;
+ mp_prec_t dif;
+ mp_size_t difn, k;
+ int shift;
+
+ /* copy c (shifted) into a */
+
+ dif = aq2 - diff_exp;
+ /* dif is the number of bits of c which overlap with a' */
+
+ difn = (dif-1)/BITS_PER_MP_LIMB + 1;
+ /* only the highest difn limbs from c have to be considered */
+ if (MPFR_UNLIKELY(difn > cn))
+ {
+ /* c doesn't have enough limbs; take into account the virtual
+ zero limbs now by zeroing the least significant limbs of a' */
+ MPFR_ASSERTD(difn - cn <= an);
+ MPN_ZERO(ap, difn - cn);
+ difn = cn;
+ }
+ k = diff_exp / BITS_PER_MP_LIMB;
+
+ /* zero the most significant k limbs of a */
+ a2p = ap + (an - k);
+ MPN_ZERO(a2p, k);
+
+ shift = diff_exp % BITS_PER_MP_LIMB;
+
+ if (MPFR_LIKELY(shift))
+ {
+ MPFR_ASSERTD(a2p - difn >= ap);
+ cc = mpn_rshift(a2p - difn, cp + (cn - difn), difn, shift);
+ if (MPFR_UNLIKELY(a2p - difn > ap))
+ *(a2p - difn - 1) = cc;
+ }
+ else
+ MPN_COPY(a2p - difn, cp + (cn - difn), difn);
+
+ /* add b to a */
+ cc = MPFR_UNLIKELY(an > bn)
+ ? mpn_add_n(ap + (an - bn), ap + (an - bn), bp, bn)
+ : mpn_add_n(ap, ap, bp + (bn - an), an);
+
+ if (MPFR_UNLIKELY(cc)) /* carry */
+ {
+ if (MPFR_UNLIKELY(exp == __gmpfr_emax))
+ {
+ inex = mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
+ goto end_of_add;
+ }
+ exp++;
+ rb = (ap[0] >> sh) & 1; /* LSB(a) --> rounding bit after the shift */
+ if (MPFR_LIKELY(sh))
+ {
+ mp_limb_t mask, bb;
+
+ mask = MPFR_LIMB_MASK (sh);
+ bb = ap[0] & mask;
+ ap[0] &= (~mask) << 1;
+ if (bb == 0)
+ fb = 0;
+ else if (bb == mask)
+ fb = 1;
+ }
+ mpn_rshift(ap, ap, an, 1);
+ ap[an-1] += MPFR_LIMB_HIGHBIT;
+ if (sh && fb < 0)
+ goto rounding;
+ } /* cc */
+ } /* aq2 > diff_exp */
+
+ /* non-significant bits of a */
+ if (MPFR_LIKELY(rb < 0 && sh))
+ {
+ mp_limb_t mask, bb;
+
+ mask = MPFR_LIMB_MASK (sh);
+ bb = ap[0] & mask;
+ ap[0] &= ~mask;
+ rb = bb >> (sh - 1);
+ if (MPFR_LIKELY(sh > 1))
+ {
+ mask >>= 1;
+ bb &= mask;
+ if (bb == 0)
+ fb = 0;
+ else if (bb == mask)
+ fb = 1;
+ else
+ goto rounding;
+ }
+ }
+
+ /* determine rounding and sticky bits (and possible carry) */
+
+ difw = (mp_exp_t) an - (mp_exp_t) (diff_exp / BITS_PER_MP_LIMB);
+ /* difw is the number of limbs from b (regarded as having an infinite
+ precision) that have already been combined with c; -n if the next
+ n limbs from b won't be combined with c. */
+
+ if (MPFR_UNLIKELY(bn > an))
+ { /* there are still limbs from b that haven't been taken into account */
+ mp_size_t bk;
+
+ if (fb == 0 && difw <= 0)
+ {
+ fb = 1; /* c hasn't been taken into account ==> sticky bit != 0 */
+ goto rounding;
+ }
+
+ bk = bn - an; /* index of lowest considered limb from b, > 0 */
+ while (difw < 0)
+ { /* ulp(next limb from b) > msb(c) */
+ mp_limb_t bb;
+
+ bb = bp[--bk];
+
+ MPFR_ASSERTD(fb != 0);
+ if (fb > 0)
+ {
+ if (bb != MP_LIMB_T_MAX)
+ {
+ fb = 1; /* c hasn't been taken into account
+ ==> sticky bit != 0 */
+ goto rounding;
+ }
+ }
+ else /* fb not initialized yet */
+ {
+ if (rb < 0) /* rb not initialized yet */
+ {
+ rb = bb >> (BITS_PER_MP_LIMB - 1);
+ bb |= MPFR_LIMB_HIGHBIT;
+ }
+ fb = 1;
+ if (bb != MP_LIMB_T_MAX)
+ goto rounding;
+ }
+
+ if (bk == 0)
+ { /* b has entirely been read */
+ fb = 1; /* c hasn't been taken into account
+ ==> sticky bit != 0 */
+ goto rounding;
+ }
+
+ difw++;
+ } /* while */
+ MPFR_ASSERTD(bk > 0 && difw >= 0);
+
+ if (difw <= cn)
+ {
+ mp_size_t ck;
+ mp_limb_t cprev;
+ int difs;
+
+ ck = cn - difw;
+ difs = diff_exp % BITS_PER_MP_LIMB;
+
+ if (difs == 0 && ck == 0)
+ goto c_read;
+
+ cprev = ck == cn ? 0 : cp[ck];
+
+ if (fb < 0)
+ {
+ mp_limb_t bb, cc;
+
+ if (difs)
+ {
+ cc = cprev << (BITS_PER_MP_LIMB - difs);
+ if (--ck >= 0)
+ {
+ cprev = cp[ck];
+ cc += cprev >> difs;
+ }
+ }
+ else
+ cc = cp[--ck];
+
+ bb = bp[--bk] + cc;
+
+ if (bb < cc /* carry */
+ && (rb < 0 || (rb ^= 1) == 0)
+ && mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh))
+ {
+ if (exp == __gmpfr_emax)
+ {
+ inex = mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
+ goto end_of_add;
+ }
+ exp++;
+ ap[an-1] = MPFR_LIMB_HIGHBIT;
+ rb = 0;
+ }
+
+ if (rb < 0) /* rb not initialized yet */
+ {
+ rb = bb >> (BITS_PER_MP_LIMB - 1);
+ bb <<= 1;
+ bb |= bb >> (BITS_PER_MP_LIMB - 1);
+ }
+
+ fb = bb != 0;
+ if (fb && bb != MP_LIMB_T_MAX)
+ goto rounding;
+ } /* fb < 0 */
+
+ while (bk > 0)
+ {
+ mp_limb_t bb, cc;
+
+ if (difs)
+ {
+ if (ck < 0)
+ goto c_read;
+ cc = cprev << (BITS_PER_MP_LIMB - difs);
+ if (--ck >= 0)
+ {
+ cprev = cp[ck];
+ cc += cprev >> difs;
+ }
+ }
+ else
+ {
+ if (ck == 0)
+ goto c_read;
+ cc = cp[--ck];
+ }
+
+ bb = bp[--bk] + cc;
+ if (bb < cc) /* carry */
+ {
+ fb ^= 1;
+ if (fb)
+ goto rounding;
+ rb ^= 1;
+ if (rb == 0 && mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh))
+ {
+ if (MPFR_UNLIKELY(exp == __gmpfr_emax))
+ {
+ inex = mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
+ goto end_of_add;
+ }
+ exp++;
+ ap[an-1] = MPFR_LIMB_HIGHBIT;
+ }
+ } /* bb < cc */
+
+ if (!fb && bb != 0)
+ {
+ fb = 1;
+ goto rounding;
+ }
+ if (fb && bb != MP_LIMB_T_MAX)
+ goto rounding;
+ } /* while */
+
+ /* b has entirely been read */
+
+ if (fb || ck < 0)
+ goto rounding;
+ if (difs && cprev << (BITS_PER_MP_LIMB - difs))
+ {
+ fb = 1;
+ goto rounding;
+ }
+ while (ck)
+ {
+ if (cp[--ck])
+ {
+ fb = 1;
+ goto rounding;
+ }
+ } /* while */
+ } /* difw <= cn */
+ else
+ { /* c has entirely been read */
+ c_read:
+ if (fb < 0) /* fb not initialized yet */
+ {
+ mp_limb_t bb;
+
+ MPFR_ASSERTD(bk > 0);
+ bb = bp[--bk];
+ if (rb < 0) /* rb not initialized yet */
+ {
+ rb = bb >> (BITS_PER_MP_LIMB - 1);
+ bb &= ~MPFR_LIMB_HIGHBIT;
+ }
+ fb = bb != 0;
+ } /* fb < 0 */
+ if (fb)
+ goto rounding;
+ while (bk)
+ {
+ if (bp[--bk])
+ {
+ fb = 1;
+ goto rounding;
+ }
+ } /* while */
+ } /* difw > cn */
+ } /* bn > an */
+ else if (fb != 1) /* if fb == 1, the sticky bit is 1 (no possible carry) */
+ { /* b has entirely been read */
+ if (difw > cn)
+ { /* c has entirely been read */
+ if (rb < 0)
+ rb = 0;
+ fb = 0;
+ }
+ else if (diff_exp > aq2)
+ { /* b is followed by at least a zero bit, then by c */
+ if (rb < 0)
+ rb = 0;
+ fb = 1;
+ }
+ else
+ {
+ mp_size_t ck;
+ int difs;
+
+ MPFR_ASSERTD(difw >= 0 && cn >= difw);
+ ck = cn - difw;
+ difs = diff_exp % BITS_PER_MP_LIMB;
+
+ if (difs == 0 && ck == 0)
+ { /* c has entirely been read */
+ if (rb < 0)
+ rb = 0;
+ fb = 0;
+ }
+ else
+ {
+ mp_limb_t cc;
+
+ cc = difs ? (MPFR_ASSERTD(ck < cn),
+ cp[ck] << (BITS_PER_MP_LIMB - difs)) : cp[--ck];
+ if (rb < 0)
+ {
+ rb = cc >> (BITS_PER_MP_LIMB - 1);
+ cc &= ~MPFR_LIMB_HIGHBIT;
+ }
+ while (cc == 0)
+ {
+ if (ck == 0)
+ {
+ fb = 0;
+ goto rounding;
+ }
+ cc = cp[--ck];
+ } /* while */
+ fb = 1;
+ }
+ }
+ } /* fb != 1 */
+
+ rounding:
+ switch (rnd_mode)
+ {
+ case GMP_RNDN:
+ if (fb == 0)
+ {
+ if (rb == 0)
+ {
+ inex = 0;
+ goto set_exponent;
+ }
+ /* round to even */
+ if (ap[0] & (MPFR_LIMB_ONE << sh))
+ goto rndn_away;
+ else
+ goto rndn_zero;
+ }
+ if (rb == 0)
+ {
+ rndn_zero:
+ inex = MPFR_IS_NEG(a) ? 1 : -1;
+ goto set_exponent;
+ }
+ else
+ {
+ rndn_away:
+ inex = MPFR_IS_POS(a) ? 1 : -1;
+ goto add_one_ulp;
+ }
+
+ case GMP_RNDZ:
+ inex = rb || fb ? (MPFR_IS_NEG(a) ? 1 : -1) : 0;
+ goto set_exponent;
+
+ case GMP_RNDU:
+ inex = rb || fb;
+ if (inex && MPFR_IS_POS(a))
+ goto add_one_ulp;
+ else
+ goto set_exponent;
+
+ case GMP_RNDD:
+ inex = - (rb || fb);
+ if (inex && MPFR_IS_NEG(a))
+ goto add_one_ulp;
+ else
+ goto set_exponent;
+
+ default:
+ MPFR_ASSERTN(0);
+ inex = 0;
+ goto set_exponent;
+ }
+
+ add_one_ulp: /* add one unit in last place to a */
+ if (MPFR_UNLIKELY(mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh)))
+ {
+ if (MPFR_UNLIKELY(exp == __gmpfr_emax))
+ {
+ inex = mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
+ goto end_of_add;
+ }
+ exp++;
+ ap[an-1] = MPFR_LIMB_HIGHBIT;
+ }
+
+ set_exponent:
+ MPFR_SET_EXP (a, exp);
+
+ end_of_add:
+ MPFR_TMP_FREE(marker);
+ MPFR_RET (inex);
+}
projects / msp430-gcc.git / blobdiff