\documentclass[a4paper,11pt,twoside]{article}
\usepackage{Lweb}
\input{defs}

\begin{document}

\section{Two elementary operations}

\subsection{Intersecting two lines}

Here we implement a $d$-dimensional intersection routine for two
non-parallel lines. We do the intersection and node calculation by the
help of some matrix calculation. We know that a common point $x$ of
two lines $s_1t_1$ and $s_2t_2$ through two $d$-dimensional points
obeys the equations:
\begin{eqnarray}
x_1 & = & s_1 + \lambda_1 (t_1 - s_1) \\
x_2 & = & s_2 + \lambda_2 (t_2 - s_2) 
\end{eqnarray}
Assume $x_1 = x_2$ and rearrange with variables $\lambda_1$ and
$\lambda_2$. Then we get the $d \times 2$ cartesian matrix system
\begin{displaymath}
\begin{bmatrix}
t_{1_0} - s_{1_0} & s_{2_0} - t_{2_0} \\
\vdots & \vdots \\
t_{1_{d-1}} - s_{1_{d-1}} & s_{2_{d-1}} - t_{2_{d-1}} 
\end{bmatrix}
* 
\begin{bmatrix} \lambda_1 \\ \lambda_2 \end{bmatrix}
=
\begin{bmatrix}
  s_{2_0} - s_{1_0} \\
  \vdots \\
  s_{2_{d-1}} - s_{1_{d-1}} 
\end{bmatrix}
\end{displaymath}
The code below puts this linear system into a matrix after
multiplication of each row with its common denominator $s_{1_d}
s_{2_d} t_{1_d} t_{2_d}$. Afterwards we solve the system. If we get a
solution we calculate the intersection point $x$ and deliver
additionally $\lambda_1 = l1_{num}/l1_{den}$ and $\lambda_2 =
l2_{num}/l2_{den}$ as rational numbers for outside use.

Note the possible exceptions here. Each of the direction vectors of
the lines can be the zero vectors in case that its spanning points 
are identical. In which case we just return that result.
<<homogeneous intersection working horses>>=
/*{\Manpage{Line_line_intersectionHd}{R}{intersecting two lines}}*/

template <class R> 
class Line_line_intersectionHd {

typedef typename R::RT RT;
typedef typename R::FT FT;
typedef typename R::LA LA;
typedef typename R::Point_d Point_d;

public:
enum Intersection_result { NO, POINT, LINE };

Intersection_result operator()(
  const Point_d& s1, const Point_d& t1,
  const Point_d& s2, const Point_d& t2,
  Point_d& p, FT& l1, FT& l2) 
/*{\Mfunop returns |NO| if the lines which are represented by |s1t1|
and |s2t2| don't intersect, returns |POINT| if they intersect in a
unique point, and returns LINE if they are identical. In the |POINT|
case the point of intersection is assigned to |p|.  Then |p = s1 + l1
* t1-s1| and |p = s2 + l2 * t2-s2|. \precond none of the point pairs
is degenerate.}*/
{ 
  int d = s1.dimension(),i; 
  CGAL_assertion_msg(d==s2.dimension(),
    "intersection: dimensions disagree!"); 
  typename LA::Matrix M(d,2),S; 
  typename LA::Vector b(d), lambda(2), c; 
  RT D; 
  
  RT s1w = s1.homogeneous(d); 
  RT t1w = t1.homogeneous(d); 
  RT s2w = s2.homogeneous(d); 
  RT t2w = t2.homogeneous(d); 
  RT g1w = s1w*t1w; 
  RT g2w = s2w*t2w; 
  RT t12w = t1w*t2w; 

  /* init $d \times 2$ - matrix |M| and $d$ - vector |b| */
  for (i = 0; i < d; i++) { 
    M(i,0) = g2w * (t1.homogeneous(i) * s1w - s1.homogeneous(i) * t1w); 
    M(i,1) = g1w * (s2.homogeneous(i) * t2w - t2.homogeneous(i) * s2w); 
    b[i]   = t12w * (s2.homogeneous(i) * s1w - s1.homogeneous(i) * s2w); 
  }

  if (LA::linear_solver(M,b,lambda,D,S,c)) { 
    if (S.column_dimension()>0) return LINE;
    l1 = R::make_FT(lambda[0],D);
    l2 = R::make_FT(lambda[1],D);
    p = s1 + l1 * (t1 - s1);
    return POINT; 
  }
  return NO; 
}

};

@ Now we implement the cartesian $d$-dimensional intersection routine
for two non-parallel lines. We do the intersection and node
calculation by the help of some matrix calculation.  The code below
puts this linear system from above directly into a matrix. Afterwards
we solve the system. If we get a solution we calculate the
intersection point $x$ and deliver additionally $\lambda_1$ and
$\lambda_2$ for outside use.

