%OLD VERSION!!!
\documentclass[a4paper,11pt,twoside]{article}
\usepackage{Lweb}

\begin{document}
\title{Affine Transformations in d-Space\\
(class Aff\_transformationHd)}
\author{M. Seel}
\maketitle
\tableofcontents
\newpage

\section{The Manual Page of class Aff\_transformationHd}

\input{Aff_transformationHd.man}

\section{The Implementation of class Aff\_transformationHd}
<<Aff_transformationHd.h>>=
//---------------------------------------------------------------------
// file generated by notangle from noweb/Aff_transformationHd.lw
// please debug or modify noweb file
// coding: K. Mehlhorn, M. Seel
//---------------------------------------------------------------------

#ifndef CGAL_AFF_TRANSFORMATIONHD_H
#define CGAL_AFF_TRANSFORMATIONHD_H

#include <CGAL/basic.h>
#include <CGAL/aff_transformation_tags.h>
#include <CGAL/Handle_for.h>
#include <CGAL/Kernel_d/PointHd.h>
#include <CGAL/Kernel_d/VectorHd.h>
#include <CGAL/Kernel_d/DirectionHd.h>
#include <CGAL/Kernel_d/HyperplaneHd.h>

CGAL_BEGIN_NAMESPACE

<<defining Aff_transformationHd_rep>>
<<defining Aff_transformationHd>>

CGAL_END_NAMESPACE
#endif // CGAL_AFF_TRANSFORMATIONHD_H

@ The implementation file contains the following.
<<Aff_transformationHd.C>>=
//---------------------------------------------------------------------
// file generated by notangle from noweb/Aff_transformationHd.lw
// please debug or modify noweb file
// coding: K. Mehlhorn, M. Seel
//---------------------------------------------------------------------

#ifndef CGAL_AFF_TRANSFORMATIONHD_C
#define CGAL_AFF_TRANSFORMATIONHD_C
CGAL_BEGIN_NAMESPACE

<<implementing Aff_transformationHd>>

CGAL_END_NAMESPACE
#endif // CGAL_AFF_TRANSFORMATIONHD_C

@ \subsection{The Representation Class}
First we provide a representation class for Aff\_transformationHd.
<<defining Aff_transformationHd_rep>>=
template <class RT, class LA > class Aff_transformationHd;
template <class RT, class LA > class Aff_transformationHd_rep;

template <class RT, class LA>
class Aff_transformationHd_rep : public Ref_counted {
  friend class Aff_transformationHd<RT,LA>;
  typedef typename LA::Matrix Matrix;
  Matrix M_;
public:
  Aff_transformationHd_rep(int d) : M_(d+1) {}
  Aff_transformationHd_rep(const Matrix& M_init) : M_(M_init) {}
  ~Aff_transformationHd_rep() {}
}; 


@ \subsection{The Datatype}
And now for the class definition. We interleave the prototyping
and the implementation.
<<defining Aff_transformationHd>>=
/*{\Moptions outfile=Aff_transformation_d.man}*/ 
/*{\Manpage{Aff_transformation_d}{R}{Affine Transformations}{t}}*/
/*{\Msubst 
Hd<RT,LA>#_d<R>
Aff_transformationHd#Aff_transformation_d
Quotient<RT>#FT
}*/

