Merge pull request #4557 from lrineau/NewKernel_d-det-GF-CGAL-5.0

NewKernel_d: Accelerate by using Mpzf instead of Gmpq (CGAL-5.0 and later)
2020-05-27 18:27:56 +02:00 · 2020-05-27 18:27:56 +02:00 · 7173530dad
parent c88d8dc88b c15d63000f
commit 7173530dad
15 changed files with 563 additions and 49 deletions
--- 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)
--- a/Kernel_23/include/CGAL/determinant.h
+++ b/Kernel_23/include/CGAL/determinant.h
@ -202,6 +202,204 @@ determinant(
  return m012345;
 }

+template <class RT>
+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
--- a/Kernel_23/test/Kernel_23/determinant_77.cpp
+++ b/Kernel_23/test/Kernel_23/determinant_77.cpp
@ -0,0 +1,121 @@
+#include <CGAL/determinant.h>
+#include <CGAL/Random.h>
+
+#include <boost/lexical_cast.hpp>
+
+#include <iostream>
+
+using CGAL::determinant;
+
+template <typename RT>
+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<int>(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<std::array<int, 7>, 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<int>(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<int>(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;
+}
--- 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<Dim,Epick_d_help1<Dim>,Epick_d_help2<Dim> >
+  Cartesian_filter_K< \
+    Epick_d_help1<Dim>, \
+    Cartesian_base_d<Interval_nt_advanced, Dim>, \
+    Cartesian_base_d<internal::Exact_ring_selector<double>::Type, Dim>, \
+    typename Functors_without_division<Dim>::type \
+  >
 template<class Dim>
 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<Dim,Epick_d_help2<Dim>,Epick_d_help3<Dim> >
+
+template<class Dim>
+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<Dim>, \
+      Epick_d_help3<Dim>, \
      Epick_d<Dim> > >
 template<class Dim>
 struct Epick_d
--- 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<class R_,class Zero_> struct Construct_LA_vector
        }
        template<class...U>
        typename std::enable_if<Constructible_from_each<RT,U...>::value &&
-                std::is_same<Dimension_tag<sizeof...(U)>, Dimension>::value,
+                std::is_same<Dimension_tag<int(sizeof...(U))>, Dimension>::value,
          result_type>::type
        operator()(U&&...u)const{
                return typename Constructor::Values()(std::forward<U>(u)...);
@ -66,7 +66,7 @@ template<class R_,class Zero_> struct Construct_LA_vector
        //template<class...U,class=typename std::enable_if<Constructible_from_each<RT,U...>::value>::type,class=typename std::enable_if<(sizeof...(U)==static_dim+1)>::type,class=void>
        template<class...U>
        typename std::enable_if<Constructible_from_each<RT,U...>::value &&
-                std::is_same<Dimension_tag<sizeof...(U)-1>, Dimension>::value,
+                std::is_same<Dimension_tag<int(sizeof...(U)-1)>, Dimension>::value,
          result_type>::type
        operator()(U&&...u)const{
                return Apply_to_last_then_rest()(typename Constructor::Values_divide(),std::forward<U>(u)...);
--- 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<AK_>, private Store_kernel2<EK_>
+  // 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<class> struct Functors_without_division { typedef typeset<> type; };
+template<> struct Functors_without_division<Dimension_tag<1> > {
+  typedef typeset<Orientation_of_points_tag, Side_of_oriented_sphere_tag> type;
+};
+template<> struct Functors_without_division<Dimension_tag<2> > {
+  typedef typeset<Orientation_of_points_tag, Side_of_oriented_sphere_tag> type;
+};
+template<> struct Functors_without_division<Dimension_tag<3> > {
+  typedef typeset<Orientation_of_points_tag, Side_of_oriented_sphere_tag> type;
+};
+template<> struct Functors_without_division<Dimension_tag<4> > {
+  typedef typeset<Orientation_of_points_tag, Side_of_oriented_sphere_tag> type;
+};
+template<> struct Functors_without_division<Dimension_tag<5> > {
+  typedef typeset<Orientation_of_points_tag, Side_of_oriented_sphere_tag> type;
+};
+template<> struct Functors_without_division<Dimension_tag<6> > {
+  typedef typeset<Orientation_of_points_tag, Side_of_oriented_sphere_tag> 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<AK_> sak;
+    CGAL_NO_UNIQUE_ADDRESS Store_kernel<EK_> 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<AK_>(a),Store_kernel2<EK_>(b){}
-    constexpr Cartesian_filter_K(int d,AK_ const&a,EK_ const&b):Base_(d),Store_kernel<AK_>(a),Store_kernel2<EK_>(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<AK_>::reference_type AK_rt;
-    AK_rt approximate_kernel()const{return this->kernel();}
-    typedef typename Store_kernel2<EK_>::reference2_type EK_rt;
-    EK_rt exact_kernel()const{return this->kernel2();}
+    AK_rt approximate_kernel()const{return sak.kernel();}
+    typedef typename Store_kernel<EK_>::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<class T> struct Type : Get_type<Kernel_base,T> {};

-    template<class T,class D=void,class=typename Get_functor_category<Cartesian_filter_K,T>::type> struct Functor :
+    template<class T,class D=void,class=typename Get_functor_category<Cartesian_filter_K,T>::type, bool=Pred_list::template contains<T>::value> struct Functor :
      Inherit_functor<Kernel_base,T,D> {};
-    template<class T,class D> struct Functor<T,D,Predicate_tag> {
+    template<class T,class D> struct Functor<T,D,Predicate_tag,true> {
      typedef typename Get_functor<AK, T>::type AP;
      typedef typename Get_functor<EK, T>::type EP;
-      typedef Filtered_predicate2<EP,AP,C2E,C2A> type;
+      typedef Filtered_predicate2<Cartesian_filter_K,EP,AP,C2E,C2A> type;
    };
 // TODO:
 //    template<class T> struct Functor<T,No_filter_tag,Predicate_tag> :
--- 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 <class EP, class AP, class C2E, class C2A, bool Protection = true>
+template <class K, class EP, class AP, class C2E, class C2A, bool Protection = true>
 class Filtered_predicate2
 {
 //TODO: pack (at least use a tuple)
@ -65,7 +65,6 @@ public:
  Filtered_predicate2()
  {}

-  template <class K>
  Filtered_predicate2(const K& k)
    : ep(k.exact_kernel()), ap(k.approximate_kernel()), c2e(k,k.exact_kernel()), c2a(k,k.approximate_kernel())
  {}
--- 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 <Eigen/Dense>
 #include <CGAL/NewKernel_d/LA_eigen/constructors.h>
 #include <CGAL/iterator_from_indices.h>
+#include <CGAL/determinant.h>
+#include <CGAL/assertions.h>

 namespace CGAL {

@ -100,15 +102,72 @@ template<class NT_,class Dim_,class Max_dim_=Dim_> struct LA_eigen {
                return (int)v.cols();
        }

-        template<class Mat_> static NT determinant(Mat_ const&m,bool=false){
+        template<class Mat_> static NT determinant_aux [[noreturn]] (Mat_ const&, Tag_true) {
+                CGAL_error();
+        }
+        template<class Mat_> 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<class Mat_> 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<class Mat_> static typename
        Same_uncertainty_nt<CGAL::Sign, NT>::type
        sign_of_determinant(Mat_ const&m,bool=false)
        {
-                return CGAL::sign(m.determinant());
+                return CGAL::sign(LA_eigen::determinant(m));
        }

        template<class Mat_> static int rank(Mat_ const&m){
--- 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 {
 template<class...>struct typelist{};
-template<int>struct arg_i{};
-template<int>struct arg_i_begin{};
-template<int>struct arg_i_end{};
-template<int>struct arg_i_ip1_range{};
+template<std::size_t>struct arg_i{};
+template<std::size_t>struct arg_i_begin{};
+template<std::size_t>struct arg_i_end{};
+template<std::size_t>struct arg_i_ip1_range{};
 template<class,class,class,class=void>struct analyze_args;
 template<class T,class U>struct analyze_args<T,U,typelist<>> {
  typedef T creator;
@ -73,24 +73,24 @@ struct analyze_args<typelist<T...>,typelist<U...>,typelist<It,It,W...>,std::enab
 analyze_args<typelist<T...,arg_i_ip1_range<sizeof...(U)>>,typelist<U...,arg_i_begin<sizeof...(T)>,arg_i_end<sizeof...(T)>>,typelist<W...>> {};
 template<class...T> using analyze_args_for_lazy = analyze_args<typelist<>,typelist<>,typelist<T...>>;
 template<class,class>struct extract1;
-template<int i,class T>struct extract1<arg_i<i>,T>:std::tuple_element<i,T>{};
-template<int i,class T>struct extract1<arg_i_ip1_range<i>,T>{
+template<std::size_t i,class T>struct extract1<arg_i<i>,T>:std::tuple_element<i,T>{};
+template<std::size_t i,class T>struct extract1<arg_i_ip1_range<i>,T>{
  typedef std::tuple_element_t<i,T> E;
  typedef std::remove_cv_t<std::remove_reference_t<E>> It;
  typedef typename std::iterator_traits<It>::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<element_type> type;
 };
-template<int i,class...T>decltype(auto)
+template<std::size_t i,class...T>decltype(auto)
 do_extract(arg_i<i>,std::tuple<T...>const&t)
 {return std::get<i>(t);}
-template<int i,class...T>decltype(auto)
+template<std::size_t i,class...T>decltype(auto)
 do_extract(arg_i_begin<i>,std::tuple<T...>const&t)
 {return std::begin(std::get<i>(t));}
-template<int i,class...T>decltype(auto)
+template<std::size_t i,class...T>decltype(auto)
 do_extract(arg_i_end<i>,std::tuple<T...>const&t)
 {return std::end(std::get<i>(t));}
-template<int i,class...T>decltype(auto)
+template<std::size_t i,class...T>decltype(auto)
 do_extract(arg_i_ip1_range<i>,std::tuple<T...>const&t)
 {
  typedef std::tuple<T...> L;
@ -294,7 +294,7 @@ struct Lazy_cartesian :
    template<class T,class D> struct Functor<T,D,Predicate_tag> {
            typedef typename Get_functor<Approximate_kernel, T>::type FA;
            typedef typename Get_functor<Exact_kernel, T>::type FE;
-            typedef Filtered_predicate2<FE,FA,C2E,C2A> type;
+            typedef Filtered_predicate2<Lazy_cartesian,FE,FA,C2E,C2A> type;
    };
    template<class T,class D> struct Functor<T,D,Compute_tag> {
            typedef Lazy_construction2<T,Kernel> type;
--- 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 <CGAL/assertions.h>
-#include <boost/type_traits/is_empty.hpp>
+#include <type_traits>

 namespace CGAL {
 namespace internal {
 BOOST_MPL_HAS_XXX_TRAIT_DEF(Do_not_store_kernel)
-template<class T,bool=boost::is_empty<T>::value,bool=has_Do_not_store_kernel<T>::value> struct Do_not_store_kernel {
+template<class T,bool=std::is_empty<T>::value,bool=has_Do_not_store_kernel<T>::value> struct Do_not_store_kernel {
        enum { value=false };
        typedef Tag_false type;
 };
--- 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<int,class...> struct Apply_to_last_then_rest_;
+template<std::size_t,class...> struct Apply_to_last_then_rest_;

-template<int d,class F,class T,class... U>
+template<std::size_t d,class F,class T,class... U>
 struct Apply_to_last_then_rest_<d,F,T,U...> {
        typedef typename Apply_to_last_then_rest_<d-1,F,U...,T>::result_type result_type;
        inline result_type operator()(F&&f,T&&t,U&&...u)const{
--- a/NewKernel_d/include/CGAL/typeset.h
+++ b/NewKernel_d/include/CGAL/typeset.h
@ -41,11 +41,17 @@ namespace CGAL {
    template<class X> using contains = std::false_type;
    template<class X> using add = typeset<X>;
  };
+  struct typeset_all {
+    typedef typeset_all type;
+    template<class X> using contains = std::true_type;
+    template<class X> using add = typeset_all;
+  };

  template<class T1, class T2> struct typeset_union_ :
    typeset_union_<typename T1::template add<typename T2::head>::type, typename T2::tail>
  {};
  template<class T> struct typeset_union_ <T, typeset<> > : T {};
+  template<class T> struct typeset_union_ <T, typeset_all > : typeset_all {};

  template<class T1, class T2>
    struct typeset_intersection_ {
@ -55,8 +61,8 @@ namespace CGAL {
        std::conditional<T2::template contains<H>::value,
        typename U::template add<H>::type, U>::type type;
    };
-  template<class T>
-    struct typeset_intersection_<typeset<>,T> : typeset<> {};
+  template<class T> struct typeset_intersection_<typeset<>, T> : typeset<> {};
+  template<class T> struct typeset_intersection_<typeset_all, T> : T {};

  template<class T1, class T2>
    using typeset_union = typename typeset_union_<T1,T2>::type;
--- 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<a)?b:a;
+  }
+  friend Mpzf const& max BOOST_PREVENT_MACRO_SUBSTITUTION (Mpzf const&a, Mpzf const&b){
+    return (a<b)?b:a;
+  }
  private:
  static Mpzf aors(Mpzf const&a, Mpzf const&b, int bsize){
    Mpzf res=noalloc();
@ -846,6 +896,9 @@ struct Mpzf {
    }
  }

+  friend Mpzf operator+(Mpzf const&x){
+    return x;
+  }
  friend Mpzf operator-(Mpzf const&x){
    Mpzf ret = x;
    ret.size = -ret.size;
--- a/Number_types/test/Number_types/Mpzf_new.cpp
+++ b/Number_types/test/Number_types/Mpzf_new.cpp
@ -0,0 +1,40 @@
+#include <CGAL/config.h>
+
+#ifdef CGAL_USE_GMP
+# include <CGAL/Mpzf.h>
+#endif
+#ifdef CGAL_HAS_MPZF
+
+#include <iostream>
+#include <CGAL/Gmpq.h>
+#include <CGAL/Test/_test_algebraic_structure.h>
+#include <CGAL/Test/_test_real_embeddable.h>
+
+int main() {
+    {
+        typedef CGAL::Mpzf NT;
+        typedef CGAL::Integral_domain_without_division_tag Tag;
+        typedef CGAL::Tag_true Is_exact;
+
+        CGAL::test_algebraic_structure<NT,Tag, Is_exact>();
+        CGAL::test_algebraic_structure<NT,Tag, Is_exact>(NT(4),NT(6),NT(15));
+        CGAL::test_algebraic_structure<NT,Tag, Is_exact>(NT(-4),NT(6),NT(15));
+        CGAL::test_algebraic_structure<NT,Tag, Is_exact>(NT(4),NT(-6),NT(15));
+        CGAL::test_algebraic_structure<NT,Tag, Is_exact>(NT(-4),NT(-6),NT(15));
+        CGAL::test_algebraic_structure<NT,Tag, Is_exact>(NT(4),NT(6),NT(-15));
+        CGAL::test_algebraic_structure<NT,Tag, Is_exact>(NT(-4),NT(6), NT(15));
+        CGAL::test_algebraic_structure<NT,Tag, Is_exact>(NT(4),NT(-6),NT(-15));
+        CGAL::test_algebraic_structure<NT,Tag, Is_exact>(NT(-4),NT(-6),NT(-15));
+
+        CGAL::test_real_embeddable<NT>();
+    }
+
+    return 0;
+}
+
+#else
+int main()
+{
+  return 0;
+}
+#endif //CGAL_HAS_MPZF
--- 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
 // ----------------