Note the possible exceptions here. Each of the direction vectors of
the lines can be the zero vectors in case that its spanning points 
are identical. In which case we just return that result.
<<cartesian intersection working horses>>=
/*{\Manpage{Line_line_intersectionCd}{R}{intersecting two lines}}*/

template <class R> 
class Line_line_intersectionCd {

typedef typename R::FT FT;
typedef typename R::LA LA;
typedef typename R::Point_d Point_d;
typedef typename R::Line_d Line_d;

public:
enum Intersection_result { NO, POINT, LINE };

Intersection_result operator()(
  const Point_d& s1, const Point_d& t1,
  const Point_d& s2, const Point_d& t2,
  Point_d& p, FT& l1, FT& l2) 
/*{\Mfunop returns |NO| if the lines which are represented by |s1t1|
and |s2t2| don't intersect, returns |POINT| if they intersect in a
unique point, and returns LINE if they are identical. In the |POINT|
case the point of intersection is assigned to |p|.  Then |p = s1 + l1
* t1-s1| and |p = s2 + l2 * t2-s2|. \precond none of the point pairs
is degenerate.}*/
{ 
  int d = s1.dimension(),i; 
  CGAL_assertion_msg(d==s2.dimension(),
    "intersection: dimensions disagree!"); 
  typename LA::Matrix M(d,2),S; 
  typename LA::Vector b(d), lambda(2), c; 
  FT D; 

  /* init $d \times 2$ - matrix |M| and $d$ - vector |b| */
  for (i = 0; i < d; i++) { 
    M(i,0) = t1.cartesian(i) - s1.cartesian(i); 
    M(i,1) = s2.cartesian(i) - t2.cartesian(i); 
    b[i]   = s2.cartesian(i) - s1.cartesian(i); 
  }

  if (LA::linear_solver(M,b,lambda,D,S,c)) {
    if ( S.column_dimension()>0 ) return LINE;
    l1 = lambda[0]; l2 = lambda[1];
    p = s1 + l1 * (t1 - s1); 
#ifdef CGAL_CHECK_EXACTNESS
    Line_d L1(s1,t1), L2(s2,t2);
    CGAL_assertion(L1.has_on(p)&&L2.has_on(p));
#endif
    return POINT; 
  }
  return NO; 
}
};

@ \subsection{The intersection of a line and a hyperplane}

We want to calculate if a line defined by two points $s$ and $t$
intersects a hyperplane $h$. For the hyperplane intersection we first
check whether $s$ and $t$ lie in $h$. If both are contained we are
done with result |LINE|.  If only one is contained we have already
found a single point of intersection for result |POINT|.  If none is
contained there is at most one point of intersection.  Let $s$ and $t$
be distinct points on $l$. Note that the evaluation of $h$ at point
$h(p) := h_d + \sum_{i=0}^{d-1} h_i |p.cartesian(i)|$ is different
from the evaluation with its homogeneous representation
$\sum_{i=0}^{d} h_i |p.homogeneous(i)|$. Let us write the latter
sloppily as an inner product of the representing vectors $h \cdot
p$. We have $h \cdot p = h(p) * p_d$.

Now let $S = s \cdot h$ and $T = t \cdot h$.  If $h(s) = S/s_d = T/t_d
= h(t)$ then the line is parallel to the plane and there is no
intersection. The latter translates to the test $D := S t_d - T s_d =
0$? If non-zero, the point with homogeneous coordinates $(S t_0 - T
s_0, \ldots, S t_i - T s_i, \ldots, S t_d - T s_d)$ lies on $h$. For
if we interpret the tuples as vectors we get $p = St - Ts$ and the
homogeneous evaluation of $h$ on this point is $s \cdot h * t \cdot h
- t \cdot h * s \cdot h$. Note that with $p = s + \lambda (t-s)$ we
get for a valid point on the hyperplane $h(p) = 0 \Leftrightarrow h(s)
+ \lambda (h(t)-h(s)) = 0$ which gets us $\lambda = - h(s) /
(h(t)-h(s)) = S t_d / D$.
<<homogeneous intersection working horses>>=
/*{\Manpage {Line_hyperplane_intersectionHd}{R} 
{intersecting a line and a hyperplane}}*/

