\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_to_paraboloidHd function object>>=
template <class R>
struct Lift_to_paraboloidHd {
typedef typename R::Point_d Point_d;
typedef typename R::RT RT;
typedef typename R::LA LA;

Point_d operator()(const Point_d& p) const
{ 
  int d = p.dimension();
  typename LA::Vector h(d+2);
  RT D = p.homogeneous(d);
  RT sum = 0;
  for (int i = 0; i<d; i++) {
    RT hi = p.homogeneous(i);
    h[i] = hi*D;
    sum += hi*hi;
  }
  h[d] = sum;
  h[d+1] = D*D;
  return Point_d(d+1,h.begin(),h.end());
}
};

<<Project_along_d_axisHd function object>>=
template <class R>
struct Project_along_d_axisHd {
typedef typename R::Point_d Point_d;
typedef typename R::RT RT;
typedef typename R::LA LA;

Point_d operator()(const Point_d& p) const
{ int d = p.dimension();
  return Point_d(d-1, p.homogeneous_begin(),p.homogeneous_end()-2,
                 p.homogeneous(d));
}
};

<<MidpointHd function object>>=
template <class R>
struct MidpointHd {
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
satiesfies 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_{dd}p_{id}(p_{ij}p_{dd} - p_{dj}p_{id})c_j/c_d = \sum_{0 \le j <
d} (p_{ij}p_{dd} - p_{dj}p_{id})(p_{ij}p_{dd} + p_{dj}p_{id}) \] 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_sphereHd function object>>=
template <class R>
struct Center_of_sphereHd {
typedef typename R::Point_d Point_d;
typedef typename R::RT RT;
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++;
  RT pdd  = pd.homogeneous(d);
  for (int i = 0; i < d; i++) { 
    // we set up the equation for p_i
    Point_d pi = *start++; 
    RT pid = pi.homogeneous(d);
    b[i] = 0;
    for (int j = 0; j < d; j++) {
      M(i,j) = RT(2) * pdd * pid * 
               (pi.homogeneous(j)*pdd - pd.homogeneous(j)*pid);
      b[i] += (pi.homogeneous(j)*pdd - pd.homogeneous(j)*pid) *
              (pi.homogeneous(j)*pdd + pd.homogeneous(j)*pid);
    }
  }
  RT D;
  typename LA::Vector x;
  LA::linear_solver(M,b,x,D);
  return Point_d(d,x.begin(),x.end(),D);
}

}; // Center_of_sphereHd


