diff --git a/NewKernel_d/include/CGAL/Kernel_d/Kernel_d_interface.h b/NewKernel_d/include/CGAL/Kernel_d/Kernel_d_interface.h
index 9fed52a0c2c..55989e2f774 100644
--- a/NewKernel_d/include/CGAL/Kernel_d/Kernel_d_interface.h
+++ b/NewKernel_d/include/CGAL/Kernel_d/Kernel_d_interface.h
@@ -109,7 +109,7 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
 	typedef typename Get_functor<Base, Affine_rank_tag>::type Affine_rank_d;
 	typedef typename Get_functor<Base, Affinely_independent_tag>::type Affinely_independent_d;
 	typedef typename Get_functor<Base, Contained_in_linear_hull_tag>::type Contained_in_linear_hull_d;
-	//typedef typename Get_functor<Base, Contained_in_simplex_tag>::type Contained_in_simplex_d;
+	typedef typename Get_functor<Base, Contained_in_simplex_tag>::type Contained_in_simplex_d;
 	typedef typename Get_functor<Base, Has_on_positive_side_tag>::type Has_on_positive_side_d;
 	typedef typename Get_functor<Base, Linear_rank_tag>::type Linear_rank_d;
 	typedef typename Get_functor<Base, Linearly_independent_tag>::type Linearly_independent_d;
@@ -137,7 +137,7 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
 	Side_of_bounded_sphere_d side_of_bounded_sphere_d_object()const{ return Side_of_bounded_sphere_d(*this); }
 	Contained_in_affine_hull_d contained_in_affine_hull_d_object()const{ return Contained_in_affine_hull_d(*this); }
 	Contained_in_linear_hull_d contained_in_linear_hull_d_object()const{ return Contained_in_linear_hull_d(*this); }
-	//Contained_in_simplex_d contained_in_simplex_d_object()const{ return Contained_in_simplex_d(*this); }
+	Contained_in_simplex_d contained_in_simplex_d_object()const{ return Contained_in_simplex_d(*this); }
 	Construct_flat_orientation_d construct_flat_orientation_d_object()const{ return Construct_flat_orientation_d(*this); }
 	In_flat_orientation_d in_flat_orientation_d_object()const{ return In_flat_orientation_d(*this); }
 	In_flat_side_of_oriented_sphere_d in_flat_side_of_oriented_sphere_d_object()const{ return In_flat_side_of_oriented_sphere_d(*this); }
diff --git a/NewKernel_d/include/CGAL/Kernel_d/function_objects_cartesian.h b/NewKernel_d/include/CGAL/Kernel_d/function_objects_cartesian.h
index 4e7eafaf5de..ee258d18e68 100644
--- a/NewKernel_d/include/CGAL/Kernel_d/function_objects_cartesian.h
+++ b/NewKernel_d/include/CGAL/Kernel_d/function_objects_cartesian.h
@@ -345,6 +345,54 @@ template<class R_> struct Affinely_independent : private Store_kernel<R_> {
 
 CGAL_KD_DEFAULT_FUNCTOR(Affinely_independent_tag,(CartesianDKernelFunctors::Affinely_independent<K>),(Point_tag),(Point_dimension_tag,Affine_rank_tag));
 
+namespace CartesianDKernelFunctors {
+template<class R_> struct Contained_in_simplex : private Store_kernel<R_> {
+	CGAL_FUNCTOR_INIT_STORE(Contained_in_simplex)
+	typedef R_ R;
+	typedef typename Get_type<R, Point_tag>::type Point;
+	// Computing a sensible Uncertain<*> is not worth it
+	typedef typename Get_type<R, Bounded_side_tag>::type result_type;
+	typedef typename Increment_dimension<typename R::Default_ambient_dimension>::type D1;
+	typedef typename Increment_dimension<typename R::Max_ambient_dimension>::type D2;
+	typedef typename R::LA::template Rebind_dimension<D1,D2>::Other LA;
+	typedef typename LA::Dynamic_matrix Matrix;
+	typedef typename LA::Dynamic_vector DynVec;
+	typedef typename LA::Vector Vec;
+
+	template<class Iter, class P>
+	result_type operator()(Iter f, Iter e, P const&q)const{
+		typename Get_functor<R, Compute_point_cartesian_coordinate_tag>::type c(this->kernel());
+		typename Get_functor<R, Point_dimension_tag>::type pd(this->kernel());
+		int n=std::distance(f,e);
+		if (n==0) return ON_UNBOUNDED_SIDE;
+		int d=pd(q);
+		Matrix m(d+1,n);
+		DynVec a(n);
+		// FIXME: Should use the proper vector constructor (Iterator_and_last)
+		Vec b(d+1);
+		for(int j=0;j<d;++j) b[j]=c(q,j);
+		b[d]=1;
+
+		for(int i=0; ++f!=e; ++i){
+		  Point const& p = *f;
+		  for(int j=0;j<d;++j){
+		    m(i,j)=c(p,j);
+		  }
+		  m(i,d)=1;
+		}
+		if (!LA::solve(a,CGAL_MOVE(m),CGAL_MOVE(b))) return false;
+		result_type res = ON_BOUNDED_SIDE;
+		for(int i=0;i<n;++i){
+		  if (a[i]<0) return ON_UNBOUNDED_SIDE;
+		  if (a[i]==0) res = ON_BOUNDARY;
+		}
+		return res;
+	}
+};
+}
+
+CGAL_KD_DEFAULT_FUNCTOR(Contained_in_simplex_tag,(CartesianDKernelFunctors::Contained_in_simplex<K>),(Point_tag),(Point_dimension_tag,Compute_point_cartesian_coordinate_tag));
+
 #if 0
 namespace CartesianDKernelFunctors {
 template<class R_,bool=boost::is_same<typename R_::Point,typename R_::Vector>::value> struct Orientation : private Store_kernel<R_> {
diff --git a/NewKernel_d/include/CGAL/LA_eigen/LA.h b/NewKernel_d/include/CGAL/LA_eigen/LA.h
index 9fff2cb04ba..000b1db5318 100644
--- a/NewKernel_d/include/CGAL/LA_eigen/LA.h
+++ b/NewKernel_d/include/CGAL/LA_eigen/LA.h
@@ -103,6 +103,12 @@ template<class NT_,class Dim_,class Max_dim_=Dim_> struct LA_eigen {
 		return decomp.rank();
 	}
 
+	// m*a==b
+	static bool solve(Dynamic_vector&a, Dynamic_matrix const&m, Vector const& b){
+		a = m.colPivHouseholderQr().solve(b);
+		return b.isApprox(m*a);
+	}
+
 	template<class Vec1,class Vec2> static Vector homogeneous_add(Vec1 const&a,Vec2 const&b){
 		//TODO: use compile-time size when available
 		int d=a.size();