template <class R> 
class Line_hyperplane_intersectionHd {
typedef typename R::RT RT;
typedef typename R::FT FT;
typedef typename R::LA LA;
typedef typename R::Point_d      Point_d;
typedef typename R::Hyperplane_d Hyperplane_d;

public:
enum Intersection_result { NO, POINT, LINE };

Intersection_result operator()(const Point_d& s, const Point_d& t,
  const Hyperplane_d& h, Point_d& p, FT& lambda) 
/*{\Mfunop returns |NO| if the line represented by |s1t1| and the
hyperplane |h| don't intersect, returns |POINT| if they intersect in a
unique point, and returns LINE if the line is part of the
hyperplane. In the |POINT| case the point of intersection is assigned
to |p|.  Then |p = s1 + lambda * t1-s1|.  \precond the point pair is
not degenerate.}*/
{ 
  CGAL_assertion_msg((h.dimension()==s.dimension() && 
                      h.dimension()==t.dimension()), 
  "Line_hyperplane_intersection_d: dimensions do not agree.");

  int d = h.dimension(),i;
  RT S(0),T(0);
  for (i=0; i<=d; ++i) {
    S += h[i]*s.homogeneous(i);
    T += h[i]*t.homogeneous(i);
  }
  bool s_contained = CGAL_NTS is_zero(S), 
       t_contained = CGAL_NTS is_zero(T);
  if (s_contained && t_contained) { p = s; return LINE; }
  if (s_contained) { p = s; return POINT; }
  if (t_contained) { p = t; return POINT; }
  // now the simple cases are done 

  RT D = S * t.homogeneous(d) - T * s.homogeneous(d);
  if (CGAL_NTS is_zero(D)) return NO;

  typename LA::Vector homog(d + 1);
  for (i = 0; i < d; ++i)
    homog[i] = S * t.homogeneous(i) - T * s.homogeneous(i); 
  homog[d] = D;
  p = Point_d(d,homog.begin(),homog.end());
  lambda = R::make_FT(S * t.homogeneous(d), D);
  return POINT;
}

};

@ Again we save on the homogenizing denominators.  We want to
calculate if a line defined by two points $s$ and $t$ intersects a
hyperplane $h$. For the hyperplane intersection we first check whether
$s$ and $t$ lie in $h$. If both are contained we are done with result
|LINE|.  If only one is contained we have already found a single point
of intersection for result |POINT|.  If none is contained there is at
most one point of intersection.  Let $s$ and $t$ be distinct points on
$l$. Note that the evaluation of $h$ at point $h(p) := h_d +
\sum_{i=0}^{d-1} h_i |p.cartesian(i)|$ is a measure for the distance
of $p$ to the hyperplane $h$.

Now let $S = h(s)$ and $T = h(t)$.  If $h(s) = h(t)$ then the line is
parallel to the plane and there is no intersection. The latter
translates to the test $D := S - T = 0$? If non-zero, the point with
homogeneous coordinates $p = \frac{S}{D} t - \frac{T}{D} s$ lies on
$h$. For the evaluation of $h$ on this point is $\frac{S}{D} h(t) -
\frac{T}{D} h(s) = (h(s)h(t) - h(t)h(s))/D = 0$.  Note that with $p =
s + \lambda (t-s)$ we get for a valid point on the hyperplane $h(p) =
0 \Leftrightarrow h(s) + \lambda (h(t)-h(s)) = 0$ which gets us
$\lambda = - h(s) / (h(t)-h(s)) = S / D$.
<<cartesian intersection working horses>>=
/*{\Manpage {Line_hyperplane_intersectionCd}{R} 
{intersecting a line and a hyperplane}}*/

template <class R> 
class Line_hyperplane_intersectionCd {

typedef typename R::FT FT;
typedef typename R::LA LA;
typedef typename R::Point_d Point_d;
typedef typename R::Hyperplane_d Hyperplane_d;

public:
enum Intersection_result { NO, POINT, LINE };

Intersection_result operator()(const Point_d& s, const Point_d& t,
  const Hyperplane_d& h, Point_d& p, FT& lambda) 
/*{\Mfunop returns |NO| if the line represented by |s1t1| and the
hyperplane |h| don't intersect, returns |POINT| if they intersect in a
unique point, and returns LINE if the line is part of the
hyperplane. In the |POINT| case the point of intersection is assigned
to |p|.  Then |p = s1 + lambda * t1-s1|.  \precond the point pair is
not degenerate.}*/
{ 
  CGAL_assertion_msg((h.dimension()==s.dimension() && 
                      h.dimension()==t.dimension()), 
  "Line_hyperplane_intersection_d: dimensions do not agree.");

  int d = h.dimension(),i;
  FT S = h.value_at(s), T = h.value_at(t);

  bool s_contained = CGAL_NTS is_zero(S),
       t_contained = CGAL_NTS is_zero(T);
  if (s_contained && t_contained) { p = s; return LINE; }
  if (s_contained) { p = s; return POINT; }
  if (t_contained) { p = t; return POINT; }
  // now the simple cases are done 

  FT D = S - T;
  if ( CGAL_NTS is_zero(D) ) return NO;

  typename LA::Vector v(d);
  for (i = 0; i < d; ++i)
    v[i] = (S * t.cartesian(i) - T * s.cartesian(i))/D;
  p = Point_d(d,v.begin(),v.end()); lambda = S/D;

#ifdef CGAL_CHECK_EXACTNESS
  Line_d l(s,t);
  CGAL_assertion(h.has_on(p)&&l.has_on(p));
#endif

  return POINT;
}

};


