New division in div.c, old one renamed in mpfr_div2. Remains to implement

exact/inexact flag. Should not be *that* hard.


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

2 files changed, 554 insertions, 158 deletions

diff --git a/div.c b/div.c
index 6e7bc7fd1..589707eb3 100644
--- a/div.c
+++ b/div.c
@@ -27,7 +27,21 @@ MA 02111-1307, USA. */
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
-/* #define DEBUG */
+// #define DEBUG
+
+#ifdef DEBUG
+void 
+affiche_mp(mp_ptr z, mp_size_t length)
+{
+  int k; 
+
+  if (length == 1) { printf("[%lu]\n", *z); return; }
+
+  printf("[%lu, ", z[length - 1]); 
+  for(k = length - 2; k >= 1; k--) { printf("%lu, ", z[k]); }
+  printf("%lu]\n", z[0]); 
+}
+#endif
 
 void
 #if __STDC__
@@ -42,17 +56,22 @@ mpfr_div (r, u, v, rnd_mode)
 {
   mp_srcptr up, vp;
   mp_ptr rp, tp, tp0, tmp, tmp2;
-  mp_size_t usize, vsize, rrsize, oldrrsize;
-  mp_size_t rsize;
-  mp_size_t sign_quotient;
-  mp_size_t prec, err;
+  mp_size_t usize, vsize, drsize, dsize;
+  mp_size_t oldrsize, rsize, sign_quotient, err;
   mp_limb_t q_limb;
   mp_exp_t rexp;
-  long k, mult, vn, t; 
+  long k, t; 
   unsigned long cc = 0, rw, nw; 
-  char can_round = 0; 
+  char can_round = 0, sh; 
   TMP_DECL (marker);
 
+
+  /**************************************************************************
+   *                                                                        *
+   *              This part of the code deals with special cases            *
+   *                                                                        *
+   **************************************************************************/
+
   if (MPFR_IS_NAN(u) || MPFR_IS_NAN(v)) { MPFR_SET_NAN(r); return; }
 
   MPFR_CLEAR_NAN(r);
@@ -81,12 +100,6 @@ mpfr_div (r, u, v, rnd_mode)
 
   MPFR_CLEAR_INF(r); /* clear Inf flag */
 
-  usize = (MPFR_PREC(u) - 1)/BITS_PER_MP_LIMB + 1; 
-  vsize = (MPFR_PREC(v) - 1)/BITS_PER_MP_LIMB + 1; 
-  sign_quotient = ((MPFR_SIGN(u) * MPFR_SIGN(v) > 0) ? 1 : -1); 
-  prec = MPFR_PREC(r);
-
-
   if (!MPFR_NOTZERO(v))
     {
       if (!MPFR_NOTZERO(u)) 
@@ -102,209 +115,282 @@ mpfr_div (r, u, v, rnd_mode)
   
   if (!MPFR_NOTZERO(u)) { MPFR_SET_ZERO(r); return; }
 
-  up = MPFR_MANT(u);
-  vp = MPFR_MANT(v);
+  /**************************************************************************
+   *                                                                        *
+   *              End of the part concerning special values.                *
+   *                                                                        *
+   **************************************************************************/
+
 
 #ifdef DEBUG
-      printf("Entering division : "); 
-      for(k = usize - 1; k >= 0; k--) { printf("%lu ", up[k]); }
-      printf(" by "); 
-      for(k = vsize - 1; k >= 0; k--) { printf("%lu ", vp[k]); }
-      printf(".\n"); 
+  printf("u = "); mpfr_out_str(stdout, 2, 0, u, GMP_RNDN); printf("\n"); 
+  printf("v = "); mpfr_out_str(stdout, 2, 0, v, GMP_RNDN); printf("\n"); 
 #endif
 
-  /* Compare the mantissas */
-  mult = mpn_cmp(up, vp, (usize > vsize ? vsize : usize)); 
-  if (mult == 0 && vsize > usize)
-    {
-      vn = vsize - usize; 
-      while (vn >= 0) if (vp[vn--]) { mult = 1; break; }
-      /* On peut diagnostiquer ici pour pas cher le cas u = v */
-    }
-  else { mult = (mult < 0 ? 1 : 0); }
+  up = MPFR_MANT(u);
+  vp = MPFR_MANT(v);
 
-  rsize = (MPFR_PREC(r) + 3)/BITS_PER_MP_LIMB + 1; 
-  rrsize = MPFR_PREC(r)/BITS_PER_MP_LIMB + 1;
-  /* Three extra bits are needed in order to get the quotient with enough
-     precision ; take one extra bit for rrsize in order to solve more 
-     easily the problem of rounding to nearest. */
+  TMP_MARK (marker);
 
-  do
-    {
-      TMP_MARK (marker);
+  usize = (MPFR_PREC(u) - 1)/BITS_PER_MP_LIMB + 1; 
+  vsize = (MPFR_PREC(v) - 1)/BITS_PER_MP_LIMB + 1; 
 
-      rexp = MPFR_EXP(u) - MPFR_EXP(v);
-      
-      err = rsize*BITS_PER_MP_LIMB; 
-      if (rsize < vsize) { err-=2; } 
-      if (rsize < usize) { err--; }
-      if (err > rrsize * BITS_PER_MP_LIMB) 
-	{ err = rrsize * BITS_PER_MP_LIMB; }
+#ifdef DEBUG
+  printf("Entering division : "); 
+  affiche_mp(up, usize); 
+  affiche_mp(vp, vsize); 
+#endif
       
-      tp0 = (mp_ptr) TMP_ALLOC ((rsize+rrsize) * BYTES_PER_MP_LIMB);
-      /* fill by zero rrsize low limbs of t */
-      MPN_ZERO(tp0, rrsize); tp = tp0 + rrsize; 
-      rp = (mp_ptr) TMP_ALLOC ((rrsize+1) * BYTES_PER_MP_LIMB); 
-
-      if (vsize >= rsize) { 
-	tmp = (mp_ptr) vp + vsize - rsize;
-      }
-      else { 
-	tmp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
-	MPN_COPY (tmp + rsize - vsize, vp, vsize);
-	MPN_ZERO (tmp, rsize - vsize); 
-      }
-
-      if (usize >= rsize) { 
-	MPN_COPY (tp, up + usize - rsize, rsize);
-      }
-      else {
-	MPN_COPY (tp + rsize - usize, up, usize); 
-	MPN_ZERO (tp, rsize - usize); 
-      }
+  /**************************************************************************
+   *                                                                        *
+   *   First try to use only part of u, v. If this is not sufficient,       *
+   *   use the full u and v, to avoid long computations eg. in the case     *
+   *   u = v.                                                               *
+   *                                                                        *
+   **************************************************************************/
+
+  dsize = (MPFR_PREC(r) + 3)/BITS_PER_MP_LIMB + 1; 
+  drsize = MPFR_PREC(r)/BITS_PER_MP_LIMB + 1 + dsize; 
+
+  /* Compute effective dividend */
+  if (vsize < dsize) { dsize = vsize; }
+  tmp = (mp_ptr) vp + vsize - dsize; 
+
+  /* Compute effective divisor. One extra bit allowed for RNDN. */
+  
+  tp0 = (mp_ptr) TMP_ALLOC(drsize * BYTES_PER_MP_LIMB); 
+  if (usize >= drsize) 
+    MPN_COPY (tp0, up + usize - drsize, drsize); 
+  else {
+    MPN_COPY (tp0 + drsize - usize, up, usize); 
+    MPN_ZERO (tp0, drsize - usize); 
+  }
 
-      /* Do the real job */
- 
+  /* Allocate limbs for quotient. */
+  rp = (mp_ptr) TMP_ALLOC ((drsize - dsize + 1) * BYTES_PER_MP_LIMB);
+      
 #ifdef DEBUG
-      printf("Dividing : "); 
-      for(k = rsize - 1; k >= 0; k--) { printf("%lu ", tp[k]); }
-      printf(" by "); 
-      for(k = rsize - 1; k >= 0; k--) { printf("%lu ", tmp[k]); }
-      printf(".\n"); 
+  printf("Dividing : "); 
+  affiche_mp(tp0, drsize); 
+  printf(" by "); 
+  affiche_mp(tmp, dsize); 
+  printf(".\n"); 
 #endif
-
+      
 #if (__GNU_MP_VERSION < 3)
-      q_limb = mpn_divrem (rp, 0, tp0, rsize+rrsize, tmp, rsize);
-      tp = tp0; /* location of remainder */
-#else /* mpn_tdiv_qr is the preferred division interface in GMP 3 */
-      tmp2 = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
-      mpn_tdiv_qr(rp, tmp2, 0, tp0, rsize+rrsize, tmp, rsize);
-      q_limb = rp[rrsize];
-      tp = tmp2; /* location of remainder */
+  q_limb = mpn_divrem (rp, 0, tp0, drsize, tmp, dsize);
+  tp = tp0;
+#else 
+  tmp2 = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB);
+  mpn_tdiv_qr(rp, tmp2, 0, tp0, drsize, tmp, dsize);
+  q_limb = rp[drsize - dsize];
+  tp = tmp2;
 #endif
+  
+  /* Estimate number of correct bits. */
+  err = (drsize - dsize) * BITS_PER_MP_LIMB; 
+  if (drsize < usize) err --; 
+  if (dsize < vsize) err -= 2; 
 
 #ifdef DEBUG
-      printf("The result is : \n"); 
-      printf("Quotient : "); 
-      for(k = rrsize - 1; k >= 0; k--) { printf("%lu ", rp[k]); }
-      printf("Remainder : "); 
-      for(k = rsize - 1; k >= 0; k--) { printf("%lu ", tp[k]); }
-      printf("(q_limb = %lu)\n", q_limb); 
+  printf("The result is : \n"); 
+  printf("Quotient : "); 
+
+  affiche_mp(rp, drsize - dsize + 1); 
+ 
+  printf("Remainder : "); 
+  affiche_mp(tp, dsize); 
+
+  printf("Number of correct bits = %lu\n", err); 
 #endif
       
-      /* msb-normalize the result */
+  /* Compute sign and exponent. Correction might occur just below. */
+  sign_quotient = ((MPFR_SIGN(u) * MPFR_SIGN(v) > 0) ? 1 : -1); 
+  rexp = MPFR_EXP(u) - MPFR_EXP(v);
+  rsize = drsize - dsize;
+  
+  /* We want to check if rounding is possible. Have to mimic normalization */
+  if (q_limb) { sh = -1; } 
+  else { count_leading_zeros(sh, rp[rsize - 1]); } 
+  
+  can_round = mpfr_can_round_raw(rp, rsize + 1, sign_quotient, err + sh + 
+				 BITS_PER_MP_LIMB, GMP_RNDN, rnd_mode, 
+				 MPFR_PREC(r) + sh + BITS_PER_MP_LIMB); 
+
+  if (!can_round && (drsize < usize || dsize < vsize)) 
+    {
+      mp_size_t ulrsize; 
+      char b = 0; 
+      mp_ptr rp2, ulorem, ulorem2; 
+
+  /**************************************************************************
+   *                                                                        *
+   *   The attempt to use only part of u and v failed. We first compute a   *
+   *   correcting term, then perform the full division.                     *
+   *   Put u = uhi + ulo, v = vhi + vlo. We have uhi = vhi * rp + tp,       *
+   *   thus u - v * rp = tp + ulo - rp*vlo, that we shall divide by v.      *
+   *                                                                        *
+   **************************************************************************/
+#ifdef DEBUG
+      printf("Using the full u and v.\n");
+#endif
       
-      if (q_limb)
-	{
-	  count_leading_zeros(k, q_limb); 
-	  mpn_rshift(rp, rp, rrsize, BITS_PER_MP_LIMB - k);
-	  rp[rrsize - 1] |= (q_limb << k);
-	  rexp += BITS_PER_MP_LIMB - k; 
-	}
-      else
-	{ 
-	  count_leading_zeros(k, rp[rrsize - 1]); 
-	  if (k) { mpn_lshift(rp, rp, rrsize, k); }
-	  rexp -= k; 
-	}
+      if (usize > drsize)
+	  ulrsize = usize - drsize + dsize; 
+	  // store the low part of u and the remainder. 
+      else ulrsize = dsize; // just store the remainder. 
 
-      can_round = (mpfr_can_round_raw(rp, rrsize, sign_quotient, err, 
-				     GMP_RNDN, rnd_mode, MPFR_PREC(r))
-	|| (usize == rsize && vsize == rsize && 
-	    mpfr_can_round_raw(rp, rrsize, sign_quotient, err, 
-			       GMP_RNDZ, rnd_mode, MPFR_PREC(r))));
+      if (vsize > dsize) ulrsize++; 
 
-      /* If we used all the limbs of both the dividend and the divisor, 
-	 then we have the correct RNDZ rounding */
+      ulorem = TMP_ALLOC(ulrsize * BYTES_PER_MP_LIMB);
+      ulorem2 = TMP_ALLOC(ulrsize * BYTES_PER_MP_LIMB);
+      ulorem[ulrsize - 1] = ulorem2 [ulrsize - 1] = 0; 
 
-      if (!can_round && (rsize < usize || rsize < vsize)) 
+      if (dsize < vsize)
 	{
+	  mp_size_t lg = vsize + rsize - dsize + 1; 
+
+	  if (rsize > vsize - dsize)
+	    mpn_mul(ulorem + ulrsize - lg, rp, rsize, vp, vsize - dsize);
+	  else
+	    mpn_mul(ulorem + ulrsize - lg, vp, vsize - dsize, rp, rsize);
+
+	  MPN_ZERO(ulorem, ulrsize - lg); 
+	}
+      else MPN_ZERO(ulorem, ulrsize); 
+
 #ifdef DEBUG
-	  printf("Increasing the precision.\n");
+      printf("vlo * q: "); 
+      affiche_mp(ulorem, ulrsize); 
+#endif 
+     
+      MPN_COPY(ulorem2 + ulrsize - 1 - dsize, tp, dsize);      
+      if (drsize < usize) 
+	MPN_COPY(ulorem2, up, usize - drsize); 
+
+#ifdef DEBUG
+      printf("(b = %d) ulo + r: ", b); 
+      affiche_mp(ulorem2, ulrsize); 
 #endif
-	  TMP_FREE(marker);
-	  /* in case we can't round at the first iteration,
-	     jump right away to the maximum precision of both
-	     operands, to avoid multiple iterations when for 
-	     example the divisor is 1.0.
-	  */
-	  if (rsize < usize) rsize = usize - 1;
-	  if (rsize < vsize) rsize = vsize - 1;
+      
+      if (mpn_cmp(ulorem2, ulorem, ulrsize) > 0)
+	mpn_sub_n(ulorem, ulorem2, ulorem, ulrsize); 
+      else 
+	{ 
+	  ulorem2[ulrsize - 1] = 
+	    mpn_add_n(ulorem2 + ulrsize - vsize - 1, 
+		      ulorem2 + ulrsize - vsize - 1, vp, vsize); 
+
+	  if (mpn_cmp(ulorem2, ulorem, ulrsize) < 0)
+	    {
+	      b = 2; 
+	      ulorem2[ulrsize - 1] += 
+		mpn_add_n(ulorem2 + ulrsize - vsize - 1, 
+			  ulorem2 + ulrsize - vsize - 1, vp, vsize); 
+	    }
+	  else b = 1; 
+
+	  mpn_sub_n(ulorem, ulorem2, ulorem, ulrsize); 
 	}
+	  
+#ifdef DEBUG
+  printf("(b = %d) ulo + r - vlo * q: ", b); 
+  affiche_mp(ulorem, ulrsize); 
+#endif
+     
+      rp2 = (mp_ptr) TMP_ALLOC((ulrsize - vsize + 1) * BYTES_PER_MP_LIMB); 
+
+#if (__GNU_MP_VERSION < 3)
+      q_limb = mpn_divrem (rp2, 0, ulorem, ulrsize, vp, vsize);
+      tp = ulorem; 
+#else 
+      tmp2 = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB);
+      mpn_tdiv_qr(rp2, tmp2, 0, ulorem, ulrsize, vp, vsize);
+      tp = tmp2; 
+#endif
+
+      if (!b)
+	mpn_add_1(rp, rp, rsize, rp2[ulrsize - vsize]);
+      else
+	mpn_sub_1(rp, rp, rsize, b); 
+
+      dsize = vsize; /* The length of the remainder, in the end of the code */
     }
-  while (!can_round && (rsize < usize || rsize < vsize) 
-	 && (rsize++) && (rrsize++));
 
-  /* ON PEUT PROBABLEMENT SE DEBROUILLER DES QUE rsize >= vsize */
-  /* MAIS IL FAUT AJOUTER LE BOUT QUI MANQUE DE usize A rsize */
-    
-  oldrrsize = rrsize;
-  rrsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1;
+  /**************************************************************************
+   *                                                                        *
+   *                       Final stuff (rounding and so.)                   *
+   *                                                                        *
+   **************************************************************************/
+
+  if (q_limb)
+    {
+      mpn_rshift(rp, rp, rsize, 1);
+      rp[rsize - 1] |= 1 << (BITS_PER_MP_LIMB - 1); rexp ++; 
+    }
+  else
+    if (sh) { mpn_lshift(rp, rp, rsize, sh); rexp -= sh; }
+  
+  oldrsize = rsize;
+  rsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1;
 
-  if (can_round) 
+  if (can_round) /* Lazy case. */
     {
       cc = mpfr_round_raw(rp, rp, err, (sign_quotient == -1 ? 1 : 0),
 			  MPFR_PREC(r), rnd_mode, NULL);
     }
   else {
-    /* Use the remainder to find out the correct rounding */
-    /* Note that at this point the division has been done */
-    /* EXACTLY. */
+    rp += oldrsize-rsize;
     if ((rnd_mode == GMP_RNDD && sign_quotient == -1) 
 	|| (rnd_mode == GMP_RNDU && sign_quotient == 1)
 	|| (rnd_mode == GMP_RNDN))
       {
 	/* We cannot round, so that the last bits of the quotient
 	   have to be zero; just look if the remainder is nonzero */
-	k = rsize - 1; 
+
+	k = dsize - 1; 
 	while (k >= 0) { if (tp[k]) break; k--; }
 	if (k >= 0) 
 	  {
 	    t = MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1); 
 	    if (t)
-	      {
-		cc = mpn_add_1(rp, rp, rrsize, 
-			       (mp_limb_t)1 << (BITS_PER_MP_LIMB - t)); 
-	      }
+	      cc = mpn_add_1(rp, rp, rsize, 
+			     (mp_limb_t)1 << (BITS_PER_MP_LIMB - t)); 
 	    else
-	      {
-		cc = mpn_add_1(rp, rp, rrsize, 1); 
-	      }
+	      cc = mpn_add_1(rp, rp, rsize, 1); 
 	  }
 	else
-	  if (rnd_mode == GMP_RNDN) /* even rounding */
-	    {
-	      rw = (MPFR_PREC(r) + 1) & (BITS_PER_MP_LIMB - 1);
-	      if (rw) { rw = BITS_PER_MP_LIMB - rw; nw = 0; } else nw = 1; 
-	      if ((rw ? (rp[nw] >> (rw + 1)) & 1 : 
-		   (rp[nw] >> (BITS_PER_MP_LIMB - 1)) & 1))
-		{
-		  cc = mpn_add_1(rp + nw, rp + nw, rrsize, 
-				 ((mp_limb_t)1) << rw); 
-		}
-	    }
+	  {
+	    if (rnd_mode == GMP_RNDN) /* even rounding */
+	      {
+		rw = (MPFR_PREC(r) + 1) & (BITS_PER_MP_LIMB - 1);
+		if (rw) { rw = BITS_PER_MP_LIMB - rw; nw = 0; } else nw = 1; 
+		if ((rw ? (rp[nw] >> (rw + 1)) & 1 : 
+		     (rp[nw] >> (BITS_PER_MP_LIMB - 1)) & 1))
+		  {
+		    cc = mpn_add_1(rp + nw, rp + nw, rsize, 
+				   ((mp_limb_t)1) << rw); 
+		  }
+	      }
 	/* cas 0111111 */
+	  }
       }
-    /* Warning: we computed the result on oldrrsize limbs, but we want it
-       on rrsize limbs only. Both can differ, especially when the target
-       precision is a multiple of the number of bits per limb, since we've
-       taken an extra bit to make rounding to nearest easier. */
-    rp += oldrrsize-rrsize;
   }
 
