songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

use custom ldexp for positive intervals

This commit is contained in:
Sébastien Loriot 2022-02-02 10:48:24 +01:00
parent 05d0cfff77
commit 07607223a2
1 changed files with 39 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

71
Number_types/include/CGAL/boost_mp.h
View File

@ -143,27 +143,28 @@ struct Algebraic_structure_traits<boost::multiprecision::detail::expression<T1,T
namespace Boost_MP_internal { namespace Boost_MP_internal {
// We use this function instead of ldexp overload from Interval_nt.h because: Interval_nt<false> shift_positive_interval( const Interval_nt<false>& intv, const int e ) {
// - see this issue: https://github.com/CGAL/cgal/issues/6004 CGAL_assertion(intv.inf() > 0.0);
// - in the original ldexp, when working with limit cases, we get CGAL_IA_MIN_DOUBLE CGAL_assertion(intv.sup() > 0.0);
// e.g. that after multiplication by scale turns into a value that is not minimal double, typedef std::numeric_limits<double> limits;
// which is expected for that limit case. So, we avoid multiplication by scale in this version
// for the limit cases and use it only for normal inf and sup cases.
// template<bool b> if (e < std::numeric_limits<double>::min_exponent)
// Interval_nt<b> my_ldexp( const Interval_nt<b>& intv, const int e ) { {
// CGAL_assertion(intv.inf() > 0.0); if ( e+std::numeric_limits<double>::digits < std::numeric_limits<double>::min_exponent)
// CGAL_assertion(intv.sup() > 0.0); return CGAL::Interval_nt<false>(0, (limits::min)());
// const double scale = std::ldexp(1.0, e); // following calls to ldexp call are exact (e+digits is larger or equal to min_exponent) is less than min_exponent
// Protect_FPU_rounding<true> P(CGAL_FE_UPWARD); return CGAL::Interval_nt<false>(std::ldexp(intv.inf(), e), std::ldexp(intv.sup(), e));
// return Interval_nt<b>( }
// CGAL_NTS is_finite(scale) ? if (e > limits::max_exponent)
// scale * intv.inf() : CGAL_IA_MAX_DOUBLE, return CGAL::Interval_nt<false>((limits::max)(), limits::infinity()); // no multiplication as intv is positive
// scale == 0.0 ? CGAL_IA_MIN_DOUBLE : scale * intv.sup());
// }
const double scale = std::ldexp(1.0, e); // ldexp call is exact (e is less than min_exponent)
return scale * intv;
}
/*
template<typename T> template<typename T>
bool is_spliiter_working(const T d) { bool is_splitter_working(const T d) {
T num = d; T num = d;
T den = 1.0; T den = 1.0;
while (std::ceil(num) != num) { while (std::ceil(num) != num) {
@ -173,40 +174,46 @@ namespace Boost_MP_internal {
if (d != num / den) return false; if (d != num / den) return false;
return true; return true;
} }
*/
// This function checks if the computed interval is correct and if it is tight. // This function checks if the computed interval is correct and if it is tight.
template<typename Type> template<typename Type>
bool are_bounds_correct( const double l, const double u, const Type& x ) { bool are_bounds_correct( const double l, const double u, const Type& x ) {
typedef std::numeric_limits<double> limits;
const double inf = std::numeric_limits<double>::infinity(); const double inf = std::numeric_limits<double>::infinity();
CGAL_assertion(u == l || u == std::nextafter(l, +inf)); if ( u!=l && (u==-inf || l==inf
|| (u==0 && l >= -(limits::min)())
|| (l==0 && u <= (limits::min)())) )
{
return x > Type((limits::max)()) ||
x < Type(-(limits::max)()) ||
(x > Type(-(limits::min)()) && x < Type((limits::min)()));
}
// CGAL_assertion(u == l || u == std::nextafter(l, +inf));
const bool are_bounds_tight = (u == l || u == std::nextafter(l, +inf)); const bool are_bounds_tight = (u == l || u == std::nextafter(l, +inf));
/*
// This check is required until a bug in double.h:split_numerator_denominator() is fixed. // This check is required until a bug in double.h:split_numerator_denominator() is fixed.
// For the moment, this splitter does not handle certain limit values. // For the moment, this splitter does not handle certain limit values.
// See also: https://github.com/CGAL/cgal/issues/5982 // See also: https://github.com/CGAL/cgal/issues/5982
// E.g. it fails for d = 2.752961027411077E-308! // E.g. it fails for d = 2.752961027411077E-308!
if (!is_spliiter_working(l) || !is_spliiter_working(u)) { if (!is_splitter_working(l) || !is_splitter_working(u)) {
return are_bounds_tight; return are_bounds_tight;
} }
CGAL_assertion(is_spliiter_working(l) && is_spliiter_working(u)); CGAL_assertion(is_spliiter_working(l) && is_spliiter_working(u));
*/
// We cannot convert inf to Type so we skip. // We cannot convert inf to Type so we skip.
if ( if (
CGAL::abs(l) == inf || CGAL::abs(l) == inf ||
CGAL::abs(u) == inf || CGAL::abs(u) == inf ||
CGAL::abs(l) == 0.0 || CGAL::abs(l) == 0.0 ||
CGAL::abs(u) == 0.0) { CGAL::abs(u) == 0.0)
{
return are_bounds_tight; 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); const Type lb(l), ub(u);
CGAL_assertion(lb <= x);
CGAL_assertion(ub >= x);
const bool are_bounds_respected = (lb <= x && x <= ub); const bool are_bounds_respected = (lb <= x && x <= ub);
return are_bounds_tight && are_bounds_respected; return are_bounds_tight && are_bounds_respected;
} }
@ -220,8 +227,8 @@ namespace Boost_MP_internal {
CGAL_assertion(pp >= 0); CGAL_assertion(pp >= 0);
const double pp_dbl = static_cast<double>(pp); const double pp_dbl = static_cast<double>(pp);
const Interval_nt<false> intv(pp_dbl, pp_dbl); const Interval_nt<false> intv(pp_dbl, pp_dbl);
Protect_FPU_rounding<true> P(CGAL_FE_UPWARD);
return CGAL::ldexp(intv, -static_cast<int>(shift)).pair(); return shift_positive_interval(intv, -static_cast<int>(shift)).pair();
} }
// This one returns 1 unit length interval. // This one returns 1 unit length interval.
@ -237,7 +244,7 @@ namespace Boost_MP_internal {
const double qq_dbl = static_cast<double>(qq); const double qq_dbl = static_cast<double>(qq);
const Interval_nt<false> intv(pp_dbl, qq_dbl); const Interval_nt<false> intv(pp_dbl, qq_dbl);
Protect_FPU_rounding<true> P(CGAL_FE_UPWARD); Protect_FPU_rounding<true> P(CGAL_FE_UPWARD);
return CGAL::ldexp(intv, -static_cast<int>(shift)).pair(); return shift_positive_interval(intv, -static_cast<int>(shift)).pair();
} }
// This is a version of to_interval that converts a rational type into a // This is a version of to_interval that converts a rational type into a