@ \section{The interface intersection pairs}

We use the engines above to write generic intersection routines.
<<line line intersection>>=
template <class R>
class Line_d_Line_d_pair {
public:
  enum Intersection_result {NO, POINT, LINE};
  typedef typename R::Point_d Point_d;
  typedef typename R::Line_d  Line_d;
  typedef typename R::FT FT;
protected:
  Line_d _l1, _l2;
  bool _known;
  Intersection_result _result;
  Point_d _ip;
public:
  Line_d_Line_d_pair() : _known(false) {}
  Line_d_Line_d_pair(const Line_d& l1, const Line_d& l2) 
    : _l1(l1), _l2(l2), _known(false) {}
  Intersection_result intersection_type();
  bool intersection(Point_d& result);
  bool intersection(Line_d& result);
};

template <class R>
Line_d_Line_d_pair<R>::Intersection_result
Line_d_Line_d_pair<R>::intersection_type()
{ 
  if (_known) return _result;
  _known = true;

  //CGAL_assertion(!_l1.is_degenerate()&&!_l2.is_degenerate());
  typedef typename R::Line_line_intersection_d Int_obj_type;
  Int_obj_type Intersect;
  FT l1,l2; 
  typename Int_obj_type::Intersection_result res = 
    Intersect(_l1.point(0),_l1.point(1),
              _l2.point(0),_l2.point(1),
              _ip,l1,l2);

  if (res == Int_obj_type::LINE)  { return _result = LINE; }
  if (res == Int_obj_type::POINT) { return _result = POINT; }
  return _result = NO; 
}

template <class R>
bool Line_d_Line_d_pair<R>::intersection(Point_d& p)
{ if (!_known) intersection_type();
  if (_result != POINT) return false;
  p = _ip; return true;
}

template <class R>
bool Line_d_Line_d_pair<R>::intersection(Line_d& l)
{ if (!_known) intersection_type();
  if (_result != LINE) return false;
  l = _l1; return true;
}

<<line ray intersection>>=
template <class R>
class Line_d_Ray_d_pair {
public:
  enum Intersection_result {NO, POINT, RAY};
  typedef typename R::Point_d Point_d;
  typedef typename R::Ray_d   Ray_d;
  typedef typename R::Line_d  Line_d;
  typedef typename R::FT FT;
protected:
  Line_d _l; Ray_d _r;
  bool _known;
  Intersection_result _result;
  Point_d _ip; 
public:
  Line_d_Ray_d_pair() : _known(false) {}
  Line_d_Ray_d_pair(const Line_d& l, const Ray_d& r) 
    : _l(l), _r(r), _known(false) {}
  Intersection_result intersection_type();
  bool intersection(Point_d& result);
  bool intersection(Ray_d& result);
};

template <class R>
Line_d_Ray_d_pair<R>::Intersection_result
Line_d_Ray_d_pair<R>::intersection_type()
{ 
  if (_known) return _result;
  _known = true;

  //CGAL_assertion(!_l.is_degenerate()&&!_r.is_degenerate());
  typedef typename R::Line_line_intersection_d Int_obj_type;
  Int_obj_type Intersect;
  FT l1,l2; 
  typename Int_obj_type::Intersection_result res = 
    Intersect(_l.point(0),_l.point(1),
              _r.point(0),_r.point(1),
              _ip,l1,l2);

  if ( res == Int_obj_type::LINE )  { return _result = RAY; }
  if ( res == Int_obj_type::POINT &&
       l2 >= FT(0) ) 
  { return _result = POINT; }
  return _result = NO; 
}

template <class R>
bool Line_d_Ray_d_pair<R>::intersection(Point_d& p)
{ if (!_known) intersection_type();
  if (_result != POINT) return false;
  p = _ip; return true;
}

template <class R>
bool Line_d_Ray_d_pair<R>::intersection(Ray_d& r)
{ if (!_known) intersection_type();
  if (_result != RAY) return false;
  r = _r; return true;
}

<<line segment intersection>>=
template <class R>
class Line_d_Segment_d_pair {
public:
  enum Intersection_result {NO, POINT, SEGMENT};
  typedef typename R::Point_d Point_d;
  typedef typename R::Segment_d Segment_d;
  typedef typename R::Line_d  Line_d;
  typedef typename R::FT FT;
protected:
  Line_d _l; Segment_d _s;
  bool _known;
  Intersection_result _result;
  Point_d _ip; 
public:
  Line_d_Segment_d_pair() : _known(false) {}
  Line_d_Segment_d_pair(const Line_d& l, const Segment_d& s) 
    : _l(l), _s(s), _known(false) {}
  Intersection_result intersection_type();
  bool intersection(Point_d& result);
  bool intersection(Segment_d& result);
};

