diff --git a/Kernel_d/doc/Kernel_d/CGAL/Epick_d.h b/Kernel_d/doc/Kernel_d/CGAL/Epick_d.h
index b9621559b8e..04b748c8bc4 100644
--- a/Kernel_d/doc/Kernel_d/CGAL/Epick_d.h
+++ b/Kernel_d/doc/Kernel_d/CGAL/Epick_d.h
@@ -73,7 +73,36 @@ Cartesian_const_iterator_d cartesian_begin()const;
 /*! returns an iterator pointing beyond the last Cartesian coordinate. */
 Cartesian_const_iterator_d cartesian_end()const;
 };
-/// @}
 
+/*! \cgalModels `Kernel_d::Center_of_sphere_d`
+ */
+struct Construct_circumcenter_d {
+/*! returns the center of the sphere defined by `A=tuple[first,last)`. The sphere is centered in the affine hull of A and passes through all the points of A. The order of the points of A does not matter.
+    \pre A is affinely independant.
+    \cgalRequires The value type of `ForwardIterator` is `Epick_d::Point_d`.
+    */
+template<typename ForwardIterator>
+Point_d operator()(ForwardIterator first, ForwardIterator last);
+};
+struct Compute_squared_radius_d {
+/*! returns the radius of the sphere defined by `A=tuple[first,last)`. The sphere is centered in the affine hull of A and passes through all the points of A. The order of the points of A does not matter.
+    \pre A is affinely independant.
+    \cgalRequires The value type of `ForwardIterator` is `Epick_d::Point_d`.
+    */
+template<class ForwardIterator>
+Point_d operator()(ForwardIterator first, ForwardIterator last);
+};
+/*! \cgalModels `Kernel_d::Side_of_bounded_sphere_d`
+ */
+struct Side_of_bounded_sphere_d {
+/*! returns the relative position of point p to the sphere defined by `A=tuple[first,last)`. The sphere is centered in the affine hull of A and passes through all the points of A. The order of the points of A does not matter.
+    \pre A is affinely independant.
+    \cgalRequires The value type of `ForwardIterator` is `Epick_d::Point_d`.
+    */
+template<class ForwardIterator>
+Bounded_side operator()(ForwardIterator first, ForwardIterator last, const Point_d&p);
+};
+Construct_circumcenter_d construct_circumcenter_d_object();
+Compute_squared_radius_d compute_squared_radius_d_object();
 }; /* end Epick_d */
 } /* end namespace CGAL */
diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Kernel_d_interface.h b/NewKernel_d/include/CGAL/NewKernel_d/Kernel_d_interface.h
index c54391ada61..43a074c2fed 100644
--- a/NewKernel_d/include/CGAL/NewKernel_d/Kernel_d_interface.h
+++ b/NewKernel_d/include/CGAL/NewKernel_d/Kernel_d_interface.h
@@ -130,6 +130,17 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
 	    return operator()(v, End_tag());
 	  }
 	};
+	struct Compute_squared_radius_d : private Store_kernel<Kernel> {
+	  typedef Kernel R_; // for the macro
+	  CGAL_FUNCTOR_INIT_STORE(Compute_squared_radius_d)
+	  typedef FT result_type;
+	  template<class S> FT operator()(CGAL_FORWARDABLE(S) s)const{
+	    return typename Get_functor<Base, Squared_radius_tag>::type(this->kernel())(CGAL_FORWARD(S,s));
+	  }
+	  template<class I> FT operator()(I b, I e)const{
+	    return typename Get_functor<Base, Squared_circumradius_tag>::type(this->kernel())(b,e);
+	  }
+	};
 	typedef typename Construct_cartesian_const_iterator_d::result_type Cartesian_const_iterator_d;
 	typedef typename Get_functor<Base, Squared_distance_tag>::type Squared_distance_d;
 	typedef typename Get_functor<Base, Squared_length_tag>::type Squared_length_d;
@@ -142,9 +153,10 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
 	typedef typename Get_functor<Base, Linear_rank_tag>::type Linear_rank_d;
 	typedef typename Get_functor<Base, Linearly_independent_tag>::type Linearly_independent_d;
 	typedef typename Get_functor<Base, Oriented_side_tag>::type Oriented_side_d;
-	typedef typename Get_functor<Base, Side_of_bounded_sphere_tag>::type Side_of_bounded_sphere_d;
+	typedef typename Get_functor<Base, Side_of_bounded_circumsphere_tag>::type Side_of_bounded_sphere_d;
 
 	typedef typename Get_functor<Base, Center_of_sphere_tag>::type Center_of_sphere_d;
