* mpfr: Update to 2001-11-16. Patch TMP handling of agm.c and sqrt.c,

        separate .c files for floor and ceil.

1 files changed, 183 insertions, 0 deletions

diff --git a/mpfr/cmp2.c b/mpfr/cmp2.c
new file mode 100644
index 000000000..725cf5867
--- /dev/null
+++ b/mpfr/cmp2.c
@@ -0,0 +1,183 @@
+/* mpfr_cmp2 -- exponent shift when subtracting two numbers.
+
+Copyright (C) 1999-2001 Free Software Foundation.
+
+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 "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#include "mpfr.h"
+#include "mpfr-impl.h"
+
+/* returns the number of cancelled bits when one subtracts abs(c) from abs(b). 
+   Assumes |b| >= |c|, which implies MPFR_EXP(b)>=MPFR_EXP(c).
+   if |b| = |c|, returns prec(b).
+
+   Assumes neither of b or c is NaN or +/- infinity.
+
+   In other terms mpfr_cmp2 (b, c) returns EXP(b) - EXP(b-c).
+*/
+mp_prec_t
+mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c)
+{
+  mp_limb_t *bp, *cp, bb, cc = 0, lastc = 0, dif, high_dif = 0;
+  mp_size_t bn, cn;
+  mp_exp_unsigned_t diff_exp;
+  mp_prec_t res = 0;
+
+  MPFR_ASSERTN(MPFR_IS_FP(b));
+  MPFR_ASSERTN(MPFR_IS_FP(c));
+
+  if (MPFR_IS_ZERO(c))
+    return 0;
+
+  MPFR_ASSERTN(MPFR_NOTZERO(b));
+
+  bp = MPFR_MANT(b);
+  cp = MPFR_MANT(c);
+
+  bn = (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
+  cn = (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB;
+
+  MPFR_ASSERTN(MPFR_EXP(b) >= MPFR_EXP(c));
+  diff_exp = (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c);
+  
+  if (diff_exp == 0) /* otherwise the shifted most significant limb of c
+			cannot match bp[bn] */
+    {
+      while (bn>=0 && cn>=0 && bp[bn] == cp[cn])
+	{
+	  bn--;
+	  cn--;
+	  res += BITS_PER_MP_LIMB;
+	}
+
+      if (bn < 0) /* b = c */
+	return MPFR_PREC(b);
+
+      if (cn < 0) /* c discards exactly the upper part of b */
+	{
+          unsigned int z;
+
+	  while (bn>=0 && bp[bn]==0)
+	    {
+	      bn--;
+	      res += BITS_PER_MP_LIMB;
+	    }
+
+	  if (bn < 0) /* b = c */
+	    return MPFR_PREC(b);
+
+	  count_leading_zeros(z, bp[bn]); /* bp[bn] <> 0 */
+	  return res + z;
+	}
+    }
+
+  /* now we have removed the identical upper limbs of b and c
+     (can happen only when diff_exp = 0): bp[bn] > cc, bn>=0, cn>=0 */
+
+  if (diff_exp < BITS_PER_MP_LIMB)
+    {
+      cc = cp[cn] >> diff_exp;
+      /* warning: a shift by BITS_PER_MP_LIMB may give wrong results */
+      if (diff_exp)
+        lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp);
+      cn--;
+    }
+  else
+    diff_exp -= BITS_PER_MP_LIMB;
+
+  dif = bp[bn--] - cc; /* necessarily dif >= 1 */
+
+  while ((cn >= 0 || lastc != 0) && (high_dif == 0) && (dif == 1))
+    { /* dif=1 implies diff_exp = 0 or 1 */
+      bb = (bn >= 0) ? bp[bn] : 0;
+      cc = lastc;
+      if (cn >= 0)
+	{
+          if (diff_exp == 0)
+            {
+              cc += cp[cn];
+            }
+          else /* diff_exp = 1 */
+            {
+              cc += cp[cn] >> 1;
+              lastc = cp[cn] << (BITS_PER_MP_LIMB - 1);
+            }
+	}
+      else
+	lastc = 0;
+      high_dif = 1 - mpn_sub_n (&dif, &bb, &cc, 1);
+      bn--;
+      cn--;
+      res += BITS_PER_MP_LIMB;
+    }
+
+  /* (cn<0 and lastc=0) or (high_dif,dif)<>(0,1) */
+
+  if (high_dif != 0) /* necessarily high_dif = 1 */
+    {
+      res--;
+      if (dif != 0)
+	return res;
+    }
+  else /* high_dif = 0 */
+    {
+      unsigned int z;
+
+      count_leading_zeros(z, dif); /* dif > 1 here */
+      res += z;
+      if (dif != (MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - z - 1)))
+	return res; /* dif is not a power of two */
+    }
+
+  /* now result is res + (low(b) < low(c)) */
+  while (bn >= 0 && (cn >= 0 || lastc != 0))
+    {
+      if (diff_exp >= BITS_PER_MP_LIMB)
+	diff_exp -= BITS_PER_MP_LIMB;
+      else
+	{
+	  cc = lastc;
+	  if (cn >= 0)
+	    {
+	      cc += cp[cn] >> diff_exp;
+	      if (diff_exp != 0)
+                lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp);
+	    }
+	  else
+	    lastc = 0;
+	  cn--;
+	}
+      if (bp[bn] != cc)
+	return res + (bp[bn] < cc);
+      bn--;
+    }
+
+  if (bn < 0)
+    {
+      if (lastc != 0)
+	return res + 1;
+      while (cn >= 0 && cp[cn] == 0)
+        cn--;
+      return res + (cn >= 0);
+    }
+
+  return res; /* remainder from c is 0 */
+}