diff --git a/Number_types/doc_tex/NumberTypeSupport_ref/MP_Float.tex b/Number_types/doc_tex/NumberTypeSupport_ref/MP_Float.tex index 99bf62fb1a7..58c74d6f608 100644 --- a/Number_types/doc_tex/NumberTypeSupport_ref/MP_Float.tex +++ b/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 diff --git a/Number_types/include/CGAL/MP_Float.h b/Number_types/include/CGAL/MP_Float.h index c6d0f09f292..adcf627fe42 100644 --- a/Number_types/include/CGAL/MP_Float.h +++ b/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 diff --git a/Number_types/test/Number_types/MP_Float.C b/Number_types/test/Number_types/MP_Float.C index ed215d309af..9dedc98f1a5 100644 --- a/Number_types/test/Number_types/MP_Float.C +++ b/Number_types/test/Number_types/MP_Float.C @@ -2,6 +2,7 @@ // Sylvain Pion. #include +#include #include #include #include @@ -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)