* mpf/ui_div.c: Use mpn_tdiv_qr. Use just one TMP_ALLOC. Use full

	divisor, since truncating can lose accuracy.

1 files changed, 50 insertions, 80 deletions

diff --git a/mpf/ui_div.c b/mpf/ui_div.c
index 74cdf97db..ba1552c25 100644
--- a/mpf/ui_div.c
+++ b/mpf/ui_div.c
@@ -1,7 +1,7 @@
 /* mpf_ui_div -- Divide an unsigned integer with a float.
 
-Copyright 1993, 1994, 1995, 1996, 2000, 2001, 2002 Free Software Foundation,
-Inc.
+Copyright 1993, 1994, 1995, 1996, 2000, 2001, 2002, 2004 Free Software
+Foundation, Inc.
 
 This file is part of the GNU MP Library.
 
@@ -20,20 +20,21 @@ along with the GNU MP 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>  /* for NULL */
 #include "gmp.h"
 #include "gmp-impl.h"
 #include "longlong.h"
 
+
 void
 mpf_ui_div (mpf_ptr r, unsigned long int u, mpf_srcptr v)
 {
   mp_srcptr vp;
-  mp_ptr rp, tp;
+  mp_ptr rp, tp, remp, new_vp;
   mp_size_t vsize;
-  mp_size_t rsize, tsize;
+  mp_size_t rsize, prospective_rsize, zeros, tsize, high_zero;
   mp_size_t sign_quotient;
   mp_size_t prec;
-  mp_limb_t q_limb;
   mp_exp_t rexp;
   TMP_DECL (marker);
 
@@ -42,10 +43,10 @@ mpf_ui_div (mpf_ptr r, unsigned long int u, mpf_srcptr v)
   vsize = ABS (vsize);
   prec = r->_mp_prec;
 
-  if (vsize == 0)
+  if (UNLIKELY (vsize == 0))
     DIVIDE_BY_ZERO;
 
-  if (u == 0)
+  if (UNLIKELY (u == 0))
     {
       r->_mp_size = 0;
       r->_mp_exp = 0;
@@ -53,94 +54,63 @@ mpf_ui_div (mpf_ptr r, unsigned long int u, mpf_srcptr v)
     }
 
   TMP_MARK (marker);
-  rexp = 1 - v->_mp_exp;
+  rexp = 1 - v->_mp_exp + 1;
 
   rp = r->_mp_d;
   vp = v->_mp_d;
 
-  if (vsize > prec)
-    {
-      vp += vsize - prec;
-      vsize = prec;
-    }
+  prospective_rsize = 1 - vsize + 1;    /* quot from using given u,v sizes */
+  rsize = prec + 1;                     /* desired quot size */
 
-  tsize = vsize + prec;
-  tp = (mp_ptr) TMP_ALLOC ((tsize + 1) * BYTES_PER_MP_LIMB);
-  MPN_ZERO (tp, tsize);
+  zeros = rsize - prospective_rsize;    /* padding u to give rsize */
+  tsize = 1 + zeros;                    /* u with zeros */
 
-  /* Normalize the divisor and the dividend.  */
-  if ((vp[vsize - 1] & GMP_NUMB_HIGHBIT) == 0)
+  if (WANT_TMP_DEBUG)
     {
-      mp_ptr tmp;
-      mp_limb_t dividend_high, dividend_low;
-      unsigned normalization_steps;
-
-      count_leading_zeros (normalization_steps, vp[vsize - 1]);
-      normalization_steps -= GMP_NAIL_BITS;
-
-      /* Shift up the divisor setting the most significant bit of
-	 the most significant limb.  Use temporary storage not to clobber
-	 the original contents of the divisor.  */
-      tmp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB);
-      mpn_lshift (tmp, vp, vsize, normalization_steps);
-      vp = tmp;
-
-      /* Shift up the dividend, possibly introducing a new most
-	 significant word.  */
-      dividend_high = (mp_limb_t) u >> (GMP_NUMB_BITS - normalization_steps);
-      dividend_low = ((mp_limb_t) u << normalization_steps) & GMP_NUMB_MASK;
-
-      tp[tsize - 1] = dividend_low;
-      if (dividend_high != 0)
-	{
-#if GMP_NAIL_BITS != 0
-	  if (dividend_high > vp[vsize - 1])
-	    {
-	      tp[tsize - 2] = dividend_low;
-	      tp[tsize - 1] = dividend_high & GMP_NUMB_MASK;
-	      tp[tsize] = dividend_high >> GMP_NUMB_BITS;
-	      tsize++;
-	      rexp += 2;
-	    }
-	  else
-#endif
-	    {
-	      tp[tsize] = dividend_high;
-	      tsize++;
-	      rexp++;
-	    }
-	}
+      /* separate alloc blocks, for malloc debugging */
+      remp = TMP_ALLOC_LIMBS (vsize);
+      tp = TMP_ALLOC_LIMBS (tsize);
+      new_vp = NULL;
+      if (rp == vp)
+        new_vp = TMP_ALLOC_LIMBS (vsize);
     }
   else
     {
-      /* The divisor is already normalized, as required.
-	 Copy it to temporary space if it overlaps with the quotient.  */
-      if (vp - rp <= tsize - vsize)
-	{
-	  mp_ptr tmp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB);
-	  MPN_COPY (tmp, vp, vsize);
-	  vp = (mp_srcptr) tmp;
-	}
-
-      tp[tsize - 1] = u & GMP_NUMB_MASK;
-#if GMP_NAIL_BITS != 0
-      if (u > GMP_NUMB_MAX)
-	{
-	  tp[tsize] = u >> GMP_NUMB_BITS;
-	  tsize++;
-	  rexp++;
-	}
-#endif
+      /* one alloc with calculated size, for efficiency */
+      remp = TMP_ALLOC_LIMBS (vsize + tsize + (rp == vp ? vsize : 0));
+      tp = remp + vsize;
+      new_vp = tp + tsize;
+    }
+
+  /* ensure divisor doesn't overlap quotient */
+  if (rp == vp)
+    {
+      MPN_COPY (new_vp, vp, vsize);
+      vp = new_vp;
     }
 
-  q_limb = mpn_divmod (rp, tp, tsize, vp, vsize);
-  rsize = tsize - vsize;
-  if (q_limb)
+  MPN_ZERO (tp, tsize-1);
+
+  tp[tsize - 1] = u & GMP_NUMB_MASK;
+#if BITS_PER_ULONG > GMP_NUMB_BITS
+  if (u > GMP_NUMB_MAX)
     {
-      rp[rsize] = q_limb;
-      rsize++;
+      /* tsize-vsize+1 == rsize, so tsize >= rsize.  rsize == prec+1 >= 2,
+         so tsize >= 2, hence there's room for 2-limb u with nails */
+      ASSERT (tsize >= 2);
+      tp[tsize - 1] = u >> GMP_NUMB_BITS;
+      tp[tsize - 2] = u & GMP_NUMB_MASK;
       rexp++;
     }
+#endif
+
+  ASSERT (tsize-vsize+1 == rsize);
+  mpn_tdiv_qr (rp, remp, (mp_size_t) 0, tp, tsize, vp, vsize);
+
+  /* strip possible zero high limb */
+  high_zero = (rp[rsize-1] == 0);
+  rsize -= high_zero;
+  rexp -= high_zero;
 
   r->_mp_size = sign_quotient >= 0 ? rsize : -rsize;
   r->_mp_exp = rexp;