path: root/src/add1sp.c

/* mpfr_add1sp -- internal function to perform a "real" addition
   All the op must have the same precision

Copyright 2004-2016 Free Software Foundation, Inc.
Contributed by the AriC and Caramba 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 */
#if MPFR_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),
               (unsigned long) MPFR_PREC (a),
               (unsigned long) MPFR_PREC (b),
               (unsigned long) 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  /* MPFR_WANT_ASSERT >= 2 */

/* 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. */
MPFR_HOT_FUNCTION_ATTR 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_ASSERTD(MPFR_IS_PURE_FP(c));
  MPFR_ASSERTD(MPFR_GET_EXP(b) >= MPFR_GET_EXP(c));

  /* Read prec and num of limbs */
  p = MPFR_GET_PREC (b);
  n = MPFR_PREC2LIMBS (p);
  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 */

      /* Fast track for faithful rounding: since b and c have the same
         precision and the same exponent, the neglected value is either
         0 or 1/2 ulp(a), thus returning a is fine. */
      if (rnd_mode == MPFR_RNDF)
        goto set_exponent;

      if ((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 ((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 (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
             Faithful: Trunc */
          if (MPFR_LIKELY(rnd_mode == MPFR_RNDF))
            goto copy_set_exponent;
          else if (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 /* 0 < d < p */
    {
      mp_limb_t mask;
      mp_limb_t bcp, bcp1; /* Cp and C'p+1 */

      /* General case: 1 <= d < p */
      cp = MPFR_TMP_LIMBS_ALLOC (n);

      /* 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) );

      /* fast track for RNDF */
      if (rnd_mode == MPFR_RNDF)
        goto clean;

      /* 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:
      /* 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_RNDF))
        goto set_exponent;
      else if (rnd_mode == MPFR_RNDN)
        {
          inexact = - (bcp1 != 0);
          if (bcp == 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));
      inexact = -(bcp != 0 || bcp1 != 0);
      if (rnd_mode == MPFR_RNDZ || inexact == 0)
        goto set_exponent;
      else
        goto add_one_ulp;
    }
  MPFR_RET_NEVER_GO_HERE();

 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));
}