+	typedef typename Get_functor<Base, Construct_circumcenter_tag>::type Construct_circumcenter_d;
 	typedef typename Get_functor<Base, Value_at_tag>::type Value_at_d;
 	typedef typename Get_functor<Base, Point_of_sphere_tag>::type Point_of_sphere_d;
 	typedef typename Get_functor<Base, Orthogonal_vector_tag>::type Orthogonal_vector_d;
@@ -196,10 +208,12 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
 	Construct_segment_d construct_segment_d_object()const{ return Construct_segment_d(*this); }
 	Construct_sphere_d construct_sphere_d_object()const{ return Construct_sphere_d(*this); }
 	Construct_hyperplane_d construct_hyperplane_d_object()const{ return Construct_hyperplane_d(*this); }
+	Compute_squared_radius_d compute_squared_radius_d_object()const{ return Compute_squared_radius_d(*this); }
 	Squared_distance_d squared_distance_d_object()const{ return Squared_distance_d(*this); }
 	Squared_length_d squared_length_d_object()const{ return Squared_length_d(*this); }
 	Scalar_product_d scalar_product_d_object()const{ return Scalar_product_d(*this); }
 	Center_of_sphere_d center_of_sphere_d_object()const{ return Center_of_sphere_d(*this); }
+	Construct_circumcenter_d construct_circumcenter_d_object()const{ return Construct_circumcenter_d(*this); }
 	Construct_direction_d construct_direction_d_object()const{ return Construct_direction_d(*this); }
 	Construct_line_d construct_line_d_object()const{ return Construct_line_d(*this); }
 	Construct_ray_d construct_ray_d_object()const{ return Construct_ray_d(*this); }
diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Types/Sphere.h b/NewKernel_d/include/CGAL/NewKernel_d/Types/Sphere.h
index 7ab10f16b99..ec565f2cc8c 100644
--- a/NewKernel_d/include/CGAL/NewKernel_d/Types/Sphere.h
+++ b/NewKernel_d/include/CGAL/NewKernel_d/Types/Sphere.h
@@ -52,39 +52,13 @@ template <class R_> struct Construct_sphere : Store_kernel<R_> {
   }
   template <class Iter>
   result_type operator()(Iter f, Iter e)const{
-    // 2*(x-y).c == x^2-y^2
-    typedef typename R_::LA LA;
-    typedef typename LA::Square_matrix Matrix;
-    typedef typename LA::Vector Vec;
-    typedef typename LA::Construct_vector CVec;
     typedef typename Get_type<R_, Point_tag>::type	Point;
-    typename Get_functor<R_, Compute_point_cartesian_coordinate_tag>::type c(this->kernel());
-    typename Get_functor<R_, Construct_ttag<Point_tag> >::type cp(this->kernel());
-    typename Get_functor<R_, Point_dimension_tag>::type pd(this->kernel());
-    typename Get_functor<R_, Squared_distance_to_origin_tag>::type sdo(this->kernel());
+    typename Get_functor<R_, Construct_circumcenter_tag>::type cc(this->kernel());
     typename Get_functor<R_, Squared_distance_tag>::type sd(this->kernel());
 
-    Point const& p0=*f;
-    int d = pd(p0);
-    FT const& n0 = sdo(p0);
-    Matrix m(d,d);
-    Vec b = typename CVec::Dimension()(d);
-    // Write the point coordinates in lines.
-    int i;
-    for(i=0; ++f!=e; ++i) {
-      Point const& p=*f;
-      for(int j=0;j<d;++j) {
-	m(i,j)=2*(c(p,j)-c(p0,j));
-	b[i] = sdo(p) - n0;
-      }
-    }
-    CGAL_assertion (i == d);
-    Vec res = typename CVec::Dimension()(d);;
-    //std::cout << "Mat: " << m << "\n Vec: " << one << std::endl;
-    LA::solve(res, CGAL_MOVE(m), CGAL_MOVE(b));
-    //std::cout << "Sol: " << res << std::endl;
-    Point center = cp(d,LA::vector_begin(res),LA::vector_end(res));
-    FT const& r2 = sd (center, p0);
+    // It should be possible to avoid copying the center by moving this code to a constructor.
+    Point center = cc(f, e);
+    FT const& r2 = sd(center, *f);
     return this->operator()(CGAL_MOVE(center), r2);
   }
 };
