From d6e39a16dc62d49c72426cff12b04be1e4efd008 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Thu, 2 Jul 2020 23:49:42 +0200 Subject: [PATCH 1/8] Special case the circumcenter computation for few points --- .../NewKernel_d/function_objects_cartesian.h | 129 +++++++++++------- 1 file changed, 78 insertions(+), 51 deletions(-) diff --git a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h index 5f534070db2..38935fd42b0 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h @@ -616,7 +616,36 @@ template struct Construct_circumcenter : Store_kernel { Point const& p0=*f; int d = pd(p0); - if (d+1 == std::distance(f,e)) + int k = static_cast(std::distance(f,e)); + if(k==1) return p0; + if(k==2){ + typename Get_functor::type mid(this->kernel()); + return mid(p0, *++f); + } + if(k==3){ + // Same equations as in the general case, but solved by hand (Cramer) + // (c-r).(p-r)=(p-r)²/2 + // (c-r).(q-r)=(q-r)²/2 + typedef typename Get_type::type Vector; + typename Get_functor::type sl(this->kernel()); + typename Get_functor::type sp(this->kernel()); + typename Get_functor::type sv(this->kernel()); + typename Get_functor::type dp(this->kernel()); + typename Get_functor::type tp(this->kernel()); + Iter f2=f; + Point const& q=*++f2; + Point const& r=*++f2; + Vector u = dp(p0, r); + Vector v = dp(q, r); + FT uv = sp(u, v); + FT u2 = sl(u); + FT v2 = sl(v); + FT den = 2 * (u2 * v2 - CGAL::square(uv)); + FT a = (u2 - uv) * v2 / den; + FT b = (v2 - uv) * u2 / den; + return tp(tp(r, sv(u, a)), sv(v, b)); + } + if (k == d+1) { // 2*(x-y).c == x^2-y^2 typedef typename LA::Square_matrix Matrix; @@ -641,61 +670,59 @@ template struct Construct_circumcenter : Store_kernel { //std::cout << "Sol: " << res << std::endl; return cp(d,LA::vector_begin(res),LA::vector_end(res)); } - else - { - /* - * Matrix P=(p1, p2, ...) (each point as a column) - * Matrix Q=2*t(p2-p1,p3-p1, ...) (each vector as a line) - * Matrix M: QP, adding a line of 1 at the top - * Vector B: (1, p2^2-p1^2, p3^2-p1^2, ...) - * Solve ML=B, the center of the sphere is PL - * - * It would likely be faster to write P then transpose, multiply, - * etc instead of doing it by hand. - */ - // TODO: check for degenerate cases? - typedef typename R_::Max_ambient_dimension D2; - typedef typename R_::LA::template Rebind_dimension::Other LAd; - typedef typename LAd::Square_matrix Matrix; - typedef typename LAd::Vector Vec; - typename Get_functor::type sp(this->kernel()); - int k=static_cast(std::distance(f,e)); - Matrix m(k,k); - Vec b(k); - Vec l(k); - int j,i=0; - // We are doing a quadratic number of *f, which can be costly with transforming_iterator. - for(Iter f2=f;f2!=e;++f2,++i){ - Point const& p2=*f2; - b(i)=m(i,i)=sdo(p2); - j=0; - for(Iter f3=f;f3!=e;++f3,++j){ - m(j,i)=m(i,j)=sp(p2,*f3); - } - } - for(i=1;i::Other LAd; + typedef typename LAd::Square_matrix Matrix; + typedef typename LAd::Vector Vec; + typename Get_functor::type sp(this->kernel()); + Matrix m(k,k); + Vec b(k); + Vec l(k); + int j,i=0; + // We are doing a quadratic number of *f, which can be costly with transforming_iterator. + for(Iter f2=f;f2!=e;++f2,++i){ + Point const& p2=*f2; + b(i)=m(i,i)=sdo(p2); j=0; - for(Iter f2=f;f2!=e;++f2,++j){ - for(i=0;i Date: Mon, 6 Jul 2020 22:54:38 +0200 Subject: [PATCH 2/8] More circumcenter tests --- .../NewKernel_d/function_objects_cartesian.h | 1 + NewKernel_d/test/NewKernel_d/Epick_d.cpp | 17 +++++++++++++++++ 2 files changed, 18 insertions(+) diff --git a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h index 38935fd42b0..b6568308226 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h @@ -617,6 +617,7 @@ template struct Construct_circumcenter : Store_kernel { Point const& p0=*f; int d = pd(p0); int k = static_cast(std::distance(f,e)); + CGAL_assume(k>=1); if(k==1) return p0; if(k==2){ typename Get_functor::type mid(this->kernel()); diff --git a/NewKernel_d/test/NewKernel_d/Epick_d.cpp b/NewKernel_d/test/NewKernel_d/Epick_d.cpp index 21c64288380..c72d5d207e6 100644 --- a/NewKernel_d/test/NewKernel_d/Epick_d.cpp +++ b/NewKernel_d/test/NewKernel_d/Epick_d.cpp @@ -712,6 +712,22 @@ void test3(){ std::ostringstream sv1; sv1 << v1; assert(sv1.str()=="3 3 2 1"); std::istringstream sv2("3 4 5 6"); sv2 >> v1; assert(v1[0]==4&&v1[1]==5); } +template +void test4(){ + typedef typename Ker::Point_d P; + typedef typename Ker::Construct_circumcenter_d CCc; + typedef typename Ker::Equal_d E; + Ker k(4); + CCc ccc Kinit(construct_circumcenter_d_object); + E ed Kinit(equal_d_object); + auto mkpt=[](auto...x){double l[]{(double)x...};return P(std::begin(l), std::end(l));}; + P tab1[]={mkpt(15,20,40,80),mkpt(10,23,36,80),mkpt(10,20,40,85),mkpt(10,15,40,80),mkpt(13,20,40,76)}; + assert(ed(ccc(tab1+0, tab1+5),mkpt(10,20,40,80))); + P tab2[]={mkpt(15,20,40,80),mkpt(13,24,40,80),mkpt(10,25,40,80),mkpt(10,20,43,84)}; + assert(ed(ccc(tab2+0, tab2+4),mkpt(10,20,40,80))); + P tab3[]={mkpt(15,20,35,80),mkpt(10,25,40,75),mkpt(13,24,37,76)}; + assert(ed(ccc(tab3+0, tab3+3),mkpt(10,20,40,80))); +} template struct CGAL::Epick_d >; template struct CGAL::Epick_d >; template struct CGAL::Epick_d; @@ -733,6 +749,7 @@ int main(){ test2>>(); test3>>(); test3>(); + test4>(); #endif } From 22e8b9cc224c7591006020f4142bec67424fe79a Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Tue, 7 Jul 2020 13:46:48 +0200 Subject: [PATCH 3/8] Rewrite general case of circumcenter --- .../NewKernel_d/function_objects_cartesian.h | 102 +++++++----------- 1 file changed, 41 insertions(+), 61 deletions(-) diff --git a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h index b6568308226..d56fe9f632d 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h @@ -603,31 +603,28 @@ namespace CartesianDKernelFunctors { template struct Construct_circumcenter : Store_kernel { CGAL_FUNCTOR_INIT_STORE(Construct_circumcenter) typedef typename Get_type::type Point; + typedef typename Get_type::type Vector; typedef Point result_type; typedef typename Get_type::type FT; template result_type operator()(Iter f, Iter e)const{ - typedef typename Get_type::type Point; - typedef typename R_::LA LA; - typename Get_functor::type c(this->kernel()); - typename Get_functor >::type cp(this->kernel()); typename Get_functor::type pd(this->kernel()); - typename Get_functor::type sdo(this->kernel()); Point const& p0=*f; int d = pd(p0); int k = static_cast(std::distance(f,e)); + // Sorted from fastest to slowest, whether the dimension is 2, 3 or 4. It may be worth checking at some point. CGAL_assume(k>=1); if(k==1) return p0; if(k==2){ typename Get_functor::type mid(this->kernel()); return mid(p0, *++f); } + // TODO: check for degenerate cases in all the following cases? if(k==3){ // Same equations as in the general case, but solved by hand (Cramer) // (c-r).(p-r)=(p-r)²/2 // (c-r).(q-r)=(q-r)²/2 - typedef typename Get_type::type Vector; typename Get_functor::type sl(this->kernel()); typename Get_functor::type sp(this->kernel()); typename Get_functor::type sv(this->kernel()); @@ -644,14 +641,19 @@ template struct Construct_circumcenter : Store_kernel { FT den = 2 * (u2 * v2 - CGAL::square(uv)); FT a = (u2 - uv) * v2 / den; FT b = (v2 - uv) * u2 / den; + // Wasteful if we only want the radius return tp(tp(r, sv(u, a)), sv(v, b)); } if (k == d+1) { // 2*(x-y).c == x^2-y^2 + typedef typename R_::LA LA; typedef typename LA::Square_matrix Matrix; typedef typename LA::Vector Vec; typedef typename LA::Construct_vector CVec; + typename Get_functor::type c(this->kernel()); + typename Get_functor >::type cp(this->kernel()); + typename Get_functor::type sdo(this->kernel()); FT const& n0 = sdo(p0); Matrix m(d,d); Vec b = typename CVec::Dimension()(d); @@ -665,70 +667,48 @@ template struct Construct_circumcenter : Store_kernel { b[i] = sdo(p) - n0; } CGAL_assertion (i == d); - Vec res = typename CVec::Dimension()(d);; - //std::cout << "Mat: " << m << "\n Vec: " << one << std::endl; - LA::solve(res, std::move(m), std::move(b)); - //std::cout << "Sol: " << res << std::endl; + //Vec res = typename CVec::Dimension()(d);; + //LA::solve(res, std::move(m), std::move(b)); + // We already assume Eigen below... + Vec res=m.partialPivLu().solve(b); return cp(d,LA::vector_begin(res),LA::vector_end(res)); } - /* The general case - * - * Matrix P=(p1, p2, ...) (each point as a column) - * Matrix Q=2*t(p2-p1,p3-p1, ...) (each vector as a line) - * Matrix M: QP, adding a line of 1 at the top - * Vector B: (1, p2^2-p1^2, p3^2-p1^2, ...) - * Solve ML=B, the center of the sphere is PL - * - * It would likely be faster to write P then transpose, multiply, - * etc instead of doing it by hand. - */ - // TODO: check for degenerate cases? - - typedef typename R_::Max_ambient_dimension D2; - typedef typename R_::LA::template Rebind_dimension::Other LAd; - typedef typename LAd::Square_matrix Matrix; - typedef typename LAd::Vector Vec; - typename Get_functor::type sp(this->kernel()); - Matrix m(k,k); - Vec b(k); - Vec l(k); - int j,i=0; - // We are doing a quadratic number of *f, which can be costly with transforming_iterator. - for(Iter f2=f;f2!=e;++f2,++i){ - Point const& p2=*f2; - b(i)=m(i,i)=sdo(p2); - j=0; - for(Iter f3=f;f3!=e;++f3,++j){ - m(j,i)=m(i,j)=sp(p2,*f3); + { + // The general case. ui=p(i+1)-p0, center-p0=c=sum ai*ui, c.ui=ui²/2, M*a=b with M symmetric + typedef typename Increment_dimension::type D2; + typedef typename R_::LA::template Rebind_dimension::Other LA; + typedef typename LA::Square_matrix Matrix; + typedef typename LA::Vector Vec; + typedef typename LA::Construct_vector CVec; + typename Get_functor::type tp(this->kernel()); + typename Get_functor::type sv(this->kernel()); + typename Get_functor::type dp(this->kernel()); + typename Get_functor::type sp(this->kernel()); + Matrix m(k-1,k-1); + Vec b = typename CVec::Dimension()(k-1); + std::vector vecs; vecs.reserve(k-1); + while(++f!=e) + vecs.emplace_back(dp(*f,p0)); + // Only need to fill the lower half + for(int i=0;i),(Point_tag),(Construct_ttag,Compute_point_cartesian_coordinate_tag,Scalar_product_tag,Squared_distance_to_origin_tag,Point_dimension_tag)); +CGAL_KD_DEFAULT_FUNCTOR(Construct_circumcenter_tag,(CartesianDKernelFunctors::Construct_circumcenter),(Point_tag,Vector_tag),(Construct_ttag,Compute_point_cartesian_coordinate_tag,Scalar_product_tag,Squared_distance_to_origin_tag,Point_dimension_tag,Translated_point_tag,Scaled_vector_tag,Difference_of_points_tag,Squared_length_tag)); namespace CartesianDKernelFunctors { template struct Squared_circumradius : Store_kernel { From e0513c243258e2f23231bb89fb439089f4c75e61 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Tue, 7 Jul 2020 16:54:04 +0200 Subject: [PATCH 4/8] Implement Eigen::NumTraits::highest Older versions of LDLT use it --- Number_types/include/CGAL/Interval_nt.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Number_types/include/CGAL/Interval_nt.h b/Number_types/include/CGAL/Interval_nt.h index ca77127d6e0..eb09e9ae3fd 100644 --- a/Number_types/include/CGAL/Interval_nt.h +++ b/Number_types/include/CGAL/Interval_nt.h @@ -1573,6 +1573,8 @@ namespace Eigen { static inline Real epsilon() { return 0; } static inline Real dummy_precision() { return 0; } + static inline Real highest() { return Real((std::numeric_limits::max)(), std::numeric_limits::infinity()); } + static inline Real lowest() { return Real(-std::numeric_limits::infinity(), std::numeric_limits::lowest()); } // Costs could depend on b. enum { From 00efb5ebef69e3de24fd24f46bdcde23e1dad11e Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Tue, 7 Jul 2020 17:51:29 +0200 Subject: [PATCH 5/8] Remove TAB --- .../CGAL/NewKernel_d/function_objects_cartesian.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h index d56fe9f632d..e1b1f6e32ae 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h @@ -689,19 +689,19 @@ template struct Construct_circumcenter : Store_kernel { Vec b = typename CVec::Dimension()(k-1); std::vector vecs; vecs.reserve(k-1); while(++f!=e) - vecs.emplace_back(dp(*f,p0)); + vecs.emplace_back(dp(*f,p0)); // Only need to fill the lower half for(int i=0;i Date: Tue, 7 Jul 2020 22:06:51 +0200 Subject: [PATCH 6/8] Rewrite Power_center_d. --- .../CGAL/NewKernel_d/Cartesian_LA_functors.h | 1 + .../CGAL/NewKernel_d/Types/Weighted_point.h | 109 +++++++++--------- NewKernel_d/test/NewKernel_d/Epick_d.cpp | 16 +++ 3 files changed, 69 insertions(+), 57 deletions(-) diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h b/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h index a379cee7ed9..b717a80bac4 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h @@ -224,6 +224,7 @@ template struct Scalar_product { }; template struct Squared_distance_to_origin_stored { + // What about weighted points, should they store sdo-w? CGAL_FUNCTOR_INIT_IGNORE(Squared_distance_to_origin_stored) typedef R_ R; typedef typename R::LA_vector LA; diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h b/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h index 8dfe875b8c8..e83af0db167 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h @@ -149,6 +149,7 @@ template struct Power_center : Store_kernel { typedef typename Get_type::type WPoint; typedef WPoint result_type; typedef typename Get_type::type Point; + typedef typename Get_type::type Vector; typedef typename Get_type::type FT; template result_type operator()(Iter f, Iter e)const{ @@ -168,11 +169,26 @@ template struct Power_center : Store_kernel { WPoint const& wp0 = *f; Point const& p0 = pdw(wp0); + FT const& w0 = pw(wp0); int d = pd(p0); int k = static_cast(std::distance(f,e)); - if (d+1 == k) - { - FT const& n0 = sdo(p0) - pw(wp0); + if (k == 1) return cwp(p0, -w0); + // TODO: check for degenerate cases? + if (k == 2) { + typename Get_functor::type dp(this->kernel()); + typename Get_functor::type sl(this->kernel()); + typename Get_functor::type tp(this->kernel()); + typename Get_functor::type sv(this->kernel()); + WPoint const& wp1 = *++f; + Point const& p1 = pdw(wp1); + FT const& w1 = pw(wp1); + Vector v01 = dp(p1, p0); + FT l01 = sl(v01); + FT coef = ((w0 - w1) / l01 + 1) / 2; + return cwp(tp(p0, sv(v01, coef)), CGAL::square(coef) * l01 - w0); + } + if (d+1 == k) { + FT const& n0 = sdo(p0) - w0; Matrix m(d,d); Vec b = typename CVec::Dimension()(d); // Write the point coordinates in lines. @@ -194,64 +210,43 @@ template struct Power_center : Store_kernel { FT const& r2 = pdp (wp0, center); return cwp(std::move(center), r2); } - else + { - /* - * Matrix P=(p1, p2, ...) (each point as a column) - * Matrix Q=2*t(p2-p1,p3-p1, ...) (each vector as a line) - * Matrix M: QP, adding a line of 1 at the top - * Vector B: (1, p2^2-p1^2, p3^2-p1^2, ...) plus weights - * Solve ML=B, the center of the sphere is PL - * - * It would likely be faster to write P then transpose, multiply, - * etc instead of doing it by hand. - */ - // TODO: check for degenerate cases? - - typedef typename R_::Max_ambient_dimension D2; - typedef typename R_::LA::template Rebind_dimension::Other LAd; - typedef typename LAd::Square_matrix Matrix; - typedef typename LAd::Vector Vec; + // The general case. ui=p(i+1)-p0, center-p0=c=sum ai*ui, c.2ui=ui²+w-w, M*a=b with M symmetric + typedef typename Increment_dimension::type D2; + typedef typename R_::LA::template Rebind_dimension::Other LA; + typedef typename LA::Square_matrix Matrix; + typedef typename LA::Vector Vec; + typedef typename LA::Construct_vector CVec; + typename Get_functor::type tp(this->kernel()); + typename Get_functor::type sv(this->kernel()); + typename Get_functor::type dp(this->kernel()); typename Get_functor::type sp(this->kernel()); - Matrix m(k,k); - Vec b(k); - Vec l(k); - int j,i=0; - for(Iter f2=f; f2!=e; ++f2,++i){ - WPoint const& wp = *f2; + typename Get_functor::type pv(this->kernel()); + typename Get_functor::type sl(this->kernel()); + + Matrix m(k-1,k-1); + Vec b = typename CVec::Dimension()(k-1); + std::vector vecs; vecs.reserve(k-1); + for(int i=0; ++f!=e; ++i) { + WPoint const& wp = *f; Point const& p = pdw(wp); - b(i) = m(i,i) = sdo(p) - pw(wp); - j=0; - for(Iter f3=f; f3!=e; ++f3,++j){ - // FIXME: scalar product of points ??? - m(j,i) = m(i,j) = sp(p,pdw(*f3)); - } + vecs.emplace_back(dp(p, p0)); + b[i] = w0 - pw(wp); } - for(i=1;i),(Weighted_point_tag),(In_flat_power_side_of_power_sphere_raw_tag,Point_drop_weight_tag,Point_weight_tag)); CGAL_KD_DEFAULT_FUNCTOR(Power_distance_tag,(CartesianDKernelFunctors::Power_distance),(Weighted_point_tag,Point_tag),(Squared_distance_tag,Point_drop_weight_tag,Point_weight_tag)); CGAL_KD_DEFAULT_FUNCTOR(Power_distance_to_point_tag,(CartesianDKernelFunctors::Power_distance_to_point),(Weighted_point_tag,Point_tag),(Squared_distance_tag,Point_drop_weight_tag,Point_weight_tag)); -CGAL_KD_DEFAULT_FUNCTOR(Power_center_tag,(CartesianDKernelFunctors::Power_center),(Weighted_point_tag,Point_tag),(Compute_point_cartesian_coordinate_tag,Construct_ttag,Construct_ttag,Point_dimension_tag,Squared_distance_to_origin_tag,Point_drop_weight_tag,Point_weight_tag,Power_distance_to_point_tag)); +CGAL_KD_DEFAULT_FUNCTOR(Power_center_tag,(CartesianDKernelFunctors::Power_center),(Weighted_point_tag,Point_tag,Vector_tag),(Compute_point_cartesian_coordinate_tag,Construct_ttag,Construct_ttag,Point_dimension_tag,Squared_distance_to_origin_tag,Point_drop_weight_tag,Point_weight_tag,Power_distance_to_point_tag,Translated_point_tag,Scaled_vector_tag,Difference_of_points_tag,Scalar_product_tag,Sum_of_vectors_tag,Squared_length_tag)); CGAL_KD_DEFAULT_FUNCTOR(Power_side_of_bounded_power_circumsphere_tag,(CartesianDKernelFunctors::Power_side_of_bounded_power_circumsphere),(Weighted_point_tag),(Power_distance_tag,Power_center_tag)); } // namespace CGAL #endif diff --git a/NewKernel_d/test/NewKernel_d/Epick_d.cpp b/NewKernel_d/test/NewKernel_d/Epick_d.cpp index c72d5d207e6..3a3b150432f 100644 --- a/NewKernel_d/test/NewKernel_d/Epick_d.cpp +++ b/NewKernel_d/test/NewKernel_d/Epick_d.cpp @@ -715,11 +715,18 @@ void test3(){ template void test4(){ typedef typename Ker::Point_d P; + typedef typename Ker::Weighted_point_d WP; typedef typename Ker::Construct_circumcenter_d CCc; typedef typename Ker::Equal_d E; + typedef typename Ker::Power_center_d PC; + typedef typename Ker::Power_distance_d PoD; + typedef typename Ker::Affine_rank_d AR; Ker k(4); CCc ccc Kinit(construct_circumcenter_d_object); E ed Kinit(equal_d_object); + PC pc Kinit(power_center_d_object); + PoD pod Kinit(power_distance_d_object); + AR ar Kinit(affine_rank_d_object); auto mkpt=[](auto...x){double l[]{(double)x...};return P(std::begin(l), std::end(l));}; P tab1[]={mkpt(15,20,40,80),mkpt(10,23,36,80),mkpt(10,20,40,85),mkpt(10,15,40,80),mkpt(13,20,40,76)}; assert(ed(ccc(tab1+0, tab1+5),mkpt(10,20,40,80))); @@ -727,6 +734,15 @@ void test4(){ assert(ed(ccc(tab2+0, tab2+4),mkpt(10,20,40,80))); P tab3[]={mkpt(15,20,35,80),mkpt(10,25,40,75),mkpt(13,24,37,76)}; assert(ed(ccc(tab3+0, tab3+3),mkpt(10,20,40,80))); + auto mkwpt=[](auto...x){double l[]{(double)x...};auto last=std::prev(std::end(l));return WP(P(std::begin(l), last),*last);}; + WP tab4[]={mkwpt(89,17,29,97,14),mkwpt(86,99,64,26,44),mkwpt(40,9,13,91,20),mkwpt(41,30,93,13,10),mkwpt(45,6,98,9,0),mkwpt(0,0,0,0,0)}; + for(int i=5;i>=1;--i){ + tab4[i]=pc(tab4+0, tab4+i); + for(int j=0;j >; template struct CGAL::Epick_d >; From 334f8ae105ccf1b335b711dc10f23bf906cb70e7 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Wed, 12 Aug 2020 07:40:09 +0200 Subject: [PATCH 7/8] Missing include --- .../include/CGAL/NewKernel_d/function_objects_cartesian.h | 1 + 1 file changed, 1 insertion(+) diff --git a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h index e1b1f6e32ae..3ca5fb3638e 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h @@ -27,6 +27,7 @@ #include #include #include +#include namespace CGAL { namespace CartesianDKernelFunctors { From 53ed991b5deef876e9a5757a3595b36e0052daac Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Thu, 27 Aug 2020 15:39:08 +0200 Subject: [PATCH 8/8] Fall back to LU instead of LDLT with old Eigen --- .../include/CGAL/NewKernel_d/Types/Weighted_point.h | 13 ++++++++++++- .../CGAL/NewKernel_d/function_objects_cartesian.h | 13 ++++++++++++- 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h b/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h index e83af0db167..6e970a5f13a 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h @@ -236,13 +236,24 @@ template struct Power_center : Store_kernel { } // Only need to fill the lower half for(int i = 0; i < k-1; ++i){ - for(int j = i; j < k-1; ++j) + for(int j = i; j < k-1; ++j){ m(j, i) = sp(vecs[i], vecs[j]); +#if ! EIGEN_VERSION_AT_LEAST(3, 3, 5) + m(i, j) = m(j, i); +#endif + } b[i] += m(i, i); b[i] /= 2; } // Assumes Eigen... +#if EIGEN_VERSION_AT_LEAST(3, 3, 5) Vec res = m.ldlt().solve(b); +#else + // Older versions of Eigen use 1/highest as tolerance, + // which we have no way to set to 0 for exact types. + // Use something slow but that should work. + Vec res = m.fullPivLu().solve(b); +#endif Vector to_center = sv(vecs[0], res[0]); for(int i=1;i struct Construct_circumcenter : Store_kernel { vecs.emplace_back(dp(*f,p0)); // Only need to fill the lower half for(int i=0;i