diff options
Diffstat (limited to 'src/add1sp.c')
| -rw-r--r-- | src/add1sp.c | 384 |
1 files changed, 384 insertions, 0 deletions
diff --git a/src/add1sp.c b/src/add1sp.c new file mode 100644 index 000000000..49e691c90 --- /dev/null +++ b/src/add1sp.c @@ -0,0 +1,384 @@ +/* mpfr_add1sp -- internal function to perform a "real" addition + All the op must have the same precision + +Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010 Free Software Foundation, Inc. +Contributed by the Arenaire and Caramel 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 3 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.LESSER. If not, see +http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., +51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ + +#define MPFR_NEED_LONGLONG_H +#include "mpfr-impl.h" + +/* Check if we have to check the result of mpfr_add1sp with mpfr_add1 */ +#ifdef WANT_ASSERT +# if WANT_ASSERT >= 2 + +int mpfr_add1sp2 (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t); +int mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) +{ + mpfr_t tmpa, tmpb, tmpc; + int inexb, inexc, inexact, inexact2; + + mpfr_init2 (tmpa, MPFR_PREC (a)); + mpfr_init2 (tmpb, MPFR_PREC (b)); + mpfr_init2 (tmpc, MPFR_PREC (c)); + + inexb = mpfr_set (tmpb, b, MPFR_RNDN); + MPFR_ASSERTN (inexb == 0); + + inexc = mpfr_set (tmpc, c, MPFR_RNDN); + MPFR_ASSERTN (inexc == 0); + + inexact2 = mpfr_add1 (tmpa, tmpb, tmpc, rnd_mode); + inexact = mpfr_add1sp2 (a, b, c, rnd_mode); + + if (mpfr_cmp (tmpa, a) || inexact != inexact2) + { + fprintf (stderr, "add1 & add1sp return different values for %s\n" + "Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nB = ", + mpfr_print_rnd_mode (rnd_mode), + MPFR_PREC (a), MPFR_PREC (b), MPFR_PREC (c)); + mpfr_fprint_binary (stderr, tmpb); + fprintf (stderr, "\nC = "); + mpfr_fprint_binary (stderr, tmpc); + fprintf (stderr, "\n\nadd1 : "); + mpfr_fprint_binary (stderr, tmpa); + fprintf (stderr, "\nadd1sp: "); + mpfr_fprint_binary (stderr, a); + fprintf (stderr, "\nInexact sp = %d | Inexact = %d\n", + inexact, inexact2); + MPFR_ASSERTN (0); + } + mpfr_clears (tmpa, tmpb, tmpc, (mpfr_ptr) 0); + return inexact; +} +# define mpfr_add1sp mpfr_add1sp2 +# endif +#endif + +/* Debugging support */ +#ifdef DEBUG +# undef DEBUG +# define DEBUG(x) (x) +#else +# define DEBUG(x) /**/ +#endif + +/* compute sign(b) * (|b| + |c|) + Returns 0 iff result is exact, + a negative value when the result is less than the exact value, + a positive value otherwise. */ +int +mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) +{ + mpfr_uexp_t d; + mpfr_prec_t p; + unsigned int sh; + mp_size_t n; + mp_limb_t *ap, *cp; + mpfr_exp_t bx; + mp_limb_t limb; + int inexact; + MPFR_TMP_DECL(marker); + + MPFR_TMP_MARK(marker); + + MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c)); + MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c)); + MPFR_ASSERTD(MPFR_GET_EXP(b) >= MPFR_GET_EXP(c)); + + /* Read prec and num of limbs */ + p = MPFR_PREC(b); + n = (p+GMP_NUMB_BITS-1)/GMP_NUMB_BITS; + MPFR_UNSIGNED_MINUS_MODULO(sh, p); + bx = MPFR_GET_EXP(b); + d = (mpfr_uexp_t) (bx - MPFR_GET_EXP(c)); + + DEBUG (printf ("New add1sp with diff=%lu\n", (unsigned long) d)); + + if (MPFR_UNLIKELY(d == 0)) + { + /* d==0 */ + DEBUG( mpfr_print_mant_binary("C= ", MPFR_MANT(c), p) ); + DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) ); + bx++; /* exp + 1 */ + ap = MPFR_MANT(a); + limb = mpn_add_n(ap, MPFR_MANT(b), MPFR_MANT(c), n); + DEBUG( mpfr_print_mant_binary("A= ", ap, p) ); + MPFR_ASSERTD(limb != 0); /* There must be a carry */ + limb = ap[0]; /* Get LSB (In fact, LSW) */ + mpn_rshift(ap, ap, n, 1); /* Shift mantissa A */ + ap[n-1] |= MPFR_LIMB_HIGHBIT; /* Set MSB */ + ap[0] &= ~MPFR_LIMB_MASK(sh); /* Clear LSB bit */ + if (MPFR_LIKELY((limb&(MPFR_LIMB_ONE<<sh)) == 0)) /* Check exact case */ + { inexact = 0; goto set_exponent; } + /* Zero: Truncate + Nearest: Even Rule => truncate or add 1 + Away: Add 1 */ + if (MPFR_LIKELY(rnd_mode==MPFR_RNDN)) + { + if (MPFR_LIKELY((ap[0]&(MPFR_LIMB_ONE<<sh))==0)) + { inexact = -1; goto set_exponent; } + else + goto add_one_ulp; + } + MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b)); + if (rnd_mode==MPFR_RNDZ) + { inexact = -1; goto set_exponent; } + else + goto add_one_ulp; + } + else if (MPFR_UNLIKELY (d >= p)) + { + if (MPFR_LIKELY (d > p)) + { + /* d > p : Copy B in A */ + /* Away: Add 1 + Nearest: Trunc + Zero: Trunc */ + if (MPFR_LIKELY (rnd_mode==MPFR_RNDN + || MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b)))) + { + copy_set_exponent: + ap = MPFR_MANT (a); + MPN_COPY (ap, MPFR_MANT(b), n); + inexact = -1; + goto set_exponent; + } + else + { + copy_add_one_ulp: + ap = MPFR_MANT(a); + MPN_COPY (ap, MPFR_MANT(b), n); + goto add_one_ulp; + } + } + else + { + /* d==p : Copy B in A */ + /* Away: Add 1 + Nearest: Even Rule if C is a power of 2, else Add 1 + Zero: Trunc */ + if (MPFR_LIKELY(rnd_mode==MPFR_RNDN)) + { + /* Check if C was a power of 2 */ + cp = MPFR_MANT(c); + if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT)) + { + mp_size_t k = n-1; + do { + k--; + } while (k>=0 && cp[k]==0); + if (MPFR_UNLIKELY(k<0)) + /* Power of 2: Even rule */ + if ((MPFR_MANT (b)[0]&(MPFR_LIMB_ONE<<sh))==0) + goto copy_set_exponent; + } + /* Not a Power of 2 */ + goto copy_add_one_ulp; + } + else if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b))) + goto copy_set_exponent; + else + goto copy_add_one_ulp; + } + } + else + { + mp_limb_t mask; + mp_limb_t bcp, bcp1; /* Cp and C'p+1 */ + + /* General case: 1 <= d < p */ + cp = (mp_limb_t*) MPFR_TMP_ALLOC(n * BYTES_PER_MP_LIMB); + + /* Shift c in temporary allocated place */ + { + mpfr_uexp_t dm; + mp_size_t m; + + dm = d % GMP_NUMB_BITS; + m = d / GMP_NUMB_BITS; + if (MPFR_UNLIKELY(dm == 0)) + { + /* dm = 0 and m > 0: Just copy */ + MPFR_ASSERTD(m!=0); + MPN_COPY(cp, MPFR_MANT(c)+m, n-m); + MPN_ZERO(cp+n-m, m); + } + else if (MPFR_LIKELY(m == 0)) + { + /* dm >=1 and m == 0: just shift */ + MPFR_ASSERTD(dm >= 1); + mpn_rshift(cp, MPFR_MANT(c), n, dm); + } + else + { + /* dm > 0 and m > 0: shift and zero */ + mpn_rshift(cp, MPFR_MANT(c)+m, n-m, dm); + MPN_ZERO(cp+n-m, m); + } + } + + DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c), p) ); + DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) ); + DEBUG( mpfr_print_mant_binary("After ", cp, p) ); + + /* Compute bcp=Cp and bcp1=C'p+1 */ + if (MPFR_LIKELY (sh > 0)) + { + /* Try to compute them from C' rather than C */ + bcp = (cp[0] & (MPFR_LIMB_ONE<<(sh-1))) ; + if (MPFR_LIKELY(cp[0]&MPFR_LIMB_MASK(sh-1))) + bcp1 = 1; + else + { + /* We can't compute C'p+1 from C'. Compute it from C */ + /* Start from bit x=p-d+sh in mantissa C + (+sh since we have already looked sh bits in C'!) */ + mpfr_prec_t x = p-d+sh-1; + if (MPFR_LIKELY(x>p)) + /* We are already looked at all the bits of c, so C'p+1 = 0*/ + bcp1 = 0; + else + { + mp_limb_t *tp = MPFR_MANT(c); + mp_size_t kx = n-1 - (x / GMP_NUMB_BITS); + mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS); + DEBUG (printf ("(First) x=%lu Kx=%ld Sx=%lu\n", + (unsigned long) x, (long) kx, + (unsigned long) sx)); + /* Looks at the last bits of limb kx (if sx=0 does nothing)*/ + if (tp[kx] & MPFR_LIMB_MASK(sx)) + bcp1 = 1; + else + { + /*kx += (sx==0);*/ + /*If sx==0, tp[kx] hasn't been checked*/ + do { + kx--; + } while (kx>=0 && tp[kx]==0); + bcp1 = (kx >= 0); + } + } + } + } + else /* sh == 0 */ + { + /* Compute Cp and C'p+1 from C with sh=0 */ + mp_limb_t *tp = MPFR_MANT(c); + /* Start from bit x=p-d in mantissa C */ + mpfr_prec_t x = p-d; + mp_size_t kx = n-1 - (x / GMP_NUMB_BITS); + mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS); + MPFR_ASSERTD(p >= d); + bcp = tp[kx] & (MPFR_LIMB_ONE<<sx); + /* Looks at the last bits of limb kx (If sx=0, does nothing)*/ + if (tp[kx]&MPFR_LIMB_MASK(sx)) + bcp1 = 1; + else + { + do { + kx--; + } while (kx>=0 && tp[kx]==0); + bcp1 = (kx>=0); + } + } + DEBUG (printf("sh=%u Cp=%lu C'p+1=%lu\n", sh, + (unsigned long) bcp, (unsigned long) bcp1)); + + /* Clean shifted C' */ + mask = ~MPFR_LIMB_MASK(sh); + cp[0] &= mask; + + /* Add the mantissa c from b in a */ + ap = MPFR_MANT(a); + limb = mpn_add_n (ap, MPFR_MANT(b), cp, n); + DEBUG( mpfr_print_mant_binary("Add= ", ap, p) ); + + /* Check for overflow */ + if (MPFR_UNLIKELY (limb)) + { + limb = ap[0] & (MPFR_LIMB_ONE<<sh); /* Get LSB */ + mpn_rshift (ap, ap, n, 1); /* Shift mantissa*/ + bx++; /* Fix exponent */ + ap[n-1] |= MPFR_LIMB_HIGHBIT; /* Set MSB */ + ap[0] &= mask; /* Clear LSB bit */ + bcp1 |= bcp; /* Recompute C'p+1 */ + bcp = limb; /* Recompute Cp */ + DEBUG (printf ("(Overflow) Cp=%lu C'p+1=%lu\n", + (unsigned long) bcp, (unsigned long) bcp1)); + DEBUG (mpfr_print_mant_binary ("Add= ", ap, p)); + } + + /* Round: + Zero: Truncate but could be exact. + Away: Add 1 if Cp or C'p+1 !=0 + Nearest: Truncate but could be exact if Cp==0 + Add 1 if C'p+1 !=0, + Even rule else */ + if (MPFR_LIKELY(rnd_mode == MPFR_RNDN)) + { + if (MPFR_LIKELY(bcp == 0)) + { inexact = MPFR_LIKELY(bcp1) ? -1 : 0; goto set_exponent; } + else if (MPFR_UNLIKELY(bcp1==0) && (ap[0]&(MPFR_LIMB_ONE<<sh))==0) + { inexact = -1; goto set_exponent; } + else + goto add_one_ulp; + } + MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b)); + if (rnd_mode == MPFR_RNDZ) + { + inexact = MPFR_LIKELY(bcp || bcp1) ? -1 : 0; + goto set_exponent; + } + else + { + if (MPFR_UNLIKELY(bcp==0 && bcp1==0)) + { inexact = 0; goto set_exponent; } + else + goto add_one_ulp; + } + } + MPFR_ASSERTN(0); + + add_one_ulp: + /* add one unit in last place to a */ + DEBUG( printf("AddOneUlp\n") ); + if (MPFR_UNLIKELY( mpn_add_1(ap, ap, n, MPFR_LIMB_ONE<<sh) )) + { + /* Case 100000x0 = 0x1111x1 + 1*/ + DEBUG( printf("Pow of 2\n") ); + bx++; + ap[n-1] = MPFR_LIMB_HIGHBIT; + } + inexact = 1; + + set_exponent: + if (MPFR_UNLIKELY(bx > __gmpfr_emax)) /* Check for overflow */ + { + DEBUG( printf("Overflow\n") ); + MPFR_TMP_FREE(marker); + MPFR_SET_SAME_SIGN(a,b); + return mpfr_overflow(a, rnd_mode, MPFR_SIGN(a)); + } + MPFR_SET_EXP (a, bx); + MPFR_SET_SAME_SIGN(a,b); + + MPFR_TMP_FREE(marker); + MPFR_RET (inexact * MPFR_INT_SIGN (a)); +} |