template <class R>
Line_d_Segment_d_pair<R>::Intersection_result
Line_d_Segment_d_pair<R>::intersection_type()
{ 
  if (_known) return _result;
  _known = true;

  //CGAL_assertion(!_l.is_degenerate());
  if ( _s.is_degenerate() ) {
    if ( _l.has_on(_s.point(0)) ) {
      _ip = _s.point(0);
      return _result = POINT;
    }
    return _result = NO; 
  }
  // _s not degenerate
  typedef typename R::Line_line_intersection_d Int_obj_type;
  Int_obj_type Intersect;
  FT l1,l2; 
  typename Int_obj_type::Intersection_result res = 
    Intersect(_l.point(0),_l.point(1),
              _s.point(0),_s.point(1),
              _ip,l1,l2);

  if ( res == Int_obj_type::LINE )  { return _result = SEGMENT; }
  if ( res == Int_obj_type::POINT &&
       FT(0) <= l2 && l2 <= FT(1) ) 
  { return _result = POINT; }
  return _result = NO; 
}

template <class R>
bool Line_d_Segment_d_pair<R>::intersection(Point_d& p)
{ if (!_known) intersection_type();
  if (_result != POINT) return false;
  p = _ip; return true;
}

template <class R>
bool Line_d_Segment_d_pair<R>::intersection(Segment_d& s)
{ if (!_known) intersection_type();
  if (_result != SEGMENT) return false;
  s = _s; return true;
}


<<ray ray intersection>>=
template <class R>
class Ray_d_Ray_d_pair {
public:
  enum Intersection_result { NO, POINT, SEGMENT, RAY };
  typedef typename R::FT FT;
  typedef typename R::Point_d Point_d;
  typedef typename R::Segment_d Segment_d;
  typedef typename R::Ray_d Ray_d;

protected:
  Ray_d _r1, _r2;
  bool _known;
  Intersection_result _result;
  Point_d _ip; Segment_d _is; Ray_d _ir;
public:
  Ray_d_Ray_d_pair() : _known(false) {}
  Ray_d_Ray_d_pair(const Ray_d& r1, const Ray_d& r2)
    : _r1(r1), _r2(r2), _known(false) {}
  Intersection_result intersection_type();
  bool intersection(Point_d& result);
  bool intersection(Segment_d& result);
  bool intersection(Ray_d& result);
};

template <class R>
Ray_d_Ray_d_pair<R>::Intersection_result
Ray_d_Ray_d_pair<R>::intersection_type()
{ 
  if (_known) return _result;
  _known = true;

  //CGAL_assertion((!_r1.is_degenerate()&&!_r2.is_degenerate());
  // none of the lines should be trivial

  typedef typename R::Line_line_intersection_d Int_obj_type;
  Int_obj_type Intersect;
  FT l1,l2; 
  typename Int_obj_type::Intersection_result res = 
    Intersect(_r1.point(0),_r1.point(1),
              _r2.point(0),_r2.point(1),_ip,l1,l2);

  if (res == Int_obj_type::LINE)   
  <<determine common points of two rays>>

  if (res == Int_obj_type::POINT && 
      FT(0) <= l1 && FT(0) <= l2) 
  { return _result = POINT; }
  return _result = NO; 
  // now not parallel
}

@ We determine the common points of a ray and a segment whose
supporting lines are equal.
<<determine common points of two rays>>=
{
  if ( _r1.direction() == _r2.direction() ) {
    if ( _r1.has_on(_r2.source()) ) _ir = _r2;
    else _ir = _r1; 
    return _result = RAY;
  }
  // now oppositely directed:
  if ( _r1.has_on(_r2.source()) ) {
    if ( _r1.source() != _r2.source() ) 
    { _is = Segment_d(_r1.source(),_r2.source()); 
      return _result = SEGMENT; }
    else
    { _ip = _r1.source(); return _result = POINT; }
  }
  return _result = NO;
}

<<ray ray intersection>>=
template <class R>
bool
Ray_d_Ray_d_pair<R>::intersection(Point_d& p)
{ if (!_known) intersection_type();
  if (_result != POINT) return false;
  p = _ip; return true;
}

template <class R>
bool
Ray_d_Ray_d_pair<R>::intersection(Segment_d& s)
{ if (!_known) intersection_type();
  if (_result != SEGMENT) return false;
  s = _is; return true;
}

template <class R>
bool
Ray_d_Ray_d_pair<R>::intersection(Ray_d& r)
{ if (!_known) intersection_type();
  if (_result != RAY) return false;
  r = _ir; return true;
}

