diff --git a/Installation/include/CGAL/config.h b/Installation/include/CGAL/config.h index a215cbbf4e5..d70989b2ffb 100644 --- a/Installation/include/CGAL/config.h +++ b/Installation/include/CGAL/config.h @@ -546,13 +546,9 @@ using std::max; #endif // Macro to specify a 'noreturn' attribute. -#if defined(__GNUG__) || __has_attribute(__noreturn__) -# define CGAL_NORETURN __attribute__ ((__noreturn__)) -#elif defined (_MSC_VER) -# define CGAL_NORETURN __declspec(noreturn) -#else -# define CGAL_NORETURN -#endif +// (This macro existed in CGAL before we switched to C++11. Let's keep +// the macro defined for backward-compatibility. That cannot harm.) +#define CGAL_NORETURN [[noreturn]] // Macro to specify [[no_unique_address]] if supported #if __has_cpp_attribute(no_unique_address) diff --git a/Kernel_23/include/CGAL/determinant.h b/Kernel_23/include/CGAL/determinant.h index fe14a2d1974..f470720d560 100644 --- a/Kernel_23/include/CGAL/determinant.h +++ b/Kernel_23/include/CGAL/determinant.h @@ -202,6 +202,204 @@ determinant( return m012345; } +template +RT +determinant( + const RT& a00, const RT& a01, const RT& a02, const RT& a03, const RT& a04, const RT& a05, const RT& a06, + const RT& a10, const RT& a11, const RT& a12, const RT& a13, const RT& a14, const RT& a15, const RT& a16, + const RT& a20, const RT& a21, const RT& a22, const RT& a23, const RT& a24, const RT& a25, const RT& a26, + const RT& a30, const RT& a31, const RT& a32, const RT& a33, const RT& a34, const RT& a35, const RT& a36, + const RT& a40, const RT& a41, const RT& a42, const RT& a43, const RT& a44, const RT& a45, const RT& a46, + const RT& a50, const RT& a51, const RT& a52, const RT& a53, const RT& a54, const RT& a55, const RT& a56, + const RT& a60, const RT& a61, const RT& a62, const RT& a63, const RT& a64, const RT& a65, const RT& a66) +{ +// First compute the det2x2 + const RT m01 = a00*a11 - a10*a01; + const RT m02 = a00*a21 - a20*a01; + const RT m03 = a00*a31 - a30*a01; + const RT m04 = a00*a41 - a40*a01; + const RT m05 = a00*a51 - a50*a01; + const RT m06 = a00*a61 - a60*a01; + + const RT m12 = a10*a21 - a20*a11; + const RT m13 = a10*a31 - a30*a11; + const RT m14 = a10*a41 - a40*a11; + const RT m15 = a10*a51 - a50*a11; + const RT m16 = a10*a61 - a60*a11; + + const RT m23 = a20*a31 - a30*a21; + const RT m24 = a20*a41 - a40*a21; + const RT m25 = a20*a51 - a50*a21; + const RT m26 = a20*a61 - a60*a21; + + const RT m34 = a30*a41 - a40*a31; + const RT m35 = a30*a51 - a50*a31; + const RT m36 = a30*a61 - a60*a31; + + const RT m45 = a40*a51 - a50*a41; + const RT m46 = a40*a61 - a60*a41; + + const RT m56 = a50*a61 - a60*a51; + +// Now compute the minors of rank 3 + const RT m012 = m01*a22 - m02*a12 + m12*a02; + const RT m013 = m01*a32 - m03*a12 + m13*a02; + const RT m014 = m01*a42 - m04*a12 + m14*a02; + const RT m015 = m01*a52 - m05*a12 + m15*a02; + const RT m016 = m01*a62 - m06*a12 + m16*a02; + + const RT m023 = m02*a32 - m03*a22 + m23*a02; + const RT m024 = m02*a42 - m04*a22 + m24*a02; + const RT m025 = m02*a52 - m05*a22 + m25*a02; + const RT m026 = m02*a62 - m06*a22 + m26*a02; + const RT m034 = m03*a42 - m04*a32 + m34*a02; + const RT m035 = m03*a52 - m05*a32 + m35*a02; + const RT m036 = m03*a62 - m06*a32 + m36*a02; + + const RT m045 = m04*a52 - m05*a42 + m45*a02; + const RT m046 = m04*a62 - m06*a42 + m46*a02; + + const RT m056 = m05*a62 - m06*a52 + m56*a02; + + + const RT m123 = m12*a32 - m13*a22 + m23*a12; + const RT m124 = m12*a42 - m14*a22 + m24*a12; + const RT m125 = m12*a52 - m15*a22 + m25*a12; + const RT m126 = m12*a62 - m16*a22 + m26*a12; + + const RT m134 = m13*a42 - m14*a32 + m34*a12; + const RT m135 = m13*a52 - m15*a32 + m35*a12; + const RT m136 = m13*a62 - m16*a32 + m36*a12; + + const RT m145 = m14*a52 - m15*a42 + m45*a12; + const RT m146 = m14*a62 - m16*a42 + m46*a12; + + const RT m156 = m15*a62 - m16*a52 + m56*a12; + + const RT m234 = m23*a42 - m24*a32 + m34*a22; + const RT m235 = m23*a52 - m25*a32 + m35*a22; + const RT m236 = m23*a62 - m26*a32 + m36*a22; + + const RT m245 = m24*a52 - m25*a42 + m45*a22; + const RT m246 = m24*a62 - m26*a42 + m46*a22; + + const RT m256 = m25*a62 - m26*a52 + m56*a22; + + const RT m345 = m34*a52 - m35*a42 + m45*a32; + const RT m346 = m34*a62 - m36*a42 + m46*a32; + + const RT m356 = m35*a62 - m36*a52 + m56*a32; + + const RT m456 = m45*a62 - m46*a52 + m56*a42; + +// Now compute the minors of rank 4 + const RT m0123 = m012*a33 - m013*a23 + m023*a13 - m123*a03; + + const RT m0124 = m012*a43 - m014*a23 + m024*a13 - m124*a03; + const RT m0125 = m012*a53 - m015*a23 + m025*a13 - m125*a03; + const RT m0126 = m012*a63 - m016*a23 + m026*a13 - m126*a03; + + const RT m0134 = m013*a43 - m014*a33 + m034*a13 - m134*a03; + const RT m0135 = m013*a53 - m015*a33 + m035*a13 - m135*a03; + const RT m0136 = m013*a63 - m016*a33 + m036*a13 - m136*a03; + + const RT m0145 = m014*a53 - m015*a43 + m045*a13 - m145*a03; + const RT m0146 = m014*a63 - m016*a43 + m046*a13 - m146*a03; + + const RT m0156 = m015*a63 - m016*a53 + m056*a13 - m156*a03; + + const RT m0234 = m023*a43 - m024*a33 + m034*a23 - m234*a03; + const RT m0235 = m023*a53 - m025*a33 + m035*a23 - m235*a03; + const RT m0236 = m023*a63 - m026*a33 + m036*a23 - m236*a03; + + const RT m0245 = m024*a53 - m025*a43 + m045*a23 - m245*a03; + const RT m0246 = m024*a63 - m026*a43 + m046*a23 - m246*a03; + + const RT m0256 = m025*a63 - m026*a53 + m056*a23 - m256*a03; + + const RT m0345 = m034*a53 - m035*a43 + m045*a33 - m345*a03; + const RT m0346 = m034*a63 - m036*a43 + m046*a33 - m346*a03; + + const RT m0356 = m035*a63 - m036*a53 + m056*a33 - m356*a03; + + const RT m0456 = m045*a63 - m046*a53 + m056*a43 - m456*a03; + + const RT m1234 = m123*a43 - m124*a33 + m134*a23 - m234*a13; + const RT m1235 = m123*a53 - m125*a33 + m135*a23 - m235*a13; + const RT m1236 = m123*a63 - m126*a33 + m136*a23 - m236*a13; + + const RT m1245 = m124*a53 - m125*a43 + m145*a23 - m245*a13; + const RT m1246 = m124*a63 - m126*a43 + m146*a23 - m246*a13; + + const RT m1256 = m125*a63 - m126*a53 + m156*a23 - m256*a13; + + const RT m1345 = m134*a53 - m135*a43 + m145*a33 - m345*a13; + const RT m1346 = m134*a63 - m136*a43 + m146*a33 - m346*a13; + + const RT m1356 = m135*a63 - m136*a53 + m156*a33 - m356*a13; + const RT m1456 = m145*a63 - m146*a53 + m156*a43 - m456*a13; + + const RT m2345 = m234*a53 - m235*a43 + m245*a33 - m345*a23; + const RT m2346 = m234*a63 - m236*a43 + m246*a33 - m346*a23; + + const RT m2356 = m235*a63 - m236*a53 + m256*a33 - m356*a23; + const RT m2456 = m245*a63 - m246*a53 + m256*a43 - m456*a23; + + const RT m3456 = m345*a63 - m346*a53 + m356*a43 - m456*a33; + + + // Now compute the minors of rank 5 + const RT m01234 = m0123*a44 - m0124*a34 + m0134*a24 - m0234*a14 + m1234*a04; + + const RT m01235 = m0123*a54 - m0125*a34 + m0135*a24 - m0235*a14 + m1235*a04; + + const RT m01236 = m0123*a64 - m0126*a34 + m0136*a24 - m0236*a14 + m1236*a04; + + const RT m01245 = m0124*a54 - m0125*a44 + m0145*a24 - m0245*a14 + m1245*a04; + const RT m01246 = m0124*a64 - m0126*a44 + m0146*a24 - m0246*a14 + m1246*a04; + + const RT m01256 = m0125*a64 - m0126*a54 + m0156*a24 - m0256*a14 + m1256*a04; + + const RT m01345 = m0134*a54 - m0135*a44 + m0145*a34 - m0345*a14 + m1345*a04; + const RT m01346 = m0134*a64 - m0136*a44 + m0146*a34 - m0346*a14 + m1346*a04; + + const RT m01356 = m0135*a64 - m0136*a54 + m0156*a34 - m0356*a14 + m1356*a04; + const RT m01456 = m0145*a64 - m0146*a54 + m0156*a44 - m0456*a14 + m1456*a04; + + const RT m02345 = m0234*a54 - m0235*a44 + m0245*a34 - m0345*a24 + m2345*a04; + const RT m02346 = m0234*a64 - m0236*a44 + m0246*a34 - m0346*a24 + m2346*a04; + + const RT m02356 = m0235*a64 - m0236*a54 + m0256*a34 - m0356*a24 + m2356*a04; + const RT m02456 = m0245*a64 - m0246*a54 + m0256*a44 - m0456*a24 + m2456*a04; + const RT m03456 = m0345*a64 - m0346*a54 + m0356*a44 - m0456*a34 + m3456*a04; + + const RT m12345 = m1234*a54 - m1235*a44 + m1245*a34 - m1345*a24 + m2345*a14; + const RT m12346 = m1234*a64 - m1236*a44 + m1246*a34 - m1346*a24 + m2346*a14; + + + const RT m12356 = m1235*a64 - m1236*a54 + m1256*a34 - m1356*a24 + m2356*a14; + const RT m12456 = m1245*a64 - m1246*a54 + m1256*a44 - m1456*a24 + m2456*a14; + const RT m13456 = m1345*a64 - m1346*a54 + m1356*a44 - m1456*a34 + m3456*a14; + + const RT m23456 = m2345*a64 - m2346*a54 + m2356*a44 - m2456*a34 + m3456*a24; + +// Now compute the minors of rank 6 + const RT m012345 = m01234*a55 - m01235*a45 + m01245*a35 - m01345*a25 + m02345*a15 - m12345*a05; + const RT m012346 = m01234*a65 - m01236*a45 + m01246*a35 - m01346*a25 + m02346*a15 - m12346*a05; + const RT m012356 = m01235*a65 - m01236*a55 + m01256*a35 - m01356*a25 + m02356*a15 - m12356*a05; + const RT m012456 = m01245*a65 - m01246*a55 + m01256*a45 - m01456*a25 + m02456*a15 - m12456*a05; + const RT m013456 = m01345*a65 - m01346*a55 + m01356*a45 - m01456*a35 + m03456*a15 - m13456*a05; + const RT m023456 = m02345*a65 - m02346*a55 + m02356*a45 - m02456*a35 + m03456*a25 - m23456*a05; + const RT m123456 = m12345*a65 - m12346*a55 + m12356*a45 - m12456*a35 + m13456*a25 - m23456*a15; + + + // Now compute the minors of rank 7 + const RT m0123456 = m012345 * a66 - m012346 * a56 + m012356 * a46 - m012456 * a36 + m013456 * a26 - m023456 * a16 + m123456 * a06; + + return m0123456; +} + + } //namespace CGAL #endif // CGAL_DETERMINANT_H diff --git a/Kernel_23/test/Kernel_23/determinant_77.cpp b/Kernel_23/test/Kernel_23/determinant_77.cpp new file mode 100644 index 00000000000..c48be9908d8 --- /dev/null +++ b/Kernel_23/test/Kernel_23/determinant_77.cpp @@ -0,0 +1,121 @@ +#include +#include + +#include + +#include + +using CGAL::determinant; + +template +RT det_77_alt(const RT& a00, const RT& a01, const RT& a02, const RT& a03, const RT& a04, const RT& a05, const RT& a06, + const RT& a10, const RT& a11, const RT& a12, const RT& a13, const RT& a14, const RT& a15, const RT& a16, + const RT& a20, const RT& a21, const RT& a22, const RT& a23, const RT& a24, const RT& a25, const RT& a26, + const RT& a30, const RT& a31, const RT& a32, const RT& a33, const RT& a34, const RT& a35, const RT& a36, + const RT& a40, const RT& a41, const RT& a42, const RT& a43, const RT& a44, const RT& a45, const RT& a46, + const RT& a50, const RT& a51, const RT& a52, const RT& a53, const RT& a54, const RT& a55, const RT& a56, + const RT& a60, const RT& a61, const RT& a62, const RT& a63, const RT& a64, const RT& a65, const RT& a66) +{ + const RT r1 = a06 * determinant( + 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, + a60, a61, a62, a63, a64, a65); + + const RT r2 = - a16 * determinant(a00, a01, a02, a03, a04, a05, + + a20, a21, a22, a23, a24, a25, + a30, a31, a32, a33, a34, a35, + a40, a41, a42, a43, a44, a45, + a50, a51, a52, a53, a54, a55, + a60, a61, a62, a63, a64, a65); + + const RT r3 = a26 * determinant(a00, a01, a02, a03, a04, a05, + a10, a11, a12, a13, a14, a15, + + a30, a31, a32, a33, a34, a35, + a40, a41, a42, a43, a44, a45, + a50, a51, a52, a53, a54, a55, + a60, a61, a62, a63, a64, a65); + + const RT r4 = - a36 * determinant(a00, a01, a02, a03, a04, a05, + a10, a11, a12, a13, a14, a15, + a20, a21, a22, a23, a24, a25, + + a40, a41, a42, a43, a44, a45, + a50, a51, a52, a53, a54, a55, + a60, a61, a62, a63, a64, a65); + + const RT r5 = a46 * determinant(a00, a01, a02, a03, a04, a05, + a10, a11, a12, a13, a14, a15, + a20, a21, a22, a23, a24, a25, + a30, a31, a32, a33, a34, a35, + + a50, a51, a52, a53, a54, a55, + a60, a61, a62, a63, a64, a65); + + const RT r6 = - a56 * determinant(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, + + a60, a61, a62, a63, a64, a65); + + const RT r7 = a66 * determinant(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 + + ); + + return r1 + r2 + r3 + r4 + r5 + r6 + r7; +} + +int main(int, char**) +{ + assert(determinant(4, 5, 1, 4, 6, 3, 1, + 4, 3, 6, 4, 2, 7, 3, + 6, 3, 3, 6, 2, 4, 5, + 1, 4, 3, 5, 5, 6 ,1, + 1, 3, 2, 7, 9, 6, 1, + 7, 6, 5, 4, 6, 2, 2, + 2, 3, 5, 7, 4, 3, 3) == 763); + + CGAL::Random rnd; + std::cout << "Seed: " << rnd.get_seed() << std::endl; + + for(int k=0; k<100; ++k) + { + std::array, 7> mat; + for(int i=0; i<7; ++i) + for(int j=0; j<7; ++j) + mat[i][j] = rnd.get_int(-18, 18); + + const int det_1 = determinant(mat[0][0], mat[0][1], mat[0][2], mat[0][3], mat[0][4], mat[0][5], mat[0][6], + mat[1][0], mat[1][1], mat[1][2], mat[1][3], mat[1][4], mat[1][5], mat[1][6], + mat[2][0], mat[2][1], mat[2][2], mat[2][3], mat[2][4], mat[2][5], mat[2][6], + mat[3][0], mat[3][1], mat[3][2], mat[3][3], mat[3][4], mat[3][5], mat[3][6], + mat[4][0], mat[4][1], mat[4][2], mat[4][3], mat[4][4], mat[4][5], mat[4][6], + mat[0][0], mat[5][1], mat[5][2], mat[5][3], mat[5][4], mat[5][5], mat[5][6], + mat[6][0], mat[6][1], mat[6][2], mat[6][3], mat[6][4], mat[6][5], mat[6][6]); + + const int det_2 = det_77_alt(mat[0][0], mat[0][1], mat[0][2], mat[0][3], mat[0][4], mat[0][5], mat[0][6], + mat[1][0], mat[1][1], mat[1][2], mat[1][3], mat[1][4], mat[1][5], mat[1][6], + mat[2][0], mat[2][1], mat[2][2], mat[2][3], mat[2][4], mat[2][5], mat[2][6], + mat[3][0], mat[3][1], mat[3][2], mat[3][3], mat[3][4], mat[3][5], mat[3][6], + mat[4][0], mat[4][1], mat[4][2], mat[4][3], mat[4][4], mat[4][5], mat[4][6], + mat[0][0], mat[5][1], mat[5][2], mat[5][3], mat[5][4], mat[5][5], mat[5][6], + mat[6][0], mat[6][1], mat[6][2], mat[6][3], mat[6][4], mat[6][5], mat[6][6]); + + std::cout << "dets: " << det_1 << " " << det_2 << std::endl; + assert(det_1 == det_2); + } + + std::cout << "Done!" << std::endl; + return EXIT_SUCCESS; +} diff --git a/NewKernel_d/include/CGAL/Epick_d.h b/NewKernel_d/include/CGAL/Epick_d.h index 54b6633ebec..39575dd034a 100644 --- a/NewKernel_d/include/CGAL/Epick_d.h +++ b/NewKernel_d/include/CGAL/Epick_d.h @@ -40,7 +40,12 @@ struct Epick_d_help1 }; #undef CGAL_BASE #define CGAL_BASE \ - Cartesian_static_filters,Epick_d_help2 > + Cartesian_filter_K< \ + Epick_d_help1, \ + Cartesian_base_d, \ + Cartesian_base_d::Type, Dim>, \ + typename Functors_without_division::type \ + > template struct Epick_d_help2 : CGAL_BASE @@ -49,10 +54,21 @@ struct Epick_d_help2 constexpr Epick_d_help2(int d):CGAL_BASE(d){} }; #undef CGAL_BASE +#define CGAL_BASE \ + Cartesian_static_filters,Epick_d_help3 > + +template +struct Epick_d_help3 +: CGAL_BASE +{ + constexpr Epick_d_help3(){} + constexpr Epick_d_help3(int d):CGAL_BASE(d){} +}; +#undef CGAL_BASE #define CGAL_BASE \ Kernel_d_interface< \ Cartesian_wrap< \ - Epick_d_help2, \ + Epick_d_help3, \ Epick_d > > template struct Epick_d 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 1bfc92a88c8..a379cee7ed9 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h @@ -58,7 +58,7 @@ template struct Construct_LA_vector } template typename std::enable_if::value && - std::is_same, Dimension>::value, + std::is_same, Dimension>::value, result_type>::type operator()(U&&...u)const{ return typename Constructor::Values()(std::forward(u)...); @@ -66,7 +66,7 @@ template struct Construct_LA_vector //template::value>::type,class=typename std::enable_if<(sizeof...(U)==static_dim+1)>::type,class=void> template typename std::enable_if::value && - std::is_same, Dimension>::value, + std::is_same, Dimension>::value, result_type>::type operator()(U&&...u)const{ return Apply_to_last_then_rest()(typename Constructor::Values_divide(),std::forward(u)...); diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_filter_K.h b/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_filter_K.h index 3496330a7ac..ed164e844da 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_filter_K.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_filter_K.h @@ -20,23 +20,49 @@ namespace CGAL { -template < typename Base_, typename AK_, typename EK_ > -struct Cartesian_filter_K : public Base_, - private Store_kernel, private Store_kernel2 + // It would be nicer to write the table in the other direction: Orientation_of_points_tag is good up to 6, Side_of_oriented_sphere_tag up to 5, etc. +template struct Functors_without_division { typedef typeset<> type; }; +template<> struct Functors_without_division > { + typedef typeset type; +}; +template<> struct Functors_without_division > { + typedef typeset type; +}; +template<> struct Functors_without_division > { + typedef typeset type; +}; +template<> struct Functors_without_division > { + typedef typeset type; +}; +template<> struct Functors_without_division > { + typedef typeset type; +}; +template<> struct Functors_without_division > { + typedef typeset type; +}; + +// FIXME: +// - Is_exact (which should be renamed to Uses_no_arithmetic) predicates should not be filtered +// - Functors_without_division should be defined near/in the actual functors + +template < typename Base_, typename AK_, typename EK_, typename Pred_list = typeset_all > +struct Cartesian_filter_K : public Base_ { + CGAL_NO_UNIQUE_ADDRESS Store_kernel sak; + CGAL_NO_UNIQUE_ADDRESS Store_kernel sek; constexpr Cartesian_filter_K(){} constexpr Cartesian_filter_K(int d):Base_(d){} //FIXME: or do we want an instance of AK and EK belonging to this kernel, //instead of a reference to external ones? - constexpr Cartesian_filter_K(AK_ const&a,EK_ const&b):Base_(),Store_kernel(a),Store_kernel2(b){} - constexpr Cartesian_filter_K(int d,AK_ const&a,EK_ const&b):Base_(d),Store_kernel(a),Store_kernel2(b){} + constexpr Cartesian_filter_K(AK_ const&a,EK_ const&b):Base_(),sak(a),sek(b){} + constexpr Cartesian_filter_K(int d,AK_ const&a,EK_ const&b):Base_(d),sak(a),sek(b){} typedef Base_ Kernel_base; typedef AK_ AK; typedef EK_ EK; typedef typename Store_kernel::reference_type AK_rt; - AK_rt approximate_kernel()const{return this->kernel();} - typedef typename Store_kernel2::reference2_type EK_rt; - EK_rt exact_kernel()const{return this->kernel2();} + AK_rt approximate_kernel()const{return sak.kernel();} + typedef typename Store_kernel::reference_type EK_rt; + EK_rt exact_kernel()const{return sek.kernel();} // MSVC is too dumb to perform the empty base optimization. typedef boost::mpl::and_< @@ -52,12 +78,12 @@ struct Cartesian_filter_K : public Base_, // TODO: only fix some types, based on some criterion? template struct Type : Get_type {}; - template::type> struct Functor : + template::type, bool=Pred_list::template contains::value> struct Functor : Inherit_functor {}; - template struct Functor { + template struct Functor { typedef typename Get_functor::type AP; typedef typename Get_functor::type EP; - typedef Filtered_predicate2 type; + typedef Filtered_predicate2 type; }; // TODO: // template struct Functor : diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Filtered_predicate2.h b/NewKernel_d/include/CGAL/NewKernel_d/Filtered_predicate2.h index a6190e42b5a..c3d2281ad96 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Filtered_predicate2.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Filtered_predicate2.h @@ -40,7 +40,7 @@ namespace CGAL { // - Some caching could be done at the Point_2 level. -template +template class Filtered_predicate2 { //TODO: pack (at least use a tuple) @@ -65,7 +65,6 @@ public: Filtered_predicate2() {} - template Filtered_predicate2(const K& k) : ep(k.exact_kernel()), ap(k.approximate_kernel()), c2e(k,k.exact_kernel()), c2a(k,k.approximate_kernel()) {} diff --git a/NewKernel_d/include/CGAL/NewKernel_d/LA_eigen/LA.h b/NewKernel_d/include/CGAL/NewKernel_d/LA_eigen/LA.h index 4dd8b5ca845..3d4afd905e5 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/LA_eigen/LA.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/LA_eigen/LA.h @@ -21,6 +21,8 @@ #include #include #include +#include +#include namespace CGAL { @@ -100,15 +102,72 @@ template struct LA_eigen { return (int)v.cols(); } - template static NT determinant(Mat_ const&m,bool=false){ + template static NT determinant_aux [[noreturn]] (Mat_ const&, Tag_true) { + CGAL_error(); + } + template static NT determinant_aux(Mat_ const& m, Tag_false) { return m.determinant(); } + // TODO: https://gitlab.com/libeigen/eigen/-/issues/1782 + // Implement a version of (sign_of_)determinant that works + // without (inexact) division in any dimension + template static NT determinant(Mat_ const&m,bool=false){ + switch(m.rows()){ + //case 0: + // return 1; + case 1: + return m(0,0); + case 2: + return CGAL::determinant( + m(0,0),m(0,1), + m(1,0),m(1,1)); + case 3: + return CGAL::determinant( + m(0,0),m(0,1),m(0,2), + m(1,0),m(1,1),m(1,2), + m(2,0),m(2,1),m(2,2)); + case 4: + return CGAL::determinant( + m(0,0),m(0,1),m(0,2),m(0,3), + m(1,0),m(1,1),m(1,2),m(1,3), + m(2,0),m(2,1),m(2,2),m(2,3), + m(3,0),m(3,1),m(3,2),m(3,3)); + case 5: + return CGAL::determinant( + m(0,0),m(0,1),m(0,2),m(0,3),m(0,4), + m(1,0),m(1,1),m(1,2),m(1,3),m(1,4), + m(2,0),m(2,1),m(2,2),m(2,3),m(2,4), + m(3,0),m(3,1),m(3,2),m(3,3),m(3,4), + m(4,0),m(4,1),m(4,2),m(4,3),m(4,4)); + case 6: + return CGAL::determinant( + m(0,0),m(0,1),m(0,2),m(0,3),m(0,4),m(0,5), + m(1,0),m(1,1),m(1,2),m(1,3),m(1,4),m(1,5), + m(2,0),m(2,1),m(2,2),m(2,3),m(2,4),m(2,5), + m(3,0),m(3,1),m(3,2),m(3,3),m(3,4),m(3,5), + m(4,0),m(4,1),m(4,2),m(4,3),m(4,4),m(4,5), + m(5,0),m(5,1),m(5,2),m(5,3),m(5,4),m(5,5)); + case 7: + return CGAL::determinant( + m(0,0),m(0,1),m(0,2),m(0,3),m(0,4),m(0,5),m(0,6), + m(1,0),m(1,1),m(1,2),m(1,3),m(1,4),m(1,5),m(1,6), + m(2,0),m(2,1),m(2,2),m(2,3),m(2,4),m(2,5),m(2,6), + m(3,0),m(3,1),m(3,2),m(3,3),m(3,4),m(3,5),m(3,6), + m(4,0),m(4,1),m(4,2),m(4,3),m(4,4),m(4,5),m(4,6), + m(5,0),m(5,1),m(5,2),m(5,3),m(5,4),m(5,5),m(5,6), + m(6,0),m(6,1),m(6,2),m(6,3),m(6,4),m(6,5),m(6,6)); + default: + return determinant_aux(m, Boolean_tag<(Mat_::MaxRowsAtCompileTime >= 1 && Mat_::MaxRowsAtCompileTime <= 7)>()); + + } + } + template static typename Same_uncertainty_nt::type sign_of_determinant(Mat_ const&m,bool=false) { - return CGAL::sign(m.determinant()); + return CGAL::sign(LA_eigen::determinant(m)); } template static int rank(Mat_ const&m){ diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Lazy_cartesian.h b/NewKernel_d/include/CGAL/NewKernel_d/Lazy_cartesian.h index d9f4112cf68..56ecd5df953 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Lazy_cartesian.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Lazy_cartesian.h @@ -56,10 +56,10 @@ namespace internal { // Whenever a construction takes iterator pairs as input, whether they point to double of Lazy objects, copy the ranges inside the lazy result so they are available for update_exact(). We analyze the input to try and guess where iterator pairs are. I would prefer if each functor had a specific signature (no overload in this layer) so we wouldn't have to guess. namespace Lazy_internal { templatestruct typelist{}; -templatestruct arg_i{}; -templatestruct arg_i_begin{}; -templatestruct arg_i_end{}; -templatestruct arg_i_ip1_range{}; +templatestruct arg_i{}; +templatestruct arg_i_begin{}; +templatestruct arg_i_end{}; +templatestruct arg_i_ip1_range{}; templatestruct analyze_args; templatestruct analyze_args> { typedef T creator; @@ -73,24 +73,24 @@ struct analyze_args,typelist,typelist,std::enab analyze_args>,typelist,arg_i_end>,typelist> {}; template using analyze_args_for_lazy = analyze_args,typelist<>,typelist>; templatestruct extract1; -templatestruct extract1,T>:std::tuple_element{}; -templatestruct extract1,T>{ +templatestruct extract1,T>:std::tuple_element{}; +templatestruct extract1,T>{ typedef std::tuple_element_t E; typedef std::remove_cv_t> It; typedef typename std::iterator_traits::value_type element_type; // TODO: find a way to use an array of the right size, at least for the most frequent constructions typedef std::vector type; }; -templatedecltype(auto) +templatedecltype(auto) do_extract(arg_i,std::tupleconst&t) {return std::get(t);} -templatedecltype(auto) +templatedecltype(auto) do_extract(arg_i_begin,std::tupleconst&t) {return std::begin(std::get(t));} -templatedecltype(auto) +templatedecltype(auto) do_extract(arg_i_end,std::tupleconst&t) {return std::end(std::get(t));} -templatedecltype(auto) +templatedecltype(auto) do_extract(arg_i_ip1_range,std::tupleconst&t) { typedef std::tuple L; @@ -294,7 +294,7 @@ struct Lazy_cartesian : template struct Functor { typedef typename Get_functor::type FA; typedef typename Get_functor::type FE; - typedef Filtered_predicate2 type; + typedef Filtered_predicate2 type; }; template struct Functor { typedef Lazy_construction2 type; diff --git a/NewKernel_d/include/CGAL/NewKernel_d/store_kernel.h b/NewKernel_d/include/CGAL/NewKernel_d/store_kernel.h index b9bfe61749f..5249dff2fcb 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/store_kernel.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/store_kernel.h @@ -13,12 +13,12 @@ #define CGAL_STORE_KERNEL_H #include -#include +#include namespace CGAL { namespace internal { BOOST_MPL_HAS_XXX_TRAIT_DEF(Do_not_store_kernel) -template::value,bool=has_Do_not_store_kernel::value> struct Do_not_store_kernel { +template::value,bool=has_Do_not_store_kernel::value> struct Do_not_store_kernel { enum { value=false }; typedef Tag_false type; }; diff --git a/NewKernel_d/include/CGAL/argument_swaps.h b/NewKernel_d/include/CGAL/argument_swaps.h index 85915098e26..a6f2cd8d794 100644 --- a/NewKernel_d/include/CGAL/argument_swaps.h +++ b/NewKernel_d/include/CGAL/argument_swaps.h @@ -18,9 +18,9 @@ namespace CGAL { namespace internal { -template struct Apply_to_last_then_rest_; +template struct Apply_to_last_then_rest_; -template +template struct Apply_to_last_then_rest_ { typedef typename Apply_to_last_then_rest_::result_type result_type; inline result_type operator()(F&&f,T&&t,U&&...u)const{ diff --git a/NewKernel_d/include/CGAL/typeset.h b/NewKernel_d/include/CGAL/typeset.h index 9c494a412a9..77c37c89c39 100644 --- a/NewKernel_d/include/CGAL/typeset.h +++ b/NewKernel_d/include/CGAL/typeset.h @@ -41,11 +41,17 @@ namespace CGAL { template using contains = std::false_type; template using add = typeset; }; + struct typeset_all { + typedef typeset_all type; + template using contains = std::true_type; + template using add = typeset_all; + }; template struct typeset_union_ : typeset_union_::type, typename T2::tail> {}; template struct typeset_union_ > : T {}; + template struct typeset_union_ : typeset_all {}; template struct typeset_intersection_ { @@ -55,8 +61,8 @@ namespace CGAL { std::conditional::value, typename U::template add::type, U>::type type; }; - template - struct typeset_intersection_,T> : typeset<> {}; + template struct typeset_intersection_, T> : typeset<> {}; + template struct typeset_intersection_ : T {}; template using typeset_union = typename typeset_union_::type; diff --git a/Number_types/include/CGAL/Mpzf.h b/Number_types/include/CGAL/Mpzf.h index d775431ebcf..2612e29b3a0 100644 --- a/Number_types/include/CGAL/Mpzf.h +++ b/Number_types/include/CGAL/Mpzf.h @@ -340,18 +340,62 @@ struct Mpzf { exp=x.exp; if(size!=0) mpn_copyi(data(),x.data(),asize); } -#if !defined(CGAL_MPZF_USE_CACHE) +#if defined(CGAL_MPZF_USE_CACHE) + Mpzf(Mpzf&& x)noexcept:size(x.size),exp(x.exp){ + auto xd = x.data(); + while(*--xd==0); + if (xd != x.cache) { + data() = x.data(); + x.init(); + } else { + init(); + if(size!=0) mpn_copyi(data(),x.data(),std::abs(size)); + } + x.size = 0; + } + Mpzf& operator=(Mpzf&& x)noexcept{ + if (this == &x) return *this; // is this needed? + size = x.size; + exp = x.exp; + auto xd = x.data(); + auto td = data(); + while(*--xd==0); + while(*--td==0); + if (xd != x.cache) { + data() = x.data(); + if (td != cache) { + pool::push(td+1); + // should we instead give it to x in case x is reused? + // x.data() = td + 1; + } + x.init(); + } else { + // In some cases data points in the middle of the buffer, reset it + data() = td + 1; + if(size!=0) mpn_copyi(data(),x.data(),std::abs(size)); + } + x.size = 0; + return *this; + } +#else Mpzf(Mpzf&& x):data_(x.data()),size(x.size),exp(x.exp){ x.init(); // yes, that's a shame... x.size = 0; x.exp = 0; } - Mpzf& operator=(Mpzf&& x){ - std::swap(size,x.size); + Mpzf& operator=(Mpzf&& x)noexcept{ + size = x.size; + // In case something tries to read it, size needs to be smaller than data + x.size = 0; exp = x.exp; std::swap(data(),x.data()); return *this; } + friend void swap(Mpzf&a, Mpzf&b)noexcept{ + std::swap(a.size, b.size); + std::swap(a.exp, b.exp); + std::swap(a.data(), b.data()); + } friend Mpzf operator-(Mpzf&& x){ Mpzf ret = std::move(x); ret.size = -ret.size; @@ -534,6 +578,12 @@ struct Mpzf { friend bool operator!=(Mpzf const&a, Mpzf const&b){ return !(a==b); } + friend Mpzf const& min BOOST_PREVENT_MACRO_SUBSTITUTION (Mpzf const&a, Mpzf const&b){ + return (b + +#ifdef CGAL_USE_GMP +# include +#endif +#ifdef CGAL_HAS_MPZF + +#include +#include +#include +#include + +int main() { + { + typedef CGAL::Mpzf NT; + typedef CGAL::Integral_domain_without_division_tag Tag; + typedef CGAL::Tag_true Is_exact; + + CGAL::test_algebraic_structure(); + CGAL::test_algebraic_structure(NT(4),NT(6),NT(15)); + CGAL::test_algebraic_structure(NT(-4),NT(6),NT(15)); + CGAL::test_algebraic_structure(NT(4),NT(-6),NT(15)); + CGAL::test_algebraic_structure(NT(-4),NT(-6),NT(15)); + CGAL::test_algebraic_structure(NT(4),NT(6),NT(-15)); + CGAL::test_algebraic_structure(NT(-4),NT(6), NT(15)); + CGAL::test_algebraic_structure(NT(4),NT(-6),NT(-15)); + CGAL::test_algebraic_structure(NT(-4),NT(-6),NT(-15)); + + CGAL::test_real_embeddable(); + } + + return 0; +} + +#else +int main() +{ + return 0; +} +#endif //CGAL_HAS_MPZF diff --git a/STL_Extension/include/CGAL/assertions.h b/STL_Extension/include/CGAL/assertions.h index 57698d51da9..caa0237fc7d 100644 --- a/STL_Extension/include/CGAL/assertions.h +++ b/STL_Extension/include/CGAL/assertions.h @@ -54,9 +54,9 @@ namespace CGAL { // ===================== // failure functions // ----------------- -CGAL_EXPORT CGAL_NORETURN void assertion_fail ( const char*, const char*, int, const char* = "") ; -CGAL_EXPORT CGAL_NORETURN void precondition_fail ( const char*, const char*, int, const char* = "") ; -CGAL_EXPORT CGAL_NORETURN void postcondition_fail ( const char*, const char*, int, const char* = "") ; +[[noreturn]] CGAL_EXPORT void assertion_fail ( const char*, const char*, int, const char* = "") ; +[[noreturn]] CGAL_EXPORT void precondition_fail ( const char*, const char*, int, const char* = "") ; +[[noreturn]] CGAL_EXPORT void postcondition_fail ( const char*, const char*, int, const char* = "") ; // warning function // ----------------