-
   if (sign_quotient * MPFR_SIGN(r) < 0) { MPFR_CHANGE_SIGN(r); } 
   MPFR_EXP(r) = rexp;
 
   if (cc) {
-    mpn_rshift(rp, rp, rrsize, 1);
-    rp[rrsize-1] |= (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1);
+    mpn_rshift(rp, rp, rsize, 1);
+    rp[rsize-1] |= (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1);
     MPFR_EXP(r)++; 
   }
     
-  rw = rrsize * BITS_PER_MP_LIMB - MPFR_PREC(r);
-  MPN_COPY(MPFR_MANT(r), rp, rrsize); 
+  rw = rsize * BITS_PER_MP_LIMB - MPFR_PREC(r);
+  MPN_COPY(MPFR_MANT(r), rp, rsize); 
   MPFR_MANT(r)[0] &= ~(((mp_limb_t)1 << rw) - 1);
-
+#ifdef DEBUG
+  printf("r = "); mpfr_out_str(stdout, 2, 0, r, GMP_RNDN); printf("\n"); 
+#endif
   TMP_FREE (marker);
 }
diff --git a/div2.c b/div2.c
new file mode 100644
index 000000000..78452f722
--- /dev/null
+++ b/div2.c
@@ -0,0 +1,310 @@
+/* mpfr_div -- divide two floating-point numbers
+
+Copyright (C) 1999 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 Library General Public License as published by
+the Free Software Foundation; either version 2 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library 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>
+#include <stdlib.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#include "mpfr.h"
+#include "mpfr-impl.h"
+
+/* #define DEBUG */
+
+void
+#if __STDC__
+mpfr_div2 (mpfr_ptr r, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode)
+#else
+mpfr_div2 (r, u, v, rnd_mode)
+     mpfr_ptr r;
+     mpfr_srcptr u;
+     mpfr_srcptr v;
+     mp_rnd_t rnd_mode;
+#endif
+{
+  mp_srcptr up, vp;
+  mp_ptr rp, tp, tp0, tmp, tmp2;
+  mp_size_t usize, vsize, rrsize, oldrrsize;
+  mp_size_t rsize;
+  mp_size_t sign_quotient;
+  mp_size_t prec, err;
+  mp_limb_t q_limb;
+  mp_exp_t rexp;
+  long k, mult, vn, t; 
+  unsigned long cc = 0, rw, nw; 
+  char can_round = 0; 
+  TMP_DECL (marker);
+
+  if (MPFR_IS_NAN(u) || MPFR_IS_NAN(v)) { MPFR_SET_NAN(r); return; }
+
+  MPFR_CLEAR_NAN(r);
+
+  if (MPFR_IS_INF(u)) 
+    { 
+      if (MPFR_IS_INF(v)) 
+	MPFR_SET_NAN(r);
+      else
+	{ 
+	  MPFR_SET_INF(r); 
+	  if (MPFR_SIGN(r) != MPFR_SIGN(u) * MPFR_SIGN(v)) 
+	    MPFR_CHANGE_SIGN(r);
+	}
+      return;
+    }
+  else 
+    if (MPFR_IS_INF(v)) 
+      {
+	MPFR_CLEAR_INF(r);
+	MPFR_SET_ZERO(r); 
+	if (MPFR_SIGN(r) != MPFR_SIGN(u) * MPFR_SIGN(v)) 
+	  MPFR_CHANGE_SIGN(r);
+	return; 
+      }
+
+  MPFR_CLEAR_INF(r); /* clear Inf flag */
+
+  usize = (MPFR_PREC(u) - 1)/BITS_PER_MP_LIMB + 1; 
+  vsize = (MPFR_PREC(v) - 1)/BITS_PER_MP_LIMB + 1; 
+  sign_quotient = ((MPFR_SIGN(u) * MPFR_SIGN(v) > 0) ? 1 : -1); 
+  prec = MPFR_PREC(r);
+
+
+  if (!MPFR_NOTZERO(v))
+    {
+      if (!MPFR_NOTZERO(u)) 
+	{ MPFR_SET_NAN(r); return; }
+      else
+	{
+	  MPFR_SET_INF(r); 
+	  if (MPFR_SIGN(r) != MPFR_SIGN(v) * MPFR_SIGN(u)) 
+	    MPFR_CHANGE_SIGN(r); 
+	  return;
+	}
+    }
+  
+  if (!MPFR_NOTZERO(u)) { MPFR_SET_ZERO(r); return; }
+
+  up = MPFR_MANT(u);
+  vp = MPFR_MANT(v);
+
+#ifdef DEBUG
+      printf("Entering division : "); 
+      for(k = usize - 1; k >= 0; k--) { printf("%lu ", up[k]); }
+      printf(" by "); 
+      for(k = vsize - 1; k >= 0; k--) { printf("%lu ", vp[k]); }
+      printf(".\n"); 
+#endif
+
+  /* Compare the mantissas */
+  mult = mpn_cmp(up, vp, (usize > vsize ? vsize : usize)); 
+  if (mult == 0 && vsize > usize)
+    {
+      vn = vsize - usize; 
+      while (vn >= 0) if (vp[vn--]) { mult = 1; break; }
+      /* On peut diagnostiquer ici pour pas cher le cas u = v */
+    }
+  else { mult = (mult < 0 ? 1 : 0); }
+
+  rsize = (MPFR_PREC(r) + 3)/BITS_PER_MP_LIMB + 1; 
+  rrsize = MPFR_PREC(r)/BITS_PER_MP_LIMB + 1;
+  /* Three extra bits are needed in order to get the quotient with enough
+     precision ; take one extra bit for rrsize in order to solve more 
+     easily the problem of rounding to nearest. */
+
+  do
+    {
+      TMP_MARK (marker);
+
+      rexp = MPFR_EXP(u) - MPFR_EXP(v);
+      
+      err = rsize*BITS_PER_MP_LIMB; 
+      if (rsize < vsize) { err-=2; } 
+      if (rsize < usize) { err--; }
+      if (err > rrsize * BITS_PER_MP_LIMB) 
+	{ err = rrsize * BITS_PER_MP_LIMB; }
+      
+      tp0 = (mp_ptr) TMP_ALLOC ((rsize+rrsize) * BYTES_PER_MP_LIMB);
+      /* fill by zero rrsize low limbs of t */
+      MPN_ZERO(tp0, rrsize); tp = tp0 + rrsize; 
+      rp = (mp_ptr) TMP_ALLOC ((rrsize+1) * BYTES_PER_MP_LIMB); 
+
+      if (vsize >= rsize) { 
+	tmp = (mp_ptr) vp + vsize - rsize;
+      }
+      else { 
+	tmp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
+	MPN_COPY (tmp + rsize - vsize, vp, vsize);
+	MPN_ZERO (tmp, rsize - vsize); 
+      }
+
+      if (usize >= rsize) { 
+	MPN_COPY (tp, up + usize - rsize, rsize);
+      }
+      else {
+	MPN_COPY (tp + rsize - usize, up, usize); 
+	MPN_ZERO (tp, rsize - usize); 
+      }
+
+      /* Do the real job */
+ 
+#ifdef DEBUG
+      printf("Dividing : "); 
+      for(k = rsize - 1; k >= 0; k--) { printf("%lu ", tp[k]); }
+      printf(" by "); 
+      for(k = rsize - 1; k >= 0; k--) { printf("%lu ", tmp[k]); }
+      printf(".\n"); 
+#endif
+
+#if (__GNU_MP_VERSION < 3)
+      q_limb = mpn_divrem (rp, 0, tp0, rsize+rrsize, tmp, rsize);
+      tp = tp0; /* location of remainder */
+#else /* mpn_tdiv_qr is the preferred division interface in GMP 3 */
+      tmp2 = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
+      mpn_tdiv_qr(rp, tmp2, 0, tp0, rsize+rrsize, tmp, rsize);
+      q_limb = rp[rrsize];
+      tp = tmp2; /* location of remainder */
+#endif
+
+#ifdef DEBUG
+      printf("The result is : \n"); 
+      printf("Quotient : "); 
+      for(k = rrsize - 1; k >= 0; k--) { printf("%lu ", rp[k]); }
+      printf("Remainder : "); 
+      for(k = rsize - 1; k >= 0; k--) { printf("%lu ", tp[k]); }
+      printf("(q_limb = %lu)\n", q_limb); 
+#endif
+      
+      /* msb-normalize the result */
+      
+      if (q_limb)
+	{
+	  count_leading_zeros(k, q_limb); 
+	  mpn_rshift(rp, rp, rrsize, BITS_PER_MP_LIMB - k);
+	  rp[rrsize - 1] |= (q_limb << k);
+	  rexp += BITS_PER_MP_LIMB - k; 
+	}
+      else
+	{ 
+	  count_leading_zeros(k, rp[rrsize - 1]); 
+	  if (k) { mpn_lshift(rp, rp, rrsize, k); }
+	  rexp -= k; 
+	}
+
+      can_round = (mpfr_can_round_raw(rp, rrsize, sign_quotient, err, 
+				     GMP_RNDN, rnd_mode, MPFR_PREC(r))
+	|| (usize == rsize && vsize == rsize && 
+	    mpfr_can_round_raw(rp, rrsize, sign_quotient, err, 
+			       GMP_RNDZ, rnd_mode, MPFR_PREC(r))));
+
+      /* If we used all the limbs of both the dividend and the divisor, 
+	 then we have the correct RNDZ rounding */
+
+      if (!can_round && (rsize < usize || rsize < vsize)) 
+	{
+#ifdef DEBUG
+	  printf("Increasing the precision.\n");
+#endif
+	  TMP_FREE(marker);
+	  /* in case we can't round at the first iteration,
+	     jump right away to the maximum precision of both
+	     operands, to avoid multiple iterations when for 
+	     example the divisor is 1.0.
+	  */
+	  if (rsize < usize) rsize = usize - 1;
+	  if (rsize < vsize) rsize = vsize - 1;
+	}
+    }
+  while (!can_round && (rsize < usize || rsize < vsize) 
+	 && (rsize++) && (rrsize++));
+
+  /* ON PEUT PROBABLEMENT SE DEBROUILLER DES QUE rsize >= vsize */
+  /* MAIS IL FAUT AJOUTER LE BOUT QUI MANQUE DE usize A rsize */
+    
+  oldrrsize = rrsize;
+  rrsize = (MPFR_PREC(r) - 1)/BITS_PER_MP_LIMB + 1;
+
+  if (can_round) 
+    {
+      cc = mpfr_round_raw(rp, rp, err, (sign_quotient == -1 ? 1 : 0),
+			  MPFR_PREC(r), rnd_mode, NULL);
+    }
+  else {
+    /* Use the remainder to find out the correct rounding */
+    /* Note that at this point the division has been done */
+    /* EXACTLY. */
+    if ((rnd_mode == GMP_RNDD && sign_quotient == -1) 
+	|| (rnd_mode == GMP_RNDU && sign_quotient == 1)
+	|| (rnd_mode == GMP_RNDN))
+      {
+	/* We cannot round, so that the last bits of the quotient
+	   have to be zero; just look if the remainder is nonzero */
+	k = rsize - 1; 
+	while (k >= 0) { if (tp[k]) break; k--; }
+	if (k >= 0) 
+	  {
+	    t = MPFR_PREC(r) & (BITS_PER_MP_LIMB - 1); 
+	    if (t)
+	      {
+		cc = mpn_add_1(rp, rp, rrsize, 
+			       (mp_limb_t)1 << (BITS_PER_MP_LIMB - t)); 
+	      }
+	    else
+	      {
+		cc = mpn_add_1(rp, rp, rrsize, 1); 
+	      }
+	  }
+	else
+	  if (rnd_mode == GMP_RNDN) /* even rounding */
+	    {
+	      rw = (MPFR_PREC(r) + 1) & (BITS_PER_MP_LIMB - 1);
+	      if (rw) { rw = BITS_PER_MP_LIMB - rw; nw = 0; } else nw = 1; 
+	      if ((rw ? (rp[nw] >> (rw + 1)) & 1 : 
+		   (rp[nw] >> (BITS_PER_MP_LIMB - 1)) & 1))
+		{
+		  cc = mpn_add_1(rp + nw, rp + nw, rrsize, 
+				 ((mp_limb_t)1) << rw); 
+		}
+	    }
+	/* cas 0111111 */
+      }
+    /* Warning: we computed the result on oldrrsize limbs, but we want it
+       on rrsize limbs only. Both can differ, especially when the target
+       precision is a multiple of the number of bits per limb, since we've
+       taken an extra bit to make rounding to nearest easier. */
+    rp += oldrrsize-rrsize;
+  }
+
+
+  if (sign_quotient * MPFR_SIGN(r) < 0) { MPFR_CHANGE_SIGN(r); } 
+  MPFR_EXP(r) = rexp;
+
+  if (cc) {
+    mpn_rshift(rp, rp, rrsize, 1);
+    rp[rrsize-1] |= (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1);
+    MPFR_EXP(r)++; 
+  }
+    
+  rw = rrsize * BITS_PER_MP_LIMB - MPFR_PREC(r);
+  MPN_COPY(MPFR_MANT(r), rp, rrsize); 
+  MPFR_MANT(r)[0] &= ~(((mp_limb_t)1 << rw) - 1);
+
+  TMP_FREE (marker);
+}