<<ray segment intersection>>=
template <class R>
class Ray_d_Segment_d_pair {
public:
  enum Intersection_result { NO, POINT, SEGMENT };
  typedef typename R::FT FT;
  typedef typename R::Point_d Point_d;
  typedef typename R::Segment_d Segment_d;
  typedef typename R::Ray_d Ray_d;

protected:
  Ray_d _r; Segment_d _s;
  bool _known;
  Intersection_result _result;
  Point_d _ip; Segment_d _is;  
public:
  Ray_d_Segment_d_pair() : _known(false) {}
  Ray_d_Segment_d_pair(const Ray_d& r, const Segment_d& s)
    : _r(r), _s(s), _known(false) {}
  Intersection_result intersection_type();
  bool intersection(Point_d& result);
  bool intersection(Segment_d& result);
};

template <class R>
Ray_d_Segment_d_pair<R>::Intersection_result
Ray_d_Segment_d_pair<R>::intersection_type()
{ 
  if (_known) return _result;
  _known = true;

  //CGAL_assertion(!_r.is_degenerate());
  if ( _s.is_degenerate() ) {
    if ( _r.has_on(_s.point(0)) ) 
    { _ip = _s.point(0); return _result = POINT; }
    return _result = NO; 
  }
  // now s is not degenerate
  typedef typename R::Line_line_intersection_d Int_obj_type;
  Int_obj_type Intersect;
  FT l1,l2; 
  typename Int_obj_type::Intersection_result res = 
    Intersect(_r.point(0),_r.point(1),
              _s.point(0),_s.point(1),_ip,l1,l2);

  if (res == Int_obj_type::LINE)   
  <<determine common points of a ray and a segment>>

  if (res == Int_obj_type::POINT && 
      FT(0) <= l1 && FT(0) <= l2 && l2 <= FT(1)) 
  { return _result = POINT; }
  return _result = NO; 
  // now not parallel
}

@ We determine the common points of a ray and a segment whose
supporting lines are equal.
<<determine common points of a ray and a segment>>=
{
  Point_d p1 = _s.point(0), p2 = _s.point(1);
  if ( _s.direction() != _r.direction() ) std::swap(p1,p2);
  // now order p2 after p1 on r underlying line 
  typename R::Position_on_line_d pos;
  pos(p1, _r.point(0),_r.point(1), l1);
  pos(p2, _r.point(0),_r.point(1), l2);
  if ( l1 < FT(0) ) p1 = _r.point(0);
  if ( l2 < FT(0) ) { return _result = NO; }
  if ( p1 == p2 ) { _ip = p1; return _result = NO; }
  _is = Segment_d(p1,p2);
  return _result = SEGMENT;
}

<<ray segment intersection>>=
template <class R>
bool
Ray_d_Segment_d_pair<R>::intersection(Point_d& p)
{ if (!_known) intersection_type();
  if (_result != POINT) return false;
  p = _ip; return true;
}

template <class R>
bool
Ray_d_Segment_d_pair<R>::intersection(Segment_d& s)
{ if (!_known) intersection_type();
  if (_result != SEGMENT) return false;
  s = _is; return true;
}


<<segment segment intersection>>=
template <class R>
class Segment_d_Segment_d_pair {
public:
  enum Intersection_result { NO, POINT, SEGMENT };
  typedef typename R::FT FT;
  typedef typename R::Point_d Point_d;
  typedef typename R::Segment_d Segment_d;

protected:
  Segment_d _s1,_s2;
  bool _known;
  Intersection_result _result;
  Point_d _ip; Segment_d _is;  
public:
  Segment_d_Segment_d_pair() : _known(false) {}
  Segment_d_Segment_d_pair(const Segment_d& s1, const Segment_d& s2)
    : _s1(s1), _s2(s2), _known(false) {}
  Intersection_result intersection_type();
  bool intersection(Point_d& result);
  bool intersection(Segment_d& result);
};

template <class R>
Segment_d_Segment_d_pair<R>::Intersection_result
Segment_d_Segment_d_pair<R>::intersection_type()
{ 
  if (_known) return _result;
  _known = true;

  if ( _s1.is_degenerate() ) {
    if ( _s2.is_degenerate() ) {
      if ( _s1.point(0) == _s2.point(0) ) 
      { _ip = _s1.point(0); return _result = POINT; }
      else { return _result = NO; }
    } else {
      if ( _s2.has_on(_s1.point(0)) ) 
      { _ip = _s1.point(0); return _result = POINT; }
      else { return _result = NO; }
    }
  }
  if ( _s2.is_degenerate() ) { 
    CGAL_assertion( !_s1.is_degenerate() );
    if ( _s1.has_on(_s2.point(0)) ) 
    { _ip = _s2.point(0); return _result = POINT; }
    else { return _result = NO; }
  }  

  // now s1,s2 not degenerate
  typedef typename R::Line_line_intersection_d Int_obj_type;
  Int_obj_type Intersect;
  FT l1,l2; 
  typename Int_obj_type::Intersection_result res = 
    Intersect(_s1.point(0),_s1.point(1),
              _s2.point(0),_s2.point(1),_ip,l1,l2);

  if (res == Int_obj_type::LINE)   
  <<determine common points of two segments>>

  if ( res == Int_obj_type::POINT && 
       FT(0) <= l1 && l1 <= FT(1) && 
       FT(0) <= l2 && l2 <= FT(1) ) 
  { return _result = POINT; }
  return _result = NO; 
}

