diff --git a/Straight_skeleton_2/include/CGAL/compute_outer_frame_margin.h b/Straight_skeleton_2/include/CGAL/compute_outer_frame_margin.h index 53e09f5d572..2a54323b534 100644 --- a/Straight_skeleton_2/include/CGAL/compute_outer_frame_margin.h +++ b/Straight_skeleton_2/include/CGAL/compute_outer_frame_margin.h @@ -95,11 +95,12 @@ boost::optional< typename Traits::FT > compute_outer_frame_margin ( ForwardPoint if ( ! lOverflow ) { - FT lDist = CGAL_SS_i::inexact_sqrt(lMaxSDist) ; + FT lDist = inexact_sqrt(lMaxSDist) ; + double approx = ceil( to_interval(lDist + ( aOffset * FT(1.05) ) ).second ); // Add a %5 gap, and ceil to get simpler values - CGAL_STSKEL_BUILDER_TRACE(4, "outer frame margin: " << ceil(lDist + ( aOffset * FT(1.05) ) ) ); - return boost::optional( ceil(lDist + ( aOffset * FT(1.05) ) ) ) ; + CGAL_STSKEL_BUILDER_TRACE(4, "outer frame margin: " << approx ); + return boost::optional ( approx ) ; } else return boost::optional(); diff --git a/Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h b/Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h index 4ec1a8e56f6..61985045684 100644 --- a/Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h +++ b/Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h @@ -18,6 +18,7 @@ #include #include #include +#include #include @@ -70,42 +71,146 @@ Trisegment_collinearity trisegment_collinearity_no_exact_constructions ( Segment return TRISEGMENT_COLLINEARITY_ALL; } -template -inline NT inexact_sqrt_implementation( NT const& n, CGAL::Null_functor /*no_sqrt*/ ) +/// + +template +typename Coercion_traits::Type +inexact_sqrt(const NT& n, CGAL::Null_functor) { typedef CGAL::Interval_nt IFT; typename IFT::Protector protector; CGAL::NT_converter to_ift; - return NT( to_double(sqrt(to_ift(n))) ); + IFT sqrt_ift = sqrt(to_ift(n)); + CGAL_STSKEL_TRAITS_TRACE("interval " << sqrt_ift.inf() << " " << sqrt_ift.sup() ) ; + CGAL_STSKEL_TRAITS_TRACE("delta " << sqrt_ift.sup() - sqrt_ift.inf() ) ; + + return NT(to_double(sqrt_ift)); } -template -inline NT inexact_sqrt_implementation( NT const& n, Sqrt sqrt_f ) +template +typename Sqrt::result_type +inexact_sqrt(const NT& nt, Sqrt sqrt) { + CGAL_STSKEL_TRAITS_TRACE("sqrt(" << typeid(NT).name() << ")"); + return sqrt(nt); +} + +template +decltype(auto) inexact_sqrt(const NT& nt) +{ + // the initial version of this function was using Algebraic_category + // for the dispatch but some ring type (like Gmpz) provides a Sqrt + // functor even if not being Field_with_sqrt. + typedef CGAL::Algebraic_structure_traits AST; + typedef typename AST::Sqrt Sqrt; + return inexact_sqrt(nt, Sqrt()); +} + +template +Quotient +inexact_sqrt(const Quotient& q) +{ + return { inexact_sqrt(q.numerator()*q.denominator()), abs(q.denominator()) }; +} + +template +Lazy_exact_nt +inexact_sqrt(const Lazy_exact_nt& lz) +{ + return inexact_sqrt(exact(lz)); +} + +// Currently the norm can be inexact even when we are in the exact pipeline of the traits +// if we use something like K = EPICK because there is no exact sqrt. +// +// @todo Ideally, we could compute how much precision is required for the sqrt +// given all possible operations that are performed in the SLS and the input values +// and compute (in the exact pipeline) an approximate sqrt with such sufficient precision. + +template +FT inexact_norm (const FT& x, const FT& y, + CGAL::Null_functor /*no_sqrt*/) +{ + CGAL_STSKEL_TRAITS_TRACE("inexact_norm(" << typeid(FT).name() << "," << typeid(FT).name() << ",Null_functor)"); + return std::hypot(CGAL::to_double(x), CGAL::to_double(y)); +} + +template +FT inexact_norm (const FT& x, const FT& y, + Sqrt sqrt_f) +{ + CGAL_STSKEL_TRAITS_TRACE("inexact_norm(" << typeid(FT).name() << "," << typeid(FT).name() << "," << typeid(Sqrt).name() << ")"); + const FT n = square(x) + square(y); return sqrt_f(n); } -template -inline NT inexact_sqrt( NT const& n ) +template +FT inexact_norm (const FT& x, const FT& y) { - typedef CGAL::Algebraic_structure_traits AST; + typedef CGAL::Algebraic_structure_traits AST; typedef typename AST::Sqrt Sqrt; - return inexact_sqrt_implementation(n,Sqrt()); + return inexact_norm(x,y, Sqrt()); } -inline Quotient inexact_sqrt( Quotient const& q ) +template +inline Lazy_exact_nt inexact_norm( Lazy_exact_nt const& lx, + Lazy_exact_nt const& ly ) { - CGAL_precondition(q > 0); - return Quotient(CGAL_SS_i::inexact_sqrt(q.numerator()*q.denominator()), q.denominator() ); + return inexact_norm( exact(lx), exact(ly) ) ; } -template -inline Lazy_exact_nt inexact_sqrt( Lazy_exact_nt const& lz) +template +inline Quotient inexact_norm( Quotient const& x, + Quotient const& y ) { - return inexact_sqrt( exact(lz) ); + CGAL_STSKEL_TRAITS_TRACE("inexact_norm(Quotient,Quotient)"); + return { inexact_norm(x.numerator()*y.denominator(), y.numerator()*x.denominator()), + abs(x.denominator()*y.denominator()) } ; } +template +inline +Interval_nt +inexact_norm (const Interval_nt& x, const Interval_nt& y) +{ + typename Interval_nt::Internal_protector P; + + CGAL_STSKEL_TRAITS_TRACE("inexact_norm(Interval_nt,Interval_nt)"); + +#if 0 // CGAL_USE_SSE2 _mm_hypot_pd is documented but does not exist...? + __m128d xx = IA_opacify128(x.simd()); + __m128d yy = IA_opacify128(y.simd()); + __m128d r = _mm_hypot_pd(xx, yy); + return Interval_nt(IA_opacify128(r)); +#else + + // @fixme is below correct? + +#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST + double i = std::nextafter(std::hypot(x.inf(), y.inf()), 0.) ; +#else + FPU_set_cw(CGAL_FE_DOWNWARD); + double i = CGAL_IA_FORCE_TO_DOUBLE(std::hypot(CGAL_IA_STOP_CPROP(x.inf()), + CGAL_IA_STOP_CPROP(y.inf()))); + FPU_set_cw(CGAL_FE_UPWARD); +#endif + + return Interval_nt(i, IA_up(std::hypot(x.sup(), y.sup()))); +#endif +} + +/// + +template +CGAL::Lazy_exact_nt ceil(const CGAL::Lazy_exact_nt& n) +{ + return { ceil(exact(n)) }; +} + +/// + + // Given an oriented 2D straight line segment 'e', computes the normalized coefficients (a,b,c) // of the supporting line, and weights them with 'aWeight'. // @@ -167,21 +272,18 @@ boost::optional< typename K::Line_2> compute_normalized_line_coeffC2( Segment_2< { FT sa = e.source().y() - e.target().y(); FT sb = e.target().x() - e.source().x(); - FT l2 = (sa*sa) + (sb*sb) ; + FT l = inexact_norm(sa, sb); - if ( CGAL_NTS is_finite(l2) ) + if ( CGAL_NTS is_finite(l) ) { - FT l = CGAL_SS_i :: inexact_sqrt(l2); - a = sa / l ; b = sb / l ; c = -e.source().x()*a - e.source().y()*b; - CGAL_STSKEL_TRAITS_TRACE("Unweighted line coefficients for line: " << s2str(e) - << "\nsa="<< n2str(sa) << "\nsb=" << n2str(sb) << "\nl2=" << n2str(l2) << "\nl=" << n2str(l) - << "\na="<< n2str(a) << "\nb=" << n2str(b) << "\nc=" << n2str(c) - ) ; + CGAL_STSKEL_TRAITS_TRACE("GENERIC line; sa="<< n2str(sa) << " sb=" << n2str(sb) + << "\nl=" << n2str(l) + << "\na="<< n2str(a) << "\nb=" << n2str(b) << "\nc=" << n2str(c) ) ; } else {