diff --git a/Algebraic_foundations/include/CGAL/Test/_test_algebraic_structure.h b/Algebraic_foundations/include/CGAL/Test/_test_algebraic_structure.h
index db74015c6758..5abb39fff255 100644
--- a/Algebraic_foundations/include/CGAL/Test/_test_algebraic_structure.h
+++ b/Algebraic_foundations/include/CGAL/Test/_test_algebraic_structure.h
@@ -338,6 +338,7 @@ void test_algebraic_structure_intern( const CGAL::Field_tag& ) {
     test_algebraic_structure_intern< AS >( CGAL::Integral_domain_tag());
      AS  a(1);
      AS  b(3);
+     AS& b_ref = b; b = b_ref; // to exercise self-copy
      AS  c = a / b;
     (void)c;  // avoid warnings for unused variables
 
diff --git a/Documentation/doc/Documentation/Preliminaries.txt b/Documentation/doc/Documentation/Preliminaries.txt
index 0fa310a54afc..e07d90e9e378 100644
--- a/Documentation/doc/Documentation/Preliminaries.txt
+++ b/Documentation/doc/Documentation/Preliminaries.txt
@@ -33,8 +33,11 @@ are:
 If the macro `CGAL_HAS_THREADS` is not defined, then \cgal assumes it can use
 any thread-unsafe code (such as static variables).  By default, this macro is not
 defined, unless `BOOST_HAS_THREADS` or `_OPENMP` is defined.  It is possible
-to force its definition on the command line, and it is possible to prevent its default
-definition by setting `CGAL_HAS_NO_THREADS` from the command line.
+to force its definition in the compiler options, and it is possible to prevent its
+default definition by defining the macro `CGAL_HAS_NO_THREADS`.
+If you are using CMake, then you can set the CMake option `CGAL_HAS_NO_THREADS` to
+`TRUE`. In addition to defining the preprocessor macro CGAL_HAS_NO_THREADS`, it will
+also avoid CMake to link with the native threads support library on your system.
 
 \section Preliminaries_cc0x C++14 Support
 
diff --git a/Filtered_kernel/include/CGAL/Lazy.h b/Filtered_kernel/include/CGAL/Lazy.h
index abb56940aa00..0b25b2e30f51 100644
--- a/Filtered_kernel/include/CGAL/Lazy.h
+++ b/Filtered_kernel/include/CGAL/Lazy.h
@@ -52,6 +52,9 @@
 #include <type_traits>
 #include <utility>
 #include <vector>
+#include <atomic>
+#include <thread>
+#include <mutex>
 
 namespace CGAL {
 
@@ -68,22 +71,12 @@ class Lazy_exact_nt;
 
 template <typename AT, typename ET, typename E2A>
 inline
-const AT&
+decltype(auto)
 approx(const Lazy<AT,ET,E2A>& l)
 {
   return l.approx();
 }
 
-// Where is this one (non-const) needed ?  Is it ?
-template <typename AT, typename ET, typename E2A>
-inline
-AT&
-approx(Lazy<AT,ET,E2A>& l)
-{
-  return l.approx();
-}
-
-
 template <typename AT, typename ET, typename E2A>
 inline
 const ET&
@@ -230,9 +223,67 @@ struct Depth_base {
 #endif
 };
 
+template<class T, bool=std::is_base_of<Handle, T>::value> struct Lazy_reset_member_1 {
+  void operator()(T& t)const{ t = T(); }
+};
+template<class T> struct Lazy_reset_member_1<T, true> {
+  void operator()(T& t)const{ t.reset(); }
+};
+template<class T>void lazy_reset_member(T&t) {
+  Lazy_reset_member_1<T>()(t);
+}
+template<class...T, std::size_t...i>void lazy_reset_member_tuple(std::tuple<T...>&t, std::index_sequence<i...>) {
+  auto ignore = [](auto&&...){};
+  ignore ( (lazy_reset_member(std::get<i>(t)), 0) ... );
+}
+template<class...T>void lazy_reset_member(std::tuple<T...>&t) {
+  lazy_reset_member_tuple(t, std::make_index_sequence<sizeof...(T)>());
+}
+
+// 0: safe default, AT is behind a pointer that can be atomically changed, and it doesn't disappear during update_exact
+// 1: use plain AT without protection
+// 2: split an interval as 2 atomic_double
+// FIXME: CGAL_USE_SSE2 is clearly not the right condition
+#ifdef CGAL_HAS_THREADS
+template<class AT>struct Lazy_rep_selector { static constexpr int value = 0; };
+# if defined CGAL_USE_SSE2 && !defined __SANITIZE_THREAD__ && !__has_feature(thread_sanitizer)
+template<bool b>struct Lazy_rep_selector<Interval_nt<b>> { static constexpr int value = 1; };
+template<bool b, int N>struct Lazy_rep_selector<std::array<Interval_nt<b>,N>> { static constexpr int value = 1; };
+ // Need some declarations, including Simple_cartesian.h would also be possible.
+ template<class>struct Simple_cartesian;
+ template<class>class Point_2;
+ template<class>class Point_3;
+template<bool b>struct Lazy_rep_selector<CGAL::Point_2<CGAL::Simple_cartesian<CGAL::Interval_nt<b>>>> { static constexpr int value = 1; };
+template<bool b>struct Lazy_rep_selector<CGAL::Point_3<CGAL::Simple_cartesian<CGAL::Interval_nt<b>>>> { static constexpr int value = 1; };
+# else
+template<bool b>struct Lazy_rep_selector<Interval_nt<b>> { static constexpr int value = 2; };
+# endif
+#else
+template<class AT>struct Lazy_rep_selector { static constexpr int value = 1; };
+#endif
+
+template<class AT>
+struct AT_wrap {
+  AT at_;
+  AT_wrap():at_(){}
+  AT_wrap(AT const& a):at_(a){}
+  AT_wrap(AT&& a):at_(std::move(a)){}
+  AT const& at()const{return at_;}
+};
+
+// TODO: avoid initializing AT for nothing
+template<class AT, class ET>
+struct AT_ET_wrap : AT_wrap<AT> {
+  ET et_;
+  AT_ET_wrap():et_(){}
+  AT_ET_wrap(ET const& e):et_(e){}
+  AT_ET_wrap(ET&& e):et_(std::move(e)){}
+  template<class A, class E>AT_ET_wrap(A&&a, E&&e):AT_wrap<AT>(std::forward<A>(a)),et_(std::forward<E>(e)){}
+  ET const& et()const{return et_;}
+};
 
 // Abstract base class for lazy numbers and lazy objects
-template <typename AT_, typename ET, typename E2A>
+template <typename AT_, typename ET, typename E2A, int=Lazy_rep_selector<AT_>::value /* 0 */>
 class Lazy_rep : public Rep, public Depth_base
 {
   Lazy_rep (const Lazy_rep&) = delete; // cannot be copied.
@@ -241,71 +292,281 @@ class Lazy_rep : public Rep, public Depth_base
 public:
 
   typedef AT_ AT;
+  typedef AT_ET_wrap<AT,ET> Indirect;
 
-  mutable AT at;
-  mutable ET *et;
+  AT_wrap<AT> at_orig{};
+  mutable std::atomic<AT_wrap<AT>*> ptr_ { &at_orig };
+  mutable std::once_flag once;
 
-  Lazy_rep ()
-    : at(), et(nullptr){}
+  Lazy_rep () {}
 
-  //move-constructor
-  Lazy_rep (Lazy_rep&& other)
-    : at(std::move(other.at)), et(other.et)
+  template<class A>
+  Lazy_rep (A&& a)
+      : at_orig(std::forward<A>(a)){}
+
+  template<class A>
+  Lazy_rep (int count, A&& a)
+    : Rep(count), at_orig(std::forward<A>(a)){}
+
+  template<class A, class E>
+  Lazy_rep (A&& a, E&& e)
+      : ptr_(new AT_ET_wrap<AT,ET>(std::forward<A>(a), std::forward<E>(e))) {}
+
+  AT const& approx() const
   {
-    other.et = nullptr;
-    this->count = std::move(other.count);
+    return ptr_.load(std::memory_order_consume)->at();
   }
 
-  //move-assignment
-  Lazy_rep& operator= (Lazy_rep&& other)
+  const ET & exact_unsafe() const
   {
-    if(this->et)
-    {
-      delete this->et;
+    CGAL_assertion(!is_lazy());
+    return static_cast<AT_ET_wrap<AT,ET>*>(ptr_.load(std::memory_order_relaxed))->et();
+  }
+
+  const ET & exact() const
+  {
+    // The test is unnecessary, only use it if benchmark says so, or in order to avoid calling Lazy_exact_Ex_Cst::update_exact() (which used to contain an assertion)
+    //if (is_lazy())
+    std::call_once(once, [this](){this->update_exact();});
+    return exact_unsafe(); // call_once already synchronized memory
+  }
+
+  template<class A>
+  void set_at(AT_ET_wrap<AT,ET>* p, A&& a) const {
+    p->at_ = std::forward<A>(a);
+  }
+  void set_at(AT_ET_wrap<AT,ET>* p) const {
+    p->at_ = E2A()(p->et());
+  }
+  void keep_at(AT_ET_wrap<AT,ET>* p) const {
+    p->at_ = at_orig.at(); // do not move!
+  }
+
+  void set_ptr(AT_ET_wrap<AT,ET>* p) const {
+    ptr_.store(p, std::memory_order_release);
+  }
+
+  // I think we should have different code for cases where there is some cleanup to do (say, a sum of 2 Lazy_exact_nt) and for cases where there isn't (a Lazy_exact_nt constructed from a double), but it may require making exact() virtual. Objects can be hidden in a tuple in Lazy_rep_n, so checking if there is something to clean requires some code. It isn't clear if we also need to restrict that to cases where update_exact doesn't touch AT. The special version would be basically: if(et==0){pet=new ET(...);if(!et.exchange(0,pet))delete pet; update at?}
+
+#ifdef CGAL_LAZY_KERNEL_DEBUG
+  void print_at_et(std::ostream& os, int level) const
+  {
+    for(int i = 0; i < level; i++){
+      os << "    ";
+    }
+    os << "Approximation: ";
+    print_at(os, approx());
+    os << std::endl;
+    if(! is_lazy()){
+      for(int i = 0; i < level; i++){
+        os << "    ";
+      }
+      os << "Exact: ";
+      print_at(os, exact_unsafe());
+      os << std::endl;
+#ifdef CGAL_LAZY_KERNEL_DEBUG_SHOW_TYPEID
+      for(int i = 0; i < level; i++){
+        os << "    ";
+      }
+      os << "  (type: " << typeid(exact_unsafe()).name() << ")" << std::endl;
+#endif // CGAL_LAZY_KERNEL_DEBUG_SHOW_TYPEID
     }
-    this->et = other.et;
-    other.et = nullptr;
-    this->at = std::move(other.at);
-    this->count = std::move(other.count);
-    return *this;
   }
 
+  virtual void print_dag(std::ostream& os, int level) const {}
+#endif
+
+  bool is_lazy() const { return ptr_.load(std::memory_order_relaxed) == &at_orig; }
+  virtual void update_exact() const = 0;
+  virtual ~Lazy_rep() {
+#if !defined __SANITIZE_THREAD__ && !__has_feature(thread_sanitizer)
+    auto* p = ptr_.load(std::memory_order_relaxed);
+    if (p != &at_orig) {
+      std::atomic_thread_fence(std::memory_order_acquire);
+      delete static_cast<Indirect*>(p);
+    }
+#else
+    auto* p = ptr_.load(std::memory_order_consume);
+    if (p != &at_orig) delete static_cast<Indirect*>(p);
+#endif
+  }
+};
+
+/* How (un)safe is this? The goal is to minimize the overhead compared to a single-thread version by making the fast path almost identical.
+ * For scalars on x86_64, the interval is aligned, so load/store instructions will not slice any double (although Intel does not explicitly guarantee it). On recent hardware, they should even be atomic, although without an official guarantee, and we don't need 128-bit atomicity anyway. The main danger is the unpredictable optimizations a compiler could apply (volatile would disable most of them, but it doesn't seem great), including replacing load/store with memcpy, where I fear some implementation/hardware combinations might slice sometimes. Making Interval_nt a pair of atomic_double would avoid this problem, but would likely incur a penalty since compilers don't optimize atomics much, and we shouldn't need to store/load all the time (TODO: benchmark).
+ * For aggregate-like types (Simple_cartesian::Point_3), it should be ok for the same reason.
+ * This is definitely NOT safe for a std::vector like a Point_d with Dynamic_dimension_tag, so it should only be enabled on a case by case basis, if at all. Storing a Point_3 piecewise with 6 atomic_double would be doable, but painful, and I didn't benchmark to check the performance. */
+template <typename AT_, typename ET, typename E2A>
+class Lazy_rep<AT_, ET, E2A, 1> : public Rep, public Depth_base
+{
+  Lazy_rep (const Lazy_rep&) = delete; // cannot be copied.
+  Lazy_rep& operator= (const Lazy_rep&) = delete; // cannot be copied.
+
+public:
+
+  typedef AT_ AT;
+  typedef ET Indirect;
+
+  mutable AT at;
+  mutable std::atomic<ET*> ptr_ { nullptr };
+#ifdef CGAL_HAS_THREADS
+  mutable std::once_flag once;
+#endif
+
+  Lazy_rep () {}
+
   template<class A>
   Lazy_rep (A&& a)
-      : at(std::forward<A>(a)), et(nullptr){}
+      : at(std::forward<A>(a)) {}
 
   template<class A>
   Lazy_rep (int count, A&& a)
-    : Rep(count), at(std::forward<A>(a)), et(nullptr){}
+    : Rep(count), at(std::forward<A>(a)){}
 
   template<class A, class E>
   Lazy_rep (A&& a, E&& e)
-      : at(std::forward<A>(a)), et(new ET(std::forward<E>(e))) {}
+      : at(std::forward<A>(a)), ptr_(new ET(std::forward<E>(e))) {}
 
-  const AT& approx() const
+  AT const& approx() const
   {
-      return at;
+    return at;
+  }
+
+  template<class A>
+  void set_at(ET*, A&& a) const {
+    at = std::forward<A>(a);
   }
 
-  AT& approx()
+  void set_at(ET* p) const {
+    set_at(p, E2A()(*p));
+  }
+  void keep_at(ET*) const { }
+
+  const ET & exact_unsafe() const
   {
-      return at;
+    return *ptr_.load(std::memory_order_relaxed);
   }
 
   const ET & exact() const
   {
-    if (et==nullptr)
-      update_exact();
-    return *et;
+#ifdef CGAL_HAS_THREADS
+    // The test is unnecessary, only use it if benchmark says so, or in order to avoid calling Lazy_exact_Ex_Cst::update_exact() (which used to contain an assertion)
+    //if (is_lazy())
+    std::call_once(once, [this](){this->update_exact();});
+#else
+    if (is_lazy())
+      this->update_exact();
+#endif
+    return exact_unsafe(); // call_once already synchronized memory
+  }
+
+  void set_ptr(ET* p) const {
+    ptr_.store(p, std::memory_order_release);
+  }
+
+#ifdef CGAL_LAZY_KERNEL_DEBUG
+  void print_at_et(std::ostream& os, int level) const
+  {
+    for(int i = 0; i < level; i++){
+      os << "    ";
+    }
+    os << "Approximation: ";
+    print_at(os, approx());
+    os << std::endl;
+    if(! is_lazy()){
+      for(int i = 0; i < level; i++){
+        os << "    ";
+      }
+      os << "Exact: ";
+      print_at(os, exact_unsafe());
+      os << std::endl;
+#ifdef CGAL_LAZY_KERNEL_DEBUG_SHOW_TYPEID
+      for(int i = 0; i < level; i++){
+        os << "    ";
+      }
+      os << "  (type: " << typeid(exact_unsafe()).name() << ")" << std::endl;
+#endif // CGAL_LAZY_KERNEL_DEBUG_SHOW_TYPEID
+    }
+  }
+
+  virtual void print_dag(std::ostream& os, int level) const {}
+#endif
+
+  bool is_lazy() const { return ptr_.load(std::memory_order_relaxed) == nullptr; }
+  virtual void update_exact() const = 0;
+  virtual ~Lazy_rep() {
+#if !defined __SANITIZE_THREAD__ && !__has_feature(thread_sanitizer)
+    auto* p = ptr_.load(std::memory_order_relaxed);
+    if (p != nullptr) {
+      std::atomic_thread_fence(std::memory_order_acquire);
+      delete p;
+    }
+#else
+    auto* p = ptr_.load(std::memory_order_consume);
+    if (p != nullptr) delete p;
+#endif
+  }
+};
+
+// do we need to (forward) declare Interval_nt?
+template <bool b, typename ET, typename E2A>
+class Lazy_rep<Interval_nt<b>, ET, E2A, 2> : public Rep, public Depth_base
+{
+  Lazy_rep (const Lazy_rep&) = delete; // cannot be copied.
+  Lazy_rep& operator= (const Lazy_rep&) = delete; // cannot be copied.
+
+public:
+
+  typedef Interval_nt<b> AT;
+  typedef ET Indirect;
+
+  mutable std::atomic<double> x, y; // -inf, +sup
+  mutable std::atomic<ET*> ptr_ { nullptr };
+  mutable std::once_flag once;
+
+  Lazy_rep () {}
+
+  Lazy_rep (AT a)
+      : x(-a.inf()), y(a.sup()) {}
+
+  template<class E>
+  Lazy_rep (AT a, E&& e)
+      : x(-a.inf()), y(a.sup()), ptr_(new ET(std::forward<E>(e))) {}
+
+  AT approx() const
+  {
+    return AT(-x.load(std::memory_order_relaxed), y.load(std::memory_order_relaxed));
+  }
+
+  void set_at(ET*, AT a) const {
+    x.store(-a.inf(), std::memory_order_relaxed);
+    y.store(a.sup(), std::memory_order_relaxed);
+  }
+
+  void set_at(ET* p) const {
+    set_at(p, E2A()(*p));
+  }
+  void keep_at(ET*) const { }
+
+  const ET & exact_unsafe() const
+  {
+    return *ptr_.load(std::memory_order_relaxed);
   }
 
-  ET & exact()
+  const ET & exact() const
   {
-    if (et==nullptr)
-      update_exact();
-    return *et;
+    // The test is unnecessary, only use it if benchmark says so, or in order to avoid calling Lazy_exact_Ex_Cst::update_exact() (which used to contain an assertion)
+    //if (is_lazy())
+    std::call_once(once, [this](){this->update_exact();});
+    return exact_unsafe(); // call_once already synchronized memory
+  }
+
+  void set_ptr(ET* p) const {
+    ptr_.store(p, std::memory_order_release);
   }
 
+  // I think we should have different code for cases where there is some cleanup to do (say, a sum of 2 Lazy_exact_nt) and for cases where there isn't (a Lazy_exact_nt constructed from a double). Objects can be hidden in a tuple in Lazy_rep_n, so checking if there is something to clean requires some code. It isn't clear if we also need to restrict that to cases where update_exact doesn't touch AT. The special version would be basically: if(et==0){pet=new ET(...);if(!et.exchange(0,pet))delete pet; update at?}
+
 #ifdef CGAL_LAZY_KERNEL_DEBUG
   void print_at_et(std::ostream& os, int level) const
   {
@@ -313,20 +574,20 @@ class Lazy_rep : public Rep, public Depth_base
       os << "    ";
     }
     os << "Approximation: ";
-    print_at(os, at);
+    print_at(os, approx());
     os << std::endl;
     if(! is_lazy()){
       for(int i = 0; i < level; i++){
         os << "    ";
       }
       os << "Exact: ";
-      print_at(os, *et);
+      print_at(os, exact_unsafe());
       os << std::endl;
 #ifdef CGAL_LAZY_KERNEL_DEBUG_SHOW_TYPEID
       for(int i = 0; i < level; i++){
         os << "    ";
       }
-      os << "  (type: " << typeid(*et).name() << ")" << std::endl;
+      os << "  (type: " << typeid(exact_unsafe()).name() << ")" << std::endl;
 #endif // CGAL_LAZY_KERNEL_DEBUG_SHOW_TYPEID
     }
   }
@@ -334,16 +595,28 @@ class Lazy_rep : public Rep, public Depth_base
   virtual void print_dag(std::ostream& os, int level) const {}
 #endif
 
-  bool is_lazy() const { return et == nullptr; }
+  bool is_lazy() const { return ptr_.load(std::memory_order_relaxed) == nullptr; }
   virtual void update_exact() const = 0;
-  virtual ~Lazy_rep() { delete et; }
+  virtual ~Lazy_rep() {
+#if !defined __SANITIZE_THREAD__ && !__has_feature(thread_sanitizer)
+    auto* p = ptr_.load(std::memory_order_relaxed);
+    if (p != nullptr) {
+      std::atomic_thread_fence(std::memory_order_acquire);
+      delete p;
+    }
+#else
+    auto* p = ptr_.load(std::memory_order_consume);
+    if (p != nullptr) delete p;
+#endif
+  }
 };
 
 
-template<typename AT, typename ET, typename AC, typename EC, typename E2A, typename...L>
-class Lazy_rep_n :
+template<typename AT, typename ET, typename AC, typename EC, typename E2A, bool noprune, typename...L>
+class Lazy_rep_n final :
   public Lazy_rep< AT, ET, E2A >, private EC
 {
+  typedef Lazy_rep< AT, ET, E2A > Base;
   // Lazy_rep_0 does not inherit from EC or take a parameter AC. It has different constructors.
   static_assert(sizeof...(L)>0, "Use Lazy_rep_0 instead");
   template <class Ei, class Ai, class E2Ai, class Ki> friend class Lazy_kernel_base;
@@ -351,9 +624,11 @@ class Lazy_rep_n :
   const EC& ec() const { return *this; }
   template<std::size_t...I>
   void update_exact_helper(std::index_sequence<I...>) const {
-    this->et = new ET(ec()( CGAL::exact( std::get<I>(l) ) ... ) );
-    this->at = E2A()(*(this->et));
-    l = std::tuple<L...>{};
+    auto* p = new typename Base::Indirect(ec()( CGAL::exact( std::get<I>(l) ) ... ) );
+    this->set_at(p);
+    this->set_ptr(p);
+    if(!noprune || is_currently_single_threaded())
+      lazy_reset_member(l);
   }
   public:
   void update_exact() const {
@@ -381,7 +656,7 @@ class Lazy_rep_n :
   }
   public:
   void print_dag(std::ostream& os, int level) const {
-    print_dag_helper(os, level, std::make_index_sequence<sizeof...L>{});
+    print_dag_helper(os, level, std::make_index_sequence<sizeof...(L)>{});
   }
 #endif
 };
@@ -400,9 +675,11 @@ class Lazy_rep_optional_n :
 
   template<std::size_t...I>
   void update_exact_helper(std::index_sequence<I...>) const {
-    this->et = new ET( * ec()( CGAL::exact( std::get<I>(l) ) ... ) );
-    this->at = E2A()(*(this->et));
-    l = std::tuple<L...>{};
+    typedef Lazy_rep< AT, ET, E2A > Base;
+    auto* p = new typename Base::Indirect( * ec()( CGAL::exact( std::get<I>(l) ) ... ) );
+    this->set_at(p);
+    this->set_ptr(p);
+    lazy_reset_member(l);
   }
   public:
 
@@ -443,7 +720,7 @@ class Lazy_rep_optional_n :
 
   public:
   void print_dag(std::ostream& os, int level) const {
-    print_dag_helper(os, level, std::make_index_sequence<sizeof...L>{});
+    print_dag_helper(os, level, std::make_index_sequence<sizeof...(L)>{});
   }
 #endif
 };
@@ -452,7 +729,7 @@ class Lazy_rep_optional_n :
 // The rep for the leaf node
 
 template <typename AT, typename ET, typename E2A>
-class Lazy_rep_0 : public Lazy_rep<AT, ET, E2A>
+class Lazy_rep_0 final : public Lazy_rep<AT, ET, E2A>
 {
 
   typedef Lazy_rep<AT, ET, E2A> Base;
@@ -461,12 +738,23 @@ class Lazy_rep_0 : public Lazy_rep<AT, ET, E2A>
   void
   update_exact() const
   {
-    this->et = new ET();
+#ifdef CGAL_HAS_THREADS
+    // Unless we add is_lazy before call_once in Lazy_rep. This test is
+    // necessary because this class can be used either for default
+    // construction, or to store a non-lazy exact value, and only the first one
+    // should have a non-empty update_exact.
+    // An alternative would be to add in the constructors taking an ET: std::call_once(this->once, [](){});
+    if(!this->is_lazy()) return;
+#endif
+    auto* p = new typename Base::Indirect();
+    this->set_ptr(p);
   }
 
   Lazy_rep_0()
     : Lazy_rep<AT,ET, E2A>() {}
 
+  // TODO: the case where the exact value is provided at construction should
+  // actually use a different class from the lazy default construction.
   template<class A, class E>
   Lazy_rep_0(A&& a, E&& e)
     : Lazy_rep<AT,ET,E2A>(std::forward<A>(a), std::forward<E>(e)) {}
@@ -509,7 +797,7 @@ struct Approx_converter
   //typedef Converter  Number_type_converter;
 
   template < typename T >
-  const typename T::AT&
+  decltype(auto)
   operator()(const T&t) const
   { return t.approx(); }
 
@@ -534,7 +822,7 @@ struct Exact_converter
   //typedef Converter  Number_type_converter;
 
   template < typename T >
-  const typename T::ET&
+  decltype(auto)
   operator()(const T&t) const
   { return t.exact(); }
 
@@ -556,7 +844,7 @@ struct Exact_converter
 
 
 template <typename AC, typename EC, typename E2A, typename L1>
-class Lazy_rep_with_vector_1
+class Lazy_rep_with_vector_1 final
   : public Lazy_rep<std::vector<Object>, std::vector<Object>, E2A>
   , private EC
 {
@@ -573,22 +861,21 @@ class Lazy_rep_with_vector_1
   void
   update_exact() const
   {
-// TODO : This looks really unfinished...
+    auto* p = new typename Base::Indirect();
+    // TODO : This looks really unfinished...
     std::vector<Object> vec;
-    this->et = new ET();
     //this->et->reserve(this->at.size());
-    ec()(CGAL::exact(l1_), std::back_inserter(*(this->et)));
-    if(this->et==nullptr)
-    E2A()(*(this->et));
-    this->at = E2A()(*(this->et));
+    ec()(CGAL::exact(l1_), std::back_inserter(p->et_));
+    this->set_at(p);
+    this->set_ptr(p);
     // Prune lazy tree
-    l1_ = L1();
+    lazy_reset_member(l1_);
   }
 
   Lazy_rep_with_vector_1(const AC& ac, const EC& /*ec*/, const L1& l1)
     : l1_(l1)
   {
-    ac(CGAL::approx(l1), std::back_inserter(this->at));
+    ac(CGAL::approx(l1), std::back_inserter(this->at_orig.at_));
   }
 
 #ifdef CGAL_LAZY_KERNEL_DEBUG
@@ -596,7 +883,7 @@ class Lazy_rep_with_vector_1
   print_dag(std::ostream& os, int level) const
   {
     this->print_at_et(os, level);
-    os << "A Lazy_rep_with_vector_1 of size " <<  this->at.size() << std::endl;
+    os << "A Lazy_rep_with_vector_1 of size " <<  this->approx().size() << std::endl;
     if(this->is_lazy()){
       CGAL::msg(os, level, "DAG with one child node:");
       CGAL::print_dag(l1_, os, level+1);
@@ -608,7 +895,7 @@ class Lazy_rep_with_vector_1
 
 
 template <typename AC, typename EC, typename E2A, typename L1, typename L2>
-class Lazy_rep_with_vector_2
+class Lazy_rep_with_vector_2 final
   : public Lazy_rep<std::vector<Object>, std::vector<Object>, E2A>
   , private EC
 {
@@ -626,19 +913,20 @@ class Lazy_rep_with_vector_2
   void
   update_exact() const
   {
-    this->et = new ET();
-    this->et->reserve(this->at.size());
-    ec()(CGAL::exact(l1_), CGAL::exact(l2_), std::back_inserter(*(this->et)));
-    this->at = E2A()(*(this->et));
+    auto* p = new typename Base::Indirect();
+    p->et_.reserve(this->at_orig.at().size());
+    ec()(CGAL::exact(l1_), CGAL::exact(l2_), std::back_inserter(p->et_));
+    this->set_at(p);
+    this->set_ptr(p);
     // Prune lazy tree
-    l1_ = L1();
-    l2_ = L2();
+    lazy_reset_member(l1_);
+    lazy_reset_member(l2_);
   }
 
   Lazy_rep_with_vector_2(const AC& ac, const EC& /*ec*/, const L1& l1, const L2& l2)
     : l1_(l1), l2_(l2)
   {
-    ac(CGAL::approx(l1), CGAL::approx(l2), std::back_inserter(this->at));
+    ac(CGAL::approx(l1), CGAL::approx(l2), std::back_inserter(this->at_orig.at_));
   }
 
 #ifdef CGAL_LAZY_KERNEL_DEBUG
@@ -646,7 +934,7 @@ class Lazy_rep_with_vector_2
   print_dag(std::ostream& os, int level) const
   {
     this->print_at_et(os, level);
-    os << "A Lazy_rep_with_vector_2 of size " <<  this->at.size() << std::endl;
+    os << "A Lazy_rep_with_vector_2 of size " <<  this->approx().size() << std::endl;
     if(this->is_lazy()){
       CGAL::msg(os, level, "DAG with two child nodes:");
       CGAL::print_dag(l1_, os, level+1);
@@ -658,7 +946,7 @@ class Lazy_rep_with_vector_2
 
 
 template <typename AC, typename EC, typename E2A, typename L1, typename L2, typename R1>
-class Lazy_rep_2_1
+class Lazy_rep_2_1 final
   : public Lazy_rep<typename R1::AT, typename R1::ET, E2A>
   , private EC
 {
@@ -676,18 +964,19 @@ class Lazy_rep_2_1
   void
   update_exact() const
   {
-    this->et = new ET();
-    ec()(CGAL::exact(l1_), CGAL::exact(l2_), *(this->et));
-    this->at = E2A()(*(this->et));
+    auto* p = new typename Base::Indirect();
+    ec()(CGAL::exact(l1_), CGAL::exact(l2_), p->et_);
+    this->set_at(p);
+    this->set_ptr(p);
     // Prune lazy tree
-    l1_ = L1();
-    l2_ = L2();
+    lazy_reset_member(l1_);
+    lazy_reset_member(l2_);
   }
 
   Lazy_rep_2_1(const AC& ac, const EC& /*ec*/, const L1& l1, const L2& l2)
     : Lazy_rep<AT,ET,E2A>(), l1_(l1), l2_(l2)
   {
-    ac(CGAL::approx(l1), CGAL::approx(l2), this->at);
+    ac(CGAL::approx(l1), CGAL::approx(l2), this->at_orig.at_);
   }
 
 #ifdef CGAL_LAZY_KERNEL_DEBUG
@@ -710,7 +999,7 @@ class Lazy_rep_2_1
 // The following rep class stores two non-const reference parameters of type R1 and R2
 
 template <typename AC, typename EC, typename E2A, typename L1, typename L2, typename R1, typename R2>
-class Lazy_rep_2_2
+class Lazy_rep_2_2 final
   : public Lazy_rep<std::pair<typename R1::AT,typename R2::AT>, std::pair<typename R1::ET, typename R2::ET>, E2A>
   , private EC
 {
@@ -728,18 +1017,19 @@ class Lazy_rep_2_2
   void
   update_exact() const
   {
-    this->et = new ET();
-    ec()(CGAL::exact(l1_), CGAL::exact(l2_), this->et->first, this->et->second );
-    this->at = E2A()(*(this->et));
+    auto* p = new typename Base::Indirect();
+    ec()(CGAL::exact(l1_), CGAL::exact(l2_), p->et_.first, p->et_.second );
+    this->set_at(p);
+    this->set_ptr(p);
     // Prune lazy tree
-    l1_ = L1();
-    l2_ = L2();
+    lazy_reset_member(l1_);
+    lazy_reset_member(l2_);
   }
 
   Lazy_rep_2_2(const AC& ac, const EC& /*ec*/, const L1& l1, const L2& l2)
     : Lazy_rep<AT,ET,E2A>(), l1_(l1), l2_(l2)
   {
-    ac(CGAL::approx(l1), CGAL::approx(l2), this->at.first, this->at.second);
+    ac(CGAL::approx(l1), CGAL::approx(l2), this->at_orig.at_.first, this->at_orig.at_.second);
   }
 
 #ifdef CGAL_LAZY_KERNEL_DEBUG
@@ -830,18 +1120,12 @@ public :
   friend void swap(Lazy& a, Lazy& b) noexcept
   { swap(static_cast<Handle&>(a), static_cast<Handle&>(b)); }
 
-  const AT& approx() const
+  decltype(auto) approx() const
   { return ptr()->approx(); }
 
   const ET& exact() const
   { return ptr()->exact(); }
 
-  AT& approx()
-  { return ptr()->approx(); }
-
-  ET& exact()
-  { return ptr()->exact(); }
-
   unsigned depth() const
   {
     return ptr()->depth();
@@ -927,7 +1211,8 @@ struct Lazy_construction_optional_for_polygonal_envelope
       CGAL_STATIC_THREAD_LOCAL_VARIABLE_0(LazyPointRep, rep);
 
       const typename AK::Point_3 ap = *oap;
-      rep = LazyPointRep(2,ap, ec, l1, l2, l3);
+      // rep = LazyPointRep(2,ap, ec, l1, l2, l3);
+      rep.~LazyPointRep(); new (&rep) LazyPointRep(2, ap, ec, l1, l2, l3);
       typename LK::Point_3 lp(&rep);
       return boost::make_optional(lp);
 
@@ -965,7 +1250,8 @@ struct Lazy_construction_optional_for_polygonal_envelope
 
       CGAL_STATIC_THREAD_LOCAL_VARIABLE_0(LazyPointRep, rep);
       const typename AK::Point_3 ap = *oap;
-      rep = LazyPointRep(2, ap, ec, l1, l2);
+      // rep = LazyPointRep(2, ap, ec, l1, l2);
+      rep.~LazyPointRep(); new (&rep) LazyPointRep(2, ap, ec, l1, l2);
       typename LK::Point_3 lp(&rep);
       return boost::make_optional(lp);
 
@@ -1011,7 +1297,7 @@ struct Lazy_construction_nt {
     CGAL_BRANCH_PROFILER(std::string(" failures/calls to   : ") + std::string(CGAL_PRETTY_FUNCTION), tmp);
     Protect_FPU_rounding<Protection> P;
     try {
-      return new Lazy_rep_n<AT, ET, AC, EC, To_interval<ET>, L... >(ac, ec, l...);
+      return new Lazy_rep_n<AT, ET, AC, EC, To_interval<ET>, false, L... >(ac, ec, l...);
     } catch (Uncertain_conversion_exception&) {
       CGAL_BRANCH_PROFILER_BRANCH(tmp);
       Protect_FPU_rounding<!Protection> P2(CGAL_FE_TONEAREST);
@@ -1300,8 +1586,8 @@ struct Lazy_functor_2_2
       typedef Lazy<std::pair<typename R1::AT, typename R2::AT>, std::pair<typename R1::ET, typename R2::ET>, E2A> Lazy_pair;
       Lazy_pair lv(new Lazy_rep_2_2<AC, EC, E2A, L1, L2, R1, R2>(ac, ec, l1, l2));
       // lv->approx() is a std::pair<R1::AT, R2::AT>;
-      r1 = R1(Handle_1(new Lazy_rep_n<void, void, First<std::pair<typename R1::AT, typename R2::AT> >, First<std::pair<typename R1::ET, typename R2::ET> >, E2A, Lazy_pair>(First<std::pair<typename R1::AT, typename R2::AT> >(), First<std::pair<typename R1::ET, typename R2::ET> >(), lv)));
-      r2 = R2(Handle_2(new Lazy_rep_n<void, void, Second<std::pair<typename R1::AT, typename R2::AT> >, Second<std::pair<typename R1::ET, typename R2::ET> >, E2A, Lazy_pair>(Second<std::pair<typename R1::AT, typename R2::AT> >(), Second<std::pair<typename R1::ET, typename R2::ET> >(), lv)));
+      r1 = R1(Handle_1(new Lazy_rep_n<void, void, First<std::pair<typename R1::AT, typename R2::AT> >, First<std::pair<typename R1::ET, typename R2::ET> >, E2A, false, Lazy_pair>(First<std::pair<typename R1::AT, typename R2::AT> >(), First<std::pair<typename R1::ET, typename R2::ET> >(), lv)));
+      r2 = R2(Handle_2(new Lazy_rep_n<void, void, Second<std::pair<typename R1::AT, typename R2::AT> >, Second<std::pair<typename R1::ET, typename R2::ET> >, E2A, false, Lazy_pair>(Second<std::pair<typename R1::AT, typename R2::AT> >(), Second<std::pair<typename R1::ET, typename R2::ET> >(), lv)));
     } catch (Uncertain_conversion_exception&) {
       CGAL_BRANCH_PROFILER_BRANCH(tmp);
       Protect_FPU_rounding<!Protection> P2(CGAL_FE_TONEAREST);
@@ -1348,7 +1634,7 @@ struct Lazy_intersect_with_iterators
 // FIXME : I'm not sure how this work...
 #define CGAL_Kernel_obj(X) if (object_cast<typename AK::X>(& (lv.approx()[i]))) { \
           *it++ = make_object(typename LK::X(new Lazy_rep_n<typename AK::X, typename EK::X, Ith<typename AK::X>, \
-                                                                      Ith<typename EK::X>, E2A, Lazy_vector> \
+                                                                      Ith<typename EK::X>, E2A, false, Lazy_vector> \
                                                  (Ith<typename AK::X>(i), Ith<typename EK::X>(i), lv))); \
           continue; \
         }
@@ -1419,14 +1705,14 @@ struct Lazy_construction_object
     CGAL_BRANCH_PROFILER(std::string(" failures/calls to   : ") + std::string(CGAL_PRETTY_FUNCTION), tmp);
     Protect_FPU_rounding<Protection> P;
     try {
-      Lazy_object lo(new Lazy_rep_n<result_type, result_type, AC, EC, E2A, L1>(ac, ec, l1));
+      Lazy_object lo(new Lazy_rep_n<result_type, result_type, AC, EC, E2A, false, L1>(ac, ec, l1));
 
       if(lo.approx().is_empty())
         return Object();
 
 #define CGAL_Kernel_obj(X) \
       if (object_cast<typename AK::X>(& (lo.approx()))) { \
-        typedef Lazy_rep_n< typename AK::X, typename EK::X, Object_cast<typename AK::X>, Object_cast<typename EK::X>, E2A, Lazy_object> Lcr; \
+        typedef Lazy_rep_n< typename AK::X, typename EK::X, Object_cast<typename AK::X>, Object_cast<typename EK::X>, E2A, false, Lazy_object> Lcr; \
         Lcr * lcr = new Lcr(Object_cast<typename AK::X>(), Object_cast<typename EK::X>(), lo); \
         return make_object(typename LK::X(lcr)); \
       }
@@ -1452,14 +1738,14 @@ struct Lazy_construction_object
     CGAL_BRANCH_PROFILER(std::string(" failures/calls to   : ") + std::string(CGAL_PRETTY_FUNCTION), tmp);
     Protect_FPU_rounding<Protection> P;
     try {
-      Lazy_object lo(new Lazy_rep_n<result_type, result_type, AC, EC, E2A, L1, L2>(ac, ec, l1, l2));
+      Lazy_object lo(new Lazy_rep_n<result_type, result_type, AC, EC, E2A, false, L1, L2>(ac, ec, l1, l2));
 
       if(lo.approx().is_empty())
         return Object();
 
 #define CGAL_Kernel_obj(X) \
       if (object_cast<typename AK::X>(& (lo.approx()))) { \
-        typedef Lazy_rep_n<typename AK::X, typename EK::X, Object_cast<typename AK::X>, Object_cast<typename EK::X>, E2A, Lazy_object> Lcr; \
+        typedef Lazy_rep_n<typename AK::X, typename EK::X, Object_cast<typename AK::X>, Object_cast<typename EK::X>, E2A, false, Lazy_object> Lcr; \
         Lcr * lcr = new Lcr(Object_cast<typename AK::X>(), Object_cast<typename EK::X>(), lo); \
         return make_object(typename LK::X(lcr)); \
       }
@@ -1476,7 +1762,7 @@ struct Lazy_construction_object
           V.resize(v_ptr->size());                           \
           for (unsigned int i = 0; i < v_ptr->size(); i++) {               \
             V[i] = typename LK::X(new Lazy_rep_n<typename AK::X, typename EK::X, Ith_for_intersection<typename AK::X>, \
-                                                 Ith_for_intersection<typename EK::X>, E2A, Lazy_object> \
+                                                 Ith_for_intersection<typename EK::X>, E2A, false, Lazy_object> \
                                   (Ith_for_intersection<typename AK::X>(i), Ith_for_intersection<typename EK::X>(i), lo)); \
           }                                                           \
           return make_object(V);                                      \
@@ -1506,14 +1792,14 @@ CGAL_Kernel_obj(Point_3)
     CGAL_BRANCH_PROFILER(std::string(" failures/calls to   : ") + std::string(CGAL_PRETTY_FUNCTION), tmp);
     Protect_FPU_rounding<Protection> P;
     try {
-      Lazy_object lo(new Lazy_rep_n<result_type, result_type, AC, EC, E2A, L1, L2, L3>(ac, ec, l1, l2, l3));
+      Lazy_object lo(new Lazy_rep_n<result_type, result_type, AC, EC, E2A, false, L1, L2, L3>(ac, ec, l1, l2, l3));
 
       if(lo.approx().is_empty())
         return Object();
 
 #define CGAL_Kernel_obj(X) \
       if (object_cast<typename AK::X>(& (lo.approx()))) { \
-        typedef Lazy_rep_n<typename AK::X, typename EK::X, Object_cast<typename AK::X>, Object_cast<typename EK::X>, E2A, Lazy_object> Lcr; \
+        typedef Lazy_rep_n<typename AK::X, typename EK::X, Object_cast<typename AK::X>, Object_cast<typename EK::X>, E2A, false, Lazy_object> Lcr; \
         Lcr * lcr = new Lcr(Object_cast<typename AK::X>(), Object_cast<typename EK::X>(), lo); \
         return make_object(typename LK::X(lcr)); \
       }
@@ -1584,7 +1870,7 @@ struct Fill_lazy_variant_visitor_2 : boost::static_visitor<> {
     typedef typename Type_mapper<AKT, AK, EK>::type EKT;
     typedef typename Type_mapper<AKT, AK, LK>::type LKT;
 
-    typedef Lazy_rep_n<AKT, EKT, Variant_cast<AKT>, Variant_cast<EKT>, typename LK::E2A, Origin> Lcr;
+    typedef Lazy_rep_n<AKT, EKT, Variant_cast<AKT>, Variant_cast<EKT>, typename LK::E2A, false, Origin> Lcr;
     Lcr * lcr = new Lcr(Variant_cast<AKT>(), Variant_cast<EKT>(), *o);
 
     *r = LKT(lcr);
@@ -1600,7 +1886,7 @@ struct Fill_lazy_variant_visitor_2 : boost::static_visitor<> {
     V.resize(t.size());
     for (unsigned int i = 0; i < t.size(); i++) {
       V[i] = LKT(new Lazy_rep_n<AKT, EKT, Ith_for_intersection<AKT>,
-                 Ith_for_intersection<EKT>, typename LK::E2A, Origin>
+                 Ith_for_intersection<EKT>, typename LK::E2A, false, Origin>
                  (Ith_for_intersection<AKT>(i), Ith_for_intersection<EKT>(i), *o));
     }
 
@@ -1683,7 +1969,7 @@ struct Lazy_construction_variant {
     Protect_FPU_rounding<Protection> P;
 
     try {
-      Lazy<AT, ET, E2A> lazy(new Lazy_rep_n<AT, ET, AC, EC, E2A, L1, L2>(AC(), EC(), l1, l2));
+      Lazy<AT, ET, E2A> lazy(new Lazy_rep_n<AT, ET, AC, EC, E2A, false, L1, L2>(AC(), EC(), l1, l2));
 
       // the approximate result requires the trait with types from the AK
       AT approx_v = lazy.approx();
@@ -1736,7 +2022,7 @@ struct Lazy_construction_variant {
     Protect_FPU_rounding<Protection> P;
 
     try {
-      Lazy<AT, ET, E2A> lazy(new Lazy_rep_n<AT, ET, AC, EC, E2A, L1, L2, L3>(AC(), EC(), l1, l2, l3));
+      Lazy<AT, ET, E2A> lazy(new Lazy_rep_n<AT, ET, AC, EC, E2A, false, L1, L2, L3>(AC(), EC(), l1, l2, l3));
 
       // the approximate result requires the trait with types from the AK
       AT approx_v = lazy.approx();
@@ -1775,6 +2061,9 @@ template<typename LK, typename AC, typename EC, typename E2A = Default,
          bool has_result_type = internal::has_result_type<AC>::value && internal::has_result_type<EC>::value >
 struct Lazy_construction;
 
+template<class AK, class AC> struct Disable_lazy_pruning { static const bool value = false; };
+template<class AK> struct Disable_lazy_pruning<AK, typename AK::Construct_weighted_point_2> { static const bool value = true; };
+template<class AK> struct Disable_lazy_pruning<AK, typename AK::Construct_weighted_point_3> { static const bool value = true; };
 
 // we have a result type, low effort
 template<typename LK, typename AC, typename EC, typename E2A_>
@@ -1792,6 +2081,8 @@ struct Lazy_construction<LK, AC, EC, E2A_, true> {
 
   typedef typename Type_mapper<AT, AK, LK>::type result_type;
 
+  static const bool noprune = Disable_lazy_pruning<AK, AC>::value;
+
   CGAL_NO_UNIQUE_ADDRESS AC ac;
   CGAL_NO_UNIQUE_ADDRESS EC ec;
 
@@ -1803,7 +2094,7 @@ struct Lazy_construction<LK, AC, EC, E2A_, true> {
     CGAL_BRANCH_PROFILER(std::string(" failures/calls to   : ") + std::string(CGAL_PRETTY_FUNCTION), tmp); \
     Protect_FPU_rounding<Protection> P;                                 \
     try {                                                               \
-      return result_type( Handle(new Lazy_rep_n<AT, ET, AC, EC, E2A, BOOST_PP_ENUM_PARAMS(n, L)>(ac, ec, BOOST_PP_ENUM_PARAMS(n, l)))); \
+      return result_type( Handle(new Lazy_rep_n<AT, ET, AC, EC, E2A, noprune, BOOST_PP_ENUM_PARAMS(n, L)>(ac, ec, BOOST_PP_ENUM_PARAMS(n, l)))); \
     } catch (Uncertain_conversion_exception&) {                          \
       CGAL_BRANCH_PROFILER_BRANCH(tmp);                                 \
       Protect_FPU_rounding<!Protection> P2(CGAL_FE_TONEAREST);          \
@@ -1819,7 +2110,7 @@ struct Lazy_construction<LK, AC, EC, E2A_, true> {
   operator()() const
   {
     typedef Lazy<AT, ET, E2A> Handle;
-    return result_type( Handle(new Lazy_rep_0<AT,ET,E2A>()) );
+    return result_type( Handle() );
   }
 
 #undef CGAL_CONSTRUCTION_OPERATOR
@@ -1842,6 +2133,8 @@ struct Lazy_construction<LK, AC, EC, E2A_, false>
     // you are on your own
   };
 
+  static const bool noprune = Disable_lazy_pruning<AK, AC>::value;
+
   CGAL_NO_UNIQUE_ADDRESS AC ac;
   CGAL_NO_UNIQUE_ADDRESS EC ec;
 
@@ -1861,7 +2154,7 @@ struct Lazy_construction<LK, AC, EC, E2A_, false>
     CGAL_BRANCH_PROFILER(std::string(" failures/calls to   : ") + std::string(CGAL_PRETTY_FUNCTION), tmp); \
     Protect_FPU_rounding<Protection> P;                                   \
     try {                                                                 \
-      return result_type( Handle(new Lazy_rep_n<AT, ET, AC, EC, E2A, BOOST_PP_ENUM_PARAMS(n, L)>(ac, ec, BOOST_PP_ENUM_PARAMS(n, l)))); \
+      return result_type( Handle(new Lazy_rep_n<AT, ET, AC, EC, E2A, noprune, BOOST_PP_ENUM_PARAMS(n, L)>(ac, ec, BOOST_PP_ENUM_PARAMS(n, l)))); \
     } catch (Uncertain_conversion_exception&) {                          \
       CGAL_BRANCH_PROFILER_BRANCH(tmp);                                 \
       Protect_FPU_rounding<!Protection> P2(CGAL_FE_TONEAREST);          \
@@ -1881,7 +2174,7 @@ struct Lazy_construction<LK, AC, EC, E2A_, false>
     typedef Lazy<AT, ET, E2A> Handle;
     typedef typename Type_mapper<AT, AK, LK>::type result_type;
 
-    return result_type( Handle(new Lazy_rep_0<AT,ET,E2A>()) );
+    return result_type( Handle() );
   }
 };
 
diff --git a/Filtered_kernel/include/CGAL/Lazy_kernel.h b/Filtered_kernel/include/CGAL/Lazy_kernel.h
index b8e5c3562ca1..71629b11826c 100644
--- a/Filtered_kernel/include/CGAL/Lazy_kernel.h
+++ b/Filtered_kernel/include/CGAL/Lazy_kernel.h
@@ -331,6 +331,7 @@ class Lazy_kernel_base
                          typename Approximate_kernel::Construct_weighted_point_2,
                          typename Exact_kernel::Construct_weighted_point_2,
                          E2A_,
+                         true,
                          Return_base_tag,
                          Point_2,
                          FT
@@ -341,6 +342,7 @@ class Lazy_kernel_base
                          typename Approximate_kernel::Construct_weighted_point_2,
                          typename Exact_kernel::Construct_weighted_point_2,
                          E2A_,
+                         true,
                          Return_base_tag,
                          Point_2,
                          int
@@ -352,6 +354,8 @@ class Lazy_kernel_base
       if(tn == typeid(LR).name()){
         LR * lr = static_cast<LR*>(p.ptr());
         if(lr->is_lazy()){
+          // Another thread could reset lr->l before this line, so we disable reset for Construct_weighted_point_2 in MT-mode.
+          // We could also always disable reset for Construct_weighted_point_2 and return lr->l here even if update_exact has run.
           return std::get<2>(lr->l);
         }
       }else{
@@ -383,6 +387,7 @@ class Lazy_kernel_base
                          typename Approximate_kernel::Construct_weighted_point_3,
                          typename Exact_kernel::Construct_weighted_point_3,
                          E2A_,
+                         true,
                          Return_base_tag,
                          Point_3,
                          FT
@@ -393,6 +398,7 @@ class Lazy_kernel_base
                          typename Approximate_kernel::Construct_weighted_point_3,
                          typename Exact_kernel::Construct_weighted_point_3,
                          E2A_,
+                         true,
                          Return_base_tag,
                          Point_3,
                          int
@@ -442,6 +448,7 @@ class Lazy_kernel_base
                          typename Approximate_kernel::Construct_weighted_point_2,
                          typename Exact_kernel::Construct_weighted_point_2,
                          E2A_,
+                         true,
                          Return_base_tag,
                          Point_2,
                          FT
@@ -452,6 +459,7 @@ class Lazy_kernel_base
                          typename Approximate_kernel::Construct_weighted_point_2,
                          typename Exact_kernel::Construct_weighted_point_2,
                          E2A_,
+                         true,
                          Return_base_tag,
                          Point_2,
                          int
@@ -501,6 +509,7 @@ class Lazy_kernel_base
                          typename Approximate_kernel::Construct_weighted_point_3,
                          typename Exact_kernel::Construct_weighted_point_3,
                          E2A_,
+                         true,
                          Return_base_tag,
                          Point_3,
                          FT
@@ -511,6 +520,7 @@ class Lazy_kernel_base
                          typename Approximate_kernel::Construct_weighted_point_3,
                          typename Exact_kernel::Construct_weighted_point_3,
                          E2A_,
+                         true,
                          Return_base_tag,
                          Point_3,
                          int
diff --git a/Installation/cmake/modules/CGAL_SetupCGALDependencies.cmake b/Installation/cmake/modules/CGAL_SetupCGALDependencies.cmake
index 0745228def8d..cb16e161fc7a 100644
--- a/Installation/cmake/modules/CGAL_SetupCGALDependencies.cmake
+++ b/Installation/cmake/modules/CGAL_SetupCGALDependencies.cmake
@@ -102,6 +102,17 @@ function(CGAL_setup_CGAL_dependencies target)
   target_include_directories(${target} INTERFACE
     $<INSTALL_INTERFACE:include>)
 
+  # Make CGAL depend on threads-support (for Epeck and Epeck_d)
+  if(CGAL_HAS_NO_THREADS)
+    target_compile_definitions(${target} INTERFACE CGAL_HAS_NO_THREADS)
+  else()
+    if(NOT TARGET Threads::Threads)
+      set(THREADS_PREFER_PTHREAD_FLAG ON)
+      find_package(Threads REQUIRED)
+    endif()
+    target_link_libraries(${target} INTERFACE Threads::Threads)
+  endif()
+
   # Now setup compilation flags
   if(MSVC)
     target_compile_options(${target} INTERFACE
diff --git a/Installation/include/CGAL/config.h b/Installation/include/CGAL/config.h
index ccca6b86cb63..02399a888557 100644
--- a/Installation/include/CGAL/config.h
+++ b/Installation/include/CGAL/config.h
@@ -590,6 +590,16 @@ using std::max;
 #define CGAL_CAN_USE_CXX11_ATOMIC
 #endif
 
+#ifndef CGAL_HAS_THREADS
+  namespace CGAL { inline bool is_currently_single_threaded(){ return true; } }
+#elif __has_include(<sys/single_threaded.h>)
+#  include <sys/single_threaded.h>
+  namespace CGAL { inline bool is_currently_single_threaded(){ return ::__libc_single_threaded; } }
+#else
+  /* This is the conservative default */
+  namespace CGAL { inline bool is_currently_single_threaded(){ return false; } }
+#endif
+
 // Support for LEDA with threads
 //   Not that, if CGAL_HAS_THREADS is defined, and you want to use LEDA,
 //   you must link with a version of LEDA libraries that support threads.
diff --git a/Kernel_23/include/CGAL/Circle_2.h b/Kernel_23/include/CGAL/Circle_2.h
index 8abf760e8d83..d4f12d902d2f 100644
--- a/Kernel_23/include/CGAL/Circle_2.h
+++ b/Kernel_23/include/CGAL/Circle_2.h
@@ -62,6 +62,9 @@ class Circle_2 : public R_::Kernel_base::Circle_2
   Circle_2(const RCircle_2& t)
     : RCircle_2(t) {}
 
+  Circle_2(RCircle_2&& t)
+    : RCircle_2(std::move(t)) {}
+
   Circle_2(const Point_2 &center, const FT &squared_radius,
            const Orientation &orientation)
     : RCircle_2(typename R::Construct_circle_2()(Return_base_tag(), center, squared_radius, orientation)) {}
diff --git a/Kernel_23/include/CGAL/Circle_3.h b/Kernel_23/include/CGAL/Circle_3.h
index 0b42941a0644..fc8a37927504 100644
--- a/Kernel_23/include/CGAL/Circle_3.h
+++ b/Kernel_23/include/CGAL/Circle_3.h
@@ -93,6 +93,9 @@ template <class R_>
   Circle_3(const Rep& r)
     : Rep(r) {}
 
+  Circle_3(Rep&& r)
+    : Rep(std::move(r)) {}
+
   decltype(auto)
   diametral_sphere() const
   {
diff --git a/Kernel_23/include/CGAL/Direction_2.h b/Kernel_23/include/CGAL/Direction_2.h
index e5b87eaad4d7..ddc28d98015b 100644
--- a/Kernel_23/include/CGAL/Direction_2.h
+++ b/Kernel_23/include/CGAL/Direction_2.h
@@ -65,6 +65,9 @@ class Direction_2 : public R_::Kernel_base::Direction_2
   Direction_2(const RDirection_2& d)
     : RDirection_2(d) {}
 
+  Direction_2(RDirection_2&& d)
+    : RDirection_2(std::move(d)) {}
+
   explicit Direction_2(const Vector_2& v)
     : RDirection_2(typename R::Construct_direction_2()(Return_base_tag(), v)) {}
 
diff --git a/Kernel_23/include/CGAL/Direction_3.h b/Kernel_23/include/CGAL/Direction_3.h
index d7268f796453..eb11081406d1 100644
--- a/Kernel_23/include/CGAL/Direction_3.h
+++ b/Kernel_23/include/CGAL/Direction_3.h
@@ -64,6 +64,9 @@ class Direction_3 : public R_::Kernel_base::Direction_3
   Direction_3(const Rep& d)
     : Rep(d) {}
 
+  Direction_3(Rep&& d)
+    : Rep(std::move(d)) {}
+
   explicit Direction_3(const Vector_3& v)
     : Rep(typename R::Construct_direction_3()(Return_base_tag(), v)) {}
 
diff --git a/Kernel_23/include/CGAL/Line_2.h b/Kernel_23/include/CGAL/Line_2.h
index 9fe90a502310..fa5dab19f10e 100644
--- a/Kernel_23/include/CGAL/Line_2.h
+++ b/Kernel_23/include/CGAL/Line_2.h
@@ -66,6 +66,9 @@ class Line_2 : public R_::Kernel_base::Line_2
   Line_2(const RLine_2& l)  // conversion impl -> interface class
     : RLine_2(l) {}
 
+  Line_2(RLine_2&& l)
+    : RLine_2(std::move(l)) {}
+
   Line_2(const Point_2 &p, const Point_2 &q)
     : RLine_2(typename R::Construct_line_2()(Return_base_tag(), p,q)) {}
 
diff --git a/Kernel_23/include/CGAL/Line_3.h b/Kernel_23/include/CGAL/Line_3.h
index c17c39eb8d40..19c5e34ae9da 100644
--- a/Kernel_23/include/CGAL/Line_3.h
+++ b/Kernel_23/include/CGAL/Line_3.h
@@ -66,6 +66,9 @@ class Line_3 : public R_::Kernel_base::Line_3
   Line_3(const Rep& l)
       : Rep(l) {}
 
+  Line_3(Rep&& l)
+      : Rep(std::move(l)) {}
+
   Line_3(const Point_3 & p, const Point_3 & q)
       : Rep(typename R::Construct_line_3()(Return_base_tag(), p, q)) {}
 
diff --git a/Kernel_23/include/CGAL/Plane_3.h b/Kernel_23/include/CGAL/Plane_3.h
index c13d57d3a056..df4d378ac4c2 100644
--- a/Kernel_23/include/CGAL/Plane_3.h
+++ b/Kernel_23/include/CGAL/Plane_3.h
@@ -66,6 +66,9 @@ class Plane_3 : public R_::Kernel_base::Plane_3
   Plane_3(const Rep& p)
     : Rep(p) {}
 
+  Plane_3(Rep&& p)
+    : Rep(std::move(p)) {}
+
   Plane_3(const Point_3& p, const Point_3& q, const Point_3& r)
     : Rep(typename R::Construct_plane_3()(Return_base_tag(), p, q, r)) {}
 
diff --git a/Kernel_23/include/CGAL/Point_2.h b/Kernel_23/include/CGAL/Point_2.h
index ac7f41cd6379..cc390178219d 100644
--- a/Kernel_23/include/CGAL/Point_2.h
+++ b/Kernel_23/include/CGAL/Point_2.h
@@ -69,6 +69,10 @@ class Point_2 : public R_::Kernel_base::Point_2
     : RPoint_2(p)
   {}
 
+  Point_2(RPoint_2&& p)
+    : RPoint_2(std::move(p))
+  {}
+
   explicit
   Point_2(const Weighted_point_2& wp)
     : Rep(wp.point())
diff --git a/Kernel_23/include/CGAL/Point_3.h b/Kernel_23/include/CGAL/Point_3.h
index 7d37a288313d..9420d995a0de 100644
--- a/Kernel_23/include/CGAL/Point_3.h
+++ b/Kernel_23/include/CGAL/Point_3.h
@@ -67,6 +67,9 @@ class Point_3 : public R_::Kernel_base::Point_3
   Point_3(const Rep& p)
       : Rep(p) {}
 
+  Point_3(Rep&& p)
+      : Rep(std::move(p)) {}
+
   explicit
   Point_3(const Weighted_point_3& wp)
     : Rep(wp.point())
diff --git a/Kernel_23/include/CGAL/Ray_2.h b/Kernel_23/include/CGAL/Ray_2.h
index 161b11fc01ef..22725a5014f8 100644
--- a/Kernel_23/include/CGAL/Ray_2.h
+++ b/Kernel_23/include/CGAL/Ray_2.h
@@ -67,6 +67,9 @@ class Ray_2 : public R_::Kernel_base::Ray_2
   Ray_2(const RRay_2& r)
     : RRay_2(r) {}
 
+  Ray_2(RRay_2&& r)
+    : RRay_2(std::move(r)) {}
+
   Ray_2(const Point_2 &sp, const Point_2 &secondp)
     : RRay_2(typename R::Construct_ray_2()(Return_base_tag(), sp, secondp)) {}
 
diff --git a/Kernel_23/include/CGAL/Ray_3.h b/Kernel_23/include/CGAL/Ray_3.h
index fe0c2e056e5c..2fc9fa40b69a 100644
--- a/Kernel_23/include/CGAL/Ray_3.h
+++ b/Kernel_23/include/CGAL/Ray_3.h
@@ -63,6 +63,9 @@ class Ray_3 : public R_::Kernel_base::Ray_3
   Ray_3(const Rep& r)
     : Rep(r) {}
 
+  Ray_3(Rep&& r)
+    : Rep(std::move(r)) {}
+
   Ray_3(const Point_3& sp, const Point_3& secondp)
     : Rep(typename R::Construct_ray_3()(Return_base_tag(), sp, secondp)) {}
 
diff --git a/Kernel_23/include/CGAL/Segment_2.h b/Kernel_23/include/CGAL/Segment_2.h
index 72773d62789f..a40f769ecf6a 100644
--- a/Kernel_23/include/CGAL/Segment_2.h
+++ b/Kernel_23/include/CGAL/Segment_2.h
@@ -66,6 +66,9 @@ class Segment_2 : public R_::Kernel_base::Segment_2
   Segment_2(const RSegment_2& s)
     : RSegment_2(s) {}
 
+  Segment_2(RSegment_2&& s)
+    : RSegment_2(std::move(s)) {}
+
   Segment_2(const Point_2 &sp, const Point_2 &ep)
     :  RSegment_2(typename R::Construct_segment_2()(Return_base_tag(), sp,ep)) {}
 
diff --git a/Kernel_23/include/CGAL/Segment_3.h b/Kernel_23/include/CGAL/Segment_3.h
index 1274e3d78689..bacd40156438 100644
--- a/Kernel_23/include/CGAL/Segment_3.h
+++ b/Kernel_23/include/CGAL/Segment_3.h
@@ -65,6 +65,9 @@ class Segment_3 : public R_::Kernel_base::Segment_3
   Segment_3(const Rep& s)
       : Rep(s) {}
 
+  Segment_3(Rep&& s)
+      : Rep(std::move(s)) {}
+
   Segment_3(const Point_3& sp, const Point_3& ep)
     : Rep(typename R::Construct_segment_3()(Return_base_tag(), sp, ep)) {}
 
diff --git a/Kernel_23/include/CGAL/Sphere_3.h b/Kernel_23/include/CGAL/Sphere_3.h
index 5f2987c2cc61..128bf1da9c56 100644
--- a/Kernel_23/include/CGAL/Sphere_3.h
+++ b/Kernel_23/include/CGAL/Sphere_3.h
@@ -62,6 +62,9 @@ class Sphere_3 : public R_::Kernel_base::Sphere_3
   Sphere_3(const Rep& s)
    : Rep(s) {}
 
+  Sphere_3(Rep&& s)
+   : Rep(std::move(s)) {}
+
   Sphere_3(const Point_3_& p, const FT& sq_rad,
            const Orientation& o = COUNTERCLOCKWISE)
    : Rep(typename R::Construct_sphere_3()(Return_base_tag(), p, sq_rad, o)) {}
diff --git a/Kernel_23/include/CGAL/Tetrahedron_3.h b/Kernel_23/include/CGAL/Tetrahedron_3.h
index edbe34c85ce5..a402ca0e9088 100644
--- a/Kernel_23/include/CGAL/Tetrahedron_3.h
+++ b/Kernel_23/include/CGAL/Tetrahedron_3.h
@@ -58,6 +58,9 @@ class Tetrahedron_3 : public R_::Kernel_base::Tetrahedron_3
   Tetrahedron_3(const Rep& t)
       : Rep(t) {}
 
+  Tetrahedron_3(Rep&& t)
+      : Rep(std::move(t)) {}
+
   Tetrahedron_3(const Point_3& p, const Point_3& q,
                 const Point_3& r, const Point_3& s)
     : Rep(typename R::Construct_tetrahedron_3()(Return_base_tag(), p, q, r, s)) {}
diff --git a/Kernel_23/include/CGAL/Vector_2.h b/Kernel_23/include/CGAL/Vector_2.h
index c86a92dfca4d..0db8248c65b3 100644
--- a/Kernel_23/include/CGAL/Vector_2.h
+++ b/Kernel_23/include/CGAL/Vector_2.h
@@ -70,6 +70,9 @@ class Vector_2 : public R_::Kernel_base::Vector_2
   Vector_2(const RVector_2& v)
       : RVector_2(v) {}
 
+  Vector_2(RVector_2&& v)
+      : RVector_2(std::move(v)) {}
+
   Vector_2(const Point_2& a, const Point_2& b)
       : RVector_2(typename R::Construct_vector_2()(Return_base_tag(), a, b)) {}
 
diff --git a/Kernel_23/include/CGAL/Vector_3.h b/Kernel_23/include/CGAL/Vector_3.h
index ae59f074c988..d349cd54bb63 100644
--- a/Kernel_23/include/CGAL/Vector_3.h
+++ b/Kernel_23/include/CGAL/Vector_3.h
@@ -70,6 +70,9 @@ class Vector_3 : public R_::Kernel_base::Vector_3
   Vector_3(const Rep& v)
       : Rep(v) {}
 
+  Vector_3(Rep&& v)
+      : Rep(std::move(v)) {}
+
   Vector_3(const Point_3& a, const Point_3& b)
     : Rep(typename R::Construct_vector_3()(Return_base_tag(), a, b)) {}
 
diff --git a/Kernel_23/include/CGAL/Weighted_point_2.h b/Kernel_23/include/CGAL/Weighted_point_2.h
index 3c090e5a00d3..bf2accb8b372 100644
--- a/Kernel_23/include/CGAL/Weighted_point_2.h
+++ b/Kernel_23/include/CGAL/Weighted_point_2.h
@@ -69,6 +69,9 @@ class Weighted_point_2 : public R_::Kernel_base::Weighted_point_2
   Weighted_point_2(const Rep& p)
       : Rep(p) {}
 
+  Weighted_point_2(Rep&& p)
+      : Rep(std::move(p)) {}
+
   explicit
   Weighted_point_2(const Point_2& p)
     : Rep(typename R::Construct_weighted_point_2()(Return_base_tag(), p, 0))
diff --git a/Kernel_23/include/CGAL/Weighted_point_3.h b/Kernel_23/include/CGAL/Weighted_point_3.h
index 66853dea7d49..094fefa26334 100644
--- a/Kernel_23/include/CGAL/Weighted_point_3.h
+++ b/Kernel_23/include/CGAL/Weighted_point_3.h
@@ -71,6 +71,9 @@ class Weighted_point_3 : public R_::Kernel_base::Weighted_point_3
   Weighted_point_3(const Rep& p)
       : Rep(p) {}
 
+  Weighted_point_3(Rep&& p)
+      : Rep(std::move(p)) {}
+
   explicit
   Weighted_point_3(const Point_3& p)
     : Rep(typename R::Construct_weighted_point_3()(Return_base_tag(), p, 0))
diff --git a/NewKernel_d/include/CGAL/Epeck_d.h b/NewKernel_d/include/CGAL/Epeck_d.h
index 62466490326f..eade53339e83 100644
--- a/NewKernel_d/include/CGAL/Epeck_d.h
+++ b/NewKernel_d/include/CGAL/Epeck_d.h
@@ -22,7 +22,7 @@
 #include <CGAL/internal/Exact_type_selector.h>
 #include <CGAL/NewKernel_d/Types/Weighted_point.h>
 
-// TODO: add static filters somewhere
+// TODO: In Kernel_23, Epeck predicates first see if they can convert their arguments to Epick types exactly, and in that case use the Epick predicates (including static filters, Mpzf, etc).
 namespace CGAL {
 #define CGAL_KA Cartesian_base_d<Interval_nt_advanced,Dim>
 #define CGAL_KE Cartesian_base_d<internal::Exact_field_selector<double>::Type, Dim>
diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Lazy_cartesian.h b/NewKernel_d/include/CGAL/NewKernel_d/Lazy_cartesian.h
index 2045f8c87c38..51ad838b0b80 100644
--- a/NewKernel_d/include/CGAL/NewKernel_d/Lazy_cartesian.h
+++ b/NewKernel_d/include/CGAL/NewKernel_d/Lazy_cartesian.h
@@ -109,6 +109,7 @@ template<typename AT, typename ET, typename AC, typename EC, typename E2A, typen
 class Lazy_rep_XXX :
   public Lazy_rep< AT, ET, E2A >, private EC
 {
+  typedef Lazy_rep< AT, ET, E2A > Base;
   // `default_construct<T>()` is the same as `T{}`. But, this is a
   // workaround to a MSVC-2015 bug (fixed in MSVC-2017): its parser
   // seemed confused by `T{}` somewhere below.
@@ -129,9 +130,10 @@ class Lazy_rep_XXX :
   const EC& ec() const { return *this; }
   template<class...T>
   void update_exact_helper(Lazy_internal::typelist<T...>) const {
-    this->et = new ET(ec()( CGAL::exact( Lazy_internal::do_extract(T{},l) ) ... ) );
-    this->at = E2A()(*(this->et));
-    l = LL(); // There should be a nicer way to clear. Destruction for instance. With this->et as a witness of whether l has already been destructed.
+    auto* p = new typename Base::Indirect(ec()( CGAL::exact( Lazy_internal::do_extract(T{},l) ) ... ) );
+    this->set_at(p);
+    this->set_ptr(p);
+    lazy_reset_member(l);
   }
   public:
   void update_exact() const {
diff --git a/Number_types/include/CGAL/Lazy_exact_nt.h b/Number_types/include/CGAL/Lazy_exact_nt.h
index 1a511325de9e..3d62ac6a7793 100644
--- a/Number_types/include/CGAL/Lazy_exact_nt.h
+++ b/Number_types/include/CGAL/Lazy_exact_nt.h
@@ -113,6 +113,10 @@ struct Lazy_exact_nt_rep : public Lazy_exact_nt<ET>::Self_rep
   Lazy_exact_nt_rep (const Interval_nt<false> & i)
       : Base(i) {}
 
+  template<class T>
+  Lazy_exact_nt_rep (const Interval_nt<false> & i, T&& e)
+      : Base(i, std::forward<T>(e)) {}
+
 #ifdef CGAL_LAZY_KERNEL_DEBUG
   void
   print_dag(std::ostream& os, int level) const
@@ -123,23 +127,33 @@ struct Lazy_exact_nt_rep : public Lazy_exact_nt<ET>::Self_rep
 };
 
 // int constant
+// Unused. This would make even more sense for a double constant but may not be worth the trouble.
+// Could be recycled as a partial specialization of Lazy_exact_Cst.
 template <typename ET>
-struct Lazy_exact_Int_Cst : public Lazy_exact_nt_rep<ET>
+struct Lazy_exact_Int_Cst final : public Lazy_exact_nt_rep<ET>
 {
   Lazy_exact_Int_Cst (int i)
       : Lazy_exact_nt_rep<ET>(double(i)) {}
 
-  void update_exact() const { this->et = new ET((int)this->approx().inf()); }
+  void update_exact() const {
+    auto* pet = new typename Lazy_exact_nt_rep<ET>::Indirect((int)this->approx().sup());
+    this->keep_at(pet);
+    this->set_ptr(pet);
+  }
 };
 
 // double constant
 template <typename ET, typename X>
-struct Lazy_exact_Cst : public Lazy_exact_nt_rep<ET>
+struct Lazy_exact_Cst final : public Lazy_exact_nt_rep<ET>
 {
   Lazy_exact_Cst (X x)
       : Lazy_exact_nt_rep<ET>(x), cste(x) {}
 
-  void update_exact() const { this->et = new ET(cste); }
+  void update_exact() const {
+    auto* pet = new typename Lazy_exact_nt_rep<ET>::Indirect(cste);
+    this->keep_at(pet);
+    this->set_ptr(pet);
+  }
 
   private:
   X cste;
@@ -147,25 +161,22 @@ struct Lazy_exact_Cst : public Lazy_exact_nt_rep<ET>
 
 // Exact constant
 template <typename ET>
-struct Lazy_exact_Ex_Cst : public Lazy_exact_nt_rep<ET>
+struct Lazy_exact_Ex_Cst final : public Lazy_exact_nt_rep<ET>
 {
-  Lazy_exact_Ex_Cst (const ET & e)
-      : Lazy_exact_nt_rep<ET>(CGAL_NTS to_interval(e))
-  {
-    this->et = new ET(e);
-  }
-  Lazy_exact_Ex_Cst (ET&& e)
-      : Lazy_exact_nt_rep<ET>(CGAL_NTS to_interval(e))
-  {
-    this->et = new ET(std::move(e));
-  }
+  template<class T>
+  Lazy_exact_Ex_Cst (T&& e)
+      : Lazy_exact_nt_rep<ET>(CGAL_NTS to_interval(e), std::forward<T>(e))
+  { }
 
-  void update_exact() const { CGAL_error(); }
+  void update_exact() const {
+    // Currently we do not check is_lazy() before calling call_once, so this is called.
+    // CGAL_error();
+  }
 };
 
 // Construction from a Lazy_exact_nt<ET1> (which keeps the lazyness).
 template <typename ET, typename ET1>
-class Lazy_lazy_exact_Cst : public Lazy_exact_nt_rep<ET>
+class Lazy_lazy_exact_Cst final : public Lazy_exact_nt_rep<ET>
 {
   mutable Lazy_exact_nt<ET1> l;
 
@@ -179,12 +190,13 @@ class Lazy_lazy_exact_Cst : public Lazy_exact_nt_rep<ET>
 
   void update_exact() const
   {
-    this->et = new ET(l.exact());
-    this->at = l.approx();
-    prune_dag();
+    auto* pet = new typename Lazy_exact_nt_rep<ET>::Indirect(l.exact());
+    this->set_at(pet, l.approx());
+    this->set_ptr(pet);
+    this->prune_dag();
   }
 
-  void prune_dag() const { l = Lazy_exact_nt<ET1>(); }
+  void prune_dag() const { l.reset(); }
 };
 
 
@@ -203,7 +215,7 @@ struct Lazy_exact_unary : public Lazy_exact_nt_rep<ET>
     this->set_depth(op1.depth() + 1);
   }
 
-  void prune_dag() const { op1 = Lazy_exact_nt<ET>(); }
+  void prune_dag() const { op1.reset(); }
 
 #ifdef CGAL_LAZY_KERNEL_DEBUG
   void
@@ -234,8 +246,8 @@ struct Lazy_exact_binary : public Lazy_exact_nt_rep<ET>
 
   void prune_dag() const
   {
-    op1 = Lazy_exact_nt<ET1>();
-    op2 = Lazy_exact_nt<ET2>();
+    op1.reset();
+    op2.reset();
   }
 
 #ifdef CGAL_LAZY_KERNEL_DEBUG
@@ -256,10 +268,12 @@ struct Lazy_exact_binary : public Lazy_exact_nt_rep<ET>
 // function objects plus, minus, multiplies, divides...).  But it would require
 // a template template parameter, and GCC 2.95 seems to crash easily with them.
 
+// Instead of having nodes only for simple operations, we should use expression templates to build nodes for arbitrary expressions, a*d-b*c should be a single node, that stores a tuple (a,b,c,d) and knows how to update_exact.
+
 // Macro for unary operations
 #define CGAL_LAZY_UNARY_OP(OP, NAME)                                     \
 template <typename ET>                                                   \
-struct NAME : public Lazy_exact_unary<ET>                                \
+struct NAME final : public Lazy_exact_unary<ET>                          \
 {                                                                        \
   typedef typename Lazy_exact_unary<ET>::AT::Protector P;                \
   NAME (const Lazy_exact_nt<ET> &a)                                      \
@@ -267,11 +281,12 @@ struct NAME : public Lazy_exact_unary<ET>                                \
                                                                          \
   void update_exact() const                                              \
   {                                                                      \
-    this->et = new ET(OP(this->op1.exact()));                            \
-    if (!this->approx().is_point())                                      \
-      this->at = CGAL_NTS to_interval(*(this->et));                      \
+    auto* pet = new typename Lazy_exact_nt_rep<ET>::Indirect(OP(this->op1.exact())); \
+    if (!this->approx().is_point())                                               \
+      this->set_at(pet);                                                 \
+    this->set_ptr(pet);                                                  \
     this->prune_dag();                                                   \
-   }                                                                     \
+  }                                                                      \
 };
 
 CGAL_LAZY_UNARY_OP(opposite,  Lazy_exact_Opp)
@@ -282,19 +297,20 @@ CGAL_LAZY_UNARY_OP(CGAL_NTS sqrt,   Lazy_exact_Sqrt)
 // A macro for +, -, * and /
 #define CGAL_LAZY_BINARY_OP(OP, NAME)                                    \
 template <typename ET, typename ET1 = ET, typename ET2 = ET>             \
-struct NAME : public Lazy_exact_binary<ET, ET1, ET2>                     \
+struct NAME final : public Lazy_exact_binary<ET, ET1, ET2>               \
 {                                                                        \
-  typedef typename Lazy_exact_binary<ET,ET1,ET2>::AT::Protector P;         \
+  typedef typename Lazy_exact_binary<ET,ET1,ET2>::AT::Protector P;       \
   NAME (const Lazy_exact_nt<ET1> &a, const Lazy_exact_nt<ET2> &b)        \
     : Lazy_exact_binary<ET, ET1, ET2>((P(), a.approx() OP b.approx()), a, b) {} \
                                                                          \
   void update_exact() const                                              \
   {                                                                      \
-    this->et = new ET(this->op1.exact() OP this->op2.exact());           \
-    if (!this->approx().is_point())                                      \
-      this->at = CGAL_NTS to_interval(*(this->et));                      \
+    auto* pet = new typename Lazy_exact_nt_rep<ET>::Indirect(this->op1.exact() OP this->op2.exact()); \
+    if (!this->approx().is_point())                                               \
+      this->set_at(pet);                                                 \
+    this->set_ptr(pet);                                                  \
     this->prune_dag();                                                   \
-   }                                                                     \
+  }                                                                      \
 };
 
 CGAL_LAZY_BINARY_OP(+, Lazy_exact_Add)
@@ -304,32 +320,35 @@ CGAL_LAZY_BINARY_OP(/, Lazy_exact_Div)
 
 // Minimum
 template <typename ET>
-struct Lazy_exact_Min : public Lazy_exact_binary<ET>
+struct Lazy_exact_Min final : public Lazy_exact_binary<ET>
 {
   Lazy_exact_Min (const Lazy_exact_nt<ET> &a, const Lazy_exact_nt<ET> &b)
     : Lazy_exact_binary<ET>((CGAL::min)(a.approx(), b.approx()), a, b) {}
 
   void update_exact() const
   {
-    this->et = new ET((CGAL::min)(this->op1.exact(), this->op2.exact()));
+    // Should we test is_point earlier, and construct ET from double in that case? Constructing from double is not free, but if op1 or op2 is not exact yet, we may be able to skip a whole tree of exact constructions.
+    auto* pet = new typename Lazy_exact_nt_rep<ET>::Indirect((CGAL::min)(this->op1.exact(), this->op2.exact()));
     if (!this->approx().is_point())
-      this->at = CGAL_NTS to_interval(*(this->et));
+      this->set_at(pet);
+    this->set_ptr(pet);
     this->prune_dag();
   }
 };
 
 // Maximum
 template <typename ET>
-struct Lazy_exact_Max : public Lazy_exact_binary<ET>
+struct Lazy_exact_Max final : public Lazy_exact_binary<ET>
 {
   Lazy_exact_Max (const Lazy_exact_nt<ET> &a, const Lazy_exact_nt<ET> &b)
     : Lazy_exact_binary<ET>((CGAL::max)(a.approx(), b.approx()), a, b) {}
 
   void update_exact() const
   {
-    this->et = new ET((CGAL::max)(this->op1.exact(), this->op2.exact()));
+    auto* pet = new typename Lazy_exact_nt_rep<ET>::Indirect((CGAL::max)(this->op1.exact(), this->op2.exact()));
     if (!this->approx().is_point())
-      this->at = CGAL_NTS to_interval(*(this->et));
+      this->set_at(pet);
+    this->set_ptr(pet);
     this->prune_dag();
   }
 };
@@ -1297,7 +1316,7 @@ operator>> (std::istream & is, Lazy_exact_nt<ET> & a)
   ET e;
   internal::read_float_or_quotient(is, e);
   if (is)
-    a = e;
+    a = std::move(e);
   return is;
 }
 
diff --git a/Number_types/include/CGAL/Root_of_traits_specializations.h b/Number_types/include/CGAL/Root_of_traits_specializations.h
index 709d6936b644..b07aa7b3478f 100644
--- a/Number_types/include/CGAL/Root_of_traits_specializations.h
+++ b/Number_types/include/CGAL/Root_of_traits_specializations.h
@@ -52,21 +52,26 @@ struct Lazy_exact_ro2
 
     void update_exact() const
     {
+        typedef typename Base::Indirect I;
+        I* pet;
         if (old_rep)
-          this->et = new RO2(make_root_of_2(op1.exact(), op2.exact(),
-                                            op3.exact(), smaller));
+          pet = new I(make_root_of_2(op1.exact(), op2.exact(),
+                                     op3.exact(), smaller));
         else
-          this->et = new RO2(make_root_of_2(op1.exact(), op2.exact(),
-                                            op3.exact()));
+          pet = new I(make_root_of_2(op1.exact(), op2.exact(),
+                                     op3.exact()));
         if (!this->approx().is_point())
-            this->at = to_interval(*(this->et));
+            this->set_at(pet);
+        this->set_ptr(pet);
         this->prune_dag();
 
     }
 
     void prune_dag() const
     {
-        op1 = op2 = op3 = Lazy_exact_nt<ET>();
+        op1.reset();
+        op2.reset();
+        op3.reset();
     }
 };
 
diff --git a/STL_Extension/include/CGAL/Handle.h b/STL_Extension/include/CGAL/Handle.h
index 744d4392cbe3..4fa3a17286b4 100644
--- a/STL_Extension/include/CGAL/Handle.h
+++ b/STL_Extension/include/CGAL/Handle.h
@@ -18,6 +18,7 @@
 #define CGAL_HANDLE_H
 
 #include <cstddef>
+#include <atomic>
 #include <CGAL/Handle_for.h>
 #include <CGAL/assertions.h>
 
@@ -25,15 +26,14 @@ namespace CGAL {
 
 class Rep
 {
-  friend class Handle;
-protected:
-  Rep() { count = 1; }
-  Rep(int count)
-    : count(count)
-  {}
-  virtual ~Rep() {}
-
-  int count;
+    friend class Handle;
+  protected:
+    Rep(int count = 1)
+      : count(count)
+    {}
+    virtual ~Rep() {}
+
+    std::atomic_int count;
 };
 
 
@@ -50,29 +50,20 @@ class Handle
     {
       CGAL_precondition( x.PTR != static_cast<Rep*>(0) );
       PTR = x.PTR;
-      CGAL_assume (PTR->count > 0);
-      PTR->count++;
+      //CGAL_assume (PTR->count > 0);
+      incref();
     }
 
-    Handle(Handle&& x) noexcept
-      : PTR(x.PTR)
-    {
-      x.PTR = static_cast<Rep*>(0);
-    }
+    Handle(Handle&& x) noexcept : PTR(x.PTR) { x.PTR = 0; }
 
-    ~Handle()
-    {
-        if ( PTR && (--PTR->count == 0))
-            delete PTR;
-    }
+    ~Handle() { reset(); }
 
     Handle&
     operator=(const Handle& x) noexcept(!CGAL_PRECONDITIONS_ENABLED)
     {
       CGAL_precondition( x.PTR != static_cast<Rep*>(0) );
-      x.PTR->count++;
-      if ( PTR && (--PTR->count == 0))
-          delete PTR;
+      x.incref();
+      if(PTR) decref(); // not reset() in case this==&x
       PTR = x.PTR;
       return *this;
     }
@@ -80,24 +71,56 @@ class Handle
     Handle&
     operator=(Handle&& x) noexcept
     {
-      swap(*this,x);
+      swap(*this, x);
       return *this;
     }
 
     friend void swap(Handle& a, Handle& b) noexcept { std::swap(a.PTR, b.PTR); }
 
+  private:
+    void incref() const noexcept
+    {
+      if (is_currently_single_threaded()) {
+        PTR->count.store(PTR->count.load(std::memory_order_relaxed) + 1, std::memory_order_relaxed);
+      } else {
+        PTR->count.fetch_add(1, std::memory_order_relaxed);
+      }
+    }
+
+    void decref()
+    {
+      if (is_currently_single_threaded()) {
+        auto c = PTR->count.load(std::memory_order_relaxed);
+        if (c == 1)
+          delete PTR;
+        else
+          PTR->count.store(c - 1, std::memory_order_relaxed);
+      } else {
+      // TSAN does not support fences :-(
+#if !defined __SANITIZE_THREAD__ && !__has_feature(thread_sanitizer)
+        if (PTR->count.load(std::memory_order_relaxed) == 1
+            || PTR->count.fetch_sub(1, std::memory_order_release) == 1) {
+          std::atomic_thread_fence(std::memory_order_acquire);
+#else
+        if (PTR->count.fetch_sub(1, std::memory_order_acq_rel) == 1) {
+#endif
+          delete PTR;
+        }
+      }
+    }
+
+  public:
     void reset()
     {
       if (PTR)
       {
-        if (--PTR->count==0)
-          delete PTR;
+        decref();
         PTR=0;
       }
     }
 
     int
-    refs()  const noexcept { return PTR->count; }
+    refs()  const noexcept { return PTR->count.load(std::memory_order_relaxed); }
 
     Id_type id() const noexcept { return PTR - static_cast<Rep*>(0); }
 
diff --git a/STL_Extension/include/CGAL/Handle_for.h b/STL_Extension/include/CGAL/Handle_for.h
index 71ac6739e363..2a93662e69b8 100644
--- a/STL_Extension/include/CGAL/Handle_for.h
+++ b/STL_Extension/include/CGAL/Handle_for.h
@@ -26,6 +26,7 @@
 #include <CGAL/memory.h>
 #include <algorithm>
 #include <cstddef>
+#include <atomic>
 
 #if defined(BOOST_MSVC)
 #  pragma warning(push)
@@ -39,7 +40,9 @@ class Handle_for
     // Wrapper that adds the reference counter.
     struct RefCounted {
         T t;
-        unsigned int count;
+        std::atomic_uint count;
+        template <class... U>
+        RefCounted(U&&...u ) : t(std::forward<U>(u)...), count(1) {}
     };
 
 
@@ -60,25 +63,19 @@ class Handle_for
     Handle_for()
     {
         pointer p = allocator.allocate(1);
-        new (&(p->t)) element_type(); // we get the warning here
-        p->count = 1;
-        ptr_ = p;
+        ptr_ = new (p) RefCounted();
     }
 
     Handle_for(const element_type& t)
     {
         pointer p = allocator.allocate(1);
-        new (&(p->t)) element_type(t);
-        p->count = 1;
-        ptr_ = p;
+        ptr_ = new (p) RefCounted(t);
     }
 
     Handle_for(element_type && t)
     {
         pointer p = allocator.allocate(1);
-        new (&(p->t)) element_type(std::move(t));
-        p->count = 1;
-        ptr_ = p;
+        ptr_ = new (p) RefCounted(std::move(t));
     }
 
 /* I comment this one for now, since it's preventing the automatic conversions
@@ -87,9 +84,7 @@ class Handle_for
     Handle_for(const T1& t1)
     {
         pointer p = allocator.allocate(1);
-        new (&(p->t)) T(t1);
-        p->count = 1;
-        ptr_ = p;
+        ptr_ = new (p) RefCounted(t1);
     }
 */
 
@@ -97,16 +92,17 @@ class Handle_for
     Handle_for(T1 && t1, T2 && t2, Args && ... args)
     {
         pointer p = allocator.allocate(1);
-        new (&(p->t)) element_type(std::forward<T1>(t1), std::forward<T2>(t2), std::forward<Args>(args)...);
-        p->count = 1;
-        ptr_ = p;
+        ptr_ = new (p) RefCounted(std::forward<T1>(t1), std::forward<T2>(t2), std::forward<Args>(args)...);
     }
 
     Handle_for(const Handle_for& h) noexcept(!CGAL_ASSERTIONS_ENABLED)
       : ptr_(h.ptr_)
     {
-        CGAL_assume (ptr_->count > 0);
-        ++(ptr_->count);
+        // CGAL_assume (ptr_->count > 0);
+        if (is_currently_single_threaded())
+          ptr_->count.store(ptr_->count.load(std::memory_order_relaxed) + 1, std::memory_order_relaxed);
+        else
+          ptr_->count.fetch_add(1, std::memory_order_relaxed);
     }
 
     Handle_for&
@@ -151,12 +147,27 @@ class Handle_for
 
     ~Handle_for()
     {
-      try {
-        if (--(ptr_->count) == 0) {
+      if (is_currently_single_threaded()) {
+        auto c = ptr_->count.load(std::memory_order_relaxed);
+        if (c == 1) {
+          Allocator_traits::destroy(allocator, ptr_);
+          allocator.deallocate(ptr_, 1);
+        } else {
+          ptr_->count.store(c - 1, std::memory_order_relaxed);
+        }
+      } else {
+      // TSAN does not support fences :-(
+#if !defined __SANITIZE_THREAD__ && !__has_feature(thread_sanitizer)
+        if (ptr_->count.load(std::memory_order_relaxed) == 1
+            || ptr_->count.fetch_sub(1, std::memory_order_release) == 1) {
+          std::atomic_thread_fence(std::memory_order_acquire);
+#else
+        if (ptr_->count.fetch_sub(1, std::memory_order_acq_rel) == 1) {
+#endif
           Allocator_traits::destroy(allocator, ptr_);
-          allocator.deallocate( ptr_, 1);
+          allocator.deallocate(ptr_, 1);
         }
-      } catch(...) {}
+      }
     }
 
     void
@@ -192,7 +203,7 @@ class Handle_for
     bool
     is_shared() const noexcept
     {
-        return ptr_->count > 1;
+        return ptr_->count.load(std::memory_order_relaxed) > 1;
     }
 
     bool
@@ -204,7 +215,7 @@ class Handle_for
     long
     use_count() const noexcept
     {
-        return ptr_->count;
+        return ptr_->count.load(std::memory_order_relaxed);
     }
 
     void