<<determine common points of two segments>>=
{ Point_d p1 = _s1.min(), p2 = _s1.max();
  Point_d q1 = _s2.min(), q2 = _s2.max();
  Point_d s,t;
  // now order the for points along the line
  typename R::Position_on_line_d pos;
  pos(p1, q1, q2, l1); pos(p2, q1, q2, l2);
  if ( l1 < FT(0) ) {
    if ( l2 < FT(0) ) { return _result = NO; }
    else { s = q1; }
  } else { s = p1; } // l1 >= 0
  if ( l2 > FT(1) ) {
    if ( l1 > FT(1) ) { return _result = NO; }
    else { t = q2; }
  } else { t = p2; } // l2 <= 1
  if ( s == t ) { _ip = s; return _result = POINT; }
  _is = Segment_d(s,t); return _result = SEGMENT;
}

<<segment segment intersection>>=
template <class R>
bool
Segment_d_Segment_d_pair<R>::intersection(Point_d& p)
{ if (!_known) intersection_type();
  if (_result != POINT) return false;
  p = _ip; return true;
}

template <class R>
bool
Segment_d_Segment_d_pair<R>::intersection(Segment_d& s)
{ if (!_known) intersection_type();
  if (_result != SEGMENT) return false;
  s = _is; return true;
}


<<line hyperplane intersection>>=
template <class R>
class Line_d_Hyperplane_d_pair {
public:
  enum Intersection_result {NO, POINT, LINE};
  typedef typename R::Point_d      Point_d;
  typedef typename R::Line_d       Line_d;
  typedef typename R::Hyperplane_d Hyperplane_d;
  typedef typename R::FT           FT;

protected:
  Line_d _l; Hyperplane_d _h;
  bool _known;
  Intersection_result _result;
  Point_d _ip;
public:
  Line_d_Hyperplane_d_pair() : _known(false) {}
  Line_d_Hyperplane_d_pair(const Line_d& l, const Hyperplane_d& h) 
   : _l(l), _h(h), _known(false) {}

  Intersection_result intersection_type();
  bool intersection(Point_d& result);
  bool intersection(Line_d& result);
};

template <class R>
Line_d_Hyperplane_d_pair<R>::Intersection_result
Line_d_Hyperplane_d_pair<R>::intersection_type()
{ if (_known) return _result;
  _known = true;

  // CGAL_assertion_msg((!_l.is_degenerate())); 
  typedef typename R::Line_hyperplane_intersection_d Int_obj_type;
  Int_obj_type Intersect;
  FT lambda; 
  typename Int_obj_type::Intersection_result res = 
    Intersect(_l.point(0),_l.point(1),_h,_ip,lambda);

  if (res == Int_obj_type::LINE)  { return _result = LINE; }
  if (res == Int_obj_type::POINT) { return _result = POINT; }
  return _result = NO; 
}

template <class R>
bool
Line_d_Hyperplane_d_pair<R>::intersection(Point_d& p)
{ if (!_known) intersection_type();
  if (_result != POINT) return false;
  p = _ip; return true;
}

template <class R>
bool
Line_d_Hyperplane_d_pair<R>::intersection(Line_d& l)
{ if (!_known) intersection_type();
  if (_result != LINE) return false;
  l = _l; return true;
}

<<ray hyperplane intersection>>=
template <class R>
class Ray_d_Hyperplane_d_pair {
public:
  enum Intersection_result {NO, POINT, RAY};
  typedef typename R::FT FT;
  typedef typename R::Point_d Point_d;
  typedef typename R::Ray_d Ray_d;
  typedef typename R::Hyperplane_d Hyperplane_d;
protected:
  Ray_d _r; Hyperplane_d _h;
  bool                    _known;
  Intersection_result     _result;
  Point_d                 _ip;  
public:
  Ray_d_Hyperplane_d_pair() : _known(false) {}
  Ray_d_Hyperplane_d_pair(const Ray_d& r, const Hyperplane_d& h) 
   : _r(r), _h(h), _known(false) {}
  Intersection_result intersection_type();
  bool intersection(Point_d& result);
  bool intersection(Ray_d& result);
};

