From 9c041567498e2e8ba5fac5a8c7d1a1d9bae7441a Mon Sep 17 00:00:00 2001 From: Dmitry Anisimov Date: Mon, 4 Oct 2021 16:29:14 +0200 Subject: [PATCH] removed duplicated code --- Number_types/include/CGAL/Exact_integer.h | 2 + Number_types/include/CGAL/boost_mp.h | 466 ++++++++-------------- 2 files changed, 169 insertions(+), 299 deletions(-) diff --git a/Number_types/include/CGAL/Exact_integer.h b/Number_types/include/CGAL/Exact_integer.h index 9e4ea3a8f5e..5b6e0988d7b 100644 --- a/Number_types/include/CGAL/Exact_integer.h +++ b/Number_types/include/CGAL/Exact_integer.h @@ -64,6 +64,8 @@ typedef leda_integer Exact_integer; typedef BOOST_cpp_arithmetic_kernel::Integer Exact_integer; #elif defined(CGAL_USE_CORE) typedef CORE::BigInt Exact_integer; +#else +#error "ERROR: Cannot determine a BigInt type!" #endif // CGAL_USE_CORE #endif // not DOXYGEN_RUNNING diff --git a/Number_types/include/CGAL/boost_mp.h b/Number_types/include/CGAL/boost_mp.h index e13c0520e62..d153dd8a6fb 100644 --- a/Number_types/include/CGAL/boost_mp.h +++ b/Number_types/include/CGAL/boost_mp.h @@ -141,6 +141,159 @@ struct Algebraic_structure_traits + my_ldexp( const Interval_nt& intv, const int e ) const { + + CGAL_assertion(intv.inf() > 0.0); + CGAL_assertion(intv.sup() > 0.0); + const double scale = std::ldexp(1.0, e); + return Interval_nt ( + CGAL_NTS is_finite(scale) ? + scale * intv.inf() : CGAL_IA_MAX_DOUBLE, + scale == 0.0 ? CGAL_IA_MIN_DOUBLE : scale * intv.sup() ); + } + + bool are_bounds_correct( const double l, const double u, const Type& x ) const { + + const double inf = std::numeric_limits::infinity(); + CGAL_assertion(u == l || u == std::nextafter(l, +inf)); + const bool are_bounds_tight = (u == l || u == std::nextafter(l, +inf)); + + if ( + CGAL::abs(l) == inf || + CGAL::abs(u) == inf || + CGAL::abs(l) == 0.0 || + CGAL::abs(u) == 0.0) { + return are_bounds_tight; + } + CGAL_assertion(CGAL::abs(l) != inf); + CGAL_assertion(CGAL::abs(u) != inf); + CGAL_assertion(CGAL::abs(l) != 0.0); + CGAL_assertion(CGAL::abs(u) != 0.0); + + const Type lb(l), ub(u); + CGAL_assertion(lb <= x); + CGAL_assertion(ub >= x); + const bool are_bounds_respected = (lb <= x && x <= ub); + return are_bounds_tight && are_bounds_respected; + } + + template + std::pair get_0ulp_interval( const int64_t shift, const ET& p ) const { + + CGAL_assertion(p >= 0); + const uint64_t pp = static_cast(p); + CGAL_assertion(pp >= 0); + const double pp_dbl = static_cast(pp); + const Interval_nt intv(pp_dbl, pp_dbl); + return my_ldexp(intv, -static_cast(shift)).pair(); + } + + template + std::pair get_1ulp_interval( const int64_t shift, const ET& p ) const { + + CGAL_assertion(p >= 0); + const uint64_t pp = static_cast(p); + const uint64_t qq = pp + 1; + CGAL_assertion(pp >= 0); + CGAL_assertion(qq > pp); + const double pp_dbl = static_cast(pp); + const double qq_dbl = static_cast(qq); + const Interval_nt intv(pp_dbl, qq_dbl); + return my_ldexp(intv, -static_cast(shift)).pair(); + } + + template + std::pair to_interval( ET xnum, ET xden ) const { + + CGAL_assertion(!CGAL::is_zero(xden)); + CGAL_assertion_code(const Type input = x); + double l = 0.0, u = 0.0; + if (CGAL::is_zero(xnum)) { // return [0.0, 0.0] + CGAL_assertion(are_bounds_correct(l, u, input)); + return std::make_pair(l, u); + } + CGAL_assertion(!CGAL::is_zero(xnum)); + + // Handle signs. + bool change_sign = false; + const bool is_num_pos = CGAL::is_positive(xnum); + const bool is_den_pos = CGAL::is_positive(xden); + if (!is_num_pos && !is_den_pos) { + xnum = -xnum; + xden = -xden; + } else if (!is_num_pos && is_den_pos) { + change_sign = true; + xnum = -xnum; + } else if (is_num_pos && !is_den_pos) { + change_sign = true; + xden = -xden; + } + CGAL_assertion(CGAL::is_positive(xnum) && CGAL::is_positive(xden)); + + const int64_t num_dbl_digits = std::numeric_limits::digits - 1; + const int64_t msb_num = static_cast(boost::multiprecision::msb(xnum)); + const int64_t msb_den = static_cast(boost::multiprecision::msb(xden)); + const int64_t msb_diff = msb_num - msb_den; + int64_t shift = num_dbl_digits - msb_diff; + + if (shift > 0) { + CGAL_assertion(msb_diff < num_dbl_digits); + xnum <<= +shift; + } else if (shift < 0) { + CGAL_assertion(msb_diff > num_dbl_digits); + xden <<= -shift; + } + CGAL_assertion(num_dbl_digits == + static_cast(boost::multiprecision::msb(xnum)) - + static_cast(boost::multiprecision::msb(xden))); + + decltype(xnum) p, r; + boost::multiprecision::divide_qr(xnum, xden, p, r); + const int64_t p_bits = static_cast(boost::multiprecision::msb(p)); + + if (r == 0) { + std::tie(l, u) = get_0ulp_interval(shift, p); + } else { + CGAL_assertion(r > 0); + CGAL_assertion(r < xden); + if (p_bits == num_dbl_digits - 1) { // we did not reach full precision + + p <<= 1; + r <<= 1; + ++shift; + + CGAL_assertion(r > 0); + const int cmp = r.compare(xden); + if (cmp > 0) { + ++p; + std::tie(l, u) = get_1ulp_interval(shift, p); + } else if (cmp == 0) { + ++p; + std::tie(l, u) = get_0ulp_interval(shift, p); + } else { + std::tie(l, u) = get_1ulp_interval(shift, p); + } + + } else { + std::tie(l, u) = get_1ulp_interval(shift, p); + } + } + + if (change_sign) { + const double t = l; + l = -u; + u = -t; + } + + CGAL_assertion(are_bounds_correct(l, u, input)); + return std::make_pair(l, u); + } + +} // Boost_MP_internal + template struct RET_boost_mp_base : public INTERN_RET::Real_embeddable_traits_base< NT , CGAL::Tag_true > { @@ -247,155 +400,11 @@ struct RET_boost_mp > { - bool are_bounds_correct( const double l, const double u, const Type& x ) const { - - const double inf = std::numeric_limits::infinity(); - CGAL_assertion(u == l || u == std::nextafter(l, +inf)); - const bool are_bounds_tight = (u == l || u == std::nextafter(l, +inf)); - - if ( - CGAL::abs(l) == inf || - CGAL::abs(u) == inf || - CGAL::abs(l) == 0.0 || - CGAL::abs(u) == 0.0) { - return are_bounds_tight; - } - CGAL_assertion(CGAL::abs(l) != inf); - CGAL_assertion(CGAL::abs(u) != inf); - CGAL_assertion(CGAL::abs(l) != 0.0); - CGAL_assertion(CGAL::abs(u) != 0.0); - - const Type lb(l), ub(u); - CGAL_assertion(lb <= x); - CGAL_assertion(ub >= x); - const bool are_bounds_respected = (lb <= x && x <= ub); - return are_bounds_tight && are_bounds_respected; - } - - Interval_nt - my_ldexp( const Interval_nt& intv, const int e ) const { - - CGAL_assertion(intv.inf() > 0.0); - CGAL_assertion(intv.sup() > 0.0); - const double scale = std::ldexp(1.0, e); - return Interval_nt ( - CGAL_NTS is_finite(scale) ? - scale * intv.inf() : CGAL_IA_MAX_DOUBLE, - scale == 0.0 ? CGAL_IA_MIN_DOUBLE : scale * intv.sup() ); - } - - template - std::pair get_0ulp_interval( const int64_t shift, const ET& p ) const { - - CGAL_assertion(p >= 0); - const uint64_t pp = static_cast(p); - CGAL_assertion(pp >= 0); - const double pp_dbl = static_cast(pp); - const Interval_nt intv(pp_dbl, pp_dbl); - return my_ldexp(intv, -static_cast(shift)).pair(); - } - - template - std::pair get_1ulp_interval( const int64_t shift, const ET& p ) const { - - CGAL_assertion(p >= 0); - const uint64_t pp = static_cast(p); - const uint64_t qq = pp + 1; - CGAL_assertion(pp >= 0); - CGAL_assertion(qq > pp); - const double pp_dbl = static_cast(pp); - const double qq_dbl = static_cast(qq); - const Interval_nt intv(pp_dbl, qq_dbl); - return my_ldexp(intv, -static_cast(shift)).pair(); - } - std::pair operator()( const Type& x ) const { - auto xnum = boost::multiprecision::numerator(x); - auto xden = boost::multiprecision::denominator(x); - - CGAL_assertion(!CGAL::is_zero(xden)); - CGAL_assertion_code(const Type input = x); - double l = 0.0, u = 0.0; - if (CGAL::is_zero(xnum)) { // return [0.0, 0.0] - CGAL_assertion(are_bounds_correct(l, u, input)); - return std::make_pair(l, u); - } - CGAL_assertion(!CGAL::is_zero(xnum)); - - // Handle signs. - bool change_sign = false; - const bool is_num_pos = CGAL::is_positive(xnum); - const bool is_den_pos = CGAL::is_positive(xden); - if (!is_num_pos && !is_den_pos) { - xnum = -xnum; - xden = -xden; - } else if (!is_num_pos && is_den_pos) { - change_sign = true; - xnum = -xnum; - } else if (is_num_pos && !is_den_pos) { - change_sign = true; - xden = -xden; - } - CGAL_assertion(CGAL::is_positive(xnum) && CGAL::is_positive(xden)); - - const int64_t num_dbl_digits = std::numeric_limits::digits - 1; - const int64_t msb_num = static_cast(boost::multiprecision::msb(xnum)); - const int64_t msb_den = static_cast(boost::multiprecision::msb(xden)); - const int64_t msb_diff = msb_num - msb_den; - int64_t shift = num_dbl_digits - msb_diff; - - if (shift > 0) { - CGAL_assertion(msb_diff < num_dbl_digits); - xnum <<= +shift; - } else if (shift < 0) { - CGAL_assertion(msb_diff > num_dbl_digits); - xden <<= -shift; - } - CGAL_assertion(num_dbl_digits == - static_cast(boost::multiprecision::msb(xnum)) - - static_cast(boost::multiprecision::msb(xden))); - - decltype(xnum) p, r; - boost::multiprecision::divide_qr(xnum, xden, p, r); - const int64_t p_bits = static_cast(boost::multiprecision::msb(p)); - - if (r == 0) { - std::tie(l, u) = get_0ulp_interval(shift, p); - } else { - CGAL_assertion(r > 0); - CGAL_assertion(r < xden); - if (p_bits == num_dbl_digits - 1) { // we did not reach full precision - - p <<= 1; - r <<= 1; - ++shift; - - CGAL_assertion(r > 0); - const int cmp = r.compare(xden); - if (cmp > 0) { - ++p; - std::tie(l, u) = get_1ulp_interval(shift, p); - } else if (cmp == 0) { - ++p; - std::tie(l, u) = get_0ulp_interval(shift, p); - } else { - std::tie(l, u) = get_1ulp_interval(shift, p); - } - - } else { - std::tie(l, u) = get_1ulp_interval(shift, p); - } - } - - if (change_sign) { - const double t = l; - l = -u; - u = -t; - } - - CGAL_assertion(are_bounds_correct(l, u, input)); - return std::make_pair(l, u); + const auto& xnum = boost::multiprecision::numerator(x); + const auto& xden = boost::multiprecision::denominator(x); + return Boost_MP_internal::to_interval(xnum, xden); } }; }; @@ -802,31 +811,6 @@ template< > class Real_embeddable_traits< Quotient > { public: - bool are_bounds_correct( const double l, const double u, const Type& x ) const { - - const double inf = std::numeric_limits::infinity(); - CGAL_assertion(u == l || u == std::nextafter(l, +inf)); - const bool are_bounds_tight = (u == l || u == std::nextafter(l, +inf)); - - if ( - CGAL::abs(l) == inf || - CGAL::abs(u) == inf || - CGAL::abs(l) == 0.0 || - CGAL::abs(u) == 0.0) { - return are_bounds_tight; - } - CGAL_assertion(CGAL::abs(l) != inf); - CGAL_assertion(CGAL::abs(u) != inf); - CGAL_assertion(CGAL::abs(l) != 0.0); - CGAL_assertion(CGAL::abs(u) != 0.0); - - const Type lb(l), ub(u); - CGAL_assertion(lb <= x); - CGAL_assertion(ub >= x); - const bool are_bounds_respected = (lb <= x && x <= ub); - return are_bounds_tight && are_bounds_respected; - } - /* // Option 1. // Inspired by the one from the gmpzf type. @@ -865,7 +849,8 @@ template< > class Real_embeddable_traits< Quotient class Real_embeddable_traits< Quotient class Real_embeddable_traits< Quotient class Real_embeddable_traits< Quotient quot = xn / xd; - CGAL_assertion(are_bounds_correct(quot.inf(), quot.sup(), x)); + CGAL_assertion(Boost_MP_internal:: + are_bounds_correct(quot.inf(), quot.sup(), x)); return std::make_pair(quot.inf(), quot.sup()); } */ @@ -972,129 +960,9 @@ template< > class Real_embeddable_traits< Quotient - my_ldexp( const Interval_nt& intv, const int e ) const { - - CGAL_assertion(intv.inf() > 0.0); - CGAL_assertion(intv.sup() > 0.0); - const double scale = std::ldexp(1.0, e); - return Interval_nt ( - CGAL_NTS is_finite(scale) ? - scale * intv.inf() : CGAL_IA_MAX_DOUBLE, - scale == 0.0 ? CGAL_IA_MIN_DOUBLE : scale * intv.sup() ); - } - - std::pair get_0ulp_interval( - const int64_t shift, - const boost::multiprecision::cpp_int& p ) const { - - CGAL_assertion(p >= 0); - const uint64_t pp = static_cast(p); - CGAL_assertion(pp >= 0); - const double pp_dbl = static_cast(pp); - const Interval_nt intv(pp_dbl, pp_dbl); - return my_ldexp(intv, -static_cast(shift)).pair(); - } - - std::pair get_1ulp_interval( - const int64_t shift, - const boost::multiprecision::cpp_int& p ) const { - - CGAL_assertion(p >= 0); - const uint64_t pp = static_cast(p); - const uint64_t qq = pp + 1; - CGAL_assertion(pp >= 0); - CGAL_assertion(qq > pp); - const double pp_dbl = static_cast(pp); - const double qq_dbl = static_cast(qq); - const Interval_nt intv(pp_dbl, qq_dbl); - return my_ldexp(intv, -static_cast(shift)).pair(); - } - - std::pair operator()( Type x ) const { - - CGAL_assertion(!CGAL::is_zero(x.den)); - CGAL_assertion_code(const Type input = x); - double l = 0.0, u = 0.0; - if (CGAL::is_zero(x.num)) { // return [0.0, 0.0] - CGAL_assertion(are_bounds_correct(l, u, input)); - return std::make_pair(l, u); - } - CGAL_assertion(!CGAL::is_zero(x.num)); - - // Handle signs. - bool change_sign = false; - const bool is_num_pos = CGAL::is_positive(x.num); - const bool is_den_pos = CGAL::is_positive(x.den); - if (!is_num_pos && !is_den_pos) { - x.num = -x.num; - x.den = -x.den; - } else if (!is_num_pos && is_den_pos) { - change_sign = true; - x.num = -x.num; - } else if (is_num_pos && !is_den_pos) { - change_sign = true; - x.den = -x.den; - } - CGAL_assertion(CGAL::is_positive(x.num) && CGAL::is_positive(x.den)); - - const int64_t num_dbl_digits = std::numeric_limits::digits - 1; - const int64_t msb_num = static_cast(boost::multiprecision::msb(x.num)); - const int64_t msb_den = static_cast(boost::multiprecision::msb(x.den)); - const int64_t msb_diff = msb_num - msb_den; - int64_t shift = num_dbl_digits - msb_diff; - - if (shift > 0) { - CGAL_assertion(msb_diff < num_dbl_digits); - x.num <<= +shift; - } else if (shift < 0) { - CGAL_assertion(msb_diff > num_dbl_digits); - x.den <<= -shift; - } - CGAL_assertion(num_dbl_digits == - static_cast(boost::multiprecision::msb(x.num)) - - static_cast(boost::multiprecision::msb(x.den))); - - boost::multiprecision::cpp_int p, r; - boost::multiprecision::divide_qr(x.num, x.den, p, r); - const int64_t p_bits = static_cast(boost::multiprecision::msb(p)); - - if (r == 0) { - std::tie(l, u) = get_0ulp_interval(shift, p); - } else { - CGAL_assertion(r > 0); - CGAL_assertion(r < x.den); - if (p_bits == num_dbl_digits - 1) { // we did not reach full precision - - p <<= 1; - r <<= 1; - ++shift; - - CGAL_assertion(r > 0); - const int cmp = r.compare(x.den); - if (cmp > 0) { - ++p; - std::tie(l, u) = get_1ulp_interval(shift, p); - } else if (cmp == 0) { - ++p; - std::tie(l, u) = get_0ulp_interval(shift, p); - } else { - std::tie(l, u) = get_1ulp_interval(shift, p); - } - - } else { - std::tie(l, u) = get_1ulp_interval(shift, p); - } - } - - if (change_sign) { - const double t = l; - l = -u; - u = -t; - } - - CGAL_assertion(are_bounds_correct(l, u, input)); - return std::make_pair(l, u); + // Option 2. Stable one! + std::pair operator()( const Type& x ) const { + return Boost_MP_internal::to_interval(x.num, x.den); } }; };