add.c -> add.c & add1.c

sub.c -> sub.c & sub1.c + some changes.


git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@1490 280ebfd0-de03-0410-8827-d642c229c3f4

1 files changed, 0 insertions, 401 deletions

diff --git a/sub.c b/sub.c
index e53305f01..65d557143 100644
--- a/sub.c
+++ b/sub.c
@@ -20,414 +20,13 @@ 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. */
 
-/* #define DEBUG */
-
-#include <stdio.h>
 #include "gmp.h"
 #include "gmp-impl.h"
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
-#define ONE ((mp_limb_t) 1)
-
-/* signs of b and c differ, abs(b) > abs(c), 
-   diff_exp = EXP(b) - EXP(c).
-   Returns 0 iff result is exact,
-   a negative value when the result is less than the exact value,
-   a positive value otherwise.
-*/
-int
-#if __STDC__
-mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c,
-	   mp_rnd_t rnd_mode, mp_exp_unsigned_t diff_exp)
-#else
-mpfr_sub1 (a, b, c, rnd_mode, diff_exp)
-     mpfr_ptr a;
-     mpfr_srcptr b;
-     mpfr_srcptr c;
-     mp_rnd_t rnd_mode;
-     mp_exp_unsigned_t diff_exp;
-#endif
-{
-  unsigned long cancel, cancel1, sh, k;
-  long int cancel2, an, bn, cn, cn0;
-  mp_limb_t *ap, *bp, *cp, carry, bb, cc, borrow = 0;
-  int inexact = 0, shift_b, shift_c, is_exact = 1, down = 0, add_exp=0;
-  TMP_DECL(marker);
-
-#ifdef DEBUG
-  printf("\nenter mpfr_sub, rnd_mode=%s:\n", mpfr_print_rnd_mode(rnd_mode));
-  printf("b="); if (MPFR_SIGN(b)>0) putchar(' '); mpfr_print_raw(b); putchar('\n');
-  printf("c="); if (MPFR_SIGN(c)>0) putchar(' '); for (k=0; k<diff_exp; k++) putchar(' '); mpfr_print_raw(c); putchar('\n');
-  printf("PREC(a)=%u PREC(b)=%u PREC(c)=%u\n", MPFR_PREC(a), MPFR_PREC(b),
-	 MPFR_PREC(c));
-#endif
-  TMP_MARK(marker);
-  ap = MPFR_MANT(a);
-  an = 1 + (MPFR_PREC(a) - 1) / BITS_PER_MP_LIMB;
-
-  cancel = mpfr_cmp2 (b, c);
-
-  /* reserve a space to store b aligned with the result, i.e. shifted by
-     (-cancel) % BITS_PER_MP_LIMB to the right */
-  bn = 1 + (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
-  shift_b = cancel % BITS_PER_MP_LIMB;
-  if (shift_b)
-    shift_b = BITS_PER_MP_LIMB - shift_b;
-  cancel1 = (cancel + shift_b) / BITS_PER_MP_LIMB;
-  /* the high cancel1 limbs from b should not be taken into account */
-  if (shift_b == 0)
-    bp = MPFR_MANT(b); /* no need of an extra space */
-  else
-    {
-      bp = TMP_ALLOC ((bn + 1) * BYTES_PER_MP_LIMB);
-      bp[0] = mpn_rshift (bp + 1, MPFR_MANT(b), bn++, shift_b);
-    }
-
-  /* reserve a space to store c aligned with the result, i.e. shifted by
-     (diff_exp-cancel) % BITS_PER_MP_LIMB to the right */
-  cn = 1 + (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB;
-  shift_c = diff_exp - (cancel % BITS_PER_MP_LIMB);
-  shift_c = (shift_c + BITS_PER_MP_LIMB) % BITS_PER_MP_LIMB;
-  if (shift_c == 0)
-    cp = MPFR_MANT(c);
-  else
-    {
-      cp = TMP_ALLOC ((cn + 1) * BYTES_PER_MP_LIMB);
-      cp[0] = mpn_rshift (cp + 1, MPFR_MANT(c), cn++, shift_c);
-    }
-
-#ifdef DEBUG
-  printf("shift_b=%u shift_c=%u\n", shift_b, shift_c);
-#endif
-
-  /* ensure ap != bp and ap != cp */
-  if (ap == bp)
-    {
-      bp = (mp_ptr) TMP_ALLOC(bn * BYTES_PER_MP_LIMB);
-      MPN_COPY (bp, ap, bn);
-      /* ap == cp cannot occur since we would have b=c, which is detected
-	 in mpfr_add or mpfr_sub */
-    }
-  else if (ap == cp)
-    {
-      cp = (mp_ptr) TMP_ALLOC (cn * BYTES_PER_MP_LIMB);
-      MPN_COPY(cp, ap, cn);
-    }
-
-  cancel2 = (long int) (cancel + shift_c - diff_exp) / BITS_PER_MP_LIMB;
-  /* the high cancel2 limbs from b should not be taken into account */
-#ifdef DEBUG
-  printf("cancel=%u cancel1=%u cancel2=%d\n", cancel, cancel1, cancel2);
-#endif
-
-  /* adjust sign of result */
-  if (MPFR_SIGN(a)*MPFR_SIGN(b) < 0)
-    MPFR_CHANGE_SIGN(a);
-
-  /*               ap[an-1]        ap[0]
-             <----------------+-----------|---->
-             <----------PREC(a)----------><-sh->
- cancel1
- limbs        bp[bn-cancel1-1]
- <--...-----><----------------+-----------+----------->
-  cancel2
-  limbs       cp[cn-cancel2-1]                                    cancel2 >= 0
-    <--...--><----------------+----------------+---------------->
-                (-cancel2)                                        cancel2 < 0
-                   limbs      <----------------+---------------->
-  */
-
-  /* first part: put in ap[0..an-1] the value of high(b) - high(c),
-     where high(b) consists of the high an+cancel1 limbs of b,
-     and high(c) consists of the high an+cancel2 limbs of c.
-   */
-
-  /* copy high(b) into a */
-  if (an + cancel1 <= bn) /* a: <----------------+-----------|---->
-		         b: <-----------------------------------------> */
-      MPN_COPY (ap, bp + bn - (an + cancel1), an);
-  else  /* a: <----------------+-----------|---->
-       b: <-------------------------> */
-    if (cancel1 < bn) /* otherwise b does not overlap with a */
-      {
-	MPN_ZERO (ap, an + cancel1 - bn);
-	MPN_COPY (ap + an + cancel1 - bn, bp, bn - cancel1);
-      }
-    else
-      MPN_ZERO (ap, an);
-
-#ifdef DEBUG
-  printf("after copying high(b), a="); mpfr_print_raw(a); putchar('\n');
-#endif
-
-  /* subtract high(c) */
-  if (an + cancel2 > 0) /* otherwise c does not overlap with a */
-    {
-      mp_limb_t *ap2;
-
-      if (cancel2 >= 0)
-	{
-	  if (an + cancel2 <= cn) /* a: <----------------------------->
-			      c: <-----------------------------------------> */
-	    mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an);
-	  else /* a: <---------------------------->
-	      c: <-------------------------> */
-	    {
-	      ap2 = ap + an + cancel2 - cn;
-	      if (cn > cancel2)
-		mpn_sub_n (ap2, ap2, cp, cn - cancel2);
-	    }
-	}
-      else /* cancel2 < 0 */
-	{
-	  if (an + cancel2 <= cn) /* a: <----------------------------->
-			                  c: <-----------------------------> */
-	      borrow = mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an + cancel2);
-	  else /* a: <---------------------------->
-	                c: <----------------> */
-	    {
-	      ap2 = ap + an + cancel2 - cn;
-	      borrow = mpn_sub_n (ap2, ap2, cp, cn);
-	    }
-	  ap2 = ap + an + cancel2;
-	  mpn_sub_1 (ap2, ap2, -cancel2, borrow);
-	}
-    }
-
-#ifdef DEBUG
-  printf("after subtracting high(c), a="); mpfr_print_raw(a); putchar('\n');
-#endif
-
-  /* now perform rounding */
-  sh = an * BITS_PER_MP_LIMB - MPFR_PREC(a); /* last unused bits from a */
-  carry = ap[0] & ((ONE << sh) - 1);
-  ap[0] -= carry;
-
-  if (rnd_mode == GMP_RNDN)
-    {
-      if (sh)
-	{
-	  is_exact = (carry == 0);
-	  /* can decide except when carry = 2^(sh-1) [middle]
-	     or carry = 0 [truncate, but cannot decide inexact flag] */
-	  down = (carry < (ONE << (sh - 1)));
-	  if (carry > (ONE << (sh - 1)))
-	    goto add_one_ulp;
-	  else if ((0 < carry) && down)
-	    {
-	      inexact = -1; /* result if smaller than exact value */
-	      goto truncate;
-	    }
-	}
-    }
-  else /* directed rounding: set rnd_mode to RNDZ iff towards zero */
-    {
-      if (((rnd_mode == GMP_RNDD) && (MPFR_SIGN(b) > 0)) ||
-	  ((rnd_mode == GMP_RNDU) && (MPFR_SIGN(b) < 0)))
-	rnd_mode = GMP_RNDZ;
-
-      if (carry)
-	{
-	  if (rnd_mode == GMP_RNDZ)
-	    {
-	      inexact = -1;
-	      goto truncate;
-	    }
-	  else /* round away */
-	    goto add_one_ulp;
-	}
-    }
-
-  /* we have to consider the low (bn - (an+cancel1)) limbs from b,
-     and the (cn - (an+cancel2)) limbs from c. */
-  bn -= an + cancel1;
-  cn0 = cn;
-  cn -= (long int) an + cancel2;
-#ifdef DEBUG
-  printf("last %u bits from a are %lu, bn=%ld, cn=%ld\n", sh, carry, bn, cn);
-#endif
-
-  for (k=0; (bn > 0) || (cn > 0); k++)
-    {
-      bb = (bn > 0) ? bp[--bn] : 0;
-      if ((cn > 0) && (cn-- <= cn0))
-	cc = cp[cn];
-      else
-	cc = 0;
-
-#ifdef DEBUG
-      printf("k=%u bb=%lu cc=%lu down=%d\n", k, bb, cc, down);
-#endif
-      if (down == 0)
-	down = (bb < cc);
-
-      if ((rnd_mode == GMP_RNDN) && !k && !sh)
-	{
-	  mp_limb_t half = ONE << (BITS_PER_MP_LIMB - 1);
-
-	  is_exact = (bb == cc);
-
-	  /* add one ulp if bb > cc + half
-	     truncate if cc - half < bb < cc + half
-	     sub one ulp if bb < cc - half
-	  */
-
-	  if (down)
-	    {
-	      if (cc >= half)
-		cc -= half;
-	      else
-		bb += half;
-	    }
-	  else /* bb >= cc */
-	    {
-	      if (cc < half)
-		cc += half;
-	      else
-		bb -= half;
-	    }
-	}
-
-#ifdef DEBUG
-      printf("    bb=%lu cc=%lu down=%d is_exact=%d\n", bb, cc, down, is_exact);
-#endif
-      if (bb < cc)
-	{
-	  if (rnd_mode == GMP_RNDZ)
-	    goto sub_one_ulp;
-	  else if (rnd_mode != GMP_RNDN) /* round away */
-	    {
-	      inexact = 1;
-	      goto truncate;
-	    }
-	  else /* round to nearest */
-	    {
-	      if (is_exact && !sh)
-		{
-		  inexact = 0;
-		  goto truncate;
-		}
-	      else if (down && !sh)
-		goto sub_one_ulp;
-	      else
-		{
-		  inexact = (is_exact) ? 1 : -1;
-		  goto truncate;
-		}
-	    }
-	}
-      else if (bb > cc)
-	{
-	  if (rnd_mode == GMP_RNDZ)
-	    {
-	      inexact = -1;
-	      goto truncate;
-	    }
-	  else if (rnd_mode != GMP_RNDN) /* round away */
-	      goto add_one_ulp;
-	  else /* round to nearest */
-	    {
-	      if (is_exact)
-		{
-		  inexact = -1;
-		  goto truncate;
-		}
-	      else if (down)
-		{
-		  inexact = 1;
-		  goto truncate;
-		}
-	      else
-		goto add_one_ulp;
-	    }
-	}
-    }
-
-  if ((rnd_mode == GMP_RNDN) && !is_exact)
-    {
-      /* even rounding rule */
-      if ((ap[0] >> sh) & 1)
-	{
-	  if (down) goto sub_one_ulp;
-	  else goto add_one_ulp;
-	}
-      else
-	inexact = (down) ? 1 : -1;
-    }
-  else
-    inexact = 0;
-  goto truncate;
-  
- sub_one_ulp: /* add one unit in last place to a */
-  mpn_sub_1 (ap, ap, an, ONE << sh);
-  inexact = -1;
-  goto end_of_sub;
-
- add_one_ulp: /* add one unit in last place to a */
-  if (mpn_add_1 (ap, ap, an, ONE << sh)) /* result is a power of two */
-    {
-      ap[an-1] |= ONE << (BITS_PER_MP_LIMB - 1);
-      add_exp = 1;
-    }
-  inexact = 1; /* result larger than exact value */
-
- truncate:
-  if ((ap[an-1] >> (BITS_PER_MP_LIMB - 1)) == 0) /* case 1 - epsilon */
-    {
-      ap[an-1] = ONE << (BITS_PER_MP_LIMB - 1);
-      add_exp = 1;
-    }
-
- end_of_sub:
-  /* we have to set MPFR_EXP(a) to MPFR_EXP(b) - cancel + add_exp, taking
-     care of underflows/overflows in that computation, and of the allowed
-     exponent range */
-  if (cancel)
-    {
-      mp_exp_t exp_b;
-
-      cancel -= add_exp; /* still valid as unsigned long */
-      exp_b = MPFR_EXP(b); /* save it in case a equals b */
-      MPFR_EXP(a) = MPFR_EXP(b) - cancel;
-      if ((MPFR_EXP(a) > exp_b) /* underflow in type mp_exp_t */
-	  || (MPFR_EXP(a) < __mpfr_emin))
-	{
-	  TMP_FREE(marker); 
-	  return mpfr_set_underflow (a, rnd_mode, MPFR_SIGN(a));
-	}
-    }
-  else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */
-    {
-      /* in case cancel = 0, add_exp can still be 1, in case b is just
-	 below a power of two, c is very small, prec(a) < prec(b),
-	 and rnd=away or nearest */
-      if (add_exp && MPFR_EXP(b) == __mpfr_emax)
-	{
-	  TMP_FREE(marker);
-	  return mpfr_set_overflow (a, rnd_mode, MPFR_SIGN(a));
-	}
-      MPFR_EXP(a) = MPFR_EXP(b) + add_exp;
-    }
-  TMP_FREE(marker);
-#ifdef DEBUG
-  printf ("result is a="); mpfr_print_raw(a); putchar('\n');
-#endif
-  /* check that result is msb-normalized */
-  ASSERT_ALWAYS(ap[an-1] > ~ap[an-1]);
-  return inexact * MPFR_SIGN(b);
-}
-
 int
-#if __STDC__
 mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
-#else
-mpfr_sub (a, b, c, rnd_mode)
-     mpfr_ptr a;
-     mpfr_srcptr b;
-     mpfr_srcptr c;
-     mp_rnd_t rnd_mode;
-#endif
 {
   if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c))
   {