From 8f77b04a955d5a6532cf2645062d9b406fae1edd Mon Sep 17 00:00:00 2001 From: Sylvain Pion Date: Thu, 5 Aug 1999 18:29:01 +0000 Subject: [PATCH] - Merge with the static filters. - Added CGAL_assertion(Interval_nt_advanced::want_exceptions) for the dynamic. --- .../predicates/Regular_triangulation_ftC2.h | 371 ++++ .../predicates/Regular_triangulation_ftC3.h | 794 ++++++++ .../predicates/Regular_triangulation_rtH2.h | 457 +++++ .../predicates/Regular_triangulation_rtH3.h | 477 +++-- .../in_smallest_orthogonalcircle_ftC2.h | 147 ++ .../predicates/sign_of_determinant.h | 1049 ++++++++++ .../Arithmetic_filter/predicates_on_ftC2.h | 1784 +++++++++++++++++ .../Arithmetic_filter/predicates_on_ftC3.h | 1220 +++++++++++ 8 files changed, 6150 insertions(+), 149 deletions(-) diff --git a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_ftC2.h b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_ftC2.h index 16a2b1b9010..4c31a4a3062 100644 --- a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_ftC2.h +++ b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_ftC2.h @@ -26,6 +26,192 @@ CGAL_BEGIN_NAMESPACE +inline +Oriented_side +power_testC2_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pwt, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qwt, + const Static_filter_error &rx, + const Static_filter_error &ry, + const Static_filter_error &rwt, + const Static_filter_error &tx, + const Static_filter_error &ty, + const Static_filter_error &twt, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + + + + FT dpx = px - tx; + FT dpy = py - ty; + FT dpz = square(dpx) + square(dpy) - pwt + twt; + FT dqx = qx - tx; + FT dqy = qy - ty; + FT dqz = square(dqx) + square(dqy) - qwt + twt; + FT drx = rx - tx; + FT dry = ry - ty; + FT drz = square(drx) + square(dry) - rwt + twt; + + return Oriented_side(sign_of_determinant3x3_SAF(dpx, dpy, dpz, + dqx, dqy, dqz, + drx, dry, drz, + epsilon_0)); +} + +inline +Oriented_side +power_testC2_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pwt, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qwt, + const Restricted_double &rx, + const Restricted_double &ry, + const Restricted_double &rwt, + const Restricted_double &tx, + const Restricted_double &ty, + const Restricted_double &twt, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + + + + FT dpx = px - tx; + FT dpy = py - ty; + FT dpz = square(dpx) + square(dpy) - pwt + twt; + FT dqx = qx - tx; + FT dqy = qy - ty; + FT dqz = square(dqx) + square(dqy) - qwt + twt; + FT drx = rx - tx; + FT dry = ry - ty; + FT drz = square(drx) + square(dry) - rwt + twt; + + return Oriented_side(sign_of_determinant3x3_SAF(dpx, dpy, dpz, + dqx, dqy, dqz, + drx, dry, drz, + epsilon_0)); +} + +inline +Oriented_side +power_testC2( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pwt, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qwt, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry, + const Static_adaptatif_filter &rwt, + const Static_adaptatif_filter &tx, + const Static_adaptatif_filter &ty, + const Static_adaptatif_filter &twt) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pwt.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qwt.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound || + fabs(rwt.value()) > SAF_bound || + fabs(tx.value()) > SAF_bound || + fabs(ty.value()) > SAF_bound || + fabs(twt.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pwt.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qwt.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + SAF_bound = std::max(SAF_bound, fabs(rwt.value())); + SAF_bound = std::max(SAF_bound, fabs(tx.value())); + SAF_bound = std::max(SAF_bound, fabs(ty.value())); + SAF_bound = std::max(SAF_bound, fabs(twt.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) power_testC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return power_testC2_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pwt.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qwt.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + Restricted_double(rwt.value()), + Restricted_double(tx.value()), + Restricted_double(ty.value()), + Restricted_double(twt.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return power_testC2( + px.exact(), + py.exact(), + pwt.exact(), + qx.exact(), + qy.exact(), + qwt.exact(), + rx.exact(), + ry.exact(), + rwt.exact(), + tx.exact(), + ty.exact(), + twt.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -45,6 +231,7 @@ power_testC2( const Filtered_exact &ty, const Filtered_exact &twt) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -83,6 +270,189 @@ power_testC2( } } +inline +Oriented_side +power_testC2_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pwt, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qwt, + const Static_filter_error &tx, + const Static_filter_error &ty, + const Static_filter_error &twt, + double & epsilon_0, + double & epsilon_1, + double & epsilon_2, + double & epsilon_3) +{ + typedef Static_filter_error FT; + + + FT dpx = px - tx; + FT dpy = py - ty; + FT dpz = square(dpx) + square(dpy) - pwt + twt; + FT dqx = qx - tx; + FT dqy = qy - ty; + FT dqz = square(dqx) + square(dqy) - qwt + twt; + + + Comparison_result cmpx = CGAL::compare_SAF(px, qx, + epsilon_0); + if (cmpx != EQUAL) + return Oriented_side(cmpx * sign_of_determinant2x2_SAF(dpx, dpz, dqx, dqz, + epsilon_1)); + + + Comparison_result cmpy = CGAL::compare_SAF(py, qy, + epsilon_2); + return Oriented_side(cmpy * sign_of_determinant2x2_SAF(dpy, dpz, dqy, dqz, + epsilon_3)); +} + +inline +Oriented_side +power_testC2_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pwt, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qwt, + const Restricted_double &tx, + const Restricted_double &ty, + const Restricted_double &twt, + const double & epsilon_0, + const double & epsilon_1, + const double & epsilon_2, + const double & epsilon_3) +{ + typedef Restricted_double FT; + + + FT dpx = px - tx; + FT dpy = py - ty; + FT dpz = square(dpx) + square(dpy) - pwt + twt; + FT dqx = qx - tx; + FT dqy = qy - ty; + FT dqz = square(dqx) + square(dqy) - qwt + twt; + + + Comparison_result cmpx = CGAL::compare_SAF(px, qx, + epsilon_0); + if (cmpx != EQUAL) + return Oriented_side(cmpx * sign_of_determinant2x2_SAF(dpx, dpz, dqx, dqz, + epsilon_1)); + + + Comparison_result cmpy = CGAL::compare_SAF(py, qy, + epsilon_2); + return Oriented_side(cmpy * sign_of_determinant2x2_SAF(dpy, dpz, dqy, dqz, + epsilon_3)); +} + +inline +Oriented_side +power_testC2( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pwt, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qwt, + const Static_adaptatif_filter &tx, + const Static_adaptatif_filter &ty, + const Static_adaptatif_filter &twt) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + static double SAF_epsilon_2; + static double SAF_epsilon_3; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pwt.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qwt.value()) > SAF_bound || + fabs(tx.value()) > SAF_bound || + fabs(ty.value()) > SAF_bound || + fabs(twt.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pwt.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qwt.value())); + SAF_bound = std::max(SAF_bound, fabs(tx.value())); + SAF_bound = std::max(SAF_bound, fabs(ty.value())); + SAF_bound = std::max(SAF_bound, fabs(twt.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) power_testC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2, + SAF_epsilon_3); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return power_testC2_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pwt.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qwt.value()), + Restricted_double(tx.value()), + Restricted_double(ty.value()), + Restricted_double(twt.value()), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2, + SAF_epsilon_3); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return power_testC2( + px.exact(), + py.exact(), + pwt.exact(), + qx.exact(), + qy.exact(), + qwt.exact(), + tx.exact(), + ty.exact(), + twt.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -99,6 +469,7 @@ power_testC2( const Filtered_exact &ty, const Filtered_exact &twt) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { diff --git a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_ftC3.h b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_ftC3.h index d203ea9d791..db16c1b6b1f 100644 --- a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_ftC3.h +++ b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_ftC3.h @@ -26,6 +26,268 @@ CGAL_BEGIN_NAMESPACE +inline +Oriented_side +power_testC3_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pz, + const Static_filter_error &pwt, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qz, + const Static_filter_error &qwt, + const Static_filter_error &rx, + const Static_filter_error &ry, + const Static_filter_error &rz, + const Static_filter_error &rwt, + const Static_filter_error &sx, + const Static_filter_error &sy, + const Static_filter_error &sz, + const Static_filter_error &swt, + const Static_filter_error &tx, + const Static_filter_error &ty, + const Static_filter_error &tz, + const Static_filter_error &twt, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + + FT dpx = px - tx; + FT dpy = py - ty; + FT dpz = pz - tz; + FT dpt = square(dpx) + square(dpy) + square(dpz) - pwt + twt; + FT dqx = qx - tx; + FT dqy = qy - ty; + FT dqz = qz - tz; + FT dqt = square(dqx) + square(dqy) + square(dqz) - qwt + twt; + FT drx = rx - tx; + FT dry = ry - ty; + FT drz = rz - tz; + FT drt = square(drx) + square(dry) + square(drz) - rwt + twt; + FT dsx = sx - tx; + FT dsy = sy - ty; + FT dsz = sz - tz; + FT dst = square(dsx) + square(dsy) + square(dsz) - swt + twt; + + return Oriented_side( - sign_of_determinant4x4_SAF(dpx, dpy, dpz, dpt, + dqx, dqy, dqz, dqt, + drx, dry, drz, drt, + dsx, dsy, dsz, dst, + epsilon_0)); +} + +inline +Oriented_side +power_testC3_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pz, + const Restricted_double &pwt, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qz, + const Restricted_double &qwt, + const Restricted_double &rx, + const Restricted_double &ry, + const Restricted_double &rz, + const Restricted_double &rwt, + const Restricted_double &sx, + const Restricted_double &sy, + const Restricted_double &sz, + const Restricted_double &swt, + const Restricted_double &tx, + const Restricted_double &ty, + const Restricted_double &tz, + const Restricted_double &twt, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + + FT dpx = px - tx; + FT dpy = py - ty; + FT dpz = pz - tz; + FT dpt = square(dpx) + square(dpy) + square(dpz) - pwt + twt; + FT dqx = qx - tx; + FT dqy = qy - ty; + FT dqz = qz - tz; + FT dqt = square(dqx) + square(dqy) + square(dqz) - qwt + twt; + FT drx = rx - tx; + FT dry = ry - ty; + FT drz = rz - tz; + FT drt = square(drx) + square(dry) + square(drz) - rwt + twt; + FT dsx = sx - tx; + FT dsy = sy - ty; + FT dsz = sz - tz; + FT dst = square(dsx) + square(dsy) + square(dsz) - swt + twt; + + return Oriented_side( - sign_of_determinant4x4_SAF(dpx, dpy, dpz, dpt, + dqx, dqy, dqz, dqt, + drx, dry, drz, drt, + dsx, dsy, dsz, dst, + epsilon_0)); +} + +inline +Oriented_side +power_testC3( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pz, + const Static_adaptatif_filter &pwt, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qz, + const Static_adaptatif_filter &qwt, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry, + const Static_adaptatif_filter &rz, + const Static_adaptatif_filter &rwt, + const Static_adaptatif_filter &sx, + const Static_adaptatif_filter &sy, + const Static_adaptatif_filter &sz, + const Static_adaptatif_filter &swt, + const Static_adaptatif_filter &tx, + const Static_adaptatif_filter &ty, + const Static_adaptatif_filter &tz, + const Static_adaptatif_filter &twt) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pz.value()) > SAF_bound || + fabs(pwt.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qz.value()) > SAF_bound || + fabs(qwt.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound || + fabs(rz.value()) > SAF_bound || + fabs(rwt.value()) > SAF_bound || + fabs(sx.value()) > SAF_bound || + fabs(sy.value()) > SAF_bound || + fabs(sz.value()) > SAF_bound || + fabs(swt.value()) > SAF_bound || + fabs(tx.value()) > SAF_bound || + fabs(ty.value()) > SAF_bound || + fabs(tz.value()) > SAF_bound || + fabs(twt.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pz.value())); + SAF_bound = std::max(SAF_bound, fabs(pwt.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qz.value())); + SAF_bound = std::max(SAF_bound, fabs(qwt.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + SAF_bound = std::max(SAF_bound, fabs(rz.value())); + SAF_bound = std::max(SAF_bound, fabs(rwt.value())); + SAF_bound = std::max(SAF_bound, fabs(sx.value())); + SAF_bound = std::max(SAF_bound, fabs(sy.value())); + SAF_bound = std::max(SAF_bound, fabs(sz.value())); + SAF_bound = std::max(SAF_bound, fabs(swt.value())); + SAF_bound = std::max(SAF_bound, fabs(tx.value())); + SAF_bound = std::max(SAF_bound, fabs(ty.value())); + SAF_bound = std::max(SAF_bound, fabs(tz.value())); + SAF_bound = std::max(SAF_bound, fabs(twt.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) power_testC3_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return power_testC3_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pz.value()), + Restricted_double(pwt.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qz.value()), + Restricted_double(qwt.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + Restricted_double(rz.value()), + Restricted_double(rwt.value()), + Restricted_double(sx.value()), + Restricted_double(sy.value()), + Restricted_double(sz.value()), + Restricted_double(swt.value()), + Restricted_double(tx.value()), + Restricted_double(ty.value()), + Restricted_double(tz.value()), + Restricted_double(twt.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return power_testC3( + px.exact(), + py.exact(), + pz.exact(), + pwt.exact(), + qx.exact(), + qy.exact(), + qz.exact(), + qwt.exact(), + rx.exact(), + ry.exact(), + rz.exact(), + rwt.exact(), + sx.exact(), + sy.exact(), + sz.exact(), + swt.exact(), + tx.exact(), + ty.exact(), + tz.exact(), + twt.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -108,6 +370,301 @@ power_testC3( } } +inline +Oriented_side +power_testC3_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pz, + const Static_filter_error &pwt, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qz, + const Static_filter_error &qwt, + const Static_filter_error &rx, + const Static_filter_error &ry, + const Static_filter_error &rz, + const Static_filter_error &rwt, + const Static_filter_error &tx, + const Static_filter_error &ty, + const Static_filter_error &tz, + const Static_filter_error &twt, + double & epsilon_0, + double & epsilon_1, + double & epsilon_2, + double & epsilon_3, + double & epsilon_4, + double & epsilon_5) +{ + typedef Static_filter_error FT; + + + FT dpx = px - tx; + FT dpy = py - ty; + FT dpz = pz - tz; + FT dpt = square(dpx) + square(dpy) + square(dpz) - pwt + twt; + FT dqx = qx - tx; + FT dqy = qy - ty; + FT dqz = qz - tz; + FT dqt = square(dqx) + square(dqy) + square(dqz) - qwt + twt; + FT drx = rx - tx; + FT dry = ry - ty; + FT drz = rz - tz; + FT drt = square(drx) + square(dry) + square(drz) - rwt + twt; + Sign cmp; + + + cmp = sign_of_determinant3x3_SAF(dpx, dpy, dpt, + dqx, dqy, dqt, + drx, dry, drt, + epsilon_0); + if (cmp != ZERO) + return Oriented_side(cmp * sign_of_determinant2x2_SAF(px-rx, py-ry, + qx-rx, qy-ry, + epsilon_1)); + + + cmp = sign_of_determinant3x3_SAF(dpx, dpz, dpt, + dqx, dqz, dqt, + drx, drz, drt, + epsilon_2); + if (cmp != ZERO) + return Oriented_side(cmp * sign_of_determinant2x2_SAF(px-rx, pz-rz, + qx-rx, qz-rz, + epsilon_3)); + + + cmp = sign_of_determinant3x3_SAF(dpy, dpz, dpt, + dqy, dqz, dqt, + dry, drz, drt, + epsilon_4); + return Oriented_side(cmp * sign_of_determinant2x2_SAF(py-ry, pz-rz, + qy-ry, qz-rz, + epsilon_5)); +} + +inline +Oriented_side +power_testC3_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pz, + const Restricted_double &pwt, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qz, + const Restricted_double &qwt, + const Restricted_double &rx, + const Restricted_double &ry, + const Restricted_double &rz, + const Restricted_double &rwt, + const Restricted_double &tx, + const Restricted_double &ty, + const Restricted_double &tz, + const Restricted_double &twt, + const double & epsilon_0, + const double & epsilon_1, + const double & epsilon_2, + const double & epsilon_3, + const double & epsilon_4, + const double & epsilon_5) +{ + typedef Restricted_double FT; + + + FT dpx = px - tx; + FT dpy = py - ty; + FT dpz = pz - tz; + FT dpt = square(dpx) + square(dpy) + square(dpz) - pwt + twt; + FT dqx = qx - tx; + FT dqy = qy - ty; + FT dqz = qz - tz; + FT dqt = square(dqx) + square(dqy) + square(dqz) - qwt + twt; + FT drx = rx - tx; + FT dry = ry - ty; + FT drz = rz - tz; + FT drt = square(drx) + square(dry) + square(drz) - rwt + twt; + Sign cmp; + + + cmp = sign_of_determinant3x3_SAF(dpx, dpy, dpt, + dqx, dqy, dqt, + drx, dry, drt, + epsilon_0); + if (cmp != ZERO) + return Oriented_side(cmp * sign_of_determinant2x2_SAF(px-rx, py-ry, + qx-rx, qy-ry, + epsilon_1)); + + + cmp = sign_of_determinant3x3_SAF(dpx, dpz, dpt, + dqx, dqz, dqt, + drx, drz, drt, + epsilon_2); + if (cmp != ZERO) + return Oriented_side(cmp * sign_of_determinant2x2_SAF(px-rx, pz-rz, + qx-rx, qz-rz, + epsilon_3)); + + + cmp = sign_of_determinant3x3_SAF(dpy, dpz, dpt, + dqy, dqz, dqt, + dry, drz, drt, + epsilon_4); + return Oriented_side(cmp * sign_of_determinant2x2_SAF(py-ry, pz-rz, + qy-ry, qz-rz, + epsilon_5)); +} + +inline +Oriented_side +power_testC3( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pz, + const Static_adaptatif_filter &pwt, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qz, + const Static_adaptatif_filter &qwt, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry, + const Static_adaptatif_filter &rz, + const Static_adaptatif_filter &rwt, + const Static_adaptatif_filter &tx, + const Static_adaptatif_filter &ty, + const Static_adaptatif_filter &tz, + const Static_adaptatif_filter &twt) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + static double SAF_epsilon_2; + static double SAF_epsilon_3; + static double SAF_epsilon_4; + static double SAF_epsilon_5; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pz.value()) > SAF_bound || + fabs(pwt.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qz.value()) > SAF_bound || + fabs(qwt.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound || + fabs(rz.value()) > SAF_bound || + fabs(rwt.value()) > SAF_bound || + fabs(tx.value()) > SAF_bound || + fabs(ty.value()) > SAF_bound || + fabs(tz.value()) > SAF_bound || + fabs(twt.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pz.value())); + SAF_bound = std::max(SAF_bound, fabs(pwt.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qz.value())); + SAF_bound = std::max(SAF_bound, fabs(qwt.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + SAF_bound = std::max(SAF_bound, fabs(rz.value())); + SAF_bound = std::max(SAF_bound, fabs(rwt.value())); + SAF_bound = std::max(SAF_bound, fabs(tx.value())); + SAF_bound = std::max(SAF_bound, fabs(ty.value())); + SAF_bound = std::max(SAF_bound, fabs(tz.value())); + SAF_bound = std::max(SAF_bound, fabs(twt.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) power_testC3_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2, + SAF_epsilon_3, + SAF_epsilon_4, + SAF_epsilon_5); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return power_testC3_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pz.value()), + Restricted_double(pwt.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qz.value()), + Restricted_double(qwt.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + Restricted_double(rz.value()), + Restricted_double(rwt.value()), + Restricted_double(tx.value()), + Restricted_double(ty.value()), + Restricted_double(tz.value()), + Restricted_double(twt.value()), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2, + SAF_epsilon_3, + SAF_epsilon_4, + SAF_epsilon_5); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return power_testC3( + px.exact(), + py.exact(), + pz.exact(), + pwt.exact(), + qx.exact(), + qy.exact(), + qz.exact(), + qwt.exact(), + rx.exact(), + ry.exact(), + rz.exact(), + rwt.exact(), + tx.exact(), + ty.exact(), + tz.exact(), + twt.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -178,6 +735,243 @@ power_testC3( } } +inline +Oriented_side +power_testC3_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pz, + const Static_filter_error &pwt, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qz, + const Static_filter_error &qwt, + const Static_filter_error &tx, + const Static_filter_error &ty, + const Static_filter_error &tz, + const Static_filter_error &twt, + double & epsilon_0, + double & epsilon_1, + double & epsilon_2, + double & epsilon_3, + double & epsilon_4, + double & epsilon_5) +{ + typedef Static_filter_error FT; + + + FT dpx = px - tx; + FT dpy = py - ty; + FT dpz = pz - tz; + FT dpt = square(dpx) + square(dpy) + square(dpz) - pwt + twt; + FT dqx = qx - tx; + FT dqy = qy - ty; + FT dqz = qz - tz; + FT dqt = square(dqx) + square(dqy) + square(dqz) - qwt + twt; + Comparison_result cmp; + + + cmp = CGAL::compare_SAF(px, qx, + epsilon_0); + if (cmp != EQUAL) + return Oriented_side(cmp * sign_of_determinant2x2_SAF(dpx, dpt, dqx, dqt, + epsilon_1)); + + + cmp = CGAL::compare_SAF(py, qy, + epsilon_2); + if (cmp != EQUAL) + return Oriented_side(cmp * sign_of_determinant2x2_SAF(dpy, dpt, dqy, dqt, + epsilon_3)); + + + cmp = CGAL::compare_SAF(pz, qz, + epsilon_4); + return Oriented_side(cmp * sign_of_determinant2x2_SAF(dpz, dpt, dqz, dqt, + epsilon_5)); +} + +inline +Oriented_side +power_testC3_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pz, + const Restricted_double &pwt, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qz, + const Restricted_double &qwt, + const Restricted_double &tx, + const Restricted_double &ty, + const Restricted_double &tz, + const Restricted_double &twt, + const double & epsilon_0, + const double & epsilon_1, + const double & epsilon_2, + const double & epsilon_3, + const double & epsilon_4, + const double & epsilon_5) +{ + typedef Restricted_double FT; + + + FT dpx = px - tx; + FT dpy = py - ty; + FT dpz = pz - tz; + FT dpt = square(dpx) + square(dpy) + square(dpz) - pwt + twt; + FT dqx = qx - tx; + FT dqy = qy - ty; + FT dqz = qz - tz; + FT dqt = square(dqx) + square(dqy) + square(dqz) - qwt + twt; + Comparison_result cmp; + + + cmp = CGAL::compare_SAF(px, qx, + epsilon_0); + if (cmp != EQUAL) + return Oriented_side(cmp * sign_of_determinant2x2_SAF(dpx, dpt, dqx, dqt, + epsilon_1)); + + + cmp = CGAL::compare_SAF(py, qy, + epsilon_2); + if (cmp != EQUAL) + return Oriented_side(cmp * sign_of_determinant2x2_SAF(dpy, dpt, dqy, dqt, + epsilon_3)); + + + cmp = CGAL::compare_SAF(pz, qz, + epsilon_4); + return Oriented_side(cmp * sign_of_determinant2x2_SAF(dpz, dpt, dqz, dqt, + epsilon_5)); +} + +inline +Oriented_side +power_testC3( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pz, + const Static_adaptatif_filter &pwt, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qz, + const Static_adaptatif_filter &qwt, + const Static_adaptatif_filter &tx, + const Static_adaptatif_filter &ty, + const Static_adaptatif_filter &tz, + const Static_adaptatif_filter &twt) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + static double SAF_epsilon_2; + static double SAF_epsilon_3; + static double SAF_epsilon_4; + static double SAF_epsilon_5; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pz.value()) > SAF_bound || + fabs(pwt.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qz.value()) > SAF_bound || + fabs(qwt.value()) > SAF_bound || + fabs(tx.value()) > SAF_bound || + fabs(ty.value()) > SAF_bound || + fabs(tz.value()) > SAF_bound || + fabs(twt.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pz.value())); + SAF_bound = std::max(SAF_bound, fabs(pwt.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qz.value())); + SAF_bound = std::max(SAF_bound, fabs(qwt.value())); + SAF_bound = std::max(SAF_bound, fabs(tx.value())); + SAF_bound = std::max(SAF_bound, fabs(ty.value())); + SAF_bound = std::max(SAF_bound, fabs(tz.value())); + SAF_bound = std::max(SAF_bound, fabs(twt.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) power_testC3_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2, + SAF_epsilon_3, + SAF_epsilon_4, + SAF_epsilon_5); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return power_testC3_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pz.value()), + Restricted_double(pwt.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qz.value()), + Restricted_double(qwt.value()), + Restricted_double(tx.value()), + Restricted_double(ty.value()), + Restricted_double(tz.value()), + Restricted_double(twt.value()), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2, + SAF_epsilon_3, + SAF_epsilon_4, + SAF_epsilon_5); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return power_testC3( + px.exact(), + py.exact(), + pz.exact(), + pwt.exact(), + qx.exact(), + qy.exact(), + qz.exact(), + qwt.exact(), + tx.exact(), + ty.exact(), + tz.exact(), + twt.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif diff --git a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_rtH2.h b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_rtH2.h index 00041b13faf..953b63af8aa 100644 --- a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_rtH2.h +++ b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_rtH2.h @@ -26,6 +26,240 @@ CGAL_BEGIN_NAMESPACE +inline +Oriented_side +power_testH2_SAF( + const Static_filter_error &phx, + const Static_filter_error &phy, + const Static_filter_error &phw, + const Static_filter_error &pwt, + const Static_filter_error &qhx, + const Static_filter_error &qhy, + const Static_filter_error &qhw, + const Static_filter_error &qwt, + const Static_filter_error &rhx, + const Static_filter_error &rhy, + const Static_filter_error &rhw, + const Static_filter_error &rwt, + const Static_filter_error &thx, + const Static_filter_error &thy, + const Static_filter_error &thw, + const Static_filter_error &twt, + double & epsilon_0) +{ + typedef Static_filter_error RT; + + RT dphx = phx*phw; + RT dphy = phy*phw; + RT dphw = square(phw); + RT dpz = square(phx) + square(phy) - pwt*dphw; + + RT dqhx = qhx*qhw; + RT dqhy = qhy*qhw; + RT dqhw = square(qhw); + RT dqz = square(qhx) + square(qhy) - qwt*dqhw; + + RT drhx = rhx*rhw; + RT drhy = rhy*rhw; + RT drhw = square(rhw); + RT drz = square(rhx) + square(rhy) - rwt*drhw; + + RT dthx = thx*thw; + RT dthy = thy*thw; + RT dthw = square(thw); + RT dtz = square(thx) + square(thy) - twt*dthw; + + return Oriented_side(sign_of_determinant4x4_SAF(dphx, dphy, dpz, dphw, + dqhx, dqhy, dqz, dqhw, + drhx, drhy, drz, drhw, + dthx, dthy, dtz, dthw, + epsilon_0)); +} + +inline +Oriented_side +power_testH2_SAF( + const Restricted_double &phx, + const Restricted_double &phy, + const Restricted_double &phw, + const Restricted_double &pwt, + const Restricted_double &qhx, + const Restricted_double &qhy, + const Restricted_double &qhw, + const Restricted_double &qwt, + const Restricted_double &rhx, + const Restricted_double &rhy, + const Restricted_double &rhw, + const Restricted_double &rwt, + const Restricted_double &thx, + const Restricted_double &thy, + const Restricted_double &thw, + const Restricted_double &twt, + const double & epsilon_0) +{ + typedef Restricted_double RT; + + RT dphx = phx*phw; + RT dphy = phy*phw; + RT dphw = square(phw); + RT dpz = square(phx) + square(phy) - pwt*dphw; + + RT dqhx = qhx*qhw; + RT dqhy = qhy*qhw; + RT dqhw = square(qhw); + RT dqz = square(qhx) + square(qhy) - qwt*dqhw; + + RT drhx = rhx*rhw; + RT drhy = rhy*rhw; + RT drhw = square(rhw); + RT drz = square(rhx) + square(rhy) - rwt*drhw; + + RT dthx = thx*thw; + RT dthy = thy*thw; + RT dthw = square(thw); + RT dtz = square(thx) + square(thy) - twt*dthw; + + return Oriented_side(sign_of_determinant4x4_SAF(dphx, dphy, dpz, dphw, + dqhx, dqhy, dqz, dqhw, + drhx, drhy, drz, drhw, + dthx, dthy, dtz, dthw, + epsilon_0)); +} + +inline +Oriented_side +power_testH2( + const Static_adaptatif_filter &phx, + const Static_adaptatif_filter &phy, + const Static_adaptatif_filter &phw, + const Static_adaptatif_filter &pwt, + const Static_adaptatif_filter &qhx, + const Static_adaptatif_filter &qhy, + const Static_adaptatif_filter &qhw, + const Static_adaptatif_filter &qwt, + const Static_adaptatif_filter &rhx, + const Static_adaptatif_filter &rhy, + const Static_adaptatif_filter &rhw, + const Static_adaptatif_filter &rwt, + const Static_adaptatif_filter &thx, + const Static_adaptatif_filter &thy, + const Static_adaptatif_filter &thw, + const Static_adaptatif_filter &twt) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(phx.value()) > SAF_bound || + fabs(phy.value()) > SAF_bound || + fabs(phw.value()) > SAF_bound || + fabs(pwt.value()) > SAF_bound || + fabs(qhx.value()) > SAF_bound || + fabs(qhy.value()) > SAF_bound || + fabs(qhw.value()) > SAF_bound || + fabs(qwt.value()) > SAF_bound || + fabs(rhx.value()) > SAF_bound || + fabs(rhy.value()) > SAF_bound || + fabs(rhw.value()) > SAF_bound || + fabs(rwt.value()) > SAF_bound || + fabs(thx.value()) > SAF_bound || + fabs(thy.value()) > SAF_bound || + fabs(thw.value()) > SAF_bound || + fabs(twt.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(phx.value())); + SAF_bound = std::max(SAF_bound, fabs(phy.value())); + SAF_bound = std::max(SAF_bound, fabs(phw.value())); + SAF_bound = std::max(SAF_bound, fabs(pwt.value())); + SAF_bound = std::max(SAF_bound, fabs(qhx.value())); + SAF_bound = std::max(SAF_bound, fabs(qhy.value())); + SAF_bound = std::max(SAF_bound, fabs(qhw.value())); + SAF_bound = std::max(SAF_bound, fabs(qwt.value())); + SAF_bound = std::max(SAF_bound, fabs(rhx.value())); + SAF_bound = std::max(SAF_bound, fabs(rhy.value())); + SAF_bound = std::max(SAF_bound, fabs(rhw.value())); + SAF_bound = std::max(SAF_bound, fabs(rwt.value())); + SAF_bound = std::max(SAF_bound, fabs(thx.value())); + SAF_bound = std::max(SAF_bound, fabs(thy.value())); + SAF_bound = std::max(SAF_bound, fabs(thw.value())); + SAF_bound = std::max(SAF_bound, fabs(twt.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) power_testH2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return power_testH2_SAF( + Restricted_double(phx.value()), + Restricted_double(phy.value()), + Restricted_double(phw.value()), + Restricted_double(pwt.value()), + Restricted_double(qhx.value()), + Restricted_double(qhy.value()), + Restricted_double(qhw.value()), + Restricted_double(qwt.value()), + Restricted_double(rhx.value()), + Restricted_double(rhy.value()), + Restricted_double(rhw.value()), + Restricted_double(rwt.value()), + Restricted_double(thx.value()), + Restricted_double(thy.value()), + Restricted_double(thw.value()), + Restricted_double(twt.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return power_testH2( + phx.exact(), + phy.exact(), + phw.exact(), + pwt.exact(), + qhx.exact(), + qhy.exact(), + qhw.exact(), + qwt.exact(), + rhx.exact(), + rhy.exact(), + rhw.exact(), + rwt.exact(), + thx.exact(), + thy.exact(), + thw.exact(), + twt.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -49,6 +283,7 @@ power_testH2( const Filtered_exact &thw, const Filtered_exact &twt) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -95,6 +330,227 @@ power_testH2( } } +inline +Oriented_side +power_testH2_SAF( + const Static_filter_error &phx, + const Static_filter_error &phy, + const Static_filter_error &phw, + const Static_filter_error &pwt, + const Static_filter_error &qhx, + const Static_filter_error &qhy, + const Static_filter_error &qhw, + const Static_filter_error &qwt, + const Static_filter_error &thx, + const Static_filter_error &thy, + const Static_filter_error &thw, + const Static_filter_error &twt, + double & epsilon_0, + double & epsilon_1) +{ + typedef Static_filter_error RT; + + + + RT pa, qa, ta; + + if (phx * qhw != qhx * phw ) + { + pa = phx*phw; + qa = qhx*qhw; + ta = thx*thw; + } + else + { + pa = phy*phw; + qa = qhy*qhw; + ta = thy*thw; + } + + RT dphw = square(phw); + RT dpz = square(phx) + square(phy) - pwt*dphw; + + RT dqhw = square(qhw); + RT dqz = square(qhx) + square(qhy) - qwt*dqhw; + + RT dthw = square(thw); + RT dtz = square(thx) + square(thy) - twt*dthw; + + return Oriented_side(CGAL::compare_SAF(pa, qa, + epsilon_0) * + sign_of_determinant3x3_SAF(pa, dpz, dphw, + qa, dqz, dqhw, + ta, dtz, dthw, + epsilon_1)); +} + +inline +Oriented_side +power_testH2_SAF( + const Restricted_double &phx, + const Restricted_double &phy, + const Restricted_double &phw, + const Restricted_double &pwt, + const Restricted_double &qhx, + const Restricted_double &qhy, + const Restricted_double &qhw, + const Restricted_double &qwt, + const Restricted_double &thx, + const Restricted_double &thy, + const Restricted_double &thw, + const Restricted_double &twt, + const double & epsilon_0, + const double & epsilon_1) +{ + typedef Restricted_double RT; + + + + RT pa, qa, ta; + + if (phx * qhw != qhx * phw ) + { + pa = phx*phw; + qa = qhx*qhw; + ta = thx*thw; + } + else + { + pa = phy*phw; + qa = qhy*qhw; + ta = thy*thw; + } + + RT dphw = square(phw); + RT dpz = square(phx) + square(phy) - pwt*dphw; + + RT dqhw = square(qhw); + RT dqz = square(qhx) + square(qhy) - qwt*dqhw; + + RT dthw = square(thw); + RT dtz = square(thx) + square(thy) - twt*dthw; + + return Oriented_side(CGAL::compare_SAF(pa, qa, + epsilon_0) * + sign_of_determinant3x3_SAF(pa, dpz, dphw, + qa, dqz, dqhw, + ta, dtz, dthw, + epsilon_1)); +} + +inline +Oriented_side +power_testH2( + const Static_adaptatif_filter &phx, + const Static_adaptatif_filter &phy, + const Static_adaptatif_filter &phw, + const Static_adaptatif_filter &pwt, + const Static_adaptatif_filter &qhx, + const Static_adaptatif_filter &qhy, + const Static_adaptatif_filter &qhw, + const Static_adaptatif_filter &qwt, + const Static_adaptatif_filter &thx, + const Static_adaptatif_filter &thy, + const Static_adaptatif_filter &thw, + const Static_adaptatif_filter &twt) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(phx.value()) > SAF_bound || + fabs(phy.value()) > SAF_bound || + fabs(phw.value()) > SAF_bound || + fabs(pwt.value()) > SAF_bound || + fabs(qhx.value()) > SAF_bound || + fabs(qhy.value()) > SAF_bound || + fabs(qhw.value()) > SAF_bound || + fabs(qwt.value()) > SAF_bound || + fabs(thx.value()) > SAF_bound || + fabs(thy.value()) > SAF_bound || + fabs(thw.value()) > SAF_bound || + fabs(twt.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(phx.value())); + SAF_bound = std::max(SAF_bound, fabs(phy.value())); + SAF_bound = std::max(SAF_bound, fabs(phw.value())); + SAF_bound = std::max(SAF_bound, fabs(pwt.value())); + SAF_bound = std::max(SAF_bound, fabs(qhx.value())); + SAF_bound = std::max(SAF_bound, fabs(qhy.value())); + SAF_bound = std::max(SAF_bound, fabs(qhw.value())); + SAF_bound = std::max(SAF_bound, fabs(qwt.value())); + SAF_bound = std::max(SAF_bound, fabs(thx.value())); + SAF_bound = std::max(SAF_bound, fabs(thy.value())); + SAF_bound = std::max(SAF_bound, fabs(thw.value())); + SAF_bound = std::max(SAF_bound, fabs(twt.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) power_testH2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return power_testH2_SAF( + Restricted_double(phx.value()), + Restricted_double(phy.value()), + Restricted_double(phw.value()), + Restricted_double(pwt.value()), + Restricted_double(qhx.value()), + Restricted_double(qhy.value()), + Restricted_double(qhw.value()), + Restricted_double(qwt.value()), + Restricted_double(thx.value()), + Restricted_double(thy.value()), + Restricted_double(thw.value()), + Restricted_double(twt.value()), + SAF_epsilon_0, + SAF_epsilon_1); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return power_testH2( + phx.exact(), + phy.exact(), + phw.exact(), + pwt.exact(), + qhx.exact(), + qhy.exact(), + qhw.exact(), + qwt.exact(), + thx.exact(), + thy.exact(), + thw.exact(), + twt.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -114,6 +570,7 @@ power_testH2( const Filtered_exact &thw, const Filtered_exact &twt) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { diff --git a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_rtH3.h b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_rtH3.h index 1f02c74d93a..eed833c9fe3 100644 --- a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_rtH3.h +++ b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/Regular_triangulation_rtH3.h @@ -26,6 +26,334 @@ CGAL_BEGIN_NAMESPACE +inline +Oriented_side +power_testH3_SAF( + const Static_filter_error &phx, + const Static_filter_error &phy, + const Static_filter_error &phz, + const Static_filter_error &phw, + const Static_filter_error &pwt, + const Static_filter_error &qhx, + const Static_filter_error &qhy, + const Static_filter_error &qhz, + const Static_filter_error &qhw, + const Static_filter_error &qwt, + const Static_filter_error &rhx, + const Static_filter_error &rhy, + const Static_filter_error &rhz, + const Static_filter_error &rhw, + const Static_filter_error &rwt, + const Static_filter_error &shx, + const Static_filter_error ­, + const Static_filter_error &shz, + const Static_filter_error &shw, + const Static_filter_error &swt, + const Static_filter_error &thx, + const Static_filter_error &thy, + const Static_filter_error &thz, + const Static_filter_error &thw, + const Static_filter_error &twt, + double & epsilon_0) +{ + typedef Static_filter_error RT; + + RT dphx = phx*phw; + RT dphy = phy*phw; + RT dphz = phz*phw; + RT dphw = square(phw); + RT dpz = square(phx) + square(phy) + square(phz) - pwt*dphw; + + RT dqhx = qhx*qhw; + RT dqhy = qhy*qhw; + RT dqhz = qhz*qhw; + RT dqhw = square(qhw); + RT dqz = square(qhx) + square(qhy) + square(qhz) - qwt*dqhw; + + RT drhx = rhx*rhw; + RT drhy = rhy*rhw; + RT drhz = rhz*rhw; + RT drhw = square(rhw); + RT drz = square(rhx) + square(rhy) + square(rhz) - rwt*drhw; + + RT dshx = shx*shw; + RT dshy = shy*shw; + RT dshz = shz*shw; + RT dshw = square(shw); + RT dsz = square(shx) + square(shy) + square(shz) - swt*dshw; + + RT dthx = thx*thw; + RT dthy = thy*thw; + RT dthz = thz*thw; + RT dthw = square(thw); + RT dtz = square(thx) + square(thy) + square(thz) - twt*dthw; + + return Oriented_side(- sign_of_determinant5x5_SAF(dphx, dphy, dphz, dpz, dphw, + dqhx, dqhy, dqhz, dqz, dqhw, + drhx, drhy, drhz, drz, drhw, + dshx, dshy, dshz, dsz, dshw, + dthx, dthy, dthz, dtz, dthw, + epsilon_0)); +} + +inline +Oriented_side +power_testH3_SAF( + const Restricted_double &phx, + const Restricted_double &phy, + const Restricted_double &phz, + const Restricted_double &phw, + const Restricted_double &pwt, + const Restricted_double &qhx, + const Restricted_double &qhy, + const Restricted_double &qhz, + const Restricted_double &qhw, + const Restricted_double &qwt, + const Restricted_double &rhx, + const Restricted_double &rhy, + const Restricted_double &rhz, + const Restricted_double &rhw, + const Restricted_double &rwt, + const Restricted_double &shx, + const Restricted_double ­, + const Restricted_double &shz, + const Restricted_double &shw, + const Restricted_double &swt, + const Restricted_double &thx, + const Restricted_double &thy, + const Restricted_double &thz, + const Restricted_double &thw, + const Restricted_double &twt, + const double & epsilon_0) +{ + typedef Restricted_double RT; + + RT dphx = phx*phw; + RT dphy = phy*phw; + RT dphz = phz*phw; + RT dphw = square(phw); + RT dpz = square(phx) + square(phy) + square(phz) - pwt*dphw; + + RT dqhx = qhx*qhw; + RT dqhy = qhy*qhw; + RT dqhz = qhz*qhw; + RT dqhw = square(qhw); + RT dqz = square(qhx) + square(qhy) + square(qhz) - qwt*dqhw; + + RT drhx = rhx*rhw; + RT drhy = rhy*rhw; + RT drhz = rhz*rhw; + RT drhw = square(rhw); + RT drz = square(rhx) + square(rhy) + square(rhz) - rwt*drhw; + + RT dshx = shx*shw; + RT dshy = shy*shw; + RT dshz = shz*shw; + RT dshw = square(shw); + RT dsz = square(shx) + square(shy) + square(shz) - swt*dshw; + + RT dthx = thx*thw; + RT dthy = thy*thw; + RT dthz = thz*thw; + RT dthw = square(thw); + RT dtz = square(thx) + square(thy) + square(thz) - twt*dthw; + + return Oriented_side(- sign_of_determinant5x5_SAF(dphx, dphy, dphz, dpz, dphw, + dqhx, dqhy, dqhz, dqz, dqhw, + drhx, drhy, drhz, drz, drhw, + dshx, dshy, dshz, dsz, dshw, + dthx, dthy, dthz, dtz, dthw, + epsilon_0)); +} + +inline +Oriented_side +power_testH3( + const Static_adaptatif_filter &phx, + const Static_adaptatif_filter &phy, + const Static_adaptatif_filter &phz, + const Static_adaptatif_filter &phw, + const Static_adaptatif_filter &pwt, + const Static_adaptatif_filter &qhx, + const Static_adaptatif_filter &qhy, + const Static_adaptatif_filter &qhz, + const Static_adaptatif_filter &qhw, + const Static_adaptatif_filter &qwt, + const Static_adaptatif_filter &rhx, + const Static_adaptatif_filter &rhy, + const Static_adaptatif_filter &rhz, + const Static_adaptatif_filter &rhw, + const Static_adaptatif_filter &rwt, + const Static_adaptatif_filter &shx, + const Static_adaptatif_filter ­, + const Static_adaptatif_filter &shz, + const Static_adaptatif_filter &shw, + const Static_adaptatif_filter &swt, + const Static_adaptatif_filter &thx, + const Static_adaptatif_filter &thy, + const Static_adaptatif_filter &thz, + const Static_adaptatif_filter &thw, + const Static_adaptatif_filter &twt) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(phx.value()) > SAF_bound || + fabs(phy.value()) > SAF_bound || + fabs(phz.value()) > SAF_bound || + fabs(phw.value()) > SAF_bound || + fabs(pwt.value()) > SAF_bound || + fabs(qhx.value()) > SAF_bound || + fabs(qhy.value()) > SAF_bound || + fabs(qhz.value()) > SAF_bound || + fabs(qhw.value()) > SAF_bound || + fabs(qwt.value()) > SAF_bound || + fabs(rhx.value()) > SAF_bound || + fabs(rhy.value()) > SAF_bound || + fabs(rhz.value()) > SAF_bound || + fabs(rhw.value()) > SAF_bound || + fabs(rwt.value()) > SAF_bound || + fabs(shx.value()) > SAF_bound || + fabs(shy.value()) > SAF_bound || + fabs(shz.value()) > SAF_bound || + fabs(shw.value()) > SAF_bound || + fabs(swt.value()) > SAF_bound || + fabs(thx.value()) > SAF_bound || + fabs(thy.value()) > SAF_bound || + fabs(thz.value()) > SAF_bound || + fabs(thw.value()) > SAF_bound || + fabs(twt.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(phx.value())); + SAF_bound = std::max(SAF_bound, fabs(phy.value())); + SAF_bound = std::max(SAF_bound, fabs(phz.value())); + SAF_bound = std::max(SAF_bound, fabs(phw.value())); + SAF_bound = std::max(SAF_bound, fabs(pwt.value())); + SAF_bound = std::max(SAF_bound, fabs(qhx.value())); + SAF_bound = std::max(SAF_bound, fabs(qhy.value())); + SAF_bound = std::max(SAF_bound, fabs(qhz.value())); + SAF_bound = std::max(SAF_bound, fabs(qhw.value())); + SAF_bound = std::max(SAF_bound, fabs(qwt.value())); + SAF_bound = std::max(SAF_bound, fabs(rhx.value())); + SAF_bound = std::max(SAF_bound, fabs(rhy.value())); + SAF_bound = std::max(SAF_bound, fabs(rhz.value())); + SAF_bound = std::max(SAF_bound, fabs(rhw.value())); + SAF_bound = std::max(SAF_bound, fabs(rwt.value())); + SAF_bound = std::max(SAF_bound, fabs(shx.value())); + SAF_bound = std::max(SAF_bound, fabs(shy.value())); + SAF_bound = std::max(SAF_bound, fabs(shz.value())); + SAF_bound = std::max(SAF_bound, fabs(shw.value())); + SAF_bound = std::max(SAF_bound, fabs(swt.value())); + SAF_bound = std::max(SAF_bound, fabs(thx.value())); + SAF_bound = std::max(SAF_bound, fabs(thy.value())); + SAF_bound = std::max(SAF_bound, fabs(thz.value())); + SAF_bound = std::max(SAF_bound, fabs(thw.value())); + SAF_bound = std::max(SAF_bound, fabs(twt.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) power_testH3_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return power_testH3_SAF( + Restricted_double(phx.value()), + Restricted_double(phy.value()), + Restricted_double(phz.value()), + Restricted_double(phw.value()), + Restricted_double(pwt.value()), + Restricted_double(qhx.value()), + Restricted_double(qhy.value()), + Restricted_double(qhz.value()), + Restricted_double(qhw.value()), + Restricted_double(qwt.value()), + Restricted_double(rhx.value()), + Restricted_double(rhy.value()), + Restricted_double(rhz.value()), + Restricted_double(rhw.value()), + Restricted_double(rwt.value()), + Restricted_double(shx.value()), + Restricted_double(shy.value()), + Restricted_double(shz.value()), + Restricted_double(shw.value()), + Restricted_double(swt.value()), + Restricted_double(thx.value()), + Restricted_double(thy.value()), + Restricted_double(thz.value()), + Restricted_double(thw.value()), + Restricted_double(twt.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return power_testH3( + phx.exact(), + phy.exact(), + phz.exact(), + phw.exact(), + pwt.exact(), + qhx.exact(), + qhy.exact(), + qhz.exact(), + qhw.exact(), + qwt.exact(), + rhx.exact(), + rhy.exact(), + rhz.exact(), + rhw.exact(), + rwt.exact(), + shx.exact(), + shy.exact(), + shz.exact(), + shw.exact(), + swt.exact(), + thx.exact(), + thy.exact(), + thz.exact(), + thw.exact(), + twt.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -123,155 +451,6 @@ power_testH3( } } -#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION -template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > -#endif -/* */ -Oriented_side -power_testH3( - const Filtered_exact &phx, - const Filtered_exact &phy, - const Filtered_exact &phz, - const Filtered_exact &phw, - const Filtered_exact &pwt, - const Filtered_exact &qhx, - const Filtered_exact &qhy, - const Filtered_exact &qhz, - const Filtered_exact &qhw, - const Filtered_exact &qwt, - const Filtered_exact &rhx, - const Filtered_exact &rhy, - const Filtered_exact &rhz, - const Filtered_exact &rhw, - const Filtered_exact &rwt, - const Filtered_exact &thx, - const Filtered_exact &thy, - const Filtered_exact &thz, - const Filtered_exact &thw, - const Filtered_exact &twt) -{ - CGAL_assertion(Interval_nt_advanced::want_exceptions); - FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); - try - { - Oriented_side result = power_testH3( - phx.interval(), - phy.interval(), - phz.interval(), - phw.interval(), - pwt.interval(), - qhx.interval(), - qhy.interval(), - qhz.interval(), - qhw.interval(), - qwt.interval(), - rhx.interval(), - rhy.interval(), - rhz.interval(), - rhw.interval(), - rwt.interval(), - thx.interval(), - thy.interval(), - thz.interval(), - thw.interval(), - twt.interval()); - FPU_set_cw(backup); - return result; - } - catch (Interval_nt_advanced::unsafe_comparison) - { - FPU_set_cw(backup); - return power_testH3( - phx.exact(), - phy.exact(), - phz.exact(), - phw.exact(), - pwt.exact(), - qhx.exact(), - qhy.exact(), - qhz.exact(), - qhw.exact(), - qwt.exact(), - rhx.exact(), - rhy.exact(), - rhz.exact(), - rhw.exact(), - rwt.exact(), - thx.exact(), - thy.exact(), - thz.exact(), - thw.exact(), - twt.exact()); - } -} - -#ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION -template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > -#endif -/* */ -Oriented_side -power_testH3( - const Filtered_exact &phx, - const Filtered_exact &phy, - const Filtered_exact &phz, - const Filtered_exact &phw, - const Filtered_exact &pwt, - const Filtered_exact &qhx, - const Filtered_exact &qhy, - const Filtered_exact &qhz, - const Filtered_exact &qhw, - const Filtered_exact &qwt, - const Filtered_exact &thx, - const Filtered_exact &thy, - const Filtered_exact &thz, - const Filtered_exact &thw, - const Filtered_exact &twt) -{ - CGAL_assertion(Interval_nt_advanced::want_exceptions); - FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); - try - { - Oriented_side result = power_testH3( - phx.interval(), - phy.interval(), - phz.interval(), - phw.interval(), - pwt.interval(), - qhx.interval(), - qhy.interval(), - qhz.interval(), - qhw.interval(), - qwt.interval(), - thx.interval(), - thy.interval(), - thz.interval(), - thw.interval(), - twt.interval()); - FPU_set_cw(backup); - return result; - } - catch (Interval_nt_advanced::unsafe_comparison) - { - FPU_set_cw(backup); - return power_testH3( - phx.exact(), - phy.exact(), - phz.exact(), - phw.exact(), - pwt.exact(), - qhx.exact(), - qhy.exact(), - qhz.exact(), - qhw.exact(), - qwt.exact(), - thx.exact(), - thy.exact(), - thz.exact(), - thw.exact(), - twt.exact()); - } -} - CGAL_END_NAMESPACE #endif // CGAL_ARITHMETIC_FILTER_REGULAR_TRIANGULATION_RTH3_H diff --git a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/in_smallest_orthogonalcircle_ftC2.h b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/in_smallest_orthogonalcircle_ftC2.h index 8fc2f0a76ae..e99632a7f2c 100644 --- a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/in_smallest_orthogonalcircle_ftC2.h +++ b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/in_smallest_orthogonalcircle_ftC2.h @@ -26,6 +26,152 @@ CGAL_BEGIN_NAMESPACE +inline +Oriented_side +in_smallest_orthogonalcircleC2_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pw, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qw, + const Static_filter_error &tx, + const Static_filter_error &ty, + const Static_filter_error &tw, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + FT dpx = px-qx; + FT dpy = py-qy; + FT dtx = tx-qx; + FT dty = ty-qy; + FT dpz = square(dpx)+square(dpy); + + return Oriented_side (sign_SAF((square(dtx)+square(dty)-tw+qw)*dpz + -(dpz-pw+qw)*(dpx*dtx+dpy*dty), + epsilon_0)); +} + +inline +Oriented_side +in_smallest_orthogonalcircleC2_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pw, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qw, + const Restricted_double &tx, + const Restricted_double &ty, + const Restricted_double &tw, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + FT dpx = px-qx; + FT dpy = py-qy; + FT dtx = tx-qx; + FT dty = ty-qy; + FT dpz = square(dpx)+square(dpy); + + return Oriented_side (sign_SAF((square(dtx)+square(dty)-tw+qw)*dpz + -(dpz-pw+qw)*(dpx*dtx+dpy*dty), + epsilon_0)); +} + +inline +Oriented_side +in_smallest_orthogonalcircleC2( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pw, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qw, + const Static_adaptatif_filter &tx, + const Static_adaptatif_filter &ty, + const Static_adaptatif_filter &tw) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pw.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qw.value()) > SAF_bound || + fabs(tx.value()) > SAF_bound || + fabs(ty.value()) > SAF_bound || + fabs(tw.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pw.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qw.value())); + SAF_bound = std::max(SAF_bound, fabs(tx.value())); + SAF_bound = std::max(SAF_bound, fabs(ty.value())); + SAF_bound = std::max(SAF_bound, fabs(tw.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) in_smallest_orthogonalcircleC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return in_smallest_orthogonalcircleC2_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pw.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qw.value()), + Restricted_double(tx.value()), + Restricted_double(ty.value()), + Restricted_double(tw.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return in_smallest_orthogonalcircleC2( + px.exact(), + py.exact(), + pw.exact(), + qx.exact(), + qy.exact(), + qw.exact(), + tx.exact(), + ty.exact(), + tw.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -42,6 +188,7 @@ in_smallest_orthogonalcircleC2( const Filtered_exact &ty, const Filtered_exact &tw) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { diff --git a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/sign_of_determinant.h b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/sign_of_determinant.h index 5c394aa67fa..3e291bf9f0c 100644 --- a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/sign_of_determinant.h +++ b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates/sign_of_determinant.h @@ -26,6 +26,94 @@ CGAL_BEGIN_NAMESPACE +inline +Sign +sign_of_determinant2x2_SAF( + const Static_filter_error &a00, + const Static_filter_error &a01, + const Static_filter_error &a10, + const Static_filter_error &a11, + double & epsilon_0) +{ + typedef Static_filter_error FT; + return static_cast(static_cast(CGAL::compare_SAF( a00*a11, a10*a01, + epsilon_0))); } + +inline +Sign +sign_of_determinant2x2_SAF( + const Restricted_double &a00, + const Restricted_double &a01, + const Restricted_double &a10, + const Restricted_double &a11, + const double & epsilon_0) +{ + typedef Restricted_double FT; + return static_cast(static_cast(CGAL::compare_SAF( a00*a11, a10*a01, + epsilon_0))); } + +inline +Sign +sign_of_determinant2x2( + const Static_adaptatif_filter &a00, + const Static_adaptatif_filter &a01, + const Static_adaptatif_filter &a10, + const Static_adaptatif_filter &a11) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(a00.value()) > SAF_bound || + fabs(a01.value()) > SAF_bound || + fabs(a10.value()) > SAF_bound || + fabs(a11.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(a00.value())); + SAF_bound = std::max(SAF_bound, fabs(a01.value())); + SAF_bound = std::max(SAF_bound, fabs(a10.value())); + SAF_bound = std::max(SAF_bound, fabs(a11.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) sign_of_determinant2x2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return sign_of_determinant2x2_SAF( + Restricted_double(a00.value()), + Restricted_double(a01.value()), + Restricted_double(a10.value()), + Restricted_double(a11.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return sign_of_determinant2x2( + a00.exact(), + a01.exact(), + a10.exact(), + a11.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -37,6 +125,7 @@ sign_of_determinant2x2( const Filtered_exact &a10, const Filtered_exact &a11) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -59,6 +148,142 @@ sign_of_determinant2x2( } } +inline +Sign +sign_of_determinant3x3_SAF( + const Static_filter_error &a00, + const Static_filter_error &a01, + const Static_filter_error &a02, + const Static_filter_error &a10, + const Static_filter_error &a11, + const Static_filter_error &a12, + const Static_filter_error &a20, + const Static_filter_error &a21, + const Static_filter_error &a22, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return CGAL::sign_SAF(det3x3_by_formula(a00, a01, a02, + a10, a11, a12, + a20, a21, a22), + epsilon_0); +} + +inline +Sign +sign_of_determinant3x3_SAF( + const Restricted_double &a00, + const Restricted_double &a01, + const Restricted_double &a02, + const Restricted_double &a10, + const Restricted_double &a11, + const Restricted_double &a12, + const Restricted_double &a20, + const Restricted_double &a21, + const Restricted_double &a22, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return CGAL::sign_SAF(det3x3_by_formula(a00, a01, a02, + a10, a11, a12, + a20, a21, a22), + epsilon_0); +} + +inline +Sign +sign_of_determinant3x3( + const Static_adaptatif_filter &a00, + const Static_adaptatif_filter &a01, + const Static_adaptatif_filter &a02, + const Static_adaptatif_filter &a10, + const Static_adaptatif_filter &a11, + const Static_adaptatif_filter &a12, + const Static_adaptatif_filter &a20, + const Static_adaptatif_filter &a21, + const Static_adaptatif_filter &a22) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(a00.value()) > SAF_bound || + fabs(a01.value()) > SAF_bound || + fabs(a02.value()) > SAF_bound || + fabs(a10.value()) > SAF_bound || + fabs(a11.value()) > SAF_bound || + fabs(a12.value()) > SAF_bound || + fabs(a20.value()) > SAF_bound || + fabs(a21.value()) > SAF_bound || + fabs(a22.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(a00.value())); + SAF_bound = std::max(SAF_bound, fabs(a01.value())); + SAF_bound = std::max(SAF_bound, fabs(a02.value())); + SAF_bound = std::max(SAF_bound, fabs(a10.value())); + SAF_bound = std::max(SAF_bound, fabs(a11.value())); + SAF_bound = std::max(SAF_bound, fabs(a12.value())); + SAF_bound = std::max(SAF_bound, fabs(a20.value())); + SAF_bound = std::max(SAF_bound, fabs(a21.value())); + SAF_bound = std::max(SAF_bound, fabs(a22.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) sign_of_determinant3x3_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return sign_of_determinant3x3_SAF( + Restricted_double(a00.value()), + Restricted_double(a01.value()), + Restricted_double(a02.value()), + Restricted_double(a10.value()), + Restricted_double(a11.value()), + Restricted_double(a12.value()), + Restricted_double(a20.value()), + Restricted_double(a21.value()), + Restricted_double(a22.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return sign_of_determinant3x3( + a00.exact(), + a01.exact(), + a02.exact(), + a10.exact(), + a11.exact(), + a12.exact(), + a20.exact(), + a21.exact(), + a22.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -75,6 +300,7 @@ sign_of_determinant3x3( const Filtered_exact &a21, const Filtered_exact &a22) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -107,6 +333,200 @@ sign_of_determinant3x3( } } +inline +Sign +sign_of_determinant4x4_SAF( + const Static_filter_error &a00, + const Static_filter_error &a01, + const Static_filter_error &a02, + const Static_filter_error &a03, + const Static_filter_error &a10, + const Static_filter_error &a11, + const Static_filter_error &a12, + const Static_filter_error &a13, + const Static_filter_error &a20, + const Static_filter_error &a21, + const Static_filter_error &a22, + const Static_filter_error &a23, + const Static_filter_error &a30, + const Static_filter_error &a31, + const Static_filter_error &a32, + const Static_filter_error &a33, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return CGAL::sign_SAF(det4x4_by_formula(a00, a01, a02, a03, + a10, a11, a12, a13, + a20, a21, a22, a23, + a30, a31, a32, a33), + epsilon_0); +} + +inline +Sign +sign_of_determinant4x4_SAF( + const Restricted_double &a00, + const Restricted_double &a01, + const Restricted_double &a02, + const Restricted_double &a03, + const Restricted_double &a10, + const Restricted_double &a11, + const Restricted_double &a12, + const Restricted_double &a13, + const Restricted_double &a20, + const Restricted_double &a21, + const Restricted_double &a22, + const Restricted_double &a23, + const Restricted_double &a30, + const Restricted_double &a31, + const Restricted_double &a32, + const Restricted_double &a33, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return CGAL::sign_SAF(det4x4_by_formula(a00, a01, a02, a03, + a10, a11, a12, a13, + a20, a21, a22, a23, + a30, a31, a32, a33), + epsilon_0); +} + +inline +Sign +sign_of_determinant4x4( + const Static_adaptatif_filter &a00, + const Static_adaptatif_filter &a01, + const Static_adaptatif_filter &a02, + const Static_adaptatif_filter &a03, + const Static_adaptatif_filter &a10, + const Static_adaptatif_filter &a11, + const Static_adaptatif_filter &a12, + const Static_adaptatif_filter &a13, + const Static_adaptatif_filter &a20, + const Static_adaptatif_filter &a21, + const Static_adaptatif_filter &a22, + const Static_adaptatif_filter &a23, + const Static_adaptatif_filter &a30, + const Static_adaptatif_filter &a31, + const Static_adaptatif_filter &a32, + const Static_adaptatif_filter &a33) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(a00.value()) > SAF_bound || + fabs(a01.value()) > SAF_bound || + fabs(a02.value()) > SAF_bound || + fabs(a03.value()) > SAF_bound || + fabs(a10.value()) > SAF_bound || + fabs(a11.value()) > SAF_bound || + fabs(a12.value()) > SAF_bound || + fabs(a13.value()) > SAF_bound || + fabs(a20.value()) > SAF_bound || + fabs(a21.value()) > SAF_bound || + fabs(a22.value()) > SAF_bound || + fabs(a23.value()) > SAF_bound || + fabs(a30.value()) > SAF_bound || + fabs(a31.value()) > SAF_bound || + fabs(a32.value()) > SAF_bound || + fabs(a33.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(a00.value())); + SAF_bound = std::max(SAF_bound, fabs(a01.value())); + SAF_bound = std::max(SAF_bound, fabs(a02.value())); + SAF_bound = std::max(SAF_bound, fabs(a03.value())); + SAF_bound = std::max(SAF_bound, fabs(a10.value())); + SAF_bound = std::max(SAF_bound, fabs(a11.value())); + SAF_bound = std::max(SAF_bound, fabs(a12.value())); + SAF_bound = std::max(SAF_bound, fabs(a13.value())); + SAF_bound = std::max(SAF_bound, fabs(a20.value())); + SAF_bound = std::max(SAF_bound, fabs(a21.value())); + SAF_bound = std::max(SAF_bound, fabs(a22.value())); + SAF_bound = std::max(SAF_bound, fabs(a23.value())); + SAF_bound = std::max(SAF_bound, fabs(a30.value())); + SAF_bound = std::max(SAF_bound, fabs(a31.value())); + SAF_bound = std::max(SAF_bound, fabs(a32.value())); + SAF_bound = std::max(SAF_bound, fabs(a33.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) sign_of_determinant4x4_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return sign_of_determinant4x4_SAF( + Restricted_double(a00.value()), + Restricted_double(a01.value()), + Restricted_double(a02.value()), + Restricted_double(a03.value()), + Restricted_double(a10.value()), + Restricted_double(a11.value()), + Restricted_double(a12.value()), + Restricted_double(a13.value()), + Restricted_double(a20.value()), + Restricted_double(a21.value()), + Restricted_double(a22.value()), + Restricted_double(a23.value()), + Restricted_double(a30.value()), + Restricted_double(a31.value()), + Restricted_double(a32.value()), + Restricted_double(a33.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return sign_of_determinant4x4( + a00.exact(), + a01.exact(), + a02.exact(), + a03.exact(), + a10.exact(), + a11.exact(), + a12.exact(), + a13.exact(), + a20.exact(), + a21.exact(), + a22.exact(), + a23.exact(), + a30.exact(), + a31.exact(), + a32.exact(), + a33.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -130,6 +550,7 @@ sign_of_determinant4x4( const Filtered_exact &a32, const Filtered_exact &a33) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -176,6 +597,274 @@ sign_of_determinant4x4( } } +inline +Sign +sign_of_determinant5x5_SAF( + const Static_filter_error &a00, + const Static_filter_error &a01, + const Static_filter_error &a02, + const Static_filter_error &a03, + const Static_filter_error &a04, + const Static_filter_error &a10, + const Static_filter_error &a11, + const Static_filter_error &a12, + const Static_filter_error &a13, + const Static_filter_error &a14, + const Static_filter_error &a20, + const Static_filter_error &a21, + const Static_filter_error &a22, + const Static_filter_error &a23, + const Static_filter_error &a24, + const Static_filter_error &a30, + const Static_filter_error &a31, + const Static_filter_error &a32, + const Static_filter_error &a33, + const Static_filter_error &a34, + const Static_filter_error &a40, + const Static_filter_error &a41, + const Static_filter_error &a42, + const Static_filter_error &a43, + const Static_filter_error &a44, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return CGAL::sign_SAF(det5x5_by_formula(a00, a01, a02, a03, a04, + a10, a11, a12, a13, a14, + a20, a21, a22, a23, a24, + a30, a31, a32, a33, a34, + a40, a41, a42, a43, a44), + epsilon_0); +} + +inline +Sign +sign_of_determinant5x5_SAF( + const Restricted_double &a00, + const Restricted_double &a01, + const Restricted_double &a02, + const Restricted_double &a03, + const Restricted_double &a04, + const Restricted_double &a10, + const Restricted_double &a11, + const Restricted_double &a12, + const Restricted_double &a13, + const Restricted_double &a14, + const Restricted_double &a20, + const Restricted_double &a21, + const Restricted_double &a22, + const Restricted_double &a23, + const Restricted_double &a24, + const Restricted_double &a30, + const Restricted_double &a31, + const Restricted_double &a32, + const Restricted_double &a33, + const Restricted_double &a34, + const Restricted_double &a40, + const Restricted_double &a41, + const Restricted_double &a42, + const Restricted_double &a43, + const Restricted_double &a44, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return CGAL::sign_SAF(det5x5_by_formula(a00, a01, a02, a03, a04, + a10, a11, a12, a13, a14, + a20, a21, a22, a23, a24, + a30, a31, a32, a33, a34, + a40, a41, a42, a43, a44), + epsilon_0); +} + +inline +Sign +sign_of_determinant5x5( + const Static_adaptatif_filter &a00, + const Static_adaptatif_filter &a01, + const Static_adaptatif_filter &a02, + const Static_adaptatif_filter &a03, + const Static_adaptatif_filter &a04, + const Static_adaptatif_filter &a10, + const Static_adaptatif_filter &a11, + const Static_adaptatif_filter &a12, + const Static_adaptatif_filter &a13, + const Static_adaptatif_filter &a14, + const Static_adaptatif_filter &a20, + const Static_adaptatif_filter &a21, + const Static_adaptatif_filter &a22, + const Static_adaptatif_filter &a23, + const Static_adaptatif_filter &a24, + const Static_adaptatif_filter &a30, + const Static_adaptatif_filter &a31, + const Static_adaptatif_filter &a32, + const Static_adaptatif_filter &a33, + const Static_adaptatif_filter &a34, + const Static_adaptatif_filter &a40, + const Static_adaptatif_filter &a41, + const Static_adaptatif_filter &a42, + const Static_adaptatif_filter &a43, + const Static_adaptatif_filter &a44) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(a00.value()) > SAF_bound || + fabs(a01.value()) > SAF_bound || + fabs(a02.value()) > SAF_bound || + fabs(a03.value()) > SAF_bound || + fabs(a04.value()) > SAF_bound || + fabs(a10.value()) > SAF_bound || + fabs(a11.value()) > SAF_bound || + fabs(a12.value()) > SAF_bound || + fabs(a13.value()) > SAF_bound || + fabs(a14.value()) > SAF_bound || + fabs(a20.value()) > SAF_bound || + fabs(a21.value()) > SAF_bound || + fabs(a22.value()) > SAF_bound || + fabs(a23.value()) > SAF_bound || + fabs(a24.value()) > SAF_bound || + fabs(a30.value()) > SAF_bound || + fabs(a31.value()) > SAF_bound || + fabs(a32.value()) > SAF_bound || + fabs(a33.value()) > SAF_bound || + fabs(a34.value()) > SAF_bound || + fabs(a40.value()) > SAF_bound || + fabs(a41.value()) > SAF_bound || + fabs(a42.value()) > SAF_bound || + fabs(a43.value()) > SAF_bound || + fabs(a44.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(a00.value())); + SAF_bound = std::max(SAF_bound, fabs(a01.value())); + SAF_bound = std::max(SAF_bound, fabs(a02.value())); + SAF_bound = std::max(SAF_bound, fabs(a03.value())); + SAF_bound = std::max(SAF_bound, fabs(a04.value())); + SAF_bound = std::max(SAF_bound, fabs(a10.value())); + SAF_bound = std::max(SAF_bound, fabs(a11.value())); + SAF_bound = std::max(SAF_bound, fabs(a12.value())); + SAF_bound = std::max(SAF_bound, fabs(a13.value())); + SAF_bound = std::max(SAF_bound, fabs(a14.value())); + SAF_bound = std::max(SAF_bound, fabs(a20.value())); + SAF_bound = std::max(SAF_bound, fabs(a21.value())); + SAF_bound = std::max(SAF_bound, fabs(a22.value())); + SAF_bound = std::max(SAF_bound, fabs(a23.value())); + SAF_bound = std::max(SAF_bound, fabs(a24.value())); + SAF_bound = std::max(SAF_bound, fabs(a30.value())); + SAF_bound = std::max(SAF_bound, fabs(a31.value())); + SAF_bound = std::max(SAF_bound, fabs(a32.value())); + SAF_bound = std::max(SAF_bound, fabs(a33.value())); + SAF_bound = std::max(SAF_bound, fabs(a34.value())); + SAF_bound = std::max(SAF_bound, fabs(a40.value())); + SAF_bound = std::max(SAF_bound, fabs(a41.value())); + SAF_bound = std::max(SAF_bound, fabs(a42.value())); + SAF_bound = std::max(SAF_bound, fabs(a43.value())); + SAF_bound = std::max(SAF_bound, fabs(a44.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) sign_of_determinant5x5_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return sign_of_determinant5x5_SAF( + Restricted_double(a00.value()), + Restricted_double(a01.value()), + Restricted_double(a02.value()), + Restricted_double(a03.value()), + Restricted_double(a04.value()), + Restricted_double(a10.value()), + Restricted_double(a11.value()), + Restricted_double(a12.value()), + Restricted_double(a13.value()), + Restricted_double(a14.value()), + Restricted_double(a20.value()), + Restricted_double(a21.value()), + Restricted_double(a22.value()), + Restricted_double(a23.value()), + Restricted_double(a24.value()), + Restricted_double(a30.value()), + Restricted_double(a31.value()), + Restricted_double(a32.value()), + Restricted_double(a33.value()), + Restricted_double(a34.value()), + Restricted_double(a40.value()), + Restricted_double(a41.value()), + Restricted_double(a42.value()), + Restricted_double(a43.value()), + Restricted_double(a44.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return sign_of_determinant5x5( + a00.exact(), + a01.exact(), + a02.exact(), + a03.exact(), + a04.exact(), + a10.exact(), + a11.exact(), + a12.exact(), + a13.exact(), + a14.exact(), + a20.exact(), + a21.exact(), + a22.exact(), + a23.exact(), + a24.exact(), + a30.exact(), + a31.exact(), + a32.exact(), + a33.exact(), + a34.exact(), + a40.exact(), + a41.exact(), + a42.exact(), + a43.exact(), + a44.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -208,6 +897,7 @@ sign_of_determinant5x5( const Filtered_exact &a43, const Filtered_exact &a44) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -272,6 +962,364 @@ sign_of_determinant5x5( } } +inline +Sign +sign_of_determinant6x6_SAF( + const Static_filter_error &a00, + const Static_filter_error &a01, + const Static_filter_error &a02, + const Static_filter_error &a03, + const Static_filter_error &a04, + const Static_filter_error &a05, + const Static_filter_error &a10, + const Static_filter_error &a11, + const Static_filter_error &a12, + const Static_filter_error &a13, + const Static_filter_error &a14, + const Static_filter_error &a15, + const Static_filter_error &a20, + const Static_filter_error &a21, + const Static_filter_error &a22, + const Static_filter_error &a23, + const Static_filter_error &a24, + const Static_filter_error &a25, + const Static_filter_error &a30, + const Static_filter_error &a31, + const Static_filter_error &a32, + const Static_filter_error &a33, + const Static_filter_error &a34, + const Static_filter_error &a35, + const Static_filter_error &a40, + const Static_filter_error &a41, + const Static_filter_error &a42, + const Static_filter_error &a43, + const Static_filter_error &a44, + const Static_filter_error &a45, + const Static_filter_error &a50, + const Static_filter_error &a51, + const Static_filter_error &a52, + const Static_filter_error &a53, + const Static_filter_error &a54, + const Static_filter_error &a55, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return CGAL::sign_SAF(det6x6_by_formula(a00, a01, a02, a03, a04, a05, + a10, a11, a12, a13, a14, a15, + a20, a21, a22, a23, a24, a25, + a30, a31, a32, a33, a34, a35, + a40, a41, a42, a43, a44, a45, + a50, a51, a52, a53, a54, a55), + epsilon_0); +} + +inline +Sign +sign_of_determinant6x6_SAF( + const Restricted_double &a00, + const Restricted_double &a01, + const Restricted_double &a02, + const Restricted_double &a03, + const Restricted_double &a04, + const Restricted_double &a05, + const Restricted_double &a10, + const Restricted_double &a11, + const Restricted_double &a12, + const Restricted_double &a13, + const Restricted_double &a14, + const Restricted_double &a15, + const Restricted_double &a20, + const Restricted_double &a21, + const Restricted_double &a22, + const Restricted_double &a23, + const Restricted_double &a24, + const Restricted_double &a25, + const Restricted_double &a30, + const Restricted_double &a31, + const Restricted_double &a32, + const Restricted_double &a33, + const Restricted_double &a34, + const Restricted_double &a35, + const Restricted_double &a40, + const Restricted_double &a41, + const Restricted_double &a42, + const Restricted_double &a43, + const Restricted_double &a44, + const Restricted_double &a45, + const Restricted_double &a50, + const Restricted_double &a51, + const Restricted_double &a52, + const Restricted_double &a53, + const Restricted_double &a54, + const Restricted_double &a55, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return CGAL::sign_SAF(det6x6_by_formula(a00, a01, a02, a03, a04, a05, + a10, a11, a12, a13, a14, a15, + a20, a21, a22, a23, a24, a25, + a30, a31, a32, a33, a34, a35, + a40, a41, a42, a43, a44, a45, + a50, a51, a52, a53, a54, a55), + epsilon_0); +} + +inline +Sign +sign_of_determinant6x6( + const Static_adaptatif_filter &a00, + const Static_adaptatif_filter &a01, + const Static_adaptatif_filter &a02, + const Static_adaptatif_filter &a03, + const Static_adaptatif_filter &a04, + const Static_adaptatif_filter &a05, + const Static_adaptatif_filter &a10, + const Static_adaptatif_filter &a11, + const Static_adaptatif_filter &a12, + const Static_adaptatif_filter &a13, + const Static_adaptatif_filter &a14, + const Static_adaptatif_filter &a15, + const Static_adaptatif_filter &a20, + const Static_adaptatif_filter &a21, + const Static_adaptatif_filter &a22, + const Static_adaptatif_filter &a23, + const Static_adaptatif_filter &a24, + const Static_adaptatif_filter &a25, + const Static_adaptatif_filter &a30, + const Static_adaptatif_filter &a31, + const Static_adaptatif_filter &a32, + const Static_adaptatif_filter &a33, + const Static_adaptatif_filter &a34, + const Static_adaptatif_filter &a35, + const Static_adaptatif_filter &a40, + const Static_adaptatif_filter &a41, + const Static_adaptatif_filter &a42, + const Static_adaptatif_filter &a43, + const Static_adaptatif_filter &a44, + const Static_adaptatif_filter &a45, + const Static_adaptatif_filter &a50, + const Static_adaptatif_filter &a51, + const Static_adaptatif_filter &a52, + const Static_adaptatif_filter &a53, + const Static_adaptatif_filter &a54, + const Static_adaptatif_filter &a55) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(a00.value()) > SAF_bound || + fabs(a01.value()) > SAF_bound || + fabs(a02.value()) > SAF_bound || + fabs(a03.value()) > SAF_bound || + fabs(a04.value()) > SAF_bound || + fabs(a05.value()) > SAF_bound || + fabs(a10.value()) > SAF_bound || + fabs(a11.value()) > SAF_bound || + fabs(a12.value()) > SAF_bound || + fabs(a13.value()) > SAF_bound || + fabs(a14.value()) > SAF_bound || + fabs(a15.value()) > SAF_bound || + fabs(a20.value()) > SAF_bound || + fabs(a21.value()) > SAF_bound || + fabs(a22.value()) > SAF_bound || + fabs(a23.value()) > SAF_bound || + fabs(a24.value()) > SAF_bound || + fabs(a25.value()) > SAF_bound || + fabs(a30.value()) > SAF_bound || + fabs(a31.value()) > SAF_bound || + fabs(a32.value()) > SAF_bound || + fabs(a33.value()) > SAF_bound || + fabs(a34.value()) > SAF_bound || + fabs(a35.value()) > SAF_bound || + fabs(a40.value()) > SAF_bound || + fabs(a41.value()) > SAF_bound || + fabs(a42.value()) > SAF_bound || + fabs(a43.value()) > SAF_bound || + fabs(a44.value()) > SAF_bound || + fabs(a45.value()) > SAF_bound || + fabs(a50.value()) > SAF_bound || + fabs(a51.value()) > SAF_bound || + fabs(a52.value()) > SAF_bound || + fabs(a53.value()) > SAF_bound || + fabs(a54.value()) > SAF_bound || + fabs(a55.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(a00.value())); + SAF_bound = std::max(SAF_bound, fabs(a01.value())); + SAF_bound = std::max(SAF_bound, fabs(a02.value())); + SAF_bound = std::max(SAF_bound, fabs(a03.value())); + SAF_bound = std::max(SAF_bound, fabs(a04.value())); + SAF_bound = std::max(SAF_bound, fabs(a05.value())); + SAF_bound = std::max(SAF_bound, fabs(a10.value())); + SAF_bound = std::max(SAF_bound, fabs(a11.value())); + SAF_bound = std::max(SAF_bound, fabs(a12.value())); + SAF_bound = std::max(SAF_bound, fabs(a13.value())); + SAF_bound = std::max(SAF_bound, fabs(a14.value())); + SAF_bound = std::max(SAF_bound, fabs(a15.value())); + SAF_bound = std::max(SAF_bound, fabs(a20.value())); + SAF_bound = std::max(SAF_bound, fabs(a21.value())); + SAF_bound = std::max(SAF_bound, fabs(a22.value())); + SAF_bound = std::max(SAF_bound, fabs(a23.value())); + SAF_bound = std::max(SAF_bound, fabs(a24.value())); + SAF_bound = std::max(SAF_bound, fabs(a25.value())); + SAF_bound = std::max(SAF_bound, fabs(a30.value())); + SAF_bound = std::max(SAF_bound, fabs(a31.value())); + SAF_bound = std::max(SAF_bound, fabs(a32.value())); + SAF_bound = std::max(SAF_bound, fabs(a33.value())); + SAF_bound = std::max(SAF_bound, fabs(a34.value())); + SAF_bound = std::max(SAF_bound, fabs(a35.value())); + SAF_bound = std::max(SAF_bound, fabs(a40.value())); + SAF_bound = std::max(SAF_bound, fabs(a41.value())); + SAF_bound = std::max(SAF_bound, fabs(a42.value())); + SAF_bound = std::max(SAF_bound, fabs(a43.value())); + SAF_bound = std::max(SAF_bound, fabs(a44.value())); + SAF_bound = std::max(SAF_bound, fabs(a45.value())); + SAF_bound = std::max(SAF_bound, fabs(a50.value())); + SAF_bound = std::max(SAF_bound, fabs(a51.value())); + SAF_bound = std::max(SAF_bound, fabs(a52.value())); + SAF_bound = std::max(SAF_bound, fabs(a53.value())); + SAF_bound = std::max(SAF_bound, fabs(a54.value())); + SAF_bound = std::max(SAF_bound, fabs(a55.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) sign_of_determinant6x6_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return sign_of_determinant6x6_SAF( + Restricted_double(a00.value()), + Restricted_double(a01.value()), + Restricted_double(a02.value()), + Restricted_double(a03.value()), + Restricted_double(a04.value()), + Restricted_double(a05.value()), + Restricted_double(a10.value()), + Restricted_double(a11.value()), + Restricted_double(a12.value()), + Restricted_double(a13.value()), + Restricted_double(a14.value()), + Restricted_double(a15.value()), + Restricted_double(a20.value()), + Restricted_double(a21.value()), + Restricted_double(a22.value()), + Restricted_double(a23.value()), + Restricted_double(a24.value()), + Restricted_double(a25.value()), + Restricted_double(a30.value()), + Restricted_double(a31.value()), + Restricted_double(a32.value()), + Restricted_double(a33.value()), + Restricted_double(a34.value()), + Restricted_double(a35.value()), + Restricted_double(a40.value()), + Restricted_double(a41.value()), + Restricted_double(a42.value()), + Restricted_double(a43.value()), + Restricted_double(a44.value()), + Restricted_double(a45.value()), + Restricted_double(a50.value()), + Restricted_double(a51.value()), + Restricted_double(a52.value()), + Restricted_double(a53.value()), + Restricted_double(a54.value()), + Restricted_double(a55.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return sign_of_determinant6x6( + a00.exact(), + a01.exact(), + a02.exact(), + a03.exact(), + a04.exact(), + a05.exact(), + a10.exact(), + a11.exact(), + a12.exact(), + a13.exact(), + a14.exact(), + a15.exact(), + a20.exact(), + a21.exact(), + a22.exact(), + a23.exact(), + a24.exact(), + a25.exact(), + a30.exact(), + a31.exact(), + a32.exact(), + a33.exact(), + a34.exact(), + a35.exact(), + a40.exact(), + a41.exact(), + a42.exact(), + a43.exact(), + a44.exact(), + a45.exact(), + a50.exact(), + a51.exact(), + a52.exact(), + a53.exact(), + a54.exact(), + a55.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -315,6 +1363,7 @@ sign_of_determinant6x6( const Filtered_exact &a54, const Filtered_exact &a55) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { diff --git a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates_on_ftC2.h b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates_on_ftC2.h index 72d58ee54d3..42e99d9d4e8 100644 --- a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates_on_ftC2.h +++ b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates_on_ftC2.h @@ -26,6 +26,139 @@ CGAL_BEGIN_NAMESPACE +inline +Comparison_result +compare_xC2_SAF( + const Static_filter_error &px, + const Static_filter_error &l1a, + const Static_filter_error &l1b, + const Static_filter_error &l1c, + const Static_filter_error &l2a, + const Static_filter_error &l2b, + const Static_filter_error &l2c, + double & epsilon_0, + double & epsilon_1) +{ + typedef Static_filter_error FT; + + Sign sign1 = sign_of_determinant2x2_SAF(l1a, l1b, l2a, l2b, + epsilon_0); + Sign sign2 = sign_of_determinant3x3_SAF(l1a, l1b, l1c, + l2a, l2b, l2c, + -FT(1), FT(0), px, + epsilon_1); + CGAL_kernel_assertion( sign1 != 0 ); + return Comparison_result (sign1 * sign2); +} + +inline +Comparison_result +compare_xC2_SAF( + const Restricted_double &px, + const Restricted_double &l1a, + const Restricted_double &l1b, + const Restricted_double &l1c, + const Restricted_double &l2a, + const Restricted_double &l2b, + const Restricted_double &l2c, + const double & epsilon_0, + const double & epsilon_1) +{ + typedef Restricted_double FT; + + Sign sign1 = sign_of_determinant2x2_SAF(l1a, l1b, l2a, l2b, + epsilon_0); + Sign sign2 = sign_of_determinant3x3_SAF(l1a, l1b, l1c, + l2a, l2b, l2c, + -FT(1), FT(0), px, + epsilon_1); + CGAL_kernel_assertion( sign1 != 0 ); + return Comparison_result (sign1 * sign2); +} + +inline +Comparison_result +compare_xC2( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &l1a, + const Static_adaptatif_filter &l1b, + const Static_adaptatif_filter &l1c, + const Static_adaptatif_filter &l2a, + const Static_adaptatif_filter &l2b, + const Static_adaptatif_filter &l2c) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(l1a.value()) > SAF_bound || + fabs(l1b.value()) > SAF_bound || + fabs(l1c.value()) > SAF_bound || + fabs(l2a.value()) > SAF_bound || + fabs(l2b.value()) > SAF_bound || + fabs(l2c.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(l1a.value())); + SAF_bound = std::max(SAF_bound, fabs(l1b.value())); + SAF_bound = std::max(SAF_bound, fabs(l1c.value())); + SAF_bound = std::max(SAF_bound, fabs(l2a.value())); + SAF_bound = std::max(SAF_bound, fabs(l2b.value())); + SAF_bound = std::max(SAF_bound, fabs(l2c.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) compare_xC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return compare_xC2_SAF( + Restricted_double(px.value()), + Restricted_double(l1a.value()), + Restricted_double(l1b.value()), + Restricted_double(l1c.value()), + Restricted_double(l2a.value()), + Restricted_double(l2b.value()), + Restricted_double(l2c.value()), + SAF_epsilon_0, + SAF_epsilon_1); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return compare_xC2( + px.exact(), + l1a.exact(), + l1b.exact(), + l1c.exact(), + l2a.exact(), + l2b.exact(), + l2c.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -40,6 +173,7 @@ compare_xC2( const Filtered_exact &l2b, const Filtered_exact &l2c) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -68,6 +202,196 @@ compare_xC2( } } +inline +Comparison_result +compare_xC2_SAF( + const Static_filter_error &l1a, + const Static_filter_error &l1b, + const Static_filter_error &l1c, + const Static_filter_error &l2a, + const Static_filter_error &l2b, + const Static_filter_error &l2c, + const Static_filter_error &h1a, + const Static_filter_error &h1b, + const Static_filter_error &h1c, + const Static_filter_error &h2a, + const Static_filter_error &h2b, + const Static_filter_error &h2c, + double & epsilon_0, + double & epsilon_1, + double & epsilon_2) +{ + typedef Static_filter_error FT; + + Sign sign1 = sign_of_determinant2x2_SAF(l1a, l1b, l2a, l2b, + epsilon_0); + Sign sign2 = sign_of_determinant2x2_SAF(h1a, h1b, h2a, h2b, + epsilon_1); + + + FT FT0(0); + Sign sign3 = sign_of_determinant4x4_SAF(l1a, l1b, FT0, l1c, + l2a, l2b, FT0, l2c, + h1a, FT0, h1b, h1c, + h2a, FT0, h2b, h2c, + epsilon_2); + CGAL_kernel_assertion( (sign1 != 0) && (sign2 != 0) ); + return Comparison_result (- sign1 * sign2 * sign3); +} + +inline +Comparison_result +compare_xC2_SAF( + const Restricted_double &l1a, + const Restricted_double &l1b, + const Restricted_double &l1c, + const Restricted_double &l2a, + const Restricted_double &l2b, + const Restricted_double &l2c, + const Restricted_double &h1a, + const Restricted_double &h1b, + const Restricted_double &h1c, + const Restricted_double &h2a, + const Restricted_double &h2b, + const Restricted_double &h2c, + const double & epsilon_0, + const double & epsilon_1, + const double & epsilon_2) +{ + typedef Restricted_double FT; + + Sign sign1 = sign_of_determinant2x2_SAF(l1a, l1b, l2a, l2b, + epsilon_0); + Sign sign2 = sign_of_determinant2x2_SAF(h1a, h1b, h2a, h2b, + epsilon_1); + + + FT FT0(0); + Sign sign3 = sign_of_determinant4x4_SAF(l1a, l1b, FT0, l1c, + l2a, l2b, FT0, l2c, + h1a, FT0, h1b, h1c, + h2a, FT0, h2b, h2c, + epsilon_2); + CGAL_kernel_assertion( (sign1 != 0) && (sign2 != 0) ); + return Comparison_result (- sign1 * sign2 * sign3); +} + +inline +Comparison_result +compare_xC2( + const Static_adaptatif_filter &l1a, + const Static_adaptatif_filter &l1b, + const Static_adaptatif_filter &l1c, + const Static_adaptatif_filter &l2a, + const Static_adaptatif_filter &l2b, + const Static_adaptatif_filter &l2c, + const Static_adaptatif_filter &h1a, + const Static_adaptatif_filter &h1b, + const Static_adaptatif_filter &h1c, + const Static_adaptatif_filter &h2a, + const Static_adaptatif_filter &h2b, + const Static_adaptatif_filter &h2c) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + static double SAF_epsilon_2; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(l1a.value()) > SAF_bound || + fabs(l1b.value()) > SAF_bound || + fabs(l1c.value()) > SAF_bound || + fabs(l2a.value()) > SAF_bound || + fabs(l2b.value()) > SAF_bound || + fabs(l2c.value()) > SAF_bound || + fabs(h1a.value()) > SAF_bound || + fabs(h1b.value()) > SAF_bound || + fabs(h1c.value()) > SAF_bound || + fabs(h2a.value()) > SAF_bound || + fabs(h2b.value()) > SAF_bound || + fabs(h2c.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(l1a.value())); + SAF_bound = std::max(SAF_bound, fabs(l1b.value())); + SAF_bound = std::max(SAF_bound, fabs(l1c.value())); + SAF_bound = std::max(SAF_bound, fabs(l2a.value())); + SAF_bound = std::max(SAF_bound, fabs(l2b.value())); + SAF_bound = std::max(SAF_bound, fabs(l2c.value())); + SAF_bound = std::max(SAF_bound, fabs(h1a.value())); + SAF_bound = std::max(SAF_bound, fabs(h1b.value())); + SAF_bound = std::max(SAF_bound, fabs(h1c.value())); + SAF_bound = std::max(SAF_bound, fabs(h2a.value())); + SAF_bound = std::max(SAF_bound, fabs(h2b.value())); + SAF_bound = std::max(SAF_bound, fabs(h2c.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) compare_xC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return compare_xC2_SAF( + Restricted_double(l1a.value()), + Restricted_double(l1b.value()), + Restricted_double(l1c.value()), + Restricted_double(l2a.value()), + Restricted_double(l2b.value()), + Restricted_double(l2c.value()), + Restricted_double(h1a.value()), + Restricted_double(h1b.value()), + Restricted_double(h1c.value()), + Restricted_double(h2a.value()), + Restricted_double(h2b.value()), + Restricted_double(h2c.value()), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return compare_xC2( + l1a.exact(), + l1b.exact(), + l1c.exact(), + l2a.exact(), + l2b.exact(), + l2c.exact(), + h1a.exact(), + h1b.exact(), + h1c.exact(), + h2a.exact(), + h2b.exact(), + h2c.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -87,6 +411,7 @@ compare_xC2( const Filtered_exact &h2b, const Filtered_exact &h2c) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -125,6 +450,119 @@ compare_xC2( } } +inline +Comparison_result +compare_y_at_xC2_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &la, + const Static_filter_error &lb, + const Static_filter_error &lc, + double & epsilon_0, + double & epsilon_1) +{ + typedef Static_filter_error FT; + + Sign sign1 = CGAL::sign_SAF(lb, + epsilon_0); + Sign sign2 = CGAL::sign_SAF(la*px + lb*py + lc, + epsilon_1); + CGAL_kernel_assertion( sign1 != 0 ); + return Comparison_result (sign1 * sign2); +} + +inline +Comparison_result +compare_y_at_xC2_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &la, + const Restricted_double &lb, + const Restricted_double &lc, + const double & epsilon_0, + const double & epsilon_1) +{ + typedef Restricted_double FT; + + Sign sign1 = CGAL::sign_SAF(lb, + epsilon_0); + Sign sign2 = CGAL::sign_SAF(la*px + lb*py + lc, + epsilon_1); + CGAL_kernel_assertion( sign1 != 0 ); + return Comparison_result (sign1 * sign2); +} + +inline +Comparison_result +compare_y_at_xC2( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &la, + const Static_adaptatif_filter &lb, + const Static_adaptatif_filter &lc) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(la.value()) > SAF_bound || + fabs(lb.value()) > SAF_bound || + fabs(lc.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(la.value())); + SAF_bound = std::max(SAF_bound, fabs(lb.value())); + SAF_bound = std::max(SAF_bound, fabs(lc.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) compare_y_at_xC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return compare_y_at_xC2_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(la.value()), + Restricted_double(lb.value()), + Restricted_double(lc.value()), + SAF_epsilon_0, + SAF_epsilon_1); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return compare_y_at_xC2( + px.exact(), + py.exact(), + la.exact(), + lb.exact(), + lc.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -137,6 +575,7 @@ compare_y_at_xC2( const Filtered_exact &lb, const Filtered_exact &lc) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -161,6 +600,144 @@ compare_y_at_xC2( } } +inline +Comparison_result +compare_y_at_xC2_SAF( + const Static_filter_error &px, + const Static_filter_error &l1a, + const Static_filter_error &l1b, + const Static_filter_error &l1c, + const Static_filter_error &l2a, + const Static_filter_error &l2b, + const Static_filter_error &l2c, + double & epsilon_0, + double & epsilon_1, + double & epsilon_2) +{ + typedef Static_filter_error FT; + + Sign sign1 = CGAL::sign_SAF(l1b, + epsilon_0); + Sign sign2 = CGAL::sign_SAF(l2b, + epsilon_1); + Sign sign3 = sign_of_determinant2x2_SAF(l1a*px+l1c,l2a*px+l2c,l1b,l2b, + epsilon_2); + CGAL_kernel_assertion( (sign1 != 0) && (sign2 != 0) ); + return Comparison_result (- sign1 * sign2 * sign3); +} + +inline +Comparison_result +compare_y_at_xC2_SAF( + const Restricted_double &px, + const Restricted_double &l1a, + const Restricted_double &l1b, + const Restricted_double &l1c, + const Restricted_double &l2a, + const Restricted_double &l2b, + const Restricted_double &l2c, + const double & epsilon_0, + const double & epsilon_1, + const double & epsilon_2) +{ + typedef Restricted_double FT; + + Sign sign1 = CGAL::sign_SAF(l1b, + epsilon_0); + Sign sign2 = CGAL::sign_SAF(l2b, + epsilon_1); + Sign sign3 = sign_of_determinant2x2_SAF(l1a*px+l1c,l2a*px+l2c,l1b,l2b, + epsilon_2); + CGAL_kernel_assertion( (sign1 != 0) && (sign2 != 0) ); + return Comparison_result (- sign1 * sign2 * sign3); +} + +inline +Comparison_result +compare_y_at_xC2( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &l1a, + const Static_adaptatif_filter &l1b, + const Static_adaptatif_filter &l1c, + const Static_adaptatif_filter &l2a, + const Static_adaptatif_filter &l2b, + const Static_adaptatif_filter &l2c) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + static double SAF_epsilon_2; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(l1a.value()) > SAF_bound || + fabs(l1b.value()) > SAF_bound || + fabs(l1c.value()) > SAF_bound || + fabs(l2a.value()) > SAF_bound || + fabs(l2b.value()) > SAF_bound || + fabs(l2c.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(l1a.value())); + SAF_bound = std::max(SAF_bound, fabs(l1b.value())); + SAF_bound = std::max(SAF_bound, fabs(l1c.value())); + SAF_bound = std::max(SAF_bound, fabs(l2a.value())); + SAF_bound = std::max(SAF_bound, fabs(l2b.value())); + SAF_bound = std::max(SAF_bound, fabs(l2c.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) compare_y_at_xC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return compare_y_at_xC2_SAF( + Restricted_double(px.value()), + Restricted_double(l1a.value()), + Restricted_double(l1b.value()), + Restricted_double(l1c.value()), + Restricted_double(l2a.value()), + Restricted_double(l2b.value()), + Restricted_double(l2c.value()), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return compare_y_at_xC2( + px.exact(), + l1a.exact(), + l1b.exact(), + l1c.exact(), + l2a.exact(), + l2b.exact(), + l2c.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -175,6 +752,7 @@ compare_y_at_xC2( const Filtered_exact &l2b, const Filtered_exact &l2c) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -203,6 +781,165 @@ compare_y_at_xC2( } } +inline +Comparison_result +compare_y_at_xC2_SAF( + const Static_filter_error &l1a, + const Static_filter_error &l1b, + const Static_filter_error &l1c, + const Static_filter_error &l2a, + const Static_filter_error &l2b, + const Static_filter_error &l2c, + const Static_filter_error &ha, + const Static_filter_error &hb, + const Static_filter_error &hc, + double & epsilon_0, + double & epsilon_1, + double & epsilon_2, + double & epsilon_3) +{ + typedef Static_filter_error FT; + + Sign sign0 = sign_of_determinant2x2_SAF(l1a,l1b,l2a,l2b, + epsilon_0); + Sign sign1 = sign_of_determinant3x3_SAF(ha,hb,hc,l1a,l1b,l1c,l2a,l2b,l2c, + epsilon_1); + CGAL_kernel_assertion( (sign0 != ZERO) && (sign_SAF(hb, + epsilon_2) != ZERO) ); + return Comparison_result (sign0 * CGAL::sign_SAF(hb, + epsilon_3) * sign1); +} + +inline +Comparison_result +compare_y_at_xC2_SAF( + const Restricted_double &l1a, + const Restricted_double &l1b, + const Restricted_double &l1c, + const Restricted_double &l2a, + const Restricted_double &l2b, + const Restricted_double &l2c, + const Restricted_double &ha, + const Restricted_double &hb, + const Restricted_double &hc, + const double & epsilon_0, + const double & epsilon_1, + const double & epsilon_2, + const double & epsilon_3) +{ + typedef Restricted_double FT; + + Sign sign0 = sign_of_determinant2x2_SAF(l1a,l1b,l2a,l2b, + epsilon_0); + Sign sign1 = sign_of_determinant3x3_SAF(ha,hb,hc,l1a,l1b,l1c,l2a,l2b,l2c, + epsilon_1); + CGAL_kernel_assertion( (sign0 != ZERO) && (sign_SAF(hb, + epsilon_2) != ZERO) ); + return Comparison_result (sign0 * CGAL::sign_SAF(hb, + epsilon_3) * sign1); +} + +inline +Comparison_result +compare_y_at_xC2( + const Static_adaptatif_filter &l1a, + const Static_adaptatif_filter &l1b, + const Static_adaptatif_filter &l1c, + const Static_adaptatif_filter &l2a, + const Static_adaptatif_filter &l2b, + const Static_adaptatif_filter &l2c, + const Static_adaptatif_filter &ha, + const Static_adaptatif_filter &hb, + const Static_adaptatif_filter &hc) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + static double SAF_epsilon_2; + static double SAF_epsilon_3; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(l1a.value()) > SAF_bound || + fabs(l1b.value()) > SAF_bound || + fabs(l1c.value()) > SAF_bound || + fabs(l2a.value()) > SAF_bound || + fabs(l2b.value()) > SAF_bound || + fabs(l2c.value()) > SAF_bound || + fabs(ha.value()) > SAF_bound || + fabs(hb.value()) > SAF_bound || + fabs(hc.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(l1a.value())); + SAF_bound = std::max(SAF_bound, fabs(l1b.value())); + SAF_bound = std::max(SAF_bound, fabs(l1c.value())); + SAF_bound = std::max(SAF_bound, fabs(l2a.value())); + SAF_bound = std::max(SAF_bound, fabs(l2b.value())); + SAF_bound = std::max(SAF_bound, fabs(l2c.value())); + SAF_bound = std::max(SAF_bound, fabs(ha.value())); + SAF_bound = std::max(SAF_bound, fabs(hb.value())); + SAF_bound = std::max(SAF_bound, fabs(hc.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) compare_y_at_xC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2, + SAF_epsilon_3); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return compare_y_at_xC2_SAF( + Restricted_double(l1a.value()), + Restricted_double(l1b.value()), + Restricted_double(l1c.value()), + Restricted_double(l2a.value()), + Restricted_double(l2b.value()), + Restricted_double(l2c.value()), + Restricted_double(ha.value()), + Restricted_double(hb.value()), + Restricted_double(hc.value()), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2, + SAF_epsilon_3); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return compare_y_at_xC2( + l1a.exact(), + l1b.exact(), + l1c.exact(), + l2a.exact(), + l2b.exact(), + l2c.exact(), + ha.exact(), + hb.exact(), + hc.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -219,6 +956,7 @@ compare_y_at_xC2( const Filtered_exact &hb, const Filtered_exact &hc) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -251,6 +989,199 @@ compare_y_at_xC2( } } +inline +Comparison_result +compare_y_at_xC2_SAF( + const Static_filter_error &l1a, + const Static_filter_error &l1b, + const Static_filter_error &l1c, + const Static_filter_error &l2a, + const Static_filter_error &l2b, + const Static_filter_error &l2c, + const Static_filter_error &h1a, + const Static_filter_error &h1b, + const Static_filter_error &h1c, + const Static_filter_error &h2a, + const Static_filter_error &h2b, + const Static_filter_error &h2c, + double & epsilon_0, + double & epsilon_1, + double & epsilon_2, + double & epsilon_3) +{ + typedef Static_filter_error FT; + + FT FT0(0); + Sign s1 = lexicographical_sign_SAF(h1b, -h1a, + epsilon_0); + Sign s2 = lexicographical_sign_SAF(h2b, -h2a, + epsilon_1); + Sign s3 = sign_of_determinant2x2_SAF(l1a, l1b, l2a, l2b, + epsilon_2); + Sign s4 = sign_of_determinant4x4_SAF(h2a, h2b, FT0, h2c, + l1a, FT0, l1b, l1c, + l2a, FT0, l2b, l2c, + h1a, h1b, FT0, h1c, + epsilon_3); + return Comparison_result (s1 * s2 * s3 * s4); +} + +inline +Comparison_result +compare_y_at_xC2_SAF( + const Restricted_double &l1a, + const Restricted_double &l1b, + const Restricted_double &l1c, + const Restricted_double &l2a, + const Restricted_double &l2b, + const Restricted_double &l2c, + const Restricted_double &h1a, + const Restricted_double &h1b, + const Restricted_double &h1c, + const Restricted_double &h2a, + const Restricted_double &h2b, + const Restricted_double &h2c, + const double & epsilon_0, + const double & epsilon_1, + const double & epsilon_2, + const double & epsilon_3) +{ + typedef Restricted_double FT; + + FT FT0(0); + Sign s1 = lexicographical_sign_SAF(h1b, -h1a, + epsilon_0); + Sign s2 = lexicographical_sign_SAF(h2b, -h2a, + epsilon_1); + Sign s3 = sign_of_determinant2x2_SAF(l1a, l1b, l2a, l2b, + epsilon_2); + Sign s4 = sign_of_determinant4x4_SAF(h2a, h2b, FT0, h2c, + l1a, FT0, l1b, l1c, + l2a, FT0, l2b, l2c, + h1a, h1b, FT0, h1c, + epsilon_3); + return Comparison_result (s1 * s2 * s3 * s4); +} + +inline +Comparison_result +compare_y_at_xC2( + const Static_adaptatif_filter &l1a, + const Static_adaptatif_filter &l1b, + const Static_adaptatif_filter &l1c, + const Static_adaptatif_filter &l2a, + const Static_adaptatif_filter &l2b, + const Static_adaptatif_filter &l2c, + const Static_adaptatif_filter &h1a, + const Static_adaptatif_filter &h1b, + const Static_adaptatif_filter &h1c, + const Static_adaptatif_filter &h2a, + const Static_adaptatif_filter &h2b, + const Static_adaptatif_filter &h2c) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + static double SAF_epsilon_2; + static double SAF_epsilon_3; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(l1a.value()) > SAF_bound || + fabs(l1b.value()) > SAF_bound || + fabs(l1c.value()) > SAF_bound || + fabs(l2a.value()) > SAF_bound || + fabs(l2b.value()) > SAF_bound || + fabs(l2c.value()) > SAF_bound || + fabs(h1a.value()) > SAF_bound || + fabs(h1b.value()) > SAF_bound || + fabs(h1c.value()) > SAF_bound || + fabs(h2a.value()) > SAF_bound || + fabs(h2b.value()) > SAF_bound || + fabs(h2c.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(l1a.value())); + SAF_bound = std::max(SAF_bound, fabs(l1b.value())); + SAF_bound = std::max(SAF_bound, fabs(l1c.value())); + SAF_bound = std::max(SAF_bound, fabs(l2a.value())); + SAF_bound = std::max(SAF_bound, fabs(l2b.value())); + SAF_bound = std::max(SAF_bound, fabs(l2c.value())); + SAF_bound = std::max(SAF_bound, fabs(h1a.value())); + SAF_bound = std::max(SAF_bound, fabs(h1b.value())); + SAF_bound = std::max(SAF_bound, fabs(h1c.value())); + SAF_bound = std::max(SAF_bound, fabs(h2a.value())); + SAF_bound = std::max(SAF_bound, fabs(h2b.value())); + SAF_bound = std::max(SAF_bound, fabs(h2c.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) compare_y_at_xC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2, + SAF_epsilon_3); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return compare_y_at_xC2_SAF( + Restricted_double(l1a.value()), + Restricted_double(l1b.value()), + Restricted_double(l1c.value()), + Restricted_double(l2a.value()), + Restricted_double(l2b.value()), + Restricted_double(l2c.value()), + Restricted_double(h1a.value()), + Restricted_double(h1b.value()), + Restricted_double(h1c.value()), + Restricted_double(h2a.value()), + Restricted_double(h2b.value()), + Restricted_double(h2c.value()), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2, + SAF_epsilon_3); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return compare_y_at_xC2( + l1a.exact(), + l1b.exact(), + l1c.exact(), + l2a.exact(), + l2b.exact(), + l2c.exact(), + h1a.exact(), + h1b.exact(), + h1c.exact(), + h2a.exact(), + h2b.exact(), + h2c.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -270,6 +1201,7 @@ compare_y_at_xC2( const Filtered_exact &h2b, const Filtered_exact &h2c) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -308,6 +1240,98 @@ compare_y_at_xC2( } } +inline +Comparison_result +compare_deltax_deltayC2_SAF( + const Static_filter_error &px, + const Static_filter_error &qx, + const Static_filter_error &ry, + const Static_filter_error &sy, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return CGAL::compare_SAF(abs(px-qx), abs(ry-sy), + epsilon_0); +} + +inline +Comparison_result +compare_deltax_deltayC2_SAF( + const Restricted_double &px, + const Restricted_double &qx, + const Restricted_double &ry, + const Restricted_double &sy, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return CGAL::compare_SAF(abs(px-qx), abs(ry-sy), + epsilon_0); +} + +inline +Comparison_result +compare_deltax_deltayC2( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &ry, + const Static_adaptatif_filter &sy) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound || + fabs(sy.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + SAF_bound = std::max(SAF_bound, fabs(sy.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) compare_deltax_deltayC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return compare_deltax_deltayC2_SAF( + Restricted_double(px.value()), + Restricted_double(qx.value()), + Restricted_double(ry.value()), + Restricted_double(sy.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return compare_deltax_deltayC2( + px.exact(), + qx.exact(), + ry.exact(), + sy.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -319,6 +1343,7 @@ compare_deltax_deltayC2( const Filtered_exact &ry, const Filtered_exact &sy) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -341,6 +1366,114 @@ compare_deltax_deltayC2( } } +inline +Orientation +orientationC2_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &rx, + const Static_filter_error &ry, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return sign_of_determinant2x2_SAF(px-rx,py-ry,qx-rx,qy-ry, + epsilon_0); +} + +inline +Orientation +orientationC2_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &rx, + const Restricted_double &ry, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return sign_of_determinant2x2_SAF(px-rx,py-ry,qx-rx,qy-ry, + epsilon_0); +} + +inline +Orientation +orientationC2( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) orientationC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return orientationC2_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return orientationC2( + px.exact(), + py.exact(), + qx.exact(), + qy.exact(), + rx.exact(), + ry.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -354,6 +1487,7 @@ orientationC2( const Filtered_exact &rx, const Filtered_exact &ry) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -380,6 +1514,160 @@ orientationC2( } } +inline +Oriented_side +side_of_oriented_circleC2_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &rx, + const Static_filter_error &ry, + const Static_filter_error &tx, + const Static_filter_error &ty, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + + + + + + + FT ptx = px-tx; + FT pty = py-ty; + FT qtx = qx-tx; + FT qty = qy-ty; + FT rtx = rx-tx; + FT rty = ry-ty; + return Oriented_side( + sign_of_determinant3x3_SAF(ptx, pty, square(ptx) + square(pty), + qtx, qty, square(qtx) + square(qty), + rtx, rty, square(rtx) + square(rty), + epsilon_0)); +} + +inline +Oriented_side +side_of_oriented_circleC2_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &rx, + const Restricted_double &ry, + const Restricted_double &tx, + const Restricted_double &ty, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + + + + + + + FT ptx = px-tx; + FT pty = py-ty; + FT qtx = qx-tx; + FT qty = qy-ty; + FT rtx = rx-tx; + FT rty = ry-ty; + return Oriented_side( + sign_of_determinant3x3_SAF(ptx, pty, square(ptx) + square(pty), + qtx, qty, square(qtx) + square(qty), + rtx, rty, square(rtx) + square(rty), + epsilon_0)); +} + +inline +Oriented_side +side_of_oriented_circleC2( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry, + const Static_adaptatif_filter &tx, + const Static_adaptatif_filter &ty) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound || + fabs(tx.value()) > SAF_bound || + fabs(ty.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + SAF_bound = std::max(SAF_bound, fabs(tx.value())); + SAF_bound = std::max(SAF_bound, fabs(ty.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) side_of_oriented_circleC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return side_of_oriented_circleC2_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + Restricted_double(tx.value()), + Restricted_double(ty.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return side_of_oriented_circleC2( + px.exact(), + py.exact(), + qx.exact(), + qy.exact(), + rx.exact(), + ry.exact(), + tx.exact(), + ty.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -395,6 +1683,7 @@ side_of_oriented_circleC2( const Filtered_exact &tx, const Filtered_exact &ty) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -425,6 +1714,143 @@ side_of_oriented_circleC2( } } +inline +Bounded_side +side_of_bounded_circleC2_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &rx, + const Static_filter_error &ry, + const Static_filter_error &tx, + const Static_filter_error &ty, + double & epsilon_0, + double & epsilon_1) +{ + typedef Static_filter_error FT; + + Oriented_side s = side_of_oriented_circleC2_SAF(px,py,qx,qy,rx,ry,tx,ty, + epsilon_0); + Orientation o = orientationC2_SAF(px,py,qx,qy,rx,ry, + epsilon_1); + + return Bounded_side (s * o); +} + +inline +Bounded_side +side_of_bounded_circleC2_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &rx, + const Restricted_double &ry, + const Restricted_double &tx, + const Restricted_double &ty, + const double & epsilon_0, + const double & epsilon_1) +{ + typedef Restricted_double FT; + + Oriented_side s = side_of_oriented_circleC2_SAF(px,py,qx,qy,rx,ry,tx,ty, + epsilon_0); + Orientation o = orientationC2_SAF(px,py,qx,qy,rx,ry, + epsilon_1); + + return Bounded_side (s * o); +} + +inline +Bounded_side +side_of_bounded_circleC2( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry, + const Static_adaptatif_filter &tx, + const Static_adaptatif_filter &ty) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound || + fabs(tx.value()) > SAF_bound || + fabs(ty.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + SAF_bound = std::max(SAF_bound, fabs(tx.value())); + SAF_bound = std::max(SAF_bound, fabs(ty.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) side_of_bounded_circleC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return side_of_bounded_circleC2_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + Restricted_double(tx.value()), + Restricted_double(ty.value()), + SAF_epsilon_0, + SAF_epsilon_1); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return side_of_bounded_circleC2( + px.exact(), + py.exact(), + qx.exact(), + qy.exact(), + rx.exact(), + ry.exact(), + tx.exact(), + ty.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -440,6 +1866,7 @@ side_of_bounded_circleC2( const Filtered_exact &tx, const Filtered_exact &ty) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -470,6 +1897,116 @@ side_of_bounded_circleC2( } } +inline +Comparison_result +cmp_dist_to_pointC2_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &rx, + const Static_filter_error &ry, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return CGAL::compare_SAF(squared_distanceC2(px,py,qx,qy), + squared_distanceC2(px,py,rx,ry), + epsilon_0); +} + +inline +Comparison_result +cmp_dist_to_pointC2_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &rx, + const Restricted_double &ry, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return CGAL::compare_SAF(squared_distanceC2(px,py,qx,qy), + squared_distanceC2(px,py,rx,ry), + epsilon_0); +} + +inline +Comparison_result +cmp_dist_to_pointC2( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) cmp_dist_to_pointC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return cmp_dist_to_pointC2_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return cmp_dist_to_pointC2( + px.exact(), + py.exact(), + qx.exact(), + qy.exact(), + rx.exact(), + ry.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -483,6 +2020,7 @@ cmp_dist_to_pointC2( const Filtered_exact &rx, const Filtered_exact &ry) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -509,6 +2047,124 @@ cmp_dist_to_pointC2( } } +inline +Comparison_result +cmp_signed_dist_to_lineC2_SAF( + const Static_filter_error &la, + const Static_filter_error &lb, + const Static_filter_error &lc, + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &qx, + const Static_filter_error &qy, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return CGAL::compare_SAF(scaled_distance_to_directionC2(la,lb,px,py), + scaled_distance_to_directionC2(la,lb,qx,qy), + epsilon_0); +} + +inline +Comparison_result +cmp_signed_dist_to_lineC2_SAF( + const Restricted_double &la, + const Restricted_double &lb, + const Restricted_double &lc, + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &qx, + const Restricted_double &qy, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return CGAL::compare_SAF(scaled_distance_to_directionC2(la,lb,px,py), + scaled_distance_to_directionC2(la,lb,qx,qy), + epsilon_0); +} + +inline +Comparison_result +cmp_signed_dist_to_lineC2( + const Static_adaptatif_filter &la, + const Static_adaptatif_filter &lb, + const Static_adaptatif_filter &lc, + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(la.value()) > SAF_bound || + fabs(lb.value()) > SAF_bound || + fabs(lc.value()) > SAF_bound || + fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(la.value())); + SAF_bound = std::max(SAF_bound, fabs(lb.value())); + SAF_bound = std::max(SAF_bound, fabs(lc.value())); + SAF_bound = std::max(SAF_bound, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) cmp_signed_dist_to_lineC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return cmp_signed_dist_to_lineC2_SAF( + Restricted_double(la.value()), + Restricted_double(lb.value()), + Restricted_double(lc.value()), + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return cmp_signed_dist_to_lineC2( + la.exact(), + lb.exact(), + lc.exact(), + px.exact(), + py.exact(), + qx.exact(), + qy.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -523,6 +2179,7 @@ cmp_signed_dist_to_lineC2( const Filtered_exact &qx, const Filtered_exact &qy) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -551,6 +2208,132 @@ cmp_signed_dist_to_lineC2( } } +inline +Comparison_result +cmp_signed_dist_to_lineC2_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &rx, + const Static_filter_error &ry, + const Static_filter_error &sx, + const Static_filter_error &sy, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return CGAL::compare_SAF(scaled_distance_to_lineC2(px,py,qx,qy,rx,ry), + scaled_distance_to_lineC2(px,py,qx,qy,sx,sy), + epsilon_0); +} + +inline +Comparison_result +cmp_signed_dist_to_lineC2_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &rx, + const Restricted_double &ry, + const Restricted_double &sx, + const Restricted_double &sy, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return CGAL::compare_SAF(scaled_distance_to_lineC2(px,py,qx,qy,rx,ry), + scaled_distance_to_lineC2(px,py,qx,qy,sx,sy), + epsilon_0); +} + +inline +Comparison_result +cmp_signed_dist_to_lineC2( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry, + const Static_adaptatif_filter &sx, + const Static_adaptatif_filter &sy) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound || + fabs(sx.value()) > SAF_bound || + fabs(sy.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + SAF_bound = std::max(SAF_bound, fabs(sx.value())); + SAF_bound = std::max(SAF_bound, fabs(sy.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) cmp_signed_dist_to_lineC2_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return cmp_signed_dist_to_lineC2_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + Restricted_double(sx.value()), + Restricted_double(sy.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return cmp_signed_dist_to_lineC2( + px.exact(), + py.exact(), + qx.exact(), + qy.exact(), + rx.exact(), + ry.exact(), + sx.exact(), + sy.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -566,6 +2349,7 @@ cmp_signed_dist_to_lineC2( const Filtered_exact &sx, const Filtered_exact &sy) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { diff --git a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates_on_ftC3.h b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates_on_ftC3.h index 08be22bd74c..29cef336343 100644 --- a/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates_on_ftC3.h +++ b/Packages/Interval_arithmetic/include/CGAL/Arithmetic_filter/predicates_on_ftC3.h @@ -26,6 +26,168 @@ CGAL_BEGIN_NAMESPACE +inline +bool +collinearC3_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pz, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qz, + const Static_filter_error &rx, + const Static_filter_error &ry, + const Static_filter_error &rz, + double & epsilon_0, + double & epsilon_1, + double & epsilon_2) +{ + typedef Static_filter_error FT; + + FT dpx = px-rx; + FT dpy = py-ry; + FT dpz = pz-rz; + FT dqx = qx-rx; + FT dqy = qy-ry; + FT dqz = qz-rz; + return (sign_of_determinant2x2_SAF(dpx,dqx,dpy,dqy, + epsilon_0) == ZERO) + && (sign_of_determinant2x2_SAF(dpx,dqx,dpz,dqz, + epsilon_1) == ZERO) + && (sign_of_determinant2x2_SAF(dpy,dqy,dpz,dqz, + epsilon_2) == ZERO); +} + +inline +bool +collinearC3_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pz, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qz, + const Restricted_double &rx, + const Restricted_double &ry, + const Restricted_double &rz, + const double & epsilon_0, + const double & epsilon_1, + const double & epsilon_2) +{ + typedef Restricted_double FT; + + FT dpx = px-rx; + FT dpy = py-ry; + FT dpz = pz-rz; + FT dqx = qx-rx; + FT dqy = qy-ry; + FT dqz = qz-rz; + return (sign_of_determinant2x2_SAF(dpx,dqx,dpy,dqy, + epsilon_0) == ZERO) + && (sign_of_determinant2x2_SAF(dpx,dqx,dpz,dqz, + epsilon_1) == ZERO) + && (sign_of_determinant2x2_SAF(dpy,dqy,dpz,dqz, + epsilon_2) == ZERO); +} + +inline +bool +collinearC3( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pz, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qz, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry, + const Static_adaptatif_filter &rz) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + static double SAF_epsilon_2; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pz.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qz.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound || + fabs(rz.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pz.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qz.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + SAF_bound = std::max(SAF_bound, fabs(rz.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) collinearC3_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return collinearC3_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pz.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qz.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + Restricted_double(rz.value()), + SAF_epsilon_0, + SAF_epsilon_1, + SAF_epsilon_2); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return collinearC3( + px.exact(), + py.exact(), + pz.exact(), + qx.exact(), + qy.exact(), + qz.exact(), + rx.exact(), + ry.exact(), + rz.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -42,6 +204,7 @@ collinearC3( const Filtered_exact &ry, const Filtered_exact &rz) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -74,6 +237,166 @@ collinearC3( } } +inline +Orientation +orientationC3_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pz, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qz, + const Static_filter_error &rx, + const Static_filter_error &ry, + const Static_filter_error &rz, + const Static_filter_error &sx, + const Static_filter_error &sy, + const Static_filter_error &sz, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return Orientation(sign_of_determinant3x3_SAF(qx-px,rx-px,sx-px, + qy-py,ry-py,sy-py, + qz-pz,rz-pz,sz-pz, + epsilon_0)); +} + +inline +Orientation +orientationC3_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pz, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qz, + const Restricted_double &rx, + const Restricted_double &ry, + const Restricted_double &rz, + const Restricted_double &sx, + const Restricted_double &sy, + const Restricted_double &sz, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return Orientation(sign_of_determinant3x3_SAF(qx-px,rx-px,sx-px, + qy-py,ry-py,sy-py, + qz-pz,rz-pz,sz-pz, + epsilon_0)); +} + +inline +Orientation +orientationC3( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pz, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qz, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry, + const Static_adaptatif_filter &rz, + const Static_adaptatif_filter &sx, + const Static_adaptatif_filter &sy, + const Static_adaptatif_filter &sz) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pz.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qz.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound || + fabs(rz.value()) > SAF_bound || + fabs(sx.value()) > SAF_bound || + fabs(sy.value()) > SAF_bound || + fabs(sz.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pz.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qz.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + SAF_bound = std::max(SAF_bound, fabs(rz.value())); + SAF_bound = std::max(SAF_bound, fabs(sx.value())); + SAF_bound = std::max(SAF_bound, fabs(sy.value())); + SAF_bound = std::max(SAF_bound, fabs(sz.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) orientationC3_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return orientationC3_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pz.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qz.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + Restricted_double(rz.value()), + Restricted_double(sx.value()), + Restricted_double(sy.value()), + Restricted_double(sz.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return orientationC3( + px.exact(), + py.exact(), + pz.exact(), + qx.exact(), + qy.exact(), + qz.exact(), + rx.exact(), + ry.exact(), + rz.exact(), + sx.exact(), + sy.exact(), + sz.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -93,6 +416,7 @@ orientationC3( const Filtered_exact &sy, const Filtered_exact &sz) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -131,6 +455,224 @@ orientationC3( } } +inline +Oriented_side +side_of_oriented_sphereC3_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pz, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qz, + const Static_filter_error &rx, + const Static_filter_error &ry, + const Static_filter_error &rz, + const Static_filter_error &sx, + const Static_filter_error &sy, + const Static_filter_error &sz, + const Static_filter_error &tx, + const Static_filter_error &ty, + const Static_filter_error &tz, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + FT ptx = px - tx; + FT pty = py - ty; + FT ptz = pz - tz; + FT pt2 = square(ptx) + square(pty) + square(ptz); + FT qtx = qx - tx; + FT qty = qy - ty; + FT qtz = qz - tz; + FT qt2 = square(qtx) + square(qty) + square(qtz); + FT rtx = rx - tx; + FT rty = ry - ty; + FT rtz = rz - tz; + FT rt2 = square(rtx) + square(rty) + square(rtz); + FT stx = sx - tx; + FT sty = sy - ty; + FT stz = sz - tz; + FT st2 = square(stx) + square(sty) + square(stz); + return Oriented_side(sign_of_determinant4x4_SAF(ptx,pty,ptz,pt2, + rtx,rty,rtz,rt2, + qtx,qty,qtz,qt2, + stx,sty,stz,st2, + epsilon_0)); +} + +inline +Oriented_side +side_of_oriented_sphereC3_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pz, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qz, + const Restricted_double &rx, + const Restricted_double &ry, + const Restricted_double &rz, + const Restricted_double &sx, + const Restricted_double &sy, + const Restricted_double &sz, + const Restricted_double &tx, + const Restricted_double &ty, + const Restricted_double &tz, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + FT ptx = px - tx; + FT pty = py - ty; + FT ptz = pz - tz; + FT pt2 = square(ptx) + square(pty) + square(ptz); + FT qtx = qx - tx; + FT qty = qy - ty; + FT qtz = qz - tz; + FT qt2 = square(qtx) + square(qty) + square(qtz); + FT rtx = rx - tx; + FT rty = ry - ty; + FT rtz = rz - tz; + FT rt2 = square(rtx) + square(rty) + square(rtz); + FT stx = sx - tx; + FT sty = sy - ty; + FT stz = sz - tz; + FT st2 = square(stx) + square(sty) + square(stz); + return Oriented_side(sign_of_determinant4x4_SAF(ptx,pty,ptz,pt2, + rtx,rty,rtz,rt2, + qtx,qty,qtz,qt2, + stx,sty,stz,st2, + epsilon_0)); +} + +inline +Oriented_side +side_of_oriented_sphereC3( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pz, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qz, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry, + const Static_adaptatif_filter &rz, + const Static_adaptatif_filter &sx, + const Static_adaptatif_filter &sy, + const Static_adaptatif_filter &sz, + const Static_adaptatif_filter &tx, + const Static_adaptatif_filter &ty, + const Static_adaptatif_filter &tz) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pz.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qz.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound || + fabs(rz.value()) > SAF_bound || + fabs(sx.value()) > SAF_bound || + fabs(sy.value()) > SAF_bound || + fabs(sz.value()) > SAF_bound || + fabs(tx.value()) > SAF_bound || + fabs(ty.value()) > SAF_bound || + fabs(tz.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pz.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qz.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + SAF_bound = std::max(SAF_bound, fabs(rz.value())); + SAF_bound = std::max(SAF_bound, fabs(sx.value())); + SAF_bound = std::max(SAF_bound, fabs(sy.value())); + SAF_bound = std::max(SAF_bound, fabs(sz.value())); + SAF_bound = std::max(SAF_bound, fabs(tx.value())); + SAF_bound = std::max(SAF_bound, fabs(ty.value())); + SAF_bound = std::max(SAF_bound, fabs(tz.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) side_of_oriented_sphereC3_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return side_of_oriented_sphereC3_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pz.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qz.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + Restricted_double(rz.value()), + Restricted_double(sx.value()), + Restricted_double(sy.value()), + Restricted_double(sz.value()), + Restricted_double(tx.value()), + Restricted_double(ty.value()), + Restricted_double(tz.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return side_of_oriented_sphereC3( + px.exact(), + py.exact(), + pz.exact(), + qx.exact(), + qy.exact(), + qz.exact(), + rx.exact(), + ry.exact(), + rz.exact(), + sx.exact(), + sy.exact(), + sz.exact(), + tx.exact(), + ty.exact(), + tz.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -153,6 +695,7 @@ side_of_oriented_sphereC3( const Filtered_exact &ty, const Filtered_exact &tz) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -197,6 +740,211 @@ side_of_oriented_sphereC3( } } +inline +Bounded_side +side_of_bounded_sphereC3_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pz, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qz, + const Static_filter_error &rx, + const Static_filter_error &ry, + const Static_filter_error &rz, + const Static_filter_error &sx, + const Static_filter_error &sy, + const Static_filter_error &sz, + const Static_filter_error &tx, + const Static_filter_error &ty, + const Static_filter_error &tz, + double & epsilon_0, + double & epsilon_1) +{ + typedef Static_filter_error FT; + + Oriented_side s = side_of_oriented_sphereC3_SAF(px, py, pz, + qx, qy, qz, + rx, ry, rz, + sx, sy, sz, + tx, ty, tz, + epsilon_0); + Orientation o = orientationC3_SAF(px, py, pz, + qx, qy, qz, + rx, ry, rz, + sx, sy, sz, + epsilon_1); + return Bounded_side(s * o); +} + +inline +Bounded_side +side_of_bounded_sphereC3_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pz, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qz, + const Restricted_double &rx, + const Restricted_double &ry, + const Restricted_double &rz, + const Restricted_double &sx, + const Restricted_double &sy, + const Restricted_double &sz, + const Restricted_double &tx, + const Restricted_double &ty, + const Restricted_double &tz, + const double & epsilon_0, + const double & epsilon_1) +{ + typedef Restricted_double FT; + + Oriented_side s = side_of_oriented_sphereC3_SAF(px, py, pz, + qx, qy, qz, + rx, ry, rz, + sx, sy, sz, + tx, ty, tz, + epsilon_0); + Orientation o = orientationC3_SAF(px, py, pz, + qx, qy, qz, + rx, ry, rz, + sx, sy, sz, + epsilon_1); + return Bounded_side(s * o); +} + +inline +Bounded_side +side_of_bounded_sphereC3( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pz, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qz, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry, + const Static_adaptatif_filter &rz, + const Static_adaptatif_filter &sx, + const Static_adaptatif_filter &sy, + const Static_adaptatif_filter &sz, + const Static_adaptatif_filter &tx, + const Static_adaptatif_filter &ty, + const Static_adaptatif_filter &tz) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + static double SAF_epsilon_1; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pz.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qz.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound || + fabs(rz.value()) > SAF_bound || + fabs(sx.value()) > SAF_bound || + fabs(sy.value()) > SAF_bound || + fabs(sz.value()) > SAF_bound || + fabs(tx.value()) > SAF_bound || + fabs(ty.value()) > SAF_bound || + fabs(tz.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pz.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qz.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + SAF_bound = std::max(SAF_bound, fabs(rz.value())); + SAF_bound = std::max(SAF_bound, fabs(sx.value())); + SAF_bound = std::max(SAF_bound, fabs(sy.value())); + SAF_bound = std::max(SAF_bound, fabs(sz.value())); + SAF_bound = std::max(SAF_bound, fabs(tx.value())); + SAF_bound = std::max(SAF_bound, fabs(ty.value())); + SAF_bound = std::max(SAF_bound, fabs(tz.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) side_of_bounded_sphereC3_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0, + SAF_epsilon_1); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return side_of_bounded_sphereC3_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pz.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qz.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + Restricted_double(rz.value()), + Restricted_double(sx.value()), + Restricted_double(sy.value()), + Restricted_double(sz.value()), + Restricted_double(tx.value()), + Restricted_double(ty.value()), + Restricted_double(tz.value()), + SAF_epsilon_0, + SAF_epsilon_1); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return side_of_bounded_sphereC3( + px.exact(), + py.exact(), + pz.exact(), + qx.exact(), + qy.exact(), + qz.exact(), + rx.exact(), + ry.exact(), + rz.exact(), + sx.exact(), + sy.exact(), + sz.exact(), + tx.exact(), + ty.exact(), + tz.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -219,6 +967,7 @@ side_of_bounded_sphereC3( const Filtered_exact &ty, const Filtered_exact &tz) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -263,6 +1012,140 @@ side_of_bounded_sphereC3( } } +inline +Comparison_result +cmp_dist_to_pointC3_SAF( + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pz, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qz, + const Static_filter_error &rx, + const Static_filter_error &ry, + const Static_filter_error &rz, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return CGAL::compare_SAF(squared_distanceC3(px,py,pz,qx,qy,qz), + squared_distanceC3(px,py,pz,rx,ry,rz), + epsilon_0); +} + +inline +Comparison_result +cmp_dist_to_pointC3_SAF( + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pz, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qz, + const Restricted_double &rx, + const Restricted_double &ry, + const Restricted_double &rz, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return CGAL::compare_SAF(squared_distanceC3(px,py,pz,qx,qy,qz), + squared_distanceC3(px,py,pz,rx,ry,rz), + epsilon_0); +} + +inline +Comparison_result +cmp_dist_to_pointC3( + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pz, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qz, + const Static_adaptatif_filter &rx, + const Static_adaptatif_filter &ry, + const Static_adaptatif_filter &rz) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pz.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qz.value()) > SAF_bound || + fabs(rx.value()) > SAF_bound || + fabs(ry.value()) > SAF_bound || + fabs(rz.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pz.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qz.value())); + SAF_bound = std::max(SAF_bound, fabs(rx.value())); + SAF_bound = std::max(SAF_bound, fabs(ry.value())); + SAF_bound = std::max(SAF_bound, fabs(rz.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) cmp_dist_to_pointC3_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return cmp_dist_to_pointC3_SAF( + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pz.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qz.value()), + Restricted_double(rx.value()), + Restricted_double(ry.value()), + Restricted_double(rz.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return cmp_dist_to_pointC3( + px.exact(), + py.exact(), + pz.exact(), + qx.exact(), + qy.exact(), + qz.exact(), + rx.exact(), + ry.exact(), + rz.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -279,6 +1162,7 @@ cmp_dist_to_pointC3( const Filtered_exact &ry, const Filtered_exact &rz) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -311,6 +1195,148 @@ cmp_dist_to_pointC3( } } +inline +Comparison_result +cmp_signed_dist_to_planeC3_SAF( + const Static_filter_error &pa, + const Static_filter_error &pb, + const Static_filter_error &pc, + const Static_filter_error &pd, + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pz, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qz, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return CGAL::compare_SAF(scaled_distance_to_planeC3(pa,pb,pc,pd,px,py,pz), + scaled_distance_to_planeC3(pa,pb,pc,pd,qx,qy,qz), + epsilon_0); +} + +inline +Comparison_result +cmp_signed_dist_to_planeC3_SAF( + const Restricted_double &pa, + const Restricted_double &pb, + const Restricted_double &pc, + const Restricted_double &pd, + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pz, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qz, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return CGAL::compare_SAF(scaled_distance_to_planeC3(pa,pb,pc,pd,px,py,pz), + scaled_distance_to_planeC3(pa,pb,pc,pd,qx,qy,qz), + epsilon_0); +} + +inline +Comparison_result +cmp_signed_dist_to_planeC3( + const Static_adaptatif_filter &pa, + const Static_adaptatif_filter &pb, + const Static_adaptatif_filter &pc, + const Static_adaptatif_filter &pd, + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pz, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qz) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(pa.value()) > SAF_bound || + fabs(pb.value()) > SAF_bound || + fabs(pc.value()) > SAF_bound || + fabs(pd.value()) > SAF_bound || + fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pz.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qz.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(pa.value())); + SAF_bound = std::max(SAF_bound, fabs(pb.value())); + SAF_bound = std::max(SAF_bound, fabs(pc.value())); + SAF_bound = std::max(SAF_bound, fabs(pd.value())); + SAF_bound = std::max(SAF_bound, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pz.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qz.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) cmp_signed_dist_to_planeC3_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return cmp_signed_dist_to_planeC3_SAF( + Restricted_double(pa.value()), + Restricted_double(pb.value()), + Restricted_double(pc.value()), + Restricted_double(pd.value()), + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pz.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qz.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return cmp_signed_dist_to_planeC3( + pa.exact(), + pb.exact(), + pc.exact(), + pd.exact(), + px.exact(), + py.exact(), + pz.exact(), + qx.exact(), + qy.exact(), + qz.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -328,6 +1354,7 @@ cmp_signed_dist_to_planeC3( const Filtered_exact &qy, const Filtered_exact &qz) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try { @@ -362,6 +1389,198 @@ cmp_signed_dist_to_planeC3( } } +inline +Comparison_result +cmp_signed_dist_to_planeC3_SAF( + const Static_filter_error &ppx, + const Static_filter_error &ppy, + const Static_filter_error &ppz, + const Static_filter_error &pqx, + const Static_filter_error &pqy, + const Static_filter_error &pqz, + const Static_filter_error &prx, + const Static_filter_error &pry, + const Static_filter_error &prz, + const Static_filter_error &px, + const Static_filter_error &py, + const Static_filter_error &pz, + const Static_filter_error &qx, + const Static_filter_error &qy, + const Static_filter_error &qz, + double & epsilon_0) +{ + typedef Static_filter_error FT; + + return CGAL::compare_SAF( + scaled_distance_to_planeC3(ppx,ppy,ppz,pqx,pqy,pqz, + prx,pry,prz,psx,psy,psz, + px,py,pz), + scaled_distance_to_planeC3(ppx,ppy,ppz,pqx,pqy,pqz, + prx,pry,prz,psx,psy,psz, + qx,qy,qz) , + epsilon_0); +} + +inline +Comparison_result +cmp_signed_dist_to_planeC3_SAF( + const Restricted_double &ppx, + const Restricted_double &ppy, + const Restricted_double &ppz, + const Restricted_double &pqx, + const Restricted_double &pqy, + const Restricted_double &pqz, + const Restricted_double &prx, + const Restricted_double &pry, + const Restricted_double &prz, + const Restricted_double &px, + const Restricted_double &py, + const Restricted_double &pz, + const Restricted_double &qx, + const Restricted_double &qy, + const Restricted_double &qz, + const double & epsilon_0) +{ + typedef Restricted_double FT; + + return CGAL::compare_SAF( + scaled_distance_to_planeC3(ppx,ppy,ppz,pqx,pqy,pqz, + prx,pry,prz,psx,psy,psz, + px,py,pz), + scaled_distance_to_planeC3(ppx,ppy,ppz,pqx,pqy,pqz, + prx,pry,prz,psx,psy,psz, + qx,qy,qz) , + epsilon_0); +} + +inline +Comparison_result +cmp_signed_dist_to_planeC3( + const Static_adaptatif_filter &ppx, + const Static_adaptatif_filter &ppy, + const Static_adaptatif_filter &ppz, + const Static_adaptatif_filter &pqx, + const Static_adaptatif_filter &pqy, + const Static_adaptatif_filter &pqz, + const Static_adaptatif_filter &prx, + const Static_adaptatif_filter &pry, + const Static_adaptatif_filter &prz, + const Static_adaptatif_filter &px, + const Static_adaptatif_filter &py, + const Static_adaptatif_filter &pz, + const Static_adaptatif_filter &qx, + const Static_adaptatif_filter &qy, + const Static_adaptatif_filter &qz) +{ + bool re_adjusted = false; + static double SAF_bound = -1.0; + static double SAF_epsilon_0; + + // Check the bounds. All arguments must be <= SAF_bound. + if ( fabs(ppx.value()) > SAF_bound || + fabs(ppy.value()) > SAF_bound || + fabs(ppz.value()) > SAF_bound || + fabs(pqx.value()) > SAF_bound || + fabs(pqy.value()) > SAF_bound || + fabs(pqz.value()) > SAF_bound || + fabs(prx.value()) > SAF_bound || + fabs(pry.value()) > SAF_bound || + fabs(prz.value()) > SAF_bound || + fabs(px.value()) > SAF_bound || + fabs(py.value()) > SAF_bound || + fabs(pz.value()) > SAF_bound || + fabs(qx.value()) > SAF_bound || + fabs(qy.value()) > SAF_bound || + fabs(qz.value()) > SAF_bound) + { +re_adjust: + // Re-adjust SAF_bound. + SAF_bound = std::max(0.0, fabs(ppx.value())); + SAF_bound = std::max(SAF_bound, fabs(ppy.value())); + SAF_bound = std::max(SAF_bound, fabs(ppz.value())); + SAF_bound = std::max(SAF_bound, fabs(pqx.value())); + SAF_bound = std::max(SAF_bound, fabs(pqy.value())); + SAF_bound = std::max(SAF_bound, fabs(pqz.value())); + SAF_bound = std::max(SAF_bound, fabs(prx.value())); + SAF_bound = std::max(SAF_bound, fabs(pry.value())); + SAF_bound = std::max(SAF_bound, fabs(prz.value())); + SAF_bound = std::max(SAF_bound, fabs(px.value())); + SAF_bound = std::max(SAF_bound, fabs(py.value())); + SAF_bound = std::max(SAF_bound, fabs(pz.value())); + SAF_bound = std::max(SAF_bound, fabs(qx.value())); + SAF_bound = std::max(SAF_bound, fabs(qy.value())); + SAF_bound = std::max(SAF_bound, fabs(qz.value())); + + // recompute the epsilons: "just" call it over Static_filter_error. + // That's the tricky part that might not work for everything. + (void) cmp_signed_dist_to_planeC3_SAF( + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + Static_filter_error(SAF_bound), + SAF_epsilon_0); + + // TODO: We should verify that all epsilons have really been updated. + } + + try // Try the epsilon variant of the predicate. + { + return cmp_signed_dist_to_planeC3_SAF( + Restricted_double(ppx.value()), + Restricted_double(ppy.value()), + Restricted_double(ppz.value()), + Restricted_double(pqx.value()), + Restricted_double(pqy.value()), + Restricted_double(pqz.value()), + Restricted_double(prx.value()), + Restricted_double(pry.value()), + Restricted_double(prz.value()), + Restricted_double(px.value()), + Restricted_double(py.value()), + Restricted_double(pz.value()), + Restricted_double(qx.value()), + Restricted_double(qy.value()), + Restricted_double(qz.value()), + SAF_epsilon_0); + } + catch (Restricted_double::unsafe_comparison) + { + // It failed, we re-adjust once. + if (!re_adjusted) { + re_adjusted = true; + goto re_adjust; + } + // This scheme definitely fails => exact computation (filtered_exact<> ?). + return cmp_signed_dist_to_planeC3( + ppx.exact(), + ppy.exact(), + ppz.exact(), + pqx.exact(), + pqy.exact(), + pqz.exact(), + prx.exact(), + pry.exact(), + prz.exact(), + px.exact(), + py.exact(), + pz.exact(), + qx.exact(), + qy.exact(), + qz.exact()); + } +} + #ifndef CGAL_CFG_NO_EXPLICIT_TEMPLATE_FUNCTION_ARGUMENT_SPECIFICATION template < class CGAL_IA_CT, class CGAL_IA_ET, class CGAL_IA_CACHE > #endif @@ -384,6 +1603,7 @@ cmp_signed_dist_to_planeC3( const Filtered_exact &qy, const Filtered_exact &qz) { + CGAL_assertion(Interval_nt_advanced::want_exceptions); FPU_CW_t backup = FPU_get_and_set_cw(FPU_cw_up); try {