Index: gcc/mpfr/add1.c |
diff --git a/gcc/mpfr/add1.c b/gcc/mpfr/add1.c |
deleted file mode 100644 |
index f5f30308f98972dfa3f5a8bef5c43d6511218e61..0000000000000000000000000000000000000000 |
--- a/gcc/mpfr/add1.c |
+++ /dev/null |
@@ -1,541 +0,0 @@ |
-/* 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); |
-} |