template <class R>
Ray_d_Hyperplane_d_pair<R>::Intersection_result
Ray_d_Hyperplane_d_pair<R>::intersection_type()
{ 
  if (_known) return _result;
  _known = true;

  // CGAL_assertion( !_r.is_degenerate() );
  typedef typename R::Line_hyperplane_intersection_d Int_obj_type;
  Int_obj_type Intersect;
  FT lambda; 
  typename Int_obj_type::Intersection_result res = 
    Intersect(_r.point(0),_r.point(1),_h,_ip,lambda);

  if ( res == Int_obj_type::POINT && FT(0) <= lambda ) 
  { return _result = POINT; }
  if ( res == Int_obj_type::LINE )  
  { return _result = RAY; }
  return _result = NO; 
}

template <class R>
bool
Ray_d_Hyperplane_d_pair<R>::intersection(Point_d& p)
{ if (!_known) intersection_type();
  if (_result != POINT) return false;
  p = _ip; return true;
}

template <class R>
bool
Ray_d_Hyperplane_d_pair<R>::intersection(Ray_d& r)
{ if (!_known) intersection_type();
  if (_result != RAY) return false;
  r = _r; return true;
}

<<segment hyperplane intersection>>=
template <class R>
class Segment_d_Hyperplane_d_pair {
public:
  enum Intersection_result {NO, POINT, SEGMENT};
  typedef typename R::FT           FT;
  typedef typename R::Point_d Point_d;
  typedef typename R::Segment_d Segment_d;
  typedef typename R::Hyperplane_d Hyperplane_d;

protected:
  Segment_d    _s;
  Hyperplane_d _h;
  bool                    _known;
  Intersection_result     _result;
  Point_d                 _ip;  
public:
  /*{\Mcreation 4}*/
  Segment_d_Hyperplane_d_pair() : _known(false) {}
  Segment_d_Hyperplane_d_pair(
    const Segment_d& s, const Hyperplane_d& h) :
    _s(s), _h(h), _known(false) {}

  Intersection_result intersection_type();
  bool intersection(Point_d& result);
  bool intersection(Segment_d& result);
};

template <class R>
Segment_d_Hyperplane_d_pair<R>::Intersection_result
Segment_d_Hyperplane_d_pair<R>::intersection_type()
{ 
  if (_known) return _result;
  _known = true;

  if ( _s.is_degenerate() ) 
    if ( _h.has_on(_s.point(0)) ) 
    { _ip = _s.point(0); return _result = POINT; } 
    else { return _result = NO; }

  typedef typename R::Line_hyperplane_intersection_d Int_obj_type;
  Int_obj_type Intersect;
  FT lambda; 
  typename Int_obj_type::Intersection_result res = 
    Intersect(_s.point(0),_s.point(1),_h,_ip,lambda);

  if ( res == Int_obj_type::LINE )
  { return _result = SEGMENT; }
  if ( res == Int_obj_type::POINT && 
       FT(0) <= lambda && lambda <= FT(1) ) 
  { return _result = POINT; }
  return _result = NO; 
}

template <class R>
bool
Segment_d_Hyperplane_d_pair<R>::intersection(Point_d& pt)
{ if (!_known) intersection_type();
  if (_result != POINT) return false;
  pt = _ip; return true;
}

template <class R>
bool
Segment_d_Hyperplane_d_pair<R>::intersection(Segment_d& s)
{ if (!_known) intersection_type();
  if (_result != LINE) return false;
  s = _s; return true;
}


<<intersection_objectsHd.h>>=
#ifndef CGAL_INTERSECTION_OBJECTSHD_H
#define CGAL_INTERSECTION_OBJECTSHD_H

#include <CGAL/basic.h>

namespace CGAL {

<<homogeneous intersection working horses>>

} //namespace CGAL

#include <CGAL/Kernel_d/intersection_objects_d.h>

#endif //CGAL_INTERSECTION_OBJECTSHD_H

<<intersection_objectsCd.h>>=
#ifndef CGAL_INTERSECTION_OBJECTSCD_H
#define CGAL_INTERSECTION_OBJECTSCD_H

#include <CGAL/basic.h>
#undef _DEBUG
#define _DEBUG 11
#include <CGAL/Kernel_d/debug.h>

namespace CGAL {

<<cartesian intersection working horses>>

} //namespace CGAL

#include <CGAL/Kernel_d/intersection_objects_d.h>

#endif //CGAL_INTERSECTION_OBJECTSCD_H



<<intersection_objects_d.h>>=
#ifndef CGAL_INTERSECTION_OBJECTS_D_H
#define CGAL_INTERSECTION_OBJECTS_D_H

namespace CGAL {

<<line line intersection>>
<<line ray intersection>>
<<line segment intersection>>
<<ray ray intersection>>
<<ray segment intersection>>
<<segment segment intersection>>
<<line hyperplane intersection>>
<<ray hyperplane intersection>>
<<segment hyperplane intersection>>

} //namespace CGAL
#endif //CGAL_INTERSECTION_OBJECTS_D_H

@ \end{document}