diff --git a/Packages/Kernel_d/dont_submit b/Packages/Kernel_d/dont_submit
new file mode 100644
index 00000000000..85b8d4a174b
--- /dev/null
+++ b/Packages/Kernel_d/dont_submit
@@ -0,0 +1 @@
+noweb
diff --git a/Packages/Kernel_d/noweb/function_objectsCd.lw b/Packages/Kernel_d/noweb/function_objectsCd.lw
new file mode 100644
index 00000000000..6e4800fef3e
--- /dev/null
+++ b/Packages/Kernel_d/noweb/function_objectsCd.lw
@@ -0,0 +1,662 @@
+\documentclass[a4paper,11pt,twoside]{article}
+\usepackage{Lweb}
+\input{defs}
+
+\begin{document}
+\title{Predicates and Constructions on Geometric Objects}
+\author{M. Seel}
+\maketitle
+
+\section{Constructions on Points}
+
+Lift a point $p$ to the paraboloid of revolution in dimension
+$d+1$. The $d+1$st coordinate of $p'$ in $d+1$-space is the
+sum of squares of the components of $p$.
+<<Lift_to_paraboloidCd function object>>=
+template <class R>
+struct Lift_to_paraboloidCd {
+typedef typename R::Point_d Point_d;
+typedef typename R::FT FT;
+typedef typename R::LA LA;
+
+Point_d operator()(const Point_d& p) const
+{ 
+  int d = p.dimension();
+  typename LA::Vector h(d+1);
+  FT sum = 0;
+  for (int i = 0; i<d; i++) {
+    h[i] = p.cartesian(i);
+    sum += h[i]*h[i];
+  }
+  h[d] = sum;
+  return Point_d(d+1,h.begin(),h.end());
+}
+};
+
+@ We just ignore the $d$ component.
+<<Project_along_d_axisCd function object>>=
+template <class R>
+struct Project_along_d_axisCd {
+typedef typename R::Point_d Point_d;
+typedef typename R::FT FT;
+
+Point_d operator()(const Point_d& p) const
+{ return Point_d(p.dimension()-1, 
+                 p.cartesian_begin(),p.cartesian_end()-1); }
+};
+
+@ Midpoint can be trivially obtained via vector arithmetic.
+<<MidpointCd function object>>=
+template <class R>
+struct MidpointCd {
+typedef typename R::Point_d Point_d;
+Point_d operator()(const Point_d& p, const Point_d& q) const
+{ return Point_d(p + (q-p)/2); }
+};
+
+@ The center of a sphere defined by a range of points is only defined
+if the set of defining points are legal. If the defining points are
+all equal the sphere is trivial. So assume otherwise. Then the center
+$c$ is the unique point with equal distance to all the defining
+points. A point $c$ has equal distance to point $p_0$ and $p_i$ if it
+lies on the perpendicual bisector of $p_d$ and $p_i$, i.e., if it
+satisfies the plane equation $(p_i - p_d)^T c = (p_i - p_d) (p_i +
+p_d)/2$. Note that $p_i - p_d$ is the normal vector of the bisector
+hyperplane and $(p_i + p_d)/2$ is the midpoint of $p_d$ and $p_i$. The
+equation above translates into the equation 
+\[ \sum_{0 \le j < d} 2* (p_{ij} - p_{dj})c_j = 
+   \sum_{0 \le j < d} (p_{ij} - p_{dj})(p_{ij} + p_{dj}) 
+\] 
+for the homogeneous coordinates of the points and the center. We may
+tentatively assume that $c_d = 1$, solve the corresponding linear
+system, and then define the center.
+<<Center_of_circleCd function object>>=
+template <class R>
+struct Center_of_circleCd {
+typedef typename R::Point_d Point_d;
+typedef typename R::FT FT;
+typedef typename R::LA LA;
+template <class Forward_iterator>
+Point_d operator()(Forward_iterator start, Forward_iterator end) const
+{ CGAL_assertion(start!=end);
+  int d = start->dimension();
+  typename LA::Matrix M(d);
+  typename LA::Vector b(d);
+  Point_d pd = *start++;
+  for (int i = 0; i < d; i++) { 
+    // we set up the equation for p_i
+    Point_d pi = *start++;
+    b[i] = 0;
+    for (int j = 0; j < d; j++) {
+      M(i,j) = FT(2)*(pi.cartesian(j) - pd.cartesian(j));
+      b[i] += (pi.cartesian(j) - pd.cartesian(j)) *
+              (pi.cartesian(j) + pd.cartesian(j));
+    }
+  }
+  FT D;
+  typename LA::Vector x;
+  LA::linear_solver(M,b,x,D);
+  return Point_d(d,x.begin(),x.end());
+}
+
+}; // Center_of_circleCd
+
+
+@ The sqared distance operations just uses the inner product:
+<<Squared_distanceCd function object>>=
+template <class R>
+struct Squared_distanceCd {
+typedef typename R::Point_d Point_d;
+typedef typename R::Vector_d Vector_d;
+typedef typename R::FT FT;
+FT operator()(const Point_d& p, const Point_d& q) const
+{ Vector_d v = p-q; return v.squared_length(); }
+};
+
+@ Finally a predicate which allows to determine affine dependence of
+three points together with a quotient determining the position of a
+point with respect to the line.  The operation returns true iff point
+|p| lies on the line through |s| and |t|.  The operation provides also
+a quotient $\lambda$ such that $p = |s| + \lambda * |(t-s)|$ in this
+case. \precond $s \neq t$.  We just calculate the $\lambda$ in the
+equation $p = s + \lambda * (t-s)$. The calculation idea is that
+$\lambda = (p-s)/(t-s)$ (component wise) and one of the components
+$t_i-s_i \neq 0$ as they are different.
+<<Position_on_lineCd function object>>=
+template <class R>
+struct Position_on_lineCd {
+typedef typename R::Point_d Point_d;
+typedef typename R::LA LA;
+typedef typename R::FT FT;
+
+bool operator()(const Point_d& p, const Point_d& s, const Point_d& t, 
+     FT& l) const
+{ int d = p.dimension(); 
+  CGAL_assertion_msg((d==s.dimension())&&(d==t.dimension()&& d>0), 
+  "position_along_line: argument dimensions disagree.");
+  CGAL_assertion_msg((s!=t), 
+  "Position_on_line_d: line defining points are equal.");
+  FT lnum = (p.cartesian(0) - s.cartesian(0)); 
+  FT lden = (t.cartesian(0) - s.cartesian(0)); 
+  FT num(lnum), den(lden), lnum_i, lden_i;
+  for (int i = 1; i < d; i++) {  
+    lnum_i = (p.cartesian(i) - s.cartesian(i)); 
+    lden_i = (t.cartesian(i) - s.cartesian(i)); 
+    if (lnum*lden_i != lnum_i*lden) return false; 
+    if (lden_i != 0) { den = lden_i; num = lnum_i; }
+  }
+  l = num/den; return true; 
+}
+};
+
+@  \section{Predicates on Points}
+
+For an iterator range |[first,last)| we define |S = set [first,last)|
+as the ordered tuple $(|S[0]|,|S[1]|, \ldots |S[d-1]|)$ where 
+$|S[i]| = |*| |++|^{(i)}|first|$ (the element obtained by $i$ times 
+forwarding the iterator by operator |++| and then dereferencing it to
+get the value to which it points). We write |d = size [first,last)|.
+This extends the syntax of random access iterators to input iterators. 
+If we index the tuple as above then we require that 
+$|++|^{(d+1)}|first == last|$.
+
+In the following we require the Iterators to be globally 
+of value type |Point_d<R>|. Also if we are handed over an iterator
+range |[first,last)|, then all points in |S = set [first,last)|
+have the same dimension |dim|.
+
+The barycentric coordinates of a point $p \in R^d$ in a affine space
+of dimension $k$ spanned by the points $p_1, \ldots ,p_k$ are the
+rational numbers $\lambda_i$ with $\sum_{i=0}^k \lambda_i p_i = p$ and
+$\sum_{i=0}^k \lambda_i = 1$. Obviously the above conditions can be
+written as a linear system $Mx=b$, where $M$ is the matrix obtained by
+the $k$-tupel of points as column vectors and $b$ as $p$ interpreted
+as a vector.
+<<Barycentric_coordinatesCd function object>>=
+template <class R>
+struct Barycentric_coordinatesCd {
+typedef typename R::Point_d Point_d;
+typedef typename R::LA LA;
+typedef typename R::FT FT;
+
+template <class ForwardIterator, class OutputIterator>
+OutputIterator operator()(ForwardIterator first, ForwardIterator last, 
+  const Point_d& p, OutputIterator result)
+{ TUPLE_DIM_CHECK(first,last,Barycentric_coordinates_d);
+  int n = std::distance(first,last); 
+  int d = p.dimension();
+  typename R::Affine_rank_d affine_rank;
+  CGAL_assertion(affine_rank(first,last)==d);
+  std::vector< Point_d > V(first,last);
+  typename LA::Matrix M(d+1,V.size());
+  typename LA::Vector b(d+1), x;
+  for (register int i=0; i<d; ++i) {
+    for (register int j=0; j<V.size(); ++j) 
+      M(i,j)=V[j].cartesian(i);
+    b[i] = p.cartesian(i);
+  }
+  for (register int j=0; j<V.size(); ++j) 
+    M(d,j) = 1;
+  b[d]=1;
+  FT D;
+  LA::linear_solver(M,b,x,D);
+  for (int i=0; i < x.dimension(); ++result, ++i) {
+    *result= x[i];
+  }
+  return result;
+}
+};
+
+
+@ We determine the orientation of the points in the set |A = set
+[first,last)| where $A$ consists of $d + 1$ points in $d$ - space.
+This is the sign of the determinant
+  \[ \left\Lvert \begin{array}{cccc}
+  1 & 1 & 1 & 1 \\
+  A[0] & A[1] & \dots & A[d]
+  \end{array}  \right\Lvert  \]
+where |A[i]| denotes the cartesian coordinate vector of the $i$-th
+point in $A$.  \precond |size [first,last) == d+1| and
+|A[i].dimension() == d| $\forall 0 \leq i \leq d$ and the value type
+of |ForwardIterator| is |Point_d<R>|.
+
+We set up this matrix and return its determinant. Actually, it is more
+convenient to transpose it and to make the first row the last. This
+changes the sign if the number of rows is even, i.e., if $d$ is odd.
+<<OrientationCd function object>>=
+template <class R>
+struct OrientationCd { 
+typedef typename R::Point_d Point_d;
+typedef typename R::LA LA;
+
+template <class ForwardIterator>
+Orientation operator()(ForwardIterator first, ForwardIterator last)
+{ TUPLE_DIM_CHECK(first,last,Orientation_d);
+  int d = std::distance(first,last); 
+  // range contains d points of dimension d-1
+  CGAL_assertion_msg(first->dimension() == d-1,
+  "Orientation_d: needs first->dimension() + 1 many points.");
+  typename LA::Matrix M(d); // quadratic
+  for (int i = 0; i < d; ++first,++i) {
+    for (int j = 0; j < d-1; ++j) 
+      M(i,j) = first->cartesian(j);
+    M(i,d-1) = 1;
+  }
+  int row_correction = ( (d % 2 == 0) ? -1 : +1 );
+  // we invert the sign if the row number is even i.e. d is odd
+  return Orientation(row_correction * LA::sign_of_determinant(M));
+}
+};
+
+@ We determine to which side of the sphere $S$ defined by the points
+in |A = set [first,last)| the point $x$ belongs, where $A$ consists of
+$d + 1$ points in $d$ - space.  The positive side is determined by the
+positive sign of the determinant
+\[
+  \left\Lvert \begin{array}{ccccc} 
+  1 & 1 & 1 & 1 & 1\\
+  |lift(A[0])| & |lift(A[1])| & \dots & |lift(A[d])| & |lift(x)| 
+  \end{array} \right\Lvert 
+\] 
+where for a point $p$ with cartesian coordinates $p_i$ we use
+|lift(p)| to denote the $d + 1$-dimensional point with cartesian
+coordinate vector $(p_0, \ldots,p_{d-1},\sum_{0 \le i < d}p_i^2)$. If
+the points in $A$ are positively oriented then the positive side is
+the inside of the sphere and the negative side is the outside of the
+sphere. \precond value type of |ForwardIterator| is |Point_d<R>|.
+
+The lifting map |lift| maps a point |p| onto the paraboloid of
+revolution and the position of the point |x| with respect to the
+oriented sphere defined by the points in $A$ is the same as the
+orientation of the lifted points. In order to evaluate the determinant
+we multiply each column by the square of the homogenizing
+coordinate. This turns any column into $(h_{d}^2,h_0
+h_{d},\ldots,h_{d-1} h_{d},\sum_{0 \le i < d} h_i^2)$, where the
+$h_i$'s denote homogeneous coordinates. Thus we set up this matrix and
+return the sign of its determinant (it is actually simpler to set up
+the transposed matrix and so that's what we do).
+<<Side_of_oriented_sphereCd function object>>=
+template <class R>
+struct Side_of_oriented_sphereCd { 
+typedef typename R::Point_d Point_d;
+typedef typename R::LA LA;
+typedef typename R::FT FT;
+
+template <class ForwardIterator> 
+Oriented_side operator()(ForwardIterator first, ForwardIterator last, 
+                         const Point_d& x)
+{ 
+  TUPLE_DIM_CHECK(first,last,Side_of_oriented_sphere_d);
+  int d = std::distance(first,last); // |A| contains |d| points
+  CGAL_assertion_msg((d-1 == first->dimension()), 
+  "Side_of_oriented_sphere_d: needs first->dimension()+1 many input points.");
+  typename LA::Matrix M(d + 1); 
+  for (register int i = 0; i < d; ++first, ++i) { 
+    FT Sum = 0;
+    M(i,0) = 1;
+    for (register int j = 0; j < d-1; j++) { 
+      FT cj = first->cartesian(j);
+      M(i,j + 1) = cj; Sum += cj*cj; 
+    }
+    M(i,d) = Sum; 
+  }
+  FT Sum = 0; 
+  M(d,0) = 1; 
+  for (register int j = 0; j < d-1; j++) { 
+    FT hj = x.cartesian(j);
+    M(d,j + 1) = hj; Sum += hj*hj; 
+  }
+  M(d,d) = Sum;
+  return Oriented_side( - LA::sign_of_determinant(M) );
+}
+};
+
+@ We determine whether the point |p| lies |ON_BOUNDED_SIDE|, 
+|ON_BOUNDARY|, or |ON_UNBOUNDED_SIDE| of the sphere defined by 
+the points in |A = set [first,last)| where $A$ consists of $d+1$ 
+points in $d$-space.
+\precond value type of |ForwardIterator| is |Point_d<R>| and
+$|orientation(first,last)| \neq |ZERO|$.
+<<Side_of_bounded_sphereCd function object>>=
+template <class R>
+struct Side_of_bounded_sphereCd { 
+typedef typename R::Point_d Point_d;
+
+template <class ForwardIterator> 
+Bounded_side operator()(ForwardIterator first, ForwardIterator last, 
+                        const Point_d& p)
+{
+  TUPLE_DIM_CHECK(first,last,region_of_sphere);
+  typename R::Orientation_d _orientation;
+  Orientation or = _orientation(first,last);
+  CGAL_assertion_msg((or != 0), "Side_of_bounded_sphere_d: \
+  A must be full dimensional.");
+  typename R::Side_of_oriented_sphere_d _side_of_oriented_sphere;
+  Oriented_side oside = _side_of_oriented_sphere(first,last,p);
+  if (or == POSITIVE) {
+    switch (oside) {
+        case ON_POSITIVE_SIDE    :   return ON_BOUNDED_SIDE;
+        case ON_ORIENTED_BOUNDARY:   return ON_BOUNDARY;
+        case ON_NEGATIVE_SIDE    :   return ON_UNBOUNDED_SIDE;
+    }       
+  } else {
+    switch (oside) {
+        case ON_POSITIVE_SIDE    :   return ON_UNBOUNDED_SIDE;
+        case ON_ORIENTED_BOUNDARY:   return ON_BOUNDARY;
+        case ON_NEGATIVE_SIDE    :   return ON_BOUNDED_SIDE;
+    }     
+  }
+  return ON_BOUNDARY; // never reached
+}
+};
+
+
+@ We determine whether |p| is contained in the simplex spanned by the
+points in |A = set[first,last)|. |A| may consists of up to $d + 1$
+points. \precond value type of |ForwardIterator| is |Point_d<R>| and
+the points in |A| are affinely independent.
+
+A point |p| is contained in the convex hull of a set |A| of points if
+|p| is a convex combination of the points in |A|, i.e., if the system
+$\sum \lambda_i A_i = p$ has a solution with $\sum \lambda_i = 1$ and
+$\lambda_i \ge 0$. If the points in $A$ are linearly independent then
+the solution is unique (if it exists at all). We therefore proceed as
+above and then check whether the solution vector is non-negative.
+<<Contained_in_simplexCd function object>>=
+template <class R>
+struct Contained_in_simplexCd { 
+typedef typename R::Point_d Point_d;
+typedef typename R::LA LA;
+typedef typename R::FT FT;
+
+template <class ForwardIterator> 
+bool operator()(ForwardIterator first, ForwardIterator last,
+                const Point_d& p)
+{
+  TUPLE_DIM_CHECK(first,last,Contained_in_simplex_d);
+  int k = std::distance(first,last); // |A| contains |k| points
+  int d = first->dimension(); 
+  typename R::Affinely_independent_d check_independence;
+  CGAL_assertion_msg(check_independence(first,last),
+  "Contained_in_simplex_d: A not affinely independent.");
+  CGAL_assertion(d==p.dimension());
+
+  typename LA::Matrix M(d + 1,k); 
+  typename LA::Vector b(d +1);
+  for (register int j = 0; j < k; ++first, ++j) {
+    for (register int i = 0; i < d; ++i) 
+      M(i,j) = first->cartesian(i);
+    M(d,j) = 1;
+  }
+  for (register int i = 0; i < d; ++i) 
+    b[i] = p.cartesian(i);
+  b[d] = 1;
+
+  FT D; 
+  typename LA::Vector lambda; 
+  if ( LA::linear_solver(M,b,lambda,D) ) { 
+    for (int j = 0; j < k; j++) { 
+      if (lambda[j] < 0) return false; 
+    }
+    return true;
+  }
+  return false; 
+}
+};
+
+@ We determine whether $p$ is contained in the affine hull of the
+points in |A = set[first,last)|.  \precond value type of
+|ForwardIterator| is |Point_d<R>|.
+
+A point |p| is contained in the affine hull of a set |A| of points if
+|p| is an affine combination of the points in |A|, i.e., if the system
+$\sum \lambda_i A_i = p$ has a solution with $\sum \lambda_i = 1$. Set
+$\lambda_i = A_{i,d} \beta_i /p_d$ with $A_{i,d}$ being the
+homogenizing coordinate of $A_i$ and $p_d$ being the homogenizing
+coordinate of $p$.  The $i$-th column of the system for the
+$\beta_i$'s is simply the homogeneous vector of $A_i$ and the right
+hand side is simply the homogeneous vector for $p$. Thus we proceed as
+above but let $i$ run up to $d$.
+<<Contained_in_affine_hullCd function object>>=
+template <class R>
+struct Contained_in_affine_hullCd { 
+typedef typename R::Point_d Point_d;
+typedef typename R::LA LA;
+
+template <class ForwardIterator> 
+bool operator()(ForwardIterator first, ForwardIterator last,
+                const Point_d& p) 
+{
+  TUPLE_DIM_CHECK(first,last,Contained_in_affine_hull_d);
+  int k = std::distance(first,last); // |A| contains |k| points
+  int d = first->dimension(); 
+  typename LA::Matrix M(d + 1,k); 
+  typename LA::Vector b(d + 1); 
+  for (register int j = 0; j < k; ++first, ++j) {
+    for (register int i = 0; i < d; ++i) 
+      M(i,j) = first->cartesian(i);
+    M(d,j) = 1;
+  }
+  for (register int i = 0; i < d; ++i)
+    b[i] = p.cartesian(i);
+  b[d] = 1;
+  return LA::is_solvable(M,b);
+}
+};
+
+
+@ We compute the affine rank of the points in |A = set [first,last)|.
+\precond value type of |ForwardIterator| is |Point_d<R>|.
+
+The affine rank of the $k+1$ points $p_0$, $p_1$,\ldots, $p_k$ is the
+linear rank of the $k$ vectors $p_1 - p_0$, \ldots, $p_k - p_0$.
+<<Affine_rankCd function object>>=
+template <class R>
+struct Affine_rankCd { 
+typedef typename R::Point_d Point_d;
+typedef typename R::Vector_d Vector_d;
+typedef typename R::LA LA;
+
+template <class ForwardIterator> 
+int operator()(ForwardIterator first, ForwardIterator last) 
+{
+  TUPLE_DIM_CHECK(first,last,Affine_rank_d);
+  int k = std::distance(first,last); // |A| contains |k| points
+  if (k == 0) return -1;
+  if (k == 1) return 0; 
+  int d = first->dimension();
+  typename LA::Matrix M(d,--k);
+  Point_d p0 = *first; ++first; // first points to second
+  for (int j = 0; j < k; ++first, ++j) {
+    Vector_d v = *first - p0;
+    for (int i = 0; i < d; i++) 
+      M(i,j) = v.cartesian(i); 
+  }
+  return LA::rank(M);
+}
+};
+
+@ We decide whether the points in |A = set [first,last)| are affinely
+independent. \precond value type of |ForwardIterator| is |Point_d<R>|.
+A set of points is affinely independent if their affine rank is equal
+to the number of points minus 1.
+<<Affinely_independentCd function object>>=
+template <class R>
+struct Affinely_independentCd { 
+typedef typename R::Point_d Point_d;
+typedef typename R::LA LA;
+
+template <class ForwardIterator> 
+bool operator()(ForwardIterator first, ForwardIterator last) 
+{ typename R::Affine_rank_d rank;
+  int n = std::distance(first,last);
+  return rank(first,last) == n-1;
+}
+};
+
+
+@ We compare the Cartesian coordiantes of points |p1| and |p2|
+lexicographically.
+<<Compare_lexicographicallyCd function object>>=
+template <class R>
+struct Compare_lexicographicallyCd {
+typedef typename R::Point_d Point_d;
+Comparison_result operator()(const Point_d& p1, const Point_d& p2)
+{ return Point_d::cmp(p1,p2); }
+};
+
+@ \section{Predicates on Vectors}
+
+We determine whether $x$ is contained in the linear hull of the
+vectors in |A = set[first,last)|.  \precond value type of
+|ForwardIterator| is |Vector_d<R>|.
+
+A vector |x| is contained in the linear hull of a set |A| of vectors
+if |x| is a linear combination of the vectors in |A|, i.e., if the
+system $\sum \lambda_i A_i = x$ has a solution.  
+<<Contained_in_linear_hullCd function object>>=
+template <class R>
+struct Contained_in_linear_hullCd {
+typedef typename R::LA LA;
+typedef typename R::FT FT;
+typedef typename R::Vector_d Vector_d;
+
+template<class ForwardIterator>
+bool operator()(
+  ForwardIterator first, ForwardIterator last, const Vector_d& x) 
+{ TUPLE_DIM_CHECK(first,last,Contained_in_linear_hull_d);
+  int k = std::distance(first,last); // |A| contains |k| vectors
+  int d = first->dimension();
+  typename LA::Matrix M(d,k);
+  typename LA::Vector b(d); 
+  for (int i = 0; i < d; i++) { 
+     b[i] = x.cartesian(i); 
+     for (int j = 0; j < k; j++) 
+       M(i,j) = (first+j)->cartesian(i); 
+  }
+  return LA::is_solvable(M,b); 
+}
+};
+
+@ We compute the linear rank of the vectors in |A = set [first,last)|.
+\precond value type of |ForwardIterator| is |Vector_d<R>|.
+
+We set up a matrix having the cartesian coordinates of the vectors in
+A as its columns. The linear rank is the rank of this matrix. 
+<<Linear_rankCd function object>>=
+template <class R>
+struct Linear_rankCd {
+typedef typename R::LA LA;
+template <class ForwardIterator>
+int operator()(ForwardIterator first, ForwardIterator last)
+{ TUPLE_DIM_CHECK(first,last,linear_rank);
+  int k = std::distance(first,last); // k vectors
+  int d = first->dimension(); 
+  typename LA::Matrix M(d,k);
+  for (int i = 0; i < d  ; i++)
+     for (int j = 0; j < k; j++)  
+       M(i,j) = (first + j)->cartesian(i);
+  return LA::rank(M);
+}
+};
+
+@ We decide whether the vectors in $A$ are linearly independent.
+\precond value type of |ForwardIterator| is |Vector_d<R>|.
+
+A set of vectors is linearly independent if their linear rank is
+equal to the number of vectors in the set.
+<<Linearly_independentCd function object>>=
+template <class R>
+struct Linearly_independentCd {
+typedef typename R::LA LA;
+template <class ForwardIterator>
+bool operator()(ForwardIterator first, ForwardIterator last)
+{ typename R::Linear_rank_d rank;
+  return rank(first,last) == std::distance(first,last);
+}
+};
+
+@ We compute a basis of the linear space spanned by the vectors in
+|set [first,last)| and return it via an iterator range starting in
+|result|. The returned iterator marks the end of the output.  \precond
+value type of |ForwardIterator| and |OutputIterator| is
+|Vector_d<R>|. To find a basis within |A| we use the corresponding
+function |independent_columns| of our linear algebra package.
+<<Linear_baseCd function object>>=
+template <class R>
+struct Linear_baseCd {
+typedef typename R::LA LA;
+typedef typename R::FT FT;
+typedef typename R::Vector_d Vector_d;
+template <class ForwardIterator, class OutputIterator>
+OutputIterator operator()(ForwardIterator first, ForwardIterator last,
+  OutputIterator result)
+{ TUPLE_DIM_CHECK(first,last,linear_base);
+  int k = std::distance(first,last); // k vectors
+  int d = first->dimension();
+  FT denom; 
+  typename LA::Matrix M(d,k); 
+  for (int j = 0; j < k; ++first, ++j)
+    for (int i = 0; i < d; i++)
+      M(i,j) = first->cartesian(i);
+
+  std::vector<int> indcols;
+  int r = LA::independent_columns(M,indcols);
+
+  for (int l=0; l < r; l++) {
+    typename LA::Vector v = M.column(indcols[l]);
+    *result++ = Vector_d(d,v.begin(),v.end());
+  }
+  return result;
+}
+};
+
+
+<<function_objectsCd.h>>=
+//---------------------------------------------------------------------
+// file generated by notangle from noweb/function_objectsCd.lw
+// please debug or modify noweb file
+// coding: K. Mehlhorn, M. Seel
+//---------------------------------------------------------------------
+#ifndef CGAL_FUNCTION_OBJECTSCD_H
+#define CGAL_FUNCTION_OBJECTSCD_H
+
+#ifndef NOCGALINCL
+#include <CGAL/basic.h>
+#include <CGAL/enum.h>
+#endif
+
+#undef TRACE
+#undef TRACEN
+#undef TRACEV
+#define TRACE(t)  std::cerr << t
+#define TRACEN(t) std::cerr << t << std::endl
+#define TRACEV(t) std::cerr << #t << " = " << (t) << std::endl
+ 
+CGAL_BEGIN_NAMESPACE
+
+<<Lift_to_paraboloidCd function object>>
+<<Project_along_d_axisCd function object>>
+<<MidpointCd function object>>
+<<Center_of_circleCd function object>>
+<<Squared_distanceCd function object>>
+<<Position_on_lineCd function object>>
+<<Barycentric_coordinatesCd function object>>
+<<OrientationCd function object>>
+<<Side_of_oriented_sphereCd function object>>
+<<Side_of_bounded_sphereCd function object>>
+<<Contained_in_simplexCd function object>>
+<<Contained_in_affine_hullCd function object>>
+<<Affine_rankCd function object>>
+<<Affinely_independentCd function object>>
+<<Compare_lexicographicallyCd function object>>
+<<Contained_in_linear_hullCd function object>>
+<<Linear_rankCd function object>>
+<<Linearly_independentCd function object>>
+<<Linear_baseCd function object>>
+
+CGAL_END_NAMESPACE
+#endif //CGAL_FUNCTION_OBJECTSCD_H
+
+@ 
+\end{document}