path: root/add1sp.c

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

Copyright 2004 Free Software Foundation.
Contributed by the Spaces project, INRIA Lorraine.

This file is part of the MPFR Library.

The 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 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 MPFR Library; see the file COPYING.LIB.  If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */

#include <stdio.h>

#define MPFR_NEED_LONGLONG_H
#include "mpfr-impl.h"

/* Define some internal MACROS */

/*#define DEBUG*/
/*#define CHECK_AGAINST_SUB1*/

#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, mp_rnd_t rnd_mode)
{
  mp_exp_unsigned_t d;
  mp_prec_t p, sh;
  mp_size_t n;
  mp_limb_t *ap, *cp;
  mp_exp_t  bx;
  mp_limb_t limb;
  int inexact;
  TMP_DECL(marker);

  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+BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB;
  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", 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==GMP_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==GMP_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 */
	  ap = MPFR_MANT(a);
	  MPN_COPY (ap, MPFR_MANT(b), n);
	  /* Away:    Add 1
	     Nearest: Trunc
	     Zero:    Trunc */
	  if (MPFR_LIKELY(rnd_mode==GMP_RNDN))
	    { inexact = -1; goto set_exponent; }
	  MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b));
	  if (rnd_mode==GMP_RNDZ)
	    { inexact = -1; goto set_exponent; }
	  else
	    goto add_one_ulp;
	}
      else
	{
	  /* d==p : Copy B in A */
	  ap = MPFR_MANT(a);
          MPN_COPY (ap, MPFR_MANT(b), n);
          /* Away:    Add 1
             Nearest: Even Rule if C is a power of 2, else Add 1 
             Zero:    Trunc */
          if (MPFR_LIKELY(rnd_mode==GMP_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 ((ap[0]&(MPFR_LIMB_ONE<<sh))==0)
		      { inexact = -1; goto set_exponent; }
		}
	      /* Not a Power of 2 */
	      goto add_one_ulp;
	    }
          MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b));
	  if (rnd_mode==GMP_RNDZ)
	    { inexact = -1; goto set_exponent; }
	  else
            goto 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*) TMP_ALLOC(n * BYTES_PER_MP_LIMB);
      
      /* Shift c in temporary allocated place */
      {
	mp_exp_unsigned_t dm;
	mp_size_t m;

	dm = d % BITS_PER_MP_LIMB;
	m = d / BITS_PER_MP_LIMB;
	if (MPFR_UNLIKELY(dm == 0))
	  {
	    /* dm = 0 and m > 0: Just copy */
	    MPFR_ASSERTD(m!=0);
	    /* Must use INCR since src & dest may overlap and dest <= src */
	    MPN_COPY_INCR(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))
        {
          /* 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 / BITS_PER_MP_LIMB);
                  mpfr_prec_t sx = BITS_PER_MP_LIMB-1-(x%BITS_PER_MP_LIMB);
                  DEBUG( printf("(First) x=%lu Kx=%ld Sx=%lu\n", x, kx, 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 / BITS_PER_MP_LIMB);
	  mpfr_prec_t sx = BITS_PER_MP_LIMB-1-(x%BITS_PER_MP_LIMB);
	  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=%lu Cp=%lu C'p+1=%lu\n", sh, bcp, 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", bcp, 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 == GMP_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 == GMP_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") );
      TMP_FREE(marker);
      MPFR_SET_SAME_SIGN(a,b);
      return mpfr_set_overflow(a, rnd_mode, MPFR_SIGN(a));
    }
  MPFR_SET_EXP (a, bx);
  MPFR_SET_SAME_SIGN(a,b);
  
  TMP_FREE(marker);
  return inexact*MPFR_INT_SIGN(a);
}