Add divides() function, and detect cases where exact_division()

is called when the division is not exact.
Also, rescale intermediate remainders.

Number_types/doc_tex/NumberTypeSupport_ref/MP_Float.tex

@@ -61,9 +61,12 @@ or \ccc{SqrtFieldNumberType} if \ccc{CGAL_MP_FLOAT_ALLOW_INEXACT} is set.
 \ccFunction{std::istream& operator>>(std::istream& in, MP_Float& m);}
 {reads a \ccc{double} from \ccc{in}, then converts it to an \ccc{MP_Float}.}
 
+\ccFunction{bool divides(const MP_Float &a, const MP_Float &b);}
+{returns true iff \ccc{b} divides \ccc{b}.}
+
 \ccFunction{MP_Float exact_division(const MP_Float &a, const MP_Float &b);}
 {computes the result of the division of \ccc{a} by \ccc{b}.
-\ccPrecond{\ccc{b} exactly divides \ccc{a}}.}
+\ccPrecond{divides(a, b)}.}
 
 \ccFunction{MP_Float approximate_division(const MP_Float &a, const MP_Float &b);}
 {computes an approximation of the division by converting the operands to

Number_types/include/CGAL/MP_Float.h

@@ -491,36 +491,79 @@ operator>> (std::istream & is, MP_Float &b);
 // TODO : needs function "bool divides(n, d)" for validity checking, and maybe other things.
 // TODO : needs function div(), with remainder.
 
+namespace CGALi {
+
 inline // Move it to libCGAL once it's stable.
 MP_Float
-exact_division(MP_Float n, MP_Float d)
+exact_division_internal(MP_Float n, MP_Float d, bool & divides)
 {
   CGAL_assertion(d != 0);
 
+  typedef MP_Float::exponent_type  exponent_type;
   // Rescale them to have to_double() values with reasonnable exponents.
-  MP_Float::exponent_type scale_n = n.find_scale();
-  MP_Float::exponent_type scale_d = d.find_scale();
+  exponent_type scale_n = n.find_scale();
+  exponent_type scale_d = d.find_scale();
   n.rescale(scale_n);
   d.rescale(scale_d);
 
+  // A simple criteria for detecting when the division is not exact
+  // is that if it is exact, then the result must have smaller or
+  // equal bit length than "n".
+  // Which we approximate with size() and a confortable margin.
+  exponent_type max_size_if_exact = n.size() + d.size() + 200;
+
   // School division algorithm.
   MP_Float res = to_double(n) / to_double(d);
   MP_Float remainder = n - res * d;
   while ( remainder != 0 )
   {
+    // We also have to rescale here, since remainder can diminish towards 0.
+    exponent_type scale_rem = remainder.find_scale();
+    remainder.rescale(scale_rem);
+    res.rescale(scale_rem);
+    scale_n += scale_rem;
+
+    // A double approximation of the quotient
+    // (imagine school division with base ~2^53).
     double approx = to_double(remainder) / to_double(d);
     CGAL_assertion(approx != 0);
     approx = (approx + (4*approx)) - (4*approx); // chop-off the last bit.
     res += approx;
     remainder -= approx * d;
-    // TODO : add detection when the division is not exact (abort).
+    if (res.size() > max_size_if_exact)
+    {
+      divides = false;
+      return MP_Float();
+    }
   }
 
+  divides = true;
   // Scale back the result.
   res.rescale(scale_d - scale_n);
   return res;
 }
 
+} // namespace CGALi
+
+inline // Move it to libCGAL once it's stable.
+MP_Float
+exact_division(const MP_Float & n, const MP_Float & d)
+{
+  bool exact_div;
+  MP_Float res = CGALi::exact_division_internal(n, d, exact_div);
+  CGAL_assertion_msg(exact_div, "exact_division() called with operands which do not divide");
+  return res;
+}
+
+inline // Move it to libCGAL once it's stable.
+bool
+divides(const MP_Float & n, const MP_Float & d)
+{
+  bool exact_div;
+  CGALi::exact_division_internal(n, d, exact_div);
+  return exact_div;
+}
+
 CGAL_END_NAMESPACE
 
 #endif // CGAL_MP_FLOAT_H

Number_types/test/Number_types/MP_Float.C

@@ -2,6 +2,7 @@
 // Sylvain Pion.
 
 #include <CGAL/basic.h>
+#include <CGAL/exceptions.h>
 #include <iostream>
 #include <cassert>
 #include <cstdlib>
@@ -21,6 +22,8 @@ void test_exact_division()
 
   // Let's pick 2 random values, multiply them, divide and check.
 
+  MPF tmp = exact_division(MPF(0), MPF(1));
+
   for (int i = 0; i < 100000; ++i) {
     MPF d = CGAL::default_random.get_double();
     MPF n = CGAL::default_random.get_double();
@@ -35,6 +38,10 @@ void test_exact_division()
     if ((i%3) < 1) n *= CGAL::default_random.get_double();
     if ((i%3) < 2) n *= CGAL::default_random.get_double();
 
+    // Try to change the signs
+    if ((i%7)  == 0) d = -d;
+    if ((i%11) == 0) n = -n;
+
     // Scale both numerator and denominator to test overflow/underflow
     // situations.
     d.rescale(107 * ((i%19) - 10));
@@ -46,9 +53,14 @@ void test_exact_division()
     //std::cout << "n = " << n << std::endl;
     assert( exact_division(nd, n) == d);
     assert( exact_division(nd, d) == n);
+    assert( divides(nd, n) );
+    assert( divides(nd, d) );
   }
-  // std::cout << "should loop..." << std::endl;
-  // MPF loop = exact_division(MPF(1), MPF(3));
+  
+  assert( ! divides(MPF(1), MPF(3)) );
+  assert( ! divides(MPF(2), MPF(7)) );
+  // test if we're lucky :)
+  assert( ! divides(MPF(CGAL::default_random.get_double()), MPF(CGAL::default_random.get_double())) );
 }
 
 void test_equality(int i)