@ The sqared distance operations just uses the inner product:
<<Squared_distanceHd function object>>=
template <class R>
struct Squared_distanceHd {
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_lineHd function object>>=
template <class R>
struct Position_on_lineHd {
typedef typename R::Point_d Point_d;
typedef typename R::LA LA;
typedef typename R::FT FT;
typedef typename R::RT RT;

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.");
  RT lnum = (p.homogeneous(0)*s.homogeneous(d) - 
             s.homogeneous(0)*p.homogeneous(d)) * t.homogeneous(d); 
  RT lden = (t.homogeneous(0)*s.homogeneous(d) - 
             s.homogeneous(0)*t.homogeneous(d)) * p.homogeneous(d); 
  RT num(lnum), den(lden), lnum_i, lden_i; 
  for (int i = 1; i < d; i++) {  
    lnum_i = (p.homogeneous(i)*s.homogeneous(d) - 
              s.homogeneous(i)*p.homogeneous(d)) * t.homogeneous(d); 
    lden_i = (t.homogeneous(i)*s.homogeneous(d) - 
              s.homogeneous(i)*t.homogeneous(d)) * p.homogeneous(d); 
    if (lnum*lden_i != lnum_i*lden) return false; 
    if (lden_i != 0) { den = lden_i; num = lnum_i; }
  }
  l = R::make_FT(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_coordinatesHd function object>>=
template <class R>
struct Barycentric_coordinatesHd {
typedef typename R::Point_d Point_d;
typedef typename R::LA LA;
typedef typename R::RT RT;

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);
  typename LA::Matrix M(first,last);
  typename LA::Vector b(p.homogeneous_begin(),p.homogeneous_end()), x;
  RT D;
  LA::linear_solver(M,b,x,D);  
  for (int i=0; i< x.dimension(); ++result, ++i) {
    *result= R::make_FT(x[i],D); 
  }
  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>|.

Multiplying the $i$-th column by the homogenizing coordinate of
|A[i]| leaves the sign of the determinant unchanged. 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.
<<OrientationHd function object>>=
template <class R>
struct OrientationHd { 
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; ++j) 
      M(i,j) = first->homogeneous(j); 
  }
  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_sphereHd function object>>=
template <class R>
struct Side_of_oriented_sphereHd { 
typedef typename R::Point_d Point_d;
typedef typename R::LA LA;
typedef typename R::RT RT;

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 (int i = 0; i < d; ++first, ++i) { 
    RT Sum = 0;
    RT hd = first->homogeneous(d-1); 
    M(i,0) = hd*hd; 
    for (int j = 0; j < d; j++) { 
      RT hj = first->homogeneous(j); 
      M(i,j + 1) = hj * hd; 
      Sum += hj*hj; 
    }
    M(i,d) = Sum; 
  }
  RT Sum = 0; 
  RT hd = x.homogeneous(d-1); 
  M(d,0) = hd*hd; 
  for (int j = 0; j < d; j++) { 
    RT hj = x.homogeneous(j); 
    M(d,j + 1) = hj * hd; 
    Sum += hj*hj; 
  }
  M(d,d) = Sum; 
  return - 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_sphereHd function object>>=
template <class R>
struct Side_of_bounded_sphereHd { 
typedef typename R::Point_d Point_d;
typedef typename R::LA LA;
typedef typename R::RT RT;

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 o = _orientation(first,last);
  CGAL_assertion_msg((o != 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 (o == 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_simplexHd function object>>=
template <class R>
struct Contained_in_simplexHd { 
typedef typename R::Point_d Point_d;
typedef typename R::LA LA;
typedef typename R::RT RT;

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(); 
  CGAL_assertion_code(
    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(p.homogeneous_begin(),p.homogeneous_end()); 
  for (int j = 0; j < k; ++first, ++j) {
    for (int i = 0; i <= d; ++i)  
      M(i,j) = first->homogeneous(i); 
  }

  RT D; 
  typename LA::Vector lambda; 
  if ( LA::linear_solver(M,b,lambda,D) ) { 
    int s = CGAL_NTS sign(D); 
    for (int j = 0; j < k; j++) { 
      int t = CGAL_NTS sign(lambda[j]); 
      if (s * t < 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_hullHd function object>>=
template <class R>
struct Contained_in_affine_hullHd { 
typedef typename R::Point_d Point_d;
typedef typename R::LA LA;
typedef typename R::RT RT;

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(p.homogeneous_begin(),p.homogeneous_end()); 
  for (int j = 0; j < k; ++first, ++j) 
    for (int i = 0; i <= d; ++i) 
      M(i,j) = first->homogeneous(i); 
  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 points $p_0$, $p_1$,\ldots, $p_k$ is the linear
rank of vectors $p_1 - p_0$, \ldots, $p_k - p_0$.
<<Affine_rankHd function object>>=
template <class R>
struct Affine_rankHd { 
typedef typename R::Point_d Point_d;
typedef typename R::Vector_d Vector_d;
typedef typename R::LA LA;
typedef typename R::RT RT;

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.homogeneous(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_independentHd function object>>=
template <class R>
struct Affinely_independentHd { 
typedef typename R::Point_d Point_d;
typedef typename R::LA LA;
typedef typename R::RT RT;

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_lexicographicallyHd function object>>=
template <class R>
struct Compare_lexicographicallyHd {
typedef typename R::Point_d Point_d;
typedef typename R::Point_d PointD; //MSVC hack
Comparison_result operator()(const Point_d& p1, const Point_d& p2)
{ return PointD::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.  We may scale each
vector by its homogenizing coordinate and hence forget about the
homogenizing coordinates.
<<Contained_in_linear_hullHd function object>>=
template <class R>
struct Contained_in_linear_hullHd {
typedef typename R::LA LA;
typedef typename R::RT RT;
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.homogeneous(i); 
     for (int j = 0; j < k; j++) 
       M(i,j) = (first+j)->homogeneous(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.  Since
the rank of a matrix does not change under multiplication of a row
with a constant we may also use the hcoord coordinates $0$ to $d-1$ of
the vectors in A.
<<Linear_rankHd function object>>=
template <class R>
struct Linear_rankHd {
typedef typename R::LA LA;
typedef typename R::RT RT;
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)->homogeneous(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_independentHd function object>>=
template <class R>
struct Linearly_independentHd {
typedef typename R::LA LA;
typedef typename R::RT RT;
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 returns 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_baseHd function object>>=
template <class R>
struct Linear_baseHd {
typedef typename R::LA LA;
typedef typename R::RT RT;
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(); 
  RT denom; 
  typename LA::Matrix M(d,k); 
  for (int j = 0; j < k; j++) 
    for (int i = 0; i < d; i++)
      M(i,j) = (first+j)->homogeneous(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(),1);
  }
  return result; 
}
};


<<function_objectsHd.h>>=
// ======================================================================
//
// Copyright (c) 2000,2001 The CGAL Consortium
//
// This software and related documentation is part of an INTERNAL release
// of the Computational Geometry Algorithms Library (CGAL). It is not
// intended for general use.
//
// ----------------------------------------------------------------------
//
// release       : $CGAL_Revision$
// release_date  : $CGAL_Date$
//
// file          : include/CGAL/Kernel_d/function_objectsHd.h
// package       : Kernel_d
// maintainer    : Michael Seel <Michael.Seel@mpi-sb.mpg.de>
// revision      : $Id$
// revision_date : $Date$
// author(s)     : Michael Seel
// coordinator   : MPI Saarbruecken (Susan.Hert@mpi-sb.mpg.de)
//
// ======================================================================
//---------------------------------------------------------------------
// file generated by notangle from noweb/function_objectsHd.lw
// please debug or modify noweb file
// coding: K. Mehlhorn, M. Seel
//---------------------------------------------------------------------
#ifndef CGAL_FUNCTION_OBJECTSHD_H
#define CGAL_FUNCTION_OBJECTSHD_H

#include <CGAL/basic.h>
#include <CGAL/enum.h>

namespace CGAL {

<<Lift_to_paraboloidHd function object>>
<<Project_along_d_axisHd function object>>
<<MidpointHd function object>>
<<Center_of_sphereHd function object>>
<<Squared_distanceHd function object>>
<<Position_on_lineHd function object>>
<<Barycentric_coordinatesHd function object>>
<<OrientationHd function object>>
<<Side_of_oriented_sphereHd function object>>
<<Side_of_bounded_sphereHd function object>>
<<Contained_in_simplexHd function object>>
<<Contained_in_affine_hullHd function object>>
<<Affine_rankHd function object>>
<<Affinely_independentHd function object>>
<<Compare_lexicographicallyHd function object>>
<<Contained_in_linear_hullHd function object>>
<<Linear_rankHd function object>>
<<Linearly_independentHd function object>>
<<Linear_baseHd function object>>

} //namespace CGAL
#endif //CGAL_FUNCTION_OBJECTSHD_H

@ 
\end{document}