template <class _RT, class _LA>
class Aff_transformationHd : 
  public Handle_for< Aff_transformationHd_rep<_RT,_LA> > { 

  typedef Aff_transformationHd_rep<_RT,_LA> Rep;
  typedef Handle_for<Rep> Base;
  typedef Aff_transformationHd<_RT,_LA> Self;

/*{\Mdefinition 
An instance of the data type |\Mname| is an affine transformation of
$d$-dimensional space. It is specified by a square matrix
$M$ of dimension $d + 1$. All entries in the last row of |M| except
the diagonal entry must be zero; the diagonal entry must be non-zero.
A point $p$ with homogeneous coordinates $(p[0], \ldots, p[d])$ can be
transformed into the point |p.transform(A)|, where |A| is an affine
transformation created from |M| by the constructors below. }*/

public: 
/*{\Mtypes 4}*/

typedef _RT RT;
/*{\Mtypemember the ring type.}*/
typedef Quotient<_RT> FT;
/*{\Mtypemember the field type.}*/
typedef _LA LA;
/*{\Mtypemember the linear algebra layer.}*/
typedef typename _LA::Matrix Matrix;
/*{\Mtypemember the matrix type.}*/
typedef typename _LA::Vector Vector;

/*{\Mcreation 3}*/

Aff_transformationHd(int d = 0) : Base( Rep(d) ) {}
/*{\Mcreate introduces a transformation in $d$-dimensional space.}*/

Aff_transformationHd(int d, Identity_transformation) : Base( Rep(d) )
/*{\Mcreate introduces the identity transformation in $d$-dimensional 
    space.}*/
{ for (int i = 0; i <= d; ++i) ptr->M_(i,i) = RT(1); }

Aff_transformationHd(const Matrix& M) : Base( Rep(M) )
/*{\Mcreate introduces the transformation of $d$ - space specified by
matrix $M$. \precond |M| is a square matrix of dimension $d + 1$. }*/
{ CGAL_assertion_msg((M.row_dimension()==M.column_dimension()),
    "Aff_transformationHd::\
     construction: initialization matrix is not quadratic.");
}

#ifndef CGAL_SIMPLE_INTERFACE

template <typename Forward_iterator>
Aff_transformationHd(Scaling, Forward_iterator start, Forward_iterator end) :
  Base( Rep(std::distance(start,end)-1) )
/*{\Mcreate introduces the transformation of $d$-space specified by a
diagonal matrix with entries |set [start,end)| on the diagonal 
(a scaling of the space). \precond |set [start,end)| is a vector of 
dimension $d+1$.}*/
{ int i=0; while (start != end) { ptr->M_(i,i) = *start++;++i; } }

#else
#define FIXATHD(I) \
Aff_transformationHd(Scaling, I start, I end) : \
  Base( Rep(end-start-1) ) \
{ int i=0; while (start != end) { ptr->M_(i,i) = *start++;++i; } }

FIXATHD(int*)
FIXATHD(const int*)
FIXATHD(RT*)
FIXATHD(const RT*)
#undef FIXATHD
#endif

Aff_transformationHd(Translation, const VectorHd<RT,LA>& v) :
  Base( Rep(v.dimension()) )
/*{\Mcreate introduces the translation by vector $v$.}*/ 
@ The homogenous representation of the rational vector is put into the
last column of the matrix and the rest of the diagonal is filled with
the homogenizing element:
\begin{displaymath}
\left( \begin{array}{ccccc}
vec.homogeneous(vec.dimension()) & 0 & \cdots & 0 & vec.homogeneous(0) \\
0 & vec.homogeneous(vec.dimension()) & & & vec.homogeneous(1) \\ 
0 & 0 & \ddots & 0 & \vdots \\
0 & \cdots & 0& 0 & vec.homogeneous(vec.dimension()) \\
       \end{array} \right)
\end{displaymath}{
<<defining Aff_transformationHd>>=
{ int d = v.dimension();
  for (int i = 0; i < d; ++i) {
    ptr->M_(i,i) = v.homogeneous(d);
    ptr->M_(i,d) = v.homogeneous(i);
  }
  ptr->M_(d,d) = v.homogeneous(d);
}

Aff_transformationHd(int d, Scaling, const RT& num, const RT& den) 
  : Base( Rep(d) ) 
/*{\Mcreate returns a scaling by a scale factor |num/den|.}*/
@ We provide a simple scaling transformation.  We put the numerator of
the rational scaling factor into the Matrix elements $M_{0,0}$ to
$M_{d-1,d-1}$ and the denominator into the lower right corner element
$M_{d,d}$: 
\begin{displaymath} 
\left( \begin{array}{cccc} num &&& \\ 
                          & \ddots && \\ 
                          && num & \\ 
                          &&& den \\ 
       \end{array} \right) 
\end{displaymath}
<<defining Aff_transformationHd>>=
{ Matrix& M = ptr->M_;
  for (int i = 0; i < d; ++i) M(i,i) = num;
  M(d,d) = den;
}

Aff_transformationHd(int d, Rotation,  
  const RT& sin_num, const RT& cos_num, const RT& den, 
  int e1 = 0, int e2 = 1); 
/*{\Mcreate returns a planar rotation with sine and cosine values
|sin_num/den| and |cos_num/den| in the plane spanned by
the base vectors $b_{e1}$ and $b_{e2}$ in $d$-space. Thus
the default use delivers a planar rotation in the $x$-$y$
plane. \precond $|sin_num|^2 + |cos_num|^2 = |den|^2$
and $0 \leq e_1 < e_2 < d$}*/

@ To implement a rational rotation we put the standard 2d rotation
matrix into the intersection points of the $e_1$-th and $e_2$-th
columns and rows and the denominator of the sine and cosine values
into the right lower corner. This looks in the default version:
\begin{displaymath}
\left( \begin{array}{ccc}
cos\_num & -sin\_num &  \\
sin\_num &  cos\_num &  \\
& & den \\
       \end{array} \right)
\end{displaymath}
<<implementing Aff_transformationHd>>=
template <class RT, class LA>
Aff_transformationHd<RT,LA>::
Aff_transformationHd(int d, Rotation,  
  const RT& sin_num, const RT& cos_num, const RT& den, 
  int e1, int e2) : Base( Rep(d) )
{
  CGAL_assertion_msg((sin_num*sin_num + cos_num*cos_num == den*den),
    "planar_rotation: rotation parameters disobey precondition.");
  CGAL_assertion_msg((0<=e1 && e1<=e2 && e2<d), 
    "planar_rotation: base vector indices wrong.");
  Matrix& M = ptr->M_;
  for (int i=0; i<d; i++) M(i,i) = 1;
  M(e1,e1) = cos_num; M(e1,e2) = -sin_num;
  M(e2,e1) = sin_num; M(e2,e2) = cos_num;
  M(d,d) = den;
}

<<defining Aff_transformationHd>>=
Aff_transformationHd(int d, Rotation, const DirectionHd<RT,LA>& dir, 
  const RT& num, const RT& den, int e1 = 0, int e2 = 1); 
/*{\Mcreate returns a planar rotation within the plane spanned by
the base vectors $b_{e1}$ and $b_{e2}$ in $d$-space.  The rotation
parameters are given by the $2$-dimensional direction |dir|, such that
the difference between the sines and cosines of the rotation given by
|dir| and the approximated rotation are at most |num/den| each.\\
\precond |dir.dimension()==2|, |!dir.is_degenerate()| and |num < den|
is positive and $0 \leq e_1 < e_2 < d$ }*/

@ Here we implement a special rotation calculation procedure starting
from a direction |dir| and a rational error bound |num/den|. We want
to find a triple $(sin,cos,denom)$ which obeys the equality $sin^2 +
cos^2 = denom^2$ and at the same time approximates the rotation given
by direction $dir$, such that the differeces between the sines and
cosines of $dir$ and the approximation are at most $num/den$.

The code is based on the rational rotation method presented by Canny
and Ressler at the 8th SCG 1992.  The approximation is based on Farey
sequences. To check the quality of the current approximation we have
to compare a rational and a (possibly) non-rational number.  The
implementation used division and modulus operation \% (the division is
always exact, that is, it is known that there is no remainder).
<<implementing Aff_transformationHd>>=
template <class RT, class LA>
Aff_transformationHd<RT,LA>::
Aff_transformationHd(int d, Rotation, const DirectionHd<RT,LA>& dir, 
   const RT& num, const RT& den, int e1, int e2) : Base( Rep(d) )
{
  CGAL_assertion_msg((dir.dimension() == 2),
    "planar_rotation: dir has to be 2 dimensional.");
  CGAL_assertion_msg((RT(0)<=num && num < den),
    "planar_rotation: num and den have to be positive.");
  CGAL_assertion_msg((0<=e1 && e1<=e2 && e2<d), 
    "planar_rotation: base vector indices wrong.");

  // now |num/den| is a rational greater zero

  RT sin;
  RT cos;
  RT denom;
  RT dx = CGAL_NTS abs(dir.dx());
  RT dy = CGAL_NTS abs(dir.dy());
  RT sq_hypotenuse = dx*dx + dy*dy;
  RT common_part;
  RT diff_part;
  RT rhs;
  bool lower_ok;
  bool upper_ok;

  if (dy > dx) {
    RT tmp = dx;
    dx = dy;
    dy = tmp;
  }
  /* approximate |sin = dy / sqrt(sq_hypotenuse)| \\
     |if ( dy / sqrt(sq_hypotenuse) < num/den )| */
  if (dy * dy * den * den < sq_hypotenuse * num * num) {
     cos = denom = 1;
     sin = 0;
  } else {
    RT p,q,p0,q0,p1,q1;
    p0 = 0;
    q0 = p1 = q1 = 1;

    for(;;) {
      p = p0 + p1;
      q = q0 + q1;
      sin = RT(2)*p*q;
      denom = CGAL_NTS square(p) + CGAL_NTS square(q);

      // sanity check for approximation
      // |       sin/denom < dy/sqrt(hypotenuse) + num/den|
      // |   &&  sin/denom > dy/sqrt(hypotenuse) - num/den|
      // | ==    sin/denom - num/den  <   dy/sqrt(sq_hypotenuse)|
      // |   &&  sin/denom + num/den  >   dy/sqrt(sq_hypotenuse)|
      // | ==    (sqr(sin) sqr(den) + sqr(num) sqr(denom)) sq_hypotenuse - 2..|
      // |       < sqr(dy) sqr(den) sqr(denom)|
      // |   &&  (sqr(sin) sqr(den) + sqr(num) sqr(denom))sq_hypotenuse + 2..|
      // |       > sqr(dy) sqr(den) sqr(denom)|

      common_part = (CGAL_NTS square(sin)*CGAL_NTS square(den) + 
                     CGAL_NTS square(num)*CGAL_NTS square(denom))*
                     sq_hypotenuse;
      diff_part   = RT(2)*num*sin*den*denom*sq_hypotenuse;
      rhs         = CGAL_NTS square(dy)*CGAL_NTS square(den)*
                    CGAL_NTS square(denom);

      upper_ok    = (common_part - diff_part < rhs);
      lower_ok    = (common_part + diff_part > rhs);

      if ( lower_ok && upper_ok ) {
         if ( CGAL_NTS square(p)%RT(2) + CGAL_NTS square(q)%RT(2) > 1) {
             sin = p*q;
             cos = (CGAL_NTS square(q) - CGAL_NTS square(p))/RT(2);     
             // exact division
             denom = (CGAL_NTS square(p) + CGAL_NTS square(q))/RT(2);   
             // exact division
         } else {
             cos = CGAL_NTS square(q) - CGAL_NTS square(p);
         }

         break;
      } else {
          /* |if ( dy/sqrt(sq_hypotenuse) < sin/denom )| */
          if ( CGAL_NTS square(dy)*CGAL_NTS square(denom) < 
               CGAL_NTS square(sin)*sq_hypotenuse )
          { p1 = p; q1 = q; }
          else
          { p0 = p; q0 = q; }
      }
    } // for(;;)
  }

  dx = dir.dx();
  dy = dir.dy();
  if (dy > dx) {
    RT tmp = dx;
    dx = dy;
    dy = tmp;
  }
  if (dx < 0) sin = - sin;
  if (dy < 0) cos = - cos;
  Matrix& M = ptr->M_;
  for (int i=0; i<d; i++) M(i,i) = 1;
  M(e1,e1) = cos; M(e1,e2) = -sin;
  M(e2,e1) = sin; M(e2,e2) = cos;
  M(d,d) = denom;
}

<<defining Aff_transformationHd>>=
/*{\Moperations 5 3}*/

int dimension() const 
{ return ptr->M_.row_dimension()-1; }
/*{\Mop the dimension of the underlying space }*/

const Matrix& matrix() const { return ptr->M_; }
/*{\Mop returns the transformation matrix }*/

Vector operator()(const Vector& iv) const
// transforms the ivector by a matrix multiplication
{ return matrix()*iv; }


Aff_transformationHd<RT,LA> inverse() const
/*{\Mop returns the inverse transformation.
\precond |\Mvar.matrix()| is invertible.}*/
{ Aff_transformationHd<RT,LA> Inv; RT D; 
  Vector dummy;
  if ( !LA::inverse(matrix(),Inv.ptr->M_,D,dummy) ) 
  CGAL_assertion_msg(0,"Aff_transformationHd::inverse: not invertible."); 
  return Inv;
}
  
Aff_transformationHd<RT,LA>  
operator*(const Aff_transformationHd<RT,LA>& s) const
/*{\Mbinop composition of transformations. Note that transformations
are not necessarily commutative. |t*s| is the transformation
which transforms first by |t| and then by |s|.}*/
{ CGAL_assertion_msg((dimension()==s.dimension()),
  "Aff_transformationHd::operator*: dimensions disagree.");
  return Aff_transformationHd<RT,LA>(matrix()*s.matrix()); 
}

bool operator==(const Aff_transformationHd<RT,LA>& a1) const
{ if ( identical(a1) ) return true;
  return ( matrix() == a1.matrix() );
}
bool operator!=(const Aff_transformationHd<RT,LA>& a1) const
{ return !operator==(a1); }

}; // Aff_transformationHd

template <class RT, class LA>
std::ostream& operator<<(
  std::ostream& os, const Aff_transformationHd<RT,LA>& t) 
{ os << t.matrix(); return os; }

template <class RT, class LA>
std::istream& operator>>(
  std::istream& is, Aff_transformationHd<RT,LA>& t)
{ typename LA::Matrix M(t.dimension());
  is >> M; t = Aff_transformationHd<RT,LA>(M); 
  return is;
}

/*{\Mimplementation 
Affine Transformations are implemented by matrices of integers as an
item type.  All operations like creation, initialization, input and
output on a transformation $t$ take time $O(|t.dimension()|^2)$. |dimension()|
takes constant time.  The operations for inversion and composition
have the cubic costs of the used matrix operations. The space
requirement is $O(|t.dimension()|^2)$. }*/

// ----------------------------- end of file ----------------------------

@ 
\end{document}