diff --git a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h
index 86f038f50b6..78602acfade 100644
--- a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h
+++ b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h
@@ -620,6 +620,123 @@ template<class R_> struct Side_of_oriented_sphere : private Store_kernel<R_> {
 
 CGAL_KD_DEFAULT_FUNCTOR(Side_of_oriented_sphere_tag,(CartesianDKernelFunctors::Side_of_oriented_sphere<K>),(Point_tag),(Point_dimension_tag,Squared_distance_to_origin_tag,Compute_point_cartesian_coordinate_tag));
 
+namespace CartesianDKernelFunctors {
+template <class R_> struct Construct_circumcenter : Store_kernel<R_> {
+  CGAL_FUNCTOR_INIT_STORE(Construct_circumcenter)
+  typedef typename Get_type<R_, Point_tag>::type Point;
+  typedef Point result_type;
+  typedef typename Get_type<R_, FT_tag>::type FT;
+  template <class Iter>
+  result_type operator()(Iter f, Iter e)const{
+    typedef typename Get_type<R_, Point_tag>::type      Point;
+    typedef typename R_::LA LA;
+    typename Get_functor<R_, Compute_point_cartesian_coordinate_tag>::type c(this->kernel());
+    typename Get_functor<R_, Construct_ttag<Point_tag> >::type cp(this->kernel());
+    typename Get_functor<R_, Point_dimension_tag>::type pd(this->kernel());
+    typename Get_functor<R_, Squared_distance_to_origin_tag>::type sdo(this->kernel());
+
+    Point const& p0=*f;
+    int d = pd(p0);
+    if (d+1 == std::distance(f,e))
+    {
+      // 2*(x-y).c == x^2-y^2
+      typedef typename LA::Square_matrix Matrix;
+      typedef typename LA::Vector Vec;
+      typedef typename LA::Construct_vector CVec;
+      FT const& n0 = sdo(p0);
+      Matrix m(d,d);
+      Vec b = typename CVec::Dimension()(d);
+      // Write the point coordinates in lines.
+      int i;
+      for(i=0; ++f!=e; ++i) {
+	Point const& p=*f;
+	for(int j=0;j<d;++j) {
+	  m(i,j)=2*(c(p,j)-c(p0,j));
+	  b[i] = sdo(p) - n0;
+	}
+      }
+      CGAL_assertion (i == d);
+      Vec res = typename CVec::Dimension()(d);;
+      //std::cout << "Mat: " << m << "\n Vec: " << one << std::endl;
+      LA::solve(res, CGAL_MOVE(m), CGAL_MOVE(b));
+      //std::cout << "Sol: " << res << std::endl;
+      return cp(d,LA::vector_begin(res),LA::vector_end(res));
+    }
+    else
+    {
+      /*
+       * Matrix P=(p1, p2, ...) (each point as a column)
+       * Matrix Q=2*t(p2-p1,p3-p1, ...) (each vector as a line)
+       * Matrix M: QP, adding a line of 1 at the top
+       * Vector B: (1, p2^2-p1^2, p3^2-p1^2, ...)
+       * Solve ML=B, the center of the sphere is PL
+       *
+       * It would likely be faster to write P then transpose, multiply,
+       * etc instead of doing it by hand.
+       */
+      // TODO: check for degenerate cases?
+
+      typedef typename R_::Max_ambient_dimension D2;
+      typedef typename R_::LA::template Rebind_dimension<Dynamic_dimension_tag,D2>::Other LAd;
+      typedef typename LAd::Square_matrix Matrix;
+      typedef typename LAd::Vector Vec;
+      typename Get_functor<R_, Scalar_product_tag>::type sp(this->kernel());
+      int k=static_cast<int>(std::distance(f,e));
+      int d=pd(p0);
+      Matrix m(k,k);
+      Vec b(k);
+      Vec l(k);
+      int j,i=0;
+      for(Iter f2=f;f2!=e;++f2,++i){
+	b(i)=m(i,i)=sdo(*f2);
+	j=0;
+	for(Iter f3=f;f3!=e;++f3,++j){
+	  m(j,i)=m(i,j)=sp(*f2,*f3);
+	}
+      }
+      for(i=1;i<k;++i){
+	b(i)-=b(0);
+	for(j=0;j<k;++j){
+	  m(i,j)=2*(m(i,j)-m(0,j));
+	}
+      }
+      for(j=0;j<k;++j) m(0,j)=1;
+      b(0)=1;
+
+      LAd::solve(l,CGAL_MOVE(m),CGAL_MOVE(b));
+
+      typename LA::Vector center=typename LA::Construct_vector::Dimension()(d);
+      for(i=0;i<d;++i) center(i)=0;
+      j=0;
+      for(Iter f2=f;f2!=e;++f2,++j){
+	for(i=0;i<d;++i){
+	  center(i)+=l(j)*c(*f2,i);
+	}
+      }
+
+      return cp(LA::vector_begin(center),LA::vector_end(center));
+    }
+  }
+};
+}
+
+CGAL_KD_DEFAULT_FUNCTOR(Construct_circumcenter_tag,(CartesianDKernelFunctors::Construct_circumcenter<K>),(Point_tag),(Construct_ttag<Point_tag>,Compute_point_cartesian_coordinate_tag,Scalar_product_tag,Squared_distance_to_origin_tag,Point_dimension_tag));
+
+namespace CartesianDKernelFunctors {
+template <class R_> struct Squared_circumradius : Store_kernel<R_> {
+  CGAL_FUNCTOR_INIT_STORE(Squared_circumradius)
+  typedef typename Get_type<R_, FT_tag>::type result_type;
+  template <class Iter>
+  result_type operator()(Iter f, Iter e)const{
+    typename Get_functor<R_, Construct_circumcenter_tag>::type cc(this->kernel());
+    typename Get_functor<R_, Squared_distance_tag>::type sd(this->kernel());
+    return sd(cc(f, e), *f);
+  }
+};
+}
+
+CGAL_KD_DEFAULT_FUNCTOR(Squared_circumradius_tag,(CartesianDKernelFunctors::Squared_circumradius<K>),(Point_tag),(Construct_circumcenter_tag,Squared_distance_tag));
+
 namespace CartesianDKernelFunctors {
 // TODO: implement it directly, it should be at least as fast as Side_of_oriented_sphere.
 template<class R_> struct Side_of_bounded_sphere : private Store_kernel<R_> {
@@ -631,6 +748,9 @@ template<class R_> struct Side_of_bounded_sphere : private Store_kernel<R_> {
 	template<class Iter>
 	result_type operator()(Iter f, Iter const& e) const {
 	  Point const& p0 = *f++; // *--e ?
+	  typename Get_functor<R, Point_dimension_tag>::type pd(this->kernel());
+	  //FIXME: Doesn't work for non-full dimension.
+	  CGAL_assertion (std::distance(f,e) == pd(p0)+1);
 	  return operator() (f, e, p0);
 	}
 
@@ -660,6 +780,24 @@ template<class R_> struct Side_of_bounded_sphere : private Store_kernel<R_> {
 
 CGAL_KD_DEFAULT_FUNCTOR(Side_of_bounded_sphere_tag,(CartesianDKernelFunctors::Side_of_bounded_sphere<K>),(Point_tag),(Side_of_oriented_sphere_tag,Orientation_of_points_tag));
 
+namespace CartesianDKernelFunctors {
+template<class R_> struct Side_of_bounded_circumsphere : private Store_kernel<R_> {
+	CGAL_FUNCTOR_INIT_STORE(Side_of_bounded_circumsphere)
+	typedef typename Get_type<R_, Bounded_side_tag>::type result_type;
+
+	template<class Iter, class P>
+	result_type operator()(Iter f, Iter const& e, P const& p0) const {
+	  // TODO: Special case when the dimension is full.
+	  typename Get_functor<R_, Construct_circumcenter_tag>::type cc(this->kernel());
+	  typename Get_functor<R_, Compare_distance_tag>::type cd(this->kernel());
+
+	  return enum_cast<Bounded_side>(cd(cc(f, e), *f, p0));
+	}
+};
+}
+
+CGAL_KD_DEFAULT_FUNCTOR(Side_of_bounded_circumsphere_tag,(CartesianDKernelFunctors::Side_of_bounded_circumsphere<K>),(Point_tag),(Squared_distance_tag,Construct_circumcenter_tag));
+
 namespace CartesianDKernelFunctors {
 template<class R_> struct Point_to_vector : private Store_kernel<R_> {
 	CGAL_FUNCTOR_INIT_STORE(Point_to_vector)
@@ -687,13 +825,13 @@ template<class R_> struct Vector_to_point : private Store_kernel<R_> {
 	typedef typename Get_type<R, RT_tag>::type RT;
 	typedef typename Get_type<R, Vector_tag>::type Vector;
 	typedef typename Get_type<R, Point_tag>::type Point;
-	typedef typename Get_functor<R, Construct_ttag<Point_tag> >::type CV;
+	typedef typename Get_functor<R, Construct_ttag<Point_tag> >::type CP;
 	typedef typename Get_functor<R, Construct_ttag<Vector_cartesian_const_iterator_tag> >::type CI;
 	typedef Point result_type;
 	typedef Vector argument_type;
 	result_type operator()(argument_type const&v)const{
 		CI ci(this->kernel());
-		return CV(this->kernel())(ci(v,Begin_tag()),ci(v,End_tag()));
+		return CP(this->kernel())(ci(v,Begin_tag()),ci(v,End_tag()));
 	}
 };
 }
diff --git a/NewKernel_d/include/CGAL/NewKernel_d/functor_tags.h b/NewKernel_d/include/CGAL/NewKernel_d/functor_tags.h
index f4e8678acc3..a184f453ba3 100644
--- a/NewKernel_d/include/CGAL/NewKernel_d/functor_tags.h
+++ b/NewKernel_d/include/CGAL/NewKernel_d/functor_tags.h
@@ -210,6 +210,7 @@ namespace CGAL {
 	CGAL_DECL_COMPUTE(Squared_distance_to_origin);
 	CGAL_DECL_COMPUTE(Squared_length);
 	CGAL_DECL_COMPUTE(Squared_radius);
+	CGAL_DECL_COMPUTE(Squared_circumradius);
 	CGAL_DECL_COMPUTE(Scalar_product);
 	CGAL_DECL_COMPUTE(Hyperplane_translation);
 	CGAL_DECL_COMPUTE(Value_at);
@@ -264,6 +265,7 @@ namespace CGAL {
 	CGAL_DECL_CONSTRUCT(Vector_to_point,Point);
 	CGAL_DECL_CONSTRUCT(Construct_min_vertex,Point);
 	CGAL_DECL_CONSTRUCT(Construct_max_vertex,Point);
+	CGAL_DECL_CONSTRUCT(Construct_circumcenter,Point);
 #undef CGAL_DECL_CONSTRUCT
 #if 0
 #define CGAL_DECL_ITER_CONSTRUCT(X,Y) struct X##_tag {}; \
@@ -293,6 +295,7 @@ namespace CGAL {
 	CGAL_DECL_PREDICATE(Orientation_of_vectors);
 	CGAL_DECL_PREDICATE(Side_of_oriented_sphere);
 	CGAL_DECL_PREDICATE(Side_of_bounded_sphere);
+	CGAL_DECL_PREDICATE(Side_of_bounded_circumsphere);
 	CGAL_DECL_PREDICATE(Contained_in_affine_hull);
 	CGAL_DECL_PREDICATE(In_flat_orientation);
 	CGAL_DECL_PREDICATE(In_flat_side_of_oriented_sphere);
diff --git a/NewKernel_d/test/NewKernel_d/Epick_d.cpp b/NewKernel_d/test/NewKernel_d/Epick_d.cpp
index ddf5a165608..860f90670c8 100644
--- a/NewKernel_d/test/NewKernel_d/Epick_d.cpp
+++ b/NewKernel_d/test/NewKernel_d/Epick_d.cpp
@@ -81,6 +81,7 @@ void test2(){
   typedef typename K1::Construct_sphere_d CSp;
   typedef typename CGAL::Get_functor<K1, CGAL::Segment_extremity_tag>::type CSE;
   typedef typename K1::Construct_cartesian_const_iterator_d CCI;
+  typedef typename K1::Construct_circumcenter_d CCc;
   typedef typename K1::Linear_base_d LB;
   typedef typename K1::Orientation_d PO;
   typedef typename K1::Side_of_oriented_sphere_d SOS;
@@ -122,6 +123,7 @@ void test2(){
   typedef typename K1::Difference_of_points_d DP;
   typedef typename K1::Construct_min_vertex_d CmV;
   typedef typename K1::Construct_max_vertex_d CMV;
+  typedef typename K1::Compute_squared_radius_d SR;
 
   CGAL_USE_TYPE(AT);
   CGAL_USE_TYPE(D);
@@ -138,6 +140,7 @@ void test2(){
   CV cv Kinit(construct_vector_d_object);
   CCI ci Kinit(construct_cartesian_const_iterator_d_object);
   CC cc Kinit(compute_coordinate_d_object);
+  CCc ccc Kinit(construct_circumcenter_d_object);
   PO po Kinit(orientation_d_object);
   CS cs Kinit(construct_segment_d_object);
   CSp csp Kinit(construct_sphere_d_object);
@@ -181,6 +184,7 @@ void test2(){
   DP dp Kinit(difference_of_points_d_object);
   CmV cmv Kinit(construct_min_vertex_d_object);
   CMV cMv Kinit(construct_max_vertex_d_object);
+  SR sr Kinit(compute_squared_radius_d_object);
 
   CGAL_USE(bc);
   CGAL_USE(pol);
@@ -229,6 +233,12 @@ void test2(){
   assert(sos(tabn+0,tabn+3,P(3,3))==CGAL::ON_POSITIVE_SIDE);
   assert(sbs(tabp+0,tabp+3,P(3,3))==CGAL::ON_UNBOUNDED_SIDE);
   assert(sbs(tabn+0,tabn+3,P(3,3))==CGAL::ON_UNBOUNDED_SIDE);
+  assert(sbs(tabp+1,tabp+3,P(1,1))==CGAL::ON_BOUNDARY);
+  assert(ccc(tabp+1,tabp+2)==tabp[1]);
+  assert(ccc(tabn+0,tabn+2)==P(0,.5));
+  assert(sr(tabp+2,tabp+3)==0);
+  assert(sr(tabp+1,tabp+3)==.5);
+  assert(sbs(tabp+1,tabp+3,P(10,-1))==CGAL::ON_UNBOUNDED_SIDE);
   assert(sos(tabp+0,tabp+3,P(.5,.5))==CGAL::ON_POSITIVE_SIDE);
   assert(sos(tabn+0,tabn+3,P(.5,.5))==CGAL::ON_NEGATIVE_SIDE);
   assert(sbs(tabp+0,tabp+3,P(.5,.5))==CGAL::ON_BOUNDED_SIDE);
@@ -401,6 +411,7 @@ void test3(){
   typedef typename K1::Point_dimension_d PD;
   typedef typename K1::Affinely_independent_d AI;
   typedef typename K1::Scaled_vector_d SV;
+  typedef typename K1::Side_of_bounded_sphere_d SBDS;
 
   Ker k
 #if 1
@@ -432,6 +443,7 @@ void test3(){
   SD sd Kinit(squared_distance_d_object);
   PD pd Kinit(point_dimension_d_object);
   AI ai Kinit(affinely_independent_d_object);
+  SBDS sbds Kinit(side_of_bounded_sphere_d_object);
   P a; // Triangulation needs this :-(
   a=cp(2,3,4);
   assert(pd(a)==3);
@@ -459,7 +471,7 @@ void test3(){
   P tab[]={a,b,c,d,e};
   std::cout << po (&tab[0],tab+4) << ' ';
   std::cout << sos(&tab[0],tab+5) << ' ';
-  std::cout << sbs(&tab[0],tab+5) << std::endl;
+  std::cout << sbs(tab+1,tab+5,tab[0]) << std::endl;
   FO fo=cfo(&tab[0],tab+3);
   std::cout << fo;
   P x[]={cp(2,2,3),cp(2,2,0),cp(1,2,1)};
@@ -559,6 +571,12 @@ void test3(){
   std::ostringstream output;
   output << showit;
   assert(output.str()=="3 1 2 4");
+  P t1[]={cp(1,2,3),cp(3,2,1),cp(2,4,2)};
+  assert(sbds(t1+0,t1+2,cp(2,2,3.414)) == CGAL::ON_BOUNDED_SIDE);
+  assert(sbds(t1+0,t1+2,cp(1,2,3)) == CGAL::ON_BOUNDARY);
+  assert(sbds(t1+0,t1+2,cp(2,2,3.415)) == CGAL::ON_UNBOUNDED_SIDE);
+  assert(sbds(t1+0,t1+3,cp(2.1,3.5,1.9)) == CGAL::ON_BOUNDED_SIDE);
+  assert(sbds(t1+0,t1+3,cp(10,10,10)) == CGAL::ON_UNBOUNDED_SIDE);
 }
 template struct CGAL::Epick_d<CGAL::Dimension_tag<2> >;
 template struct CGAL::Epick_d<CGAL::Dimension_tag<3> >;