#ifndef GEOMETRY_HH
#define GEOMETRY_HH

// Oriented projective geometry classes
// Copyright (C) 2002  Bruce J. Bell
//
// Intro:
//   for context, see Jorge Stolfi's book or thesis --
//     book:   _Oriented Projective Geometry_
//     thesis: "Primitives for Computational Geometry"
//
//   generally, T^N is the  N-d "oriented projective space":
//     concretely, a double embedding ofthe N-d projective space P^N
//     in the usual homogeneous (N+1)-d vector space, that distinguishes
//     between positive and negative multiples of the homogeneous vector
//     in a systematic and useful way.
//
//   Plu:cker/Grassman-type coordinates are (so far...) only used for
//     3-D lines (class line_T3) [see below]
//
//
//
// General notation summary:
//   ===  means "is defined as"
//   ->  <-  and  <->  are arrows, used informally (e.g., for assignment)
//   =    in commentary is mathematical equality (not assignment)
//   ~    means "is equivalent to";  specifically, if not otherwise stated,
//          it means oriented projective equivalence:
//                  x ~ y  <->  x = s*y  [for some positive scalar s]
//   ~=   means "approximately equal", informally
//   ^    can be used as an exponentiation operator;
//          to distinguish between meet and exponential operators,
//            the meet should have spaces on either side and
//            the exponential should not have spaces on either side.
//          to distinguish between superscripts and exponentials,
//            the exponential should be immediately followed by a
//            parenthetical expression, while the superscript
//            should not be.
//          examples:
//            X^(r-s)  -> exponential (followed by "(", preceded by non-space)
//            X ^ Y  -> meet (preceded and followed by " ")
//            Z_ij^mn c_i d^j  -> superscript (followed by index index name)
//   ...  is a general indication of a point that perhaps needs more work.
//          It can also be used in code to represent unwritten portions --
//          it will cause a syntax error if compiled as code, and it can
//          also be easily searched for.
//
//
// Generic names:
//   s generally indicates a *positive* scalar
//   b generally indicates a "base" vector
//   p generally indicates a point
//   d generally indicates a dual-of-point
//   M generally indicates a projective transform
//   D generally indicates a duomorphism (mapping objects to their dual space)
//   L generally indicates a 3-D line, or linelike object
//     (e.g. antisymmetric 2-form in 3-D)
//   T generally indicates a matrix or other unconstrained order-2 object
//
//   otherwise unspecified lowercase letters generally represent
//     scalars, vectors, points, dual points, or other objects of order < 2
//   otherwise unspecified uppercase letters generally represent
//     matrices, maps, lines, or other objects of order >= 2
//
//   rep(x) is the concrete representation of geometric object x
//   sgn(s) is the sign function: 1, 0, or -1 according to the sign of s
//
// Geometric notation:
//   sign/orientation conventions are the same as in Stolfi's system
//   the "meet" operation is written:  A ^ B   [spaces necessary]
//   the "join" operation is written:  A v B   [spaces necessary]
//     thus the variable v should  be left out of commentary
//     (in code v is used for internal representation data;
//     this likewise means that function argument names should avoid v)
//
// Matrix notation:
//   T^-1 is T-inverse;  T^t is T-transpose  (similarly for vectors)
//   |T| is determinant of T (any matrix norms should be written as functions)
//   product of matrix M and row vector r = r M
//   product of matrix M and column vector c = M c
//   points are represented by homogeneous row vectors
//   duals of points are represented by their dual point, also a row vector
//   however, when making a product, duals should be transposed and used as
//     a column vector --  e.g., geometric product of a point p and dual d:
//     rep(p v d) = rep(p) rep(d)^t
//
// Tensor notation:
//   subscripts are indicated as T_ij
//   superscripts are indicated as T^ij  (try to avoid t as index...)
//   mixed indices are indicated like T_i^j or T^i_j
//   points should have subscripts and their duals superscripts
//   when binding unrelated spaces, covariant & contravariant don't make sense;
//     in such a case, the indices should all be subscripts.
//
//
// Implementation notes:
//   to take advantage of duality, operations may be written
//     in terms of "base" classes (e.g., base_T2)
//     and coerced or translated to geometric interfaces as needed
//   "basic" operation, with monomials of degree <=2, may not be exact,
//     but should be as (anti)symmetric as conveniently possible
//   "compound" operations, with monomials of degree >2, should be
//     implemented in terms of the "basic" operations.  This avoids
//     placing undue strain on the accumulator's dynamic range, but
//     may yield less precise results, and possibly perturb
//     the (anti)symmetry present in "basic" operations.
//   for this reason, "basic" operations should be preferred to
//     "compound" operations where feasible.


#include "linear.hh"
#include <math.h>


// class summary:


// basic geometric objects --

// 1-D projective objects (2-vector homogeneous space) --
class point_T1; // point on 1-D projective line (self-dual)

// 2-D projective objects (3-vector homogeneous space) --
class base_T2;    // base representation used for 2-D point, line
class point_T2;   // point in 2-D projective plane
class line_T2;    // line in 2-D projective plane

// 3-D projective objects (4-vector homogeneous space) --
class base_T3;    // base representation used for 3-D point, plane
class point_T3;   // point in 3-D projective space
class plane_T3;   // plane in 3-D projective space
class line_T3;    // line in 3-D projective space (self-dual)

// 5-D projective objects (6-vector homogeneous space)
class base_T5;    // base representation used for line_T3

// linear transforms --

// base transforms
class matrix_T2;  // transform matrix: base_T2 -> base_T2
class matrix_T3;  // transform matrix: base_T3 -> base_T3
class matrix_T5;  // transform matrix: line_T3 -> line_T3 (specialized)
// (don't do symmetric representations yet -- it's just an optimization
// -- for now, use full matrix for metric tensors and such...)

// geometric transforms, from geometric space to itself --
//   to represent M, note that rep(p), rep(d) are both row vectors:
//     rep(p v d) = rep(p) rep(d)^t
//
// define A,B such that:
//   rep(M(point)) = rep(point) A
//   rep(M(dual)) = rep(dual) B
// then space orientation consistency (see Stolfi) requires:
//   M(point) v M(dual) = M(point v dual)  ~  (point v dual)*|A|
//   ->  point A B^t dual^t  ~  |A|*(point dual^t)
//   ->  A B^t ~ |A|* unit  ->  B^t ~ |A| A^-1
//   ->  B ~ |A| (A^-1)^t === cofactor(A)
//
// inverse map of (A, B) = (A^-1, B^-1) = (B^t/|A|, A^t*|A|)
//   ->  (A^-1, B^-1) ~ (sgn|A|*B^t, sgn|A|*A^t)
//
// note that orientation-reversing maps ( |A| < 0 )
//   will reverse the sign of (many/most/all?...) scalar results
//   from basic element operations, since they (formally) map
//   the positive space onto the negative space.
//
// 3-D line mapping may be done 3 ways:
//  - select 2 distinct points on line, map them, then reconstitute line
//  - select 2 distinct planes through lie, map them, then reconstitute line
//  - expand matrix to 6x6, and map 6-element line coordinates directly
//
class gmap_T2; // geometric projective map, T2 -> T2 (point_T2, line_T2, etc)
class gmap_T3; // geometric projective map, T3 -> T3 (point_T3, plane_T3, etc)

// "Duomorphic" transforms
//   given duomorphism D, let A,B be such that:
//     rep(D(point)) = rep(point) A   [D(point) is a dual]
//     rep(D(dual))  = rep(dual) B    [D(dual) is a point]
//
//   duomorphism consistency requirement is:
//     D("vacuum" rank-0 element) = "reference space" highest-rank element
//     ->  rep(D(point ^ dual)) ~ rep(point v dual) = rep(point) rep(dual)^t
//   by definition of D:
//     D(point ^ dual) = D(point) v D(dual)
//   so:  rep(point) rep(dual)^t ~ rep(D(point ^ dual))
//       = rep(D(point) v D(dual))
//       = rep( (-1)^(rank-1) * D(dual) v D(point) )
//       = (-1)^(rank-1) * rep(D(dual)) rep(D(dual))^t
//       = (-1)^(rank-1) * rep(dual) B  A^t rep(point)^t
//     ->  B A^t  ~  (-1)^(rank-1) * unit ->  B  ~  (-1)^(rank-1) * (A^-1)^t
//
//   inverse map:
//     let D' = inverse(D)
//         rep(D'(p)) = rep(p) A'  (D'(p) is a dual)
//         rep(D'(d)) = rep(d) B'  (D'(d) is a point)
//     then p ~ rep(D'(D(p))) = rep(D(p)) B' = p A B'
//       -> B' ~ A^-1 ~ (-1)^(rank-1) * B^t
//     and d ~ rep(D'(D(d))) = rep(D(d)) A' = d B A'
//       -> A' ~ B^-1 ~ (-1)^(rank-1) * A^t
//
//   in particular, for T^2:  dim=2, rank=3 so B ~ (A^-1)^t ~ cofactor(A)/|A|
//                            and for D^-1=D', (A', B') ~ (A^t, B^t)
//                  for T^3:  dim=3, rank=4 so B ~ -(A^-1)^t ~ -cofactor(A)/|A|
//                            and for D^-1=D', (A', B') ~ (-A^t, -B^t)
//
// 3-D line mapping... (work out later)
//
class dmap_T2;  // symmetric dual map (point_T2 <-> line_T2, etc)
class dmap_T3;  // symmetric dual map (point_T3 <-> plane_T3, etc)



// class definitions:


// basic geometric objects (e.g., vectors)

class point_T1 {
  // 1-D homogeneous point
  // self-dual element
  //
  // (since there are no other interesting 1-D objects,
  // there is no need for a separate "base" class)
  //
  // besides representing the 1-D homogeneous line,
  // point_T1 is useful for:
  //     periodic values like angles,
  //     extending the dynamic range of scalars
  //     (note that, while extending dynamic range to 2*scalar.bits,
  //     it does not generally increase resolution significantly)
  //
  // usage note: 1-D homogeneous line, like other
  //   homogeneous representations, can often simplify
  //   algorithms by eliminating special cases and
  //   unnecessary singularities, and/or implicitly
  //   storing orientation flag as well as analog value

private: // internal data
  //scalar v[2];  // homogeneous representation
  vect<2> v;  // homogeneous representation

  // name convention:
  // v[0] = w
  // v[1] = x
  //
  // affine line convention:
  //   val = x/w   (conventional real value normalizes w to 1)
  //   1/val = w/x
  //   orientation = sgn(w)  [w>0 -> point on front range]
  //
  // angle convention:
  //   w = K cos(angle)
  //   x = K sin(angle)
  //     (where K is some positive real number)
  //   -> angle = atan2(x,w) = atan(x/w) + N pi = atan(val) + N pi
  //     (where N is one of {-1,0,1}, consistent with above K being positive)

public: // ctors, dtor, etc.
  point_T1(void)  {} // use default scalar initialization
  point_T1(scalar w, scalar x) { set_point(w,x); }
  // use default copy ctor
  // use default dtor
  // use default assignment operator

public: // access methods
  scalar get_val(int i) const;
  scalar operator[](int i) const { return get_val(i);  }
    // note: there is no get_ref or non-const operator[]
    // (design choice -- if you want one, you should probably rethink why...)
  bool isnull(void) const;

public: // operations returning point_T1 (update/replace *this)
  void set_point(const point_T1& a) { *this= a; }
  void set_point(scalar w, scalar x);  // *this  <- [w,x]
  void set_null(void);  // *this  <-  (0,0)
  void set_anti(const point_T1& a);  // *this  <-  -a

  // right and left duals are antipodal to each other
  void set_dual_right(const point_T1& a);  // *this  <-  right-dual(a)
  void set_dual_left(const point_T1& a);  // *this  <-  left-dual(a)

public: // operations returning scalars
  scalar dot(const point_T1& a) const;  // conventional dot product
    // horsepower:  requires 2*, 1+
  scalar det(const point_T1& a) const;
    // returns 2x2 determinant  |w  a.w|  =  |*this, a|
    //                          |x  a.x|
    // tensor notation:  epsilon^ij (*this)_i a_j
    //
    // horsepower:  requires 2*, 1+

public: // projective operations returning scalars
  // concrete justification/mnemonic in terms of Stolfi's "straight" model:
  //   If (p,q) is a positive simplex on the front span of the
  //     representation space, 0 < rep(p v q)  and  (p_w>0, q_w>0)
  //   also,  val(p) < val(q)  ->  p_x/p_w < q_x/q_w
  //   ->  p_x*q_w < q_x*p_w  ->  0 < p_w*q_x - p_x*q_w
  // Thus,  sgn( rep(p v q) ) = sgn( |rep(p),rep(q)| )
  scalar eval_point(const point_T1& a) const { return det(a); }
  scalar operator*(const point_T1& a) const { return eval_point(a); }
  scalar operator()(const point_T1& a) const { return eval_point(a); }

public: // affine and other operations
  bool isfront(void) const  { return get_val(0).ispos(); }
  bool isback(void) const  { return get_val(0).isneg(); }
  bool isideal(void) const  { return get_val(0).iszero(); }

  double fpval(void) const;  // returns x/w
  double fpinv(void) const;  // returns w/x

  double fpangle(void) const;  // returns atan2(x,w) = atan(x/w) + N pi

public: // debug operations
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;

}; // end class point_T1


class base_T2 {
  // 2-D homogeneous point or plane
  friend class matrix_T2;

private: // internal data
  //scalar v[3];  // homogeneous representation
  vect<3> v;

  // name convention:
  // v[0] = w
  // v[1] = x
  // v[2] = y

public: // ctors, dtor, set methods
  base_T2(void) {} // use default scalar initialization
  base_T2(scalar w, scalar x, scalar y) { set_base(w,x,y); }
  // use default copy ctor
  // use default dtor
  // use default assignment operator

public: // access methods
  scalar get_val(int i) const;
  scalar operator[](int i) const { return get_val(i); }
    // note: there is no get_ref or non-const operator[]
    // (design choice -- if you want one, you should probably rethink why...)
  bool isnull(void) const;

public: // operations returning base_T2 (update/replace *this)
  void set_base(const base_T2& a) { *this= a; }
  void set_base(scalar w, scalar x, scalar y);  // *this <- [w,x,y]
  void set_null(void);  // *this  <-  (0,0,0)
  void set_anti(const base_T2& a);  // *this  <-  -a

  void set_cross(const base_T2& a, const base_T2& b);  // *this  <-  cross(a,b)
    // (*this)^k = epsilon^ijk a_i b_j   [use result as dual of argument type]
    // horsepower:  requires 6*, 3+

public: // operations returning scalars
  scalar dot(const  base_T2& a) const;  // conventional dot product
    // horsepower:  requires 3*, 2+

  scalar det(const base_T2& b, const base_T2& c) const;
    // returns 3x3 determinant |w b.w c.w| = |(*this), b, c|
    //                         |x b.x c.x|
    //                         |y b.y c.y|
    // tensor notation:  epsilon^ijk (*this)_i b_j c_k
    //
    // use triple product === dot(a, cross(b,c))
    // horsepower:  requires 9*, 5+

public: // debug operations
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;
}; // end class base_T2


class point_T2: public base_T2 {
  // 2-D homogeneous point

  // affine plane convention:
  //   coords (X, Y) = (x/w, y/w)
  //   orientation = sgn(w)  [w>0 -> point on front range]

public: // ctors, etc.
  point_T2(void) {}
  point_T2(scalar w, scalar x, scalar y): base_T2(w,x,y) {}

public: // projective operations returning point_T2 (update/replace *this)
  void set_point(const point_T2& p) { *this= p; }
  void set_point(scalar w, scalar x, scalar y) { set_base(w,x,y); }
  // set_null inherited from base_T2
  void set_anti(const point_T2& p) { base_T2::set_anti((const base_T2&)p); }

  // right and left duals of 2-D elements are the same
  void set_dual(const line_T2& d)  { set_base((const base_T2&)d); }

  void set_meet_lines(const line_T2& d, const line_T2& e)
    // *this  <-  oriented intersection of two 2-D lines === d ^ e
    { set_cross((const base_T2&)d, (const base_T2&)e); }

public: // projective operations returning scalar
  scalar eval_line(const line_T2& d) const { return dot((base_T2&)d); }
    // returns rep(*this v d)
  scalar operator*(const line_T2& d) const { return eval_line(d); }
  scalar operator()(const line_T2& d) const { return eval_line(d); }

  scalar orient(const point_T2& p, const point_T2& q) const
    // returns measure of point-triangle (*this, p, q)
    { return det(p,q); }

public: // affine methods
  bool isfront(void) const  { return get_val(0).ispos(); }
  bool isback(void) const  { return get_val(0).isneg(); }
  bool isideal(void) const  { return get_val(0).iszero(); }

  double fpX(void) const;  // returns x/w
  double fpXinv(void) const;  // returns w/x
  double fpY(void) const;  // returns y/w
  double fpYinv(void) const;  // returns w/y

public: // debug operations
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;

}; // end class point_T2


class line_T2: public base_T2 {
  // 2-D homogeneous line
public: // ctors, etc.
  line_T2(void) {}
  line_T2(scalar w, scalar x, scalar y): base_T2(w,x,y) {}

public: // projective operations returning line_T2 (update/replace *this)
  void set_line(const line_T2& d) { *this= d; }
  void set_line(scalar w, scalar x, scalar y) { set_base(w,x,y); }
  // set_null inherited from base_T2
  void set_anti(const line_T2& d)  { base_T2::set_anti((const base_T2&)d); }

  // right and left duals of 2-D elements are the same
  void set_dual(const point_T2& p) { set_base((const base_T2&)p); }

  void set_join_points(const point_T2& p, const point_T2& q)
    // *this  <-  oriented line through two 2-D points === p v q
    { set_cross((const base_T2&)p,(const base_T2&)q); }

public: // projective operations returning scalar
  scalar eval_point(const point_T2& p) const { return dot(p); }
    // returns rep(*this v p)
  scalar operator*(const point_T2& p) const { return eval_point(p); }
  scalar operator()(const point_T2& p) const { return eval_point(p); }

  scalar orient(const line_T2& d, const line_T2& e) const
    // returns measure of line-triangle (*this, d, e)
    { return det(d,e); }

public: // affine methods
  void get_dir(point_T2& p_out) const {
    // p_out  <-  "direction" of line *this === oriented point-at-infinity
    // p_out  <-  *this ^ [1,0,0] === cross(*this,  [1,0,0])
    // (p_out)_k  <-  epsilon_ijk (*this)^i [1,0,0]^j
    //
    // specifically:  p_out = [0,-y, x]
    p_out.set_point(0,-get_val(2), get_val(1));
  }
  bool isideal(void) const
    { return get_val(1).iszero() && get_val(2).iszero(); }

public: // debug operations
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;

}; // end class line_T2


class base_T3 {
  // 3-D homogeneous point or plane

  friend class matrix_T3;
  friend class matrix_T5;

private: // internal data
  //scalar v[4];  // homogeneous representation
  vect<4> v;  // homogeneous representation

  // name convention:
  // v[0]= w
  // v[1]= x
  // v[2]= y
  // v[3]= z

public: // ctors, dtor, set methods
  base_T3(void) {} // use default scalar initialization
  base_T3(scalar w, scalar x, scalar y, scalar z) { set_base(w,x,y,z); }
  // use default copy ctor
  // use default dtor
  // use default assignment operator

public: // access methods
  scalar get_val(int i) const;
  scalar operator[](int i) const { return get_val(i); }
    // note: there is no get_ref or non-const operator[]
    // (design choice -- if you want one, you should probably rethink why...)
  bool isnull(void) const;

public: // base operations returning base_T3 (update/replace *this)
  void set_base(const base_T3& a) { *this= a; }
  void set_base(scalar w, scalar x, scalar y, scalar z);  // *this <- [w,x,y,z]
  void set_null(void);  // *this  <-  (0,0,0,0)
  void set_anti(const base_T3& a);  // *this  <-  -a

  // 4-D "cross product" cross(a,b,c)^m === epsilon^ijkm a_i b_j c_k
  // is a compound operation -- use nuller(a,b).map(c)
  // (if commonly needed, go ahead and implement it here...)

public:  // base operations returning scalars
  scalar dot(const base_T3& b) const;  // conventional dot product
    // horsepower:  requires 4*, 3+

  scalar det(const base_T3& b, const base_T3& c, const base_T3& d) const;
    // returns 4x4 determinant |w b.w c.w d.w| = |(*this),b,c,d|
    //                         |x b.x c.x d.x|
    //                         |y b.y c.y d.y|
    //                         |z b.z c.z d.z|
    // tensor notation:  epsilon^ijmn (*this)_i b_j c_m d_n
    //
    // use quadruple product === orient(line(*this,b), line(c,d))
    // horsepower:  uses 30*, 17+   (2 line_T3::set_dyad, 1 line_T3::orient)

public: // debug operations
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;
}; // end class base_T3


class point_T3: public base_T3 {
  // 3-D homogeneous point

  // affine space convention:
  //   coords (X, Y, Z) = (x/w, y/w, z/w)
  //   orientation = sgn(w)  [w>0 -> point on front range]

public: // ctors, etc.
  point_T3(void) {}
  point_T3(scalar w, scalar x, scalar y, scalar z): base_T3(w,x,y,z) {}

public: // projective operations returning point_T3 (update/replace *this)
  void set_point(const point_T3& p) { *this= p; }
  void set_point(scalar w, scalar x, scalar y, scalar z) { set_base(w,x,y,z); }
  // set_null inherited from base_T3
  void set_anti(const point_T3& p) { base_T3::set_anti((const base_T3&)p); }

  // right and left duals are antipodes of each other:
  void set_dual_right(const plane_T3& d)
    // *this  <-  right dual of 3-D plane d
    // such that rep(d v (*this)) > 0
    { base_T3::set_anti((const base_T3&)d); }
  void set_dual_left(const plane_T3& d)
    // *this  <-  left dual of 3-D plane d
    // such that rep((*this) v a) > 0
    { set_base((const base_T3&)d); }

public: // projective operations returning scalar
  scalar eval_plane(const plane_T3& d) const { return dot((base_T3&)d); }
    // returns rep(*this v d)
    // in 3-D, point*plane is anti_commutative:  plane*point = -point*plane
  scalar operator*(const plane_T3& a) const { return eval_plane(a); }
  scalar operator()(const plane_T3& a) const { return eval_plane(a); }

  scalar orient(const point_T3& p, const point_T3& q, const point_T3& r) const
    // returns measure of point-tetrahedron (*this,p,q,r)
    { return det(p,q,r); }

public: // affine methods
  bool isfront(void) const  { return get_val(0).ispos(); }
  bool isback(void) const  { return get_val(0).isneg(); }
  bool isideal(void) const  { return get_val(0).iszero(); }

  double fpX(void) const;  // returns x/w
  double fpXinv(void) const;  // returns w/x
  double fpY(void) const;  // returns y/w
  double fpYinv(void) const;  // returns w/y
  double fpZ(void) const;  // returns z/w
  double fpZinv(void) const;  // returns w/z

public: // debug operations
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;

}; // end class point_T3

class plane_T3: public base_T3 {
  // 3-D homogeneous plane
  //
  // rep(plane) === rep(left-dual of plane)
public: // ctors, etc.
  plane_T3(void) {}
  plane_T3(scalar w, scalar x, scalar y, scalar z): base_T3(w,x,y,z) {}

public: // projective operations returning plane_T3 (update/replace *this)
  void set_plane(const plane_T3& d) { *this= d; }
  void set_plane(scalar w, scalar x, scalar y, scalar z) {set_base(w,x,y,z); }
  // set_null inherited from base_T3
  void set_anti(const plane_T3& d) { base_T3::set_anti((const base_T3&)d); }

  // right and left duals are antipodes of each other
  void set_dual_right(const point_T3& p)
    // *this  <-  right dual of 3-D point p
    // such that rep(p v *this) > 0
    { set_base((const base_T3&)p); }
  void set_dual_left(const point_T3& p)
    // *this  <-  left dual of 3-D point p
    // such that rep(*this v p) > 0
    { base_T3::set_anti((const base_T3&)p); }

  // join of 3 points to determine 3-D plane is a compound operation
  // use 4-D "cross product" to find intersection, as above
  // (if commonly needed, go ahead and implement it...)

public: // projective operations returning scalar
  scalar eval_point(const point_T3& p) const {return -dot((const base_T3&)p); }
    // returns rep(*this v p)
    // in 3-D, plane*point is _anti_commutative:  point*plane = -plane*point
  scalar operator*(const point_T3& a) const { return eval_point(a); }
  scalar operator()(const point_T3& a) const { return eval_point(a); }

  scalar orient(const plane_T3& b, const plane_T3& c, const plane_T3& d) const
    // returns measure of plane-tetrahedron (*this,b,c,d)
    { return det(b,c,d); }

public: // affine methods
  void get_dir(line_T3& L_out) const;
    // L_out  <-  "direction" of plane *this === oriented line-at-infinity
    // implementation depends on line representation, so delegate to line_T3...
  bool isideal(void) const
  { return get_val(1).iszero() && get_val(2).iszero() && get_val(3).iszero(); }

public: // debug operations
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;

}; // end class plane_T3


class base_T5 {
  // 5-D homogeneous base element
  // (used for Plu:cker coordinates of 3-D line)
  friend class matrix_T5;

private: // internal data
  //scalar v[6]; // modified Plu:cker coordinates
  vect<6> v;  // modified Plu:cker coordinates

public: // ctors, dtor, etc.
  base_T5(void) {}
  base_T5(scalar a, scalar b, scalar c, scalar d, scalar e, scalar f)
    { set_base(a,b,c,d,e,f); }
  // use default copy ctor
  // use default dtor
  // use default assignment operator

protected: // access methods
  scalar& get_ref(int i);
public: // access methods
  scalar get_val(int i) const;
  scalar operator[](int i) const { return get_val(i); }
    // note: there is no get_ref or non-const operator[]
    // (design choice -- if you want one, you should probably rethink why...)
  bool isnull(void) const;

public: // decision-making and Monte Carlo operations
  int find_big_elt(void) const;
    // return index (0-5) of an element with largest absolute value

public: // base operations returning base_T5 (update_replace *this)
  void set_base(const base_T5& x) { *this= x; }
  void set_base(scalar a, scalar b, scalar c, scalar d, scalar e, scalar f);
  void set_null(void);  // *this  <-  (0,0,0,0,0,0)
  void set_anti(const base_T5& b);  // *this  <-  -b

public:  // base operations returning scalars
  scalar dot(const base_T5& b) const;  // conventional dot product
  // horsepower:  requires 6*, 5+

  // add determinant here if needed...

public: // debug operations
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;

}; // end class base_T5

class line_T3:  public base_T5 {
  // 3-D homogeneous line class
  // self-dual element
  //
  //   sign/order convention is neither Plu:cker's or Stolfi's;
  //   consider line as antisymmetric 2-form binding homogeneous points:
  //   scalar result is L^ij p_i q_j
  //
  //   Antisymmetry means L^ij = - L^ji  ->  6 independent nonzero elements
  //   (contex is homogeneous, so class can be derived from base_T5)
  //   Representation is upper matrix elements in lexicographic order
  //   wrt (column,row) [not wrt bitset as in Stolfi]
  //
  // (When considering L^ij as a tensor, note that both point arguments are
  //   covariant, so the rows/columns choice below is arbitrary...)
  //
  //  L = [ 0    L^01 L^02 L^03 ] === [ 0  a  b  c ]
  //      [ L^10 0    L^12 L^13 ]     [-a  0  d  e ]
  //      [ L^20 L^21 0    L^23 ]     [-b -d  0  f ]
  //      [ L^30 L^31 L^32 0    ]     [-c -e -f  0 ]
  //
  //  Since this is the only real use of Plu:cker coordinates, I feel no
  //  need to be systematic, and I find this order easy to remember
  //  (even so, this kind of ordering should extend to higher dimensions --
  //  I have no idea if it would be better or worse than other orderings...)
  //  Again, considering L^ij as a tensor, with epsilon being the 4-D
  //  fully antisymmetric tensor (AKA determinant tensor), the dualizing
  //  operation is:
  //
  //  (L*)_mn === (1/2) epsilon_ijmn L^ij = [ 0  f -e  d ]
  //                                        [-f  0  c -b ]
  //                                        [ e -c  0  a ]
  //                                        [-d  b -a  0 ]
  //
  //  As in Stolfi's ordering, the 3-D line dual operation is self-inverse:
  //    if representation of L is ( a, b, c, d, e, f),
  //       representation of L* = ( f,-e, d, c,-b, a)

  // name convention:
  // v[0] = L^01 (= a above)
  // v[1] = L^02 (= b above)
  // v[2] = L^03 (= c above)
  // v[3] = L^12 (= d above)
  // v[4] = L^13 (= e above)
  // v[5] = L^23 (= f above)

public: // private static functions
  // in keeping with notation above, first_index(a) < second_col(a)
  static int first_index(int i);
    // return base_T3 first-index of element i in line *this 
  static int second_index(int i);
    // return base_T3 second-index of element i in line *this

public: // ctors, dtors, set methods
  line_T3(void)  {} // use default scalar initialization
  line_T3(scalar a, scalar b, scalar c, scalar d, scalar e, scalar f)
    { set_line(a,b,c,d,e,f); }
  // use default copy ctor
  // use default dtor
  // use default assignment operator

public: // access methods
  // get_val() inherited from base_T5
  // operator[]() inherited from base_T5
  // isnull() inherited from base_T5


public: // operations returning line_T3 (update/replace *this)
  void set_line(scalar a, scalar b, scalar c, scalar d, scalar e, scalar f);
  void set_line(const line_T3& L) { *this= L; }
  void set_anti(const line_T3& L) { base_T5::set_anti((const base_T5&)L); }

  // for 3-D lines, dual is self-inverse
  void set_dual(const line_T3& L);  // *this  <-  dual of 3-D line L
    // such that rep(*this v L) = rep(L v *this) > 0
    // geometrically:  (*this)^ab  <-  metric^ma metric^nb epsilon_ijmn L^ij
    // alternately:  (*this)_ij  <-  epsilon_ijmn L^ij   [breaks convention]
    //
    // other notes:
    //   for all linelike elements J,K:  orient(J, K) === dot(dual(J), K)
    //   for proper lines L:  orient(L,L) = 0

  void set_dyad(const base_T3& a, const base_T3& b);
    // (*this)^ij  <-  a^i b^j - b^i a^j
    // === dyadic 2-form: antisymmetric outer product of (a,b)
    // (result tensor should bind dual of argument type)
    // horsepower:  requires 12*, 6+
  void set_nuller(const base_T3& a, const base_T3& b);
    // (*this)^mn  <-  epsilon^ijmn a_i b_j   [result nulls a and b]
    // === derive antisymmetric 2-form using 4-D determinant tensor
    // (result tensor should bind argument type)
    // horsepower:  requires 12*, 6+

public: // operations returning base_T3
  void map_vect(const base_T3& b_in, base_T3& b_out) const;
    // (b_out)^i  <-  (*this)^ij (b_in)_j   [output should be dual of arg type]
    // horsepower:  requires 12*, 8+

public: // operations that involve selecting basis-related vectors
  // based on vector index [0..3] or pivot index [0..5]
  // (for best results, pivot element should have a large absolute value
  // -- find_big_elt() is previously defined)
  void get_vect(int r, base_T3& b_out) const;
    // - select 1 row of L, determined by vector index r
    // (b_out)_i  <-  (*this)_ir   [for respective column, negate result...]
  void get_vects(int pivot, base_T3& b_out, base_T3& c_out) const;
    // - select 2 rows of L, determined by base_T5 index pivot
    //   (result should maintain line orientation)
  void set_fix_line(const line_T3& L, int pivot);
    // *this  <-  line L "regenerated" from rows:
    // - select 2 vectors using get_elts(pivot)
    // - regenerate proper line from elements
    //   (should be identity (modulo a constant factor) for proper lines)
  void set_xform_as_rows(const line_T3& L, int pivot, const matrix_T3& T);
    // *this  <-  line L transformed as rows by T:
    // - select 2 vectors using get_elts(pivot)
    // - transform selected vectors as rows by T
    // - generate transformed line from resulting vectors
  void set_xform_as_cols(const line_T3& L, int pivot, const matrix_T3& T);
    // *this  <-  line L transformed as columns by T:
    // - select 2 vectors using get_elts(pivot)
    // - transform selected vectors as columns by T
    // - generate transformed line from resulting vectors

public: // operations returning scalars
  scalar dot(const line_T3& L) const  { return base_T5::dot(L); }
    // in basis metric:  value= metric_im metric_jn (*this)^ij L^mn
    // or vector/covector:  value= (*this)^ij L_ij
  scalar orient(const line_T3& L) const;
    // 4x4 determinant from 2 lines:
    // value= epsilon_ijmn (*this)^ij L^mn   [commutative operation]
    // horsepower:  requires 6*, 5+  (same as dot product)
 
public: // projective operations returning line_T3 (replace/update *this)
  // geometric operations represent line_T3 as 6 independent coordinates
  //   [in (representation basis)^2] of an antisymmetric tensor

  // set_dual(const line_T3& L) is declared above

  void set_join_points(const point_T3& p, const point_T3& q)
    // *this  <-  oriented line through 2 points === p v q
    { set_nuller((const base_T3&)p,(const base_T3&)q); }
  void set_meet_planes(const plane_T3& d, const plane_T3& e)
    // *this  <-  oriented line intersection of 2 planes === d ^ e
    { set_dyad((const base_T3&)d,(const base_T3&)e); }

public: // geometric mapping operations (use map_vect())
  // join and dual with 3-D lines are commutative
  void join_point(const point_T3& p_in, plane_T3& d_out) const
    // d_out  <-  plane through line (*this) and point (p_in)
    //            === *this v p_in
    // (d_out)^j  <-  (*this)^ij (p_in)_i
    { map_vect(p_in, d_out); }
  void meet_plane(const plane_T3& d_in, point_T3& p_out) const
    // p_out  <-  point intersection of line (*this) and plane (d_in)
    //            === *this ^ d_in
    // (p_out)_j  <-  dual(*this)_ij (d_in)^i
    { line_T3 tmp; tmp.set_dual(*this);  tmp.map_vect(d_in, p_out); }

public: // projective operations returning scalars
  // geometric evaluation with 3-D lines is commutative
  scalar eval_line(const line_T3& L) const  { return orient(L); }
    // returns orient(*this,a) = epsilon_ijmn (*this)^ij a^mn
  scalar operator*(const line_T3& L) const  { return eval_line(L); }
  scalar operator()(const line_T3& L) const  { return eval_line(L); }

public: // affine methods
  void get_dir(point_T3& p_out) const {
    // p_out  <-  "direction" of line (*this) === oriented point-at-infinity
    // p_out  <-  (*this) ^ [1,0,0,0] === dual( dual(*this) v (1,0,0,0) )
    // (p_out)_j = dual(*this)_ij [1,0,0,0]^i
    //
    // specifically:  p_out = (0,-f,e,-c)
    p_out.set_point(0, -get_val(5), get_val(4), -get_val(3));
  }
  bool isideal(void) const {
    return get_val(5).iszero() && get_val(4).iszero() && get_val(3).iszero();
  }

  void set_dir(const plane_T3& d) {
    // *this  <-  "direction" of plane (d) === oriented line-at-infinity
    // *this  <-  d ^ [1,0,0,0] = dyad( d , (1,0,0,0) )
    // (*this)^ij  <-  d^i [1,0,0,0]^j - [1,0,0,0]^i d^j
    //
    // specifically:  *this  <-  [ 0, x, y, z]
    //                           [-x, 0, 0, 0]
    //                           [-y, 0, 0, 0]
    //                           [-z, 0, 0, 0]
    set_line(d[1],d[2],d[3], 0,0,0);
  }

public: // debug operations
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;

}; // end class line_T3


// linear transforms of basic objects

class matrix_T2 {
  // general matrix operating on base_T2 objects

private: // representation
  tens2<3,3> v;  // 3x3 matrix

public: // ctors, dtors, etc.
  matrix_T2(void) {}  // use default initializer (start as null)
  // use default copy ctor
  // use default dtor
  // use default assignment operator

public: // access methods
  bool isnull(void) const;
  scalar get_val(int row, int col) const;
  void get_col(int col, base_T2& b) const;
  void get_row(int row, base_T2& b) const;

public: // operations resulting in a matrix (replace/update *this)
  void set_null(void);
  void set_unit(void);
  void set_cols(const base_T2 &w, const base_T2 &x, const base_T2 &y);
  void set_rows(const base_T2 &w, const base_T2 &x, const base_T2 &y);

  void set_neg(const matrix_T2& A);
  void negate(void);  // inplace operation
  void set_transpose(const matrix_T2& A);
  void transpose(void);  // inplace operation

  void set_product(const matrix_T2 &A, const matrix_T2 &B);
    // *this  <-  A B
    // (*this)_ij  <-  A_ik B_kj
    // horsepower:  requires 27*, 18+
  scalar set_cofactor(const matrix_T2 &A);  // returns det(A)
    // *this  <-  cofactor matrix of A
    // horsepower:  requires 18*, 9+ for cofactor, plus add'l 3*, 2+ for det(A)
    // [generate 3 cross-products:  3(6*,3+)]

public: // operations returning base_T2
  void xform_row(const base_T2& r_in, base_T2& r_out) const;
    // r_out  <-  r_in (*this)   [treat r_in as row vector]
    // (r_out)_j  <-  (r_in)_i (*this)^i_j
    // horsepower:  requires 9*, 6+
  void xform_col(const base_T2& c_in, base_T2& c_out) const;
    // c_out  <-  (*this) c_in   [treat c_in as column vector]
    // (c_out)^i  <-  (*this)^i_j (c_in)^j
    // horsepower:  requires 9*, 6+

public: // operations resulting in a scalar
  scalar det(void) const;

public: // debug methods
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;
  
}; // end class matrix_T2


class matrix_T3 {
  friend class matrix_T5;

private: // representation
  tens2<4,4> v;  // 4x4 matrix

public: // ctors, etc.
  matrix_T3(void) {}  // use default initializer (start as null)
  // use default copy ctor
  // use default dtor
  // use default assignment operator

public: // access methods
  bool isnull(void) const;
  void get_col(int col, base_T3& b) const;
  void get_row(int col, base_T3& b) const;
  scalar get_val(int col, int row) const;

public: // operations resulting in a matrix (replace/update *this)
  void set_null(void);
  void set_unit(void);
  void set_cols(const base_T3 &w, const base_T3 &x,
		const base_T3 &y, const base_T3 &z);
  void set_rows(const base_T3 &w, const base_T3 &x,
		const base_T3 &y, const base_T3 &z);

  void set_neg(const matrix_T3& A);
  void negate(void);  // inplace operation
  void set_transpose(const matrix_T3& A);
  void transpose(void);  // inplace operation

  void set_product(const matrix_T3 &A, const matrix_T3 &B);
    // *this  <-  A B
    // (*this)_ij  <- A_ik B_kj
    // horsepower:  requires 64*, 48+
  scalar set_cofactor(const matrix_T3 &A);  // returns det(A)
    // *this  <-  cofactor matrix of A
    // horsepower: requires 72*, 44+ for cofactor, plus add'l 6*, 5+ for det(A)
    // [generate 2 lines, map 4 vectors:  2(12*,6+) + 4(12*,8+)...]

public: // operations returning base_T3
  void xform_row(const base_T3& r_in, base_T3& r_out) const;
    // r_out  <-  r_in (*this)   [treat r_in as row vector]
    // (r_out)_j  <-  (r_in)_i (*this)^i_j
    // horsepower:  requires 16*, 12+
  void xform_col(const base_T3& c_in, base_T3& c_out) const;
    // c_out  <-  (*this) c_in   [treat c_in as column vector]
    // (c_out)^i  <-  (*this)^i_j (c_in)^j
    // horsepower:  requires 16*, 12+

public: // operations resulting in a scalar
  scalar det(void) const;

public: // debug methods
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;
}; // end class matrix_T3

class matrix_T5 {
  // specialized support for 3-D line coordinates:
  //
  // treat default line (whose indices bind point vectors)
  //   as row vector of 6 coordinates (in abovementioned order)
  //   result is comparable row vector (also binding point vectors)
  //
  // if line matrix is L^pq,
  //   and transform tensor represented by *this is N^rs_tu
  //   transform result L'^ij = L^mn N^ij_mn = L^mn (T^i_m T^j_n - T^j_m T^i_n)
  //
  // representation is native for line_T3::set_xform_as_row()
  // thus, expand operation is:
  //   N^ij_mn = T^i_m T^j_n - T^j_m T^i_n
  //
  // so each column of the expanded matrix is a line determined by
  //   N^ij = set_dyad(T^i, T^j)  [of corresponding columns of (matrix_T3) T]

private: // representation
  tens2<6,6> v;  // 6x6 matrix

public: // ctors, etc.
  matrix_T5(void) {}  // use default initializer (start as null?)
  // use default copy ctor
  // use default dtor
  // use default assignment operator

public: // access methods
  bool isnull(void) const;

public: // operations resulting in a matrix (replace/update *this)
  void set_null(void);
  void set_neg(const matrix_T5& A);
  void negate(void);
  void set_transpose(const matrix_T5& A);
  void transpose(void);

  void set_expand(const matrix_T3& T);  // expand from 4x4 transform matrix
    // horsepower:  requires 72*, 36+
    // (could generate cofactor(T) at same time for an additional 48*,24+
    // and determinant(T) for an additional 6*,5+ ...)

public: // operations returning base_T3 or descendents
  void xform_row(const base_T5& r_in, base_T5& r_out) const;
    // r_out  <-  r_in (*this)   [treat r_in as row vector]
    // (r_out)_j  <-  (r_in)_i (*this)^i_j
    // horsepower:  requires 36*, 30+
  void xform_col(const base_T5& c_in, base_T5& c_out) const;
    // c_out  <-  (*this) c_in   [treat c_in as column vector]
    // (c_out)^i  <-  (*this)^i_j (c_in)^j
    // horsepower:  requires 36*, 30+

public: // debug methods
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;
}; // end class matrix_T5


class gmap_T2 {
  // general projective transform T2 -> T2

private: // representation
  bool pos;  // true if det(ptrans) > 0
  bool neg;  // true if det(ptrans) < 0
  matrix_T2 ptrans;  // transforms point -> point
  matrix_T2 dtrans;  // transforms dual -> dual
  // dtrans ~ cofactor(ptrans)

public: // ctors, etc.
  gmap_T2(void):pos(false),neg(false) {}  // start as null
  gmap_T2(const matrix_T2 &A) { set_p2p(A); }
  // use default copy ctor
  // use default dtor
  // use default assignment operator

public: // access operations
  bool ispos(void) const {return pos;} // det(ptrans) > 0 if true
  bool isneg(void) const {return neg;} // det(ptrans) < 0 if true

public: // operations resulting in a gmap_T2 (replace/update *this)
  void set_null(void);
  void set_unit(void);

  void set_gmap(const gmap_T2& M);  // (*this)  <-  M
  void set_inv(const gmap_T2& M);  // (*this)  <-  M^-1

  // ... Note: I may want to revise semantics/names for inverse operations ...
  void set_p2p(const matrix_T2& A);  // specify (*this)(point) -> point:
    // set *this  <-  M such that:  rep(M(point)) ~ rep(point) A
  void set_inv_p2p(const matrix_T2& Ainv);  // specify (*this)(point) -> point:
    // set *this  <-  M such that:  rep(point) ~ rep(M(point)) Ainv
  void set_d2d(const matrix_T2& B);  // specify (*this)(line) -> line:
    // set *this  <-  M such that:  rep(M(line)) ~ rep(line) B
  void set_inv_d2d(const matrix_T2& Binv);  // specify (*this)(line) -> line:
    // set *this  <-  M such that:  rep(line) ~ rep(M(line)) Binv

public: // geometric mapping operations
  void gmap_point(const point_T2& p_in, point_T2& p_out) const;
    // p_out  <-  M(p_in)
  void gmap_inv_point(const point_T2& p_in, point_T2& p_out) const;
    // p_out  <-  M^-1(p_in)
  void gmap_line(const line_T2& d_in, line_T2& d_out) const;
    // d_out  <-  M(d_in)
  void gmap_inv_line(const line_T2& d_in, line_T2& d_out) const;
    // d_out  <-  M^-1(d_in)

public: // debug methods
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;
}; // end class gmap_T2


class gmap_T3 {
  // general projective transform T3 -> T3

private: // representation
  bool pos;
  bool neg;
  matrix_T3 ptrans;  // transforms point -> point
  matrix_T3 dtrans;  // transforms dual -> dual
    // dtrans ~ cofactor(ptrans)
  matrix_T5 ltrans;  // transforms line -> line
    // ltrans is expanded 6x6 matrix

public: // ctors, etc.
  gmap_T3(void):pos(false),neg(false) {}  // start as null
  gmap_T3(const matrix_T3 &A) { set_p2p(A); }
  // use default copy ctor
  // use default dtor
  // use default assignment operator

public: // access operations
  bool ispos(void) const {return pos;} // det(ptrans) > 0 if true
  bool isneg(void) const {return neg;} // det(ptrans) < 0 if true
  int sgn(void) const;  // sgn(det(ptrans)) [one of -1,0,1]

public: // operations resulting in a gmap_T2 (replace/update *this)
  void set_null();
  void set_unit();

  void set_gmap(const gmap_T3& M);  // (*this)  <-  M
  void set_inv(const gmap_T3& M);  // (*this)  <-  M^-1

  // ... Note:  I may want to revise semantics/names for inverse operations ...
  void set_p2p(const matrix_T3& A);  // specify (*this)(point) -> point:
    // set *this  <-  M such that:  rep(M(point)) ~ rep(point) A
  void set_inv_p2p(const matrix_T3& Ainv);  // specify (*this)(point) -> point:
    // set *this  <-  M such that:  rep(point) ~ rep(M(point)) Ainv
  void set_d2d(const matrix_T3& B);  // specify (*this)(plane) -> plane:
    // set *this  <-  M such that:  rep(M(plane)) ~ rep(plane) B
  void set_inv_d2d(const matrix_T3& binv);  // specify (*this)(plane) -> plane:
    // set *this  <-  M such that:  rep(plane) ~ rep(M(d)) Binv

public: // geometric mapping operations
  void gmap_point(const point_T3& p_in, point_T3& p_out) const;
    // p_out  <-  M(p_in)
  void gmap_inv_point(const point_T3& p_in, point_T3& p_out) const;
    // p_out  <-  M^-1(p_in)
  void gmap_plane(const plane_T3& d_in, plane_T3& d_out) const;
    // d_out  <-  M(d_in)
  void gmap_inv_plane(const plane_T3& d_in, plane_T3& d_out) const;
    // d_out  <-  M^-1(d_in)
  void gmap_line(const line_T3& L_in, line_T3& L_out) const;
    // L_out  <-  M(L_in)
  void gmap_inv_line(const line_T3& L_in, line_T3& L_out) const;
    // L_out  <-  M^-1(L_in)

public: // debug methods
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;
}; // end class gmap_T3


class dmap_T2 {
  // general duomorphism T2 -> T2
  friend class point_T2;
  friend class line_T2;

private: // representation
  matrix_T2 ptrans;  // transforms point -> dual
  matrix_T2 dtrans;  // transforms dual -> point
  // dtrans ~ ptrans^-1

public: // ctors, etc.
  dmap_T2(void) {}  // start as null
  dmap_T2(const matrix_T2 &A) { set_p2d(A); }
  // use default copy ctor
  // use default dtor
  // use default assignment operator

public: // operations resulting in a dmap_T2 (replace/update *this)
  void set_null();
  void set_unit();

  void set_dmap(const dmap_T2& D);  // (*this)  <-  D
  void set_inv(const dmap_T2& D);  // (*this)  <-  D^-1

  // ... Note:  I expect half of these are redundant ...
  void set_p2d(const matrix_T2& A); // specify (*this)(point) -> line:
    // set *this  <-  D such that:  rep(D(point)) ~ rep(point) A
  void set_inv_p2d(const matrix_T2& Ainv); // specify (*this)(point) -> line:
    // set *this  <-  D such that:  rep(point) ~ rep(D(point)) Ainv
  void set_d2p(const matrix_T2& B); // specify (*this)(line) -> point:
    // set *this  <-  D such that:  rep(D(line)) ~ rep(line) B
  void set_inv_d2p(const matrix_T2& Binv); // specify (*this)(line) -> point:
    // set *this  <-  D such that:  rep(line) ~ rep(D(line)) Binv

public: // geometric mapping operations
  void dmap_line(const line_T2& d_in, point_T2& p_out) const;
    // p_out  <-  D(d_in)
  void dmap_inv_line(const line_T2& d_in, point_T2& p_out) const;
    // p_out  <-  D^-1(d_in)
  void dmap_point(const point_T2& p_in, line_T2& d_out) const;
    // d_out  <-  D(p_in)
  void dmap_inv_point(const point_T2& p_in, line_T2& d_out) const;
    // d_out  <-  D^-1(p_in)

public: // debug methods
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;
}; // end class dmap_T2


class dmap_T3 {
  // general duomorphism T3 -> T3
  friend class point_T3;
  friend class plane_T3;
  friend class line_T3;

private: // representation
  matrix_T3 ptrans;  // transforms point -> dual
  matrix_T3 dtrans;  // transforms dual -> point
    // dtrans ~ ptrans^-1
  matrix_T5 ltrans;  // transforms line -> dual line
    // (consistent with above...)

public: // ctors, etc.
  dmap_T3(void) {}  // start as null
  dmap_T3(const matrix_T3 &A) { set_p2d(A); }
  // use default copy ctor
  // use default dtor
  // use default assignment operator

public: // operations resulting in a gmap_T3 (replace/update *this)
  void set_null();
  void set_unit();
  void set_dmap(const dmap_T3& D);
  void set_inv(const dmap_T3& D);

  // ... Note:  I expect half of these are redundant ...
  void set_p2d(const matrix_T3& A); // specify (*this)(point) -> plane:
    // set *this  <-  D such that:  rep(D(point)) ~ rep(point) A
  void set_inv_p2d(const matrix_T3& Ainv); // specify (*this)(point) -> plane:
    // set *this  <-  D such that:  rep(point) ~ rep(D(point)) Ainv
  void set_d2p(const matrix_T3& B); // specify (*this)(plane) -> point:
    // set *this  <-  D such that:  rep(D(plane)) ~ rep(plane) B
  void set_inv_d2p(const matrix_T3& Binv); // specify (*this)(plane) -> point:
    // set *this  <-  D such that:  rep(plane) ~ rep(D(plane)) Binv

public: // geometric mapping operations
  void dmap_point(const point_T3& p_in, plane_T3& d_out) const;
    // d_out  <-  D(p_in)
  void dmap_inv_point(const point_T3& p_in, plane_T3& d_out) const;
    // d_out  <-  D^-1(p_in)
  void dmap_plane(const plane_T3& d_in, point_T3& p_out) const;
    // p_out  <-  D(d_in)
  void dmap_inv_plane(const plane_T3& d_in, point_T3& p_out) const;
    // p_out  <-  D^-1(d_in)
  void dmap_line(const line_T3& L_in, line_T3& L_out) const;
    // L_out  <-  D(L_in)
  void dmap_inv_line(const line_T3& L_in, line_T3& L_out) const;
    // L_out  <-  D^-1(L_in)

public: // debug methods
  bool ok(void) const;
  void dump(FILE *f) const;
  bool analyze(FILE *f) const;
}; // end class dmap_T3




// inlined method definitions:
//   -- small, trivial, self-documenting functions may be in class declaration
//   -- small and potentially time-critical functions may be inlined here
//      as well as trivial functions that require forwards declaration
//   -- remainder should be defined in geometry.cc

// class point_T1
inline scalar point_T1::get_val(int i) const {
  assert(0<=i && i<2);
  return v[i];
}
inline bool point_T1::isnull(void) const {
  return get_val(0).iszero() && get_val(1).iszero();
}

inline void point_T1::set_point(scalar w, scalar x) {
  v[0]= w; v[1]= x;
}
inline void point_T1::set_null(void) {
  set_point(0,0);
}
inline void point_T1::set_anti(const point_T1& a) {
  // set *this to -a
  // (rotate 180 degrees)
  set_point(-a[0], -a[1]);
}
inline void point_T1::set_dual_right(const point_T1& a) {
  // set *this to right dual of a in representation basis
  // (rotate 90 degrees counterclockwise)
  // such that |a.w w| === |a, *this|  >=  0
  //           |a.x x|
  set_point(-a[1], a[0]);
}
inline void point_T1::set_dual_left(const point_T1& a) {
  // left *this to left dual of a in representation basis
  // (rotate 90 degrees clockwise)
  // such that |w a.w| === |*this, a|  >=  0
  //           |x a.x|
  set_point(a[1], -a[0]);
}
inline scalar point_T1::dot(const point_T1& a) const {
  accumulator s;
  s.mac(get_val(0),a[0]);
  s.mac(get_val(1),a[1]);
  return s.convert(2, 1); // convert sum of 2 products -> logsums = 1
}
inline scalar point_T1::det(const point_T1& a) const {
  accumulator s;
  s.mac (get_val(0), a[1]);
  s.msub(get_val(1), a[0]);
  return s.convert(2,1); // convert sum of 2 products -> logsums = 1
}



// class base_T2
inline scalar base_T2::get_val(int i) const {
  assert(0<=i && i<3);
  return v[i];
}
inline bool base_T2::isnull(void) const {
  return get_val(0).iszero() && get_val(1).iszero() && get_val(2).iszero();
}
inline void base_T2::set_base(scalar w, scalar x, scalar y)
  { v[0]= w; v[1]= x; v[2]= y; }
inline void base_T2::set_null(void) {
  set_base(0,0,0);
}
inline void base_T2::set_anti(const base_T2& a) {
  set_base(-a[0], -a[1], -a[2]);
}

inline void base_T2::set_cross(const base_T2& a, const base_T2& b) {
  // sum of 2 products -> logsums = 1
  accumulator s;
  
  s.zero();
  s.mac (a.v[1], b.v[2]);
  s.msub(a.v[2], b.v[1]);
  v[0]= s.convert(2,1);
  
  s.zero();
  s.mac (a.v[2], b.v[0]);
  s.msub(a.v[0], b.v[2]);
  v[1]= s.convert(2,1);
  
  s.zero();
  s.mac (a.v[0], b.v[1]);
  s.msub(a.v[1], b.v[0]);
  v[2]= s.convert(2,1);
}
inline scalar base_T2::dot(const  base_T2& a) const {
  // sum of 3 products -> logsums = 2
  accumulator s;
  s.mac(v[0],a.v[0]);
  s.mac(v[1],a.v[1]);
  s.mac(v[2],a.v[2]);
  return s.convert(2,2);
}
inline scalar base_T2::det(const base_T2& b, const base_T2& c) const {
  base_T2 tmp;  tmp.set_cross(b,c);
  return dot(tmp);
}



// class point_T2


// class line_T2


// class base_T3
inline scalar base_T3::get_val(int i) const {
  assert(0<=i && i<4);
  return v[i];
}
inline bool base_T3::isnull(void) const {
  return get_val(0).iszero() && get_val(1).iszero()
      && get_val(2).iszero() && get_val(3).iszero();
}
inline void base_T3::set_base(scalar w, scalar x, scalar y, scalar z)
  { v[0]=w; v[1]=x; v[2]= y; v[3]= z; }
inline void base_T3::set_null(void) {
  set_base(0,0,0,0);
}
inline void base_T3::set_anti(const base_T3& a) {
  set_base(-a[0], -a[1], -a[2], -a[3]);
}
inline scalar base_T3::dot(const base_T3& a) const {
  // sum of 4 products -> logsums = 2
  accumulator s;
  s.mac(get_val(0),a[0]);
  s.mac(get_val(1),a[1]);
  s.mac(get_val(2),a[2]);
  s.mac(get_val(3),a[3]);
  return s.convert(2,2);
}
inline scalar
base_T3::det(const base_T3& b, const base_T3& c, const base_T3& d) const {
  line_T3 p,q;
  p.set_dyad(*this,b);  q.set_dyad(c,d);
  return p.orient(q);
}


// class point_T3
/*
inline void point_T3::set_dmap(const plane_T3& d, const dmap_T3& D) {
  // *this= rep(D(d)) = rep(d) D.dtrans  [= rep(d) B]
  set_xform_row((const base_T3&)d, D.dtrans);
}
inline void point_T3::set_dmap_inv(const plane_T3& d, const dmap_T3& D) {
  // *this= rep(D^-1(d)) = - rep(d) (D.dtrans)^t  [= - rep(d) B^t]
  base_T3 tmp;
  tmp.set_xform_col((const base_T3&)d, D.dtrans);
  base_T3::set_anti(tmp);
}
*/


// class plane_T3

inline void plane_T3::get_dir(line_T3& L_out) const {
  L_out.set_dir(*this);
}

/*
inline void plane_T3::set_dmap(const point_T3& p, const dmap_T3& D) {
  // *this= rep(D(p)) = rep(p) D.ptrans  [= rep(p) A]
  set_xform_row((const base_T3&)p, D.ptrans);
}
inline void plane_T3::set_dmap_inv(const point_T3& p, const dmap_T3& D) {
  // *this= rep(D^-1(p)) = - rep(p) (D.ptrans)^t [= rep(p) A^t]
  base_T3 tmp;
  tmp.set_xform_col((const base_T3&)p, D.ptrans);
  base_T3::set_anti(tmp);
}
*/


// class base_T5
inline scalar& base_T5::get_ref(int i) {
  assert(0<=i && i<6);
  return v[i];
}
inline scalar base_T5::get_val(int i) const {
  assert(0<=i && i<6);
  return v[i];
}
inline bool base_T5::isnull(void) const {
  return get_val(0).iszero() && get_val(1).iszero() && get_val(2).iszero()
      && get_val(3).iszero() && get_val(4).iszero() && get_val(5).iszero();
}
inline void base_T5::set_base(scalar a, scalar b, scalar c,
                              scalar d, scalar e, scalar f) {
  v[0]= a; v[1]= b; v[2]= c; v[3]= d; v[4]= e; v[5]= f;
}
inline void base_T5::set_null(void) {
  set_base(0,0,0,0,0,0);
}
inline void base_T5::set_anti(const base_T5& b) {
  set_base(-b[0], -b[1], -b[2], -b[3], -b[4], -b[5]);
}

inline scalar base_T5::dot(const base_T5& a) const {
  accumulator s;
  s.mac(get_val(0),a[0]);
  s.mac(get_val(1),a[1]);
  s.mac(get_val(2),a[2]);
  s.mac(get_val(3),a[3]);
  s.mac(get_val(4),a[4]);
  s.mac(get_val(5),a[5]);
  return s.convert(2,3); // convert sum of 6 products -> logsums=3
}

// class line_T3
inline void
line_T3::set_line(scalar a, scalar b, scalar c, scalar d, scalar e, scalar f) 
  { set_base(a,b,c,d,e,f); }
inline void line_T3::set_dual(const line_T3& L)
  { set_line( L[5], -L[4],  L[3],  L[2], -L[1],  L[0]); }
inline void line_T3::set_dyad(const base_T3& a, const base_T3& b) {
  // sum of 2 products -> logsums = 1
  accumulator s;
  
  s.zero();
  s.mac (a[0], b[1]);
  s.msub(a[1], b[0]);
  get_ref(0)= s.convert(2,1);   // v[0]  <-  L^01
  
  s.zero();
  s.mac (a[0], b[2]);
  s.msub(a[2], b[0]);
  get_ref(1)= s.convert(2,1);   // v[1]  <-  L^02
  
  s.zero();
  s.mac (a[0], b[3]);
  s.msub(a[3], b[0]);
  get_ref(2)= s.convert(2,1);   // v[2]  <-  L^03
  
  s.zero();
  s.mac (a[1], b[2]);
  s.msub(a[2], b[1]);
  get_ref(3)= s.convert(2,1);   // v[3]  <-  L^12
  
  s.zero();
  s.mac (a[1], b[3]);
  s.msub(a[3], b[1]);
  get_ref(4)= s.convert(2,1);   // v[4]  <- L^13
  
  s.zero();
  s.mac (a[2], b[3]);
  s.msub(a[3], b[2]);
  get_ref(5)= s.convert(2,1);   // v[5]  <-  L^23
}
inline void line_T3::set_nuller(const base_T3& a, const base_T3& b) {
  line_T3 tmp;  tmp.set_dyad(a,b);
  set_dual(tmp);
}
inline void line_T3::map_vect(const base_T3& b_in, base_T3& b_out) const {
  // sum of 3 products -> logsums=2
  vect<4> tmp;
  accumulator s;

  s.zero();
  //          ignore b_in[0]    -- L^00 = 0
  s.mac(get_val(0),  b_in[1]);  // L^01 = (*this)[0]
  s.mac(get_val(1),  b_in[2]);  // L^02 = (*this)[1]
  s.mac(get_val(2),  b_in[3]);  // L^03 = (*this)[2]
  tmp[0]= s.convert(2,2);  // L^00 b_0 + L^01 b_1 + L^02 b_2 + L^03 b_3

  s.zero();
  s.msub(get_val(0), b_in[0]);  // L^10 = -(*this)[0]
  //          ignore b_in[1]    -- L^11 = 0
  s.mac(get_val(3),  b_in[2]);  // L^12 =  (*this)[3]
  s.mac(get_val(4),  b_in[3]);  // L^13 =  (*this)[4]
  tmp[1]= s.convert(2,2);  // L^10 b_0 + L^11 b_1 + L^12 b_2 + L^13 b_3

  s.zero();
  s.msub(get_val(1), b_in[0]);  // L^20 = -(*this)[1]
  s.msub(get_val(3), b_in[1]);  // L^21 = -(*this)[3]
  //          ignore b_in[2]    -- L^22 = 0
  s.mac(get_val(5),  b_in[3]);  // L^23 =  (*this)[5]
  tmp[2]= s.convert(2,2);  // L^20 b_0 + L^21 b_1 + L^22 b_2 + L^23 b_3

  s.zero();
  s.msub(get_val(2), b_in[0]);  // L^30 = -(*this)[2]
  s.msub(get_val(4), b_in[1]);  // L^31 = -(*this)[4]
  s.msub(get_val(5), b_in[2]);  // L^32 = -(*this)[5]
  //          ignore b_in[3]    -- L^22 = 0
  tmp[3]= s.convert(2,2);  // L^30 b_0 + L^31 b_1 + L^32 b_2 + L^33 b_3

  b_out.set_base(tmp[0],tmp[1],tmp[2],tmp[3]);
}
inline void line_T3::get_vects(int pivot, base_T3& b, base_T3& c) const {
  // to maintain line orientation, swap b,c for positive pivot
  if(get_val(pivot).ispos()) {
    get_vect(first_index(pivot), b);
    get_vect(second_index(pivot), c);
  } else {
    get_vect(first_index(pivot), c);
    get_vect(second_index(pivot), b);
  }
}
inline void line_T3::set_fix_line(const line_T3& L, int pivot) {
  base_T3 p,q;  L.get_vects(pivot, p,q);
  set_dyad(p,q);
}
inline void line_T3::set_xform_as_rows(
const line_T3& L, int pivot, const matrix_T3& T) {
  base_T3 p,q;  get_vects(pivot, p,q);

  base_T3 r,s;
  T.xform_row(p,r);
  T.xform_row(q,s);

  set_dyad(r,s);
}
inline void line_T3::set_xform_as_cols(
const line_T3& L, int pivot, const matrix_T3& T) {
  base_T3 p,q;  get_vects(pivot, p,q);

  base_T3 r,s;
  T.xform_col(p,r);
  T.xform_col(q,s);

  set_dyad(r,s);
}

inline scalar line_T3::orient(const line_T3& L) const {
  accumulator s;
  s.mac(get_val(0),L[5]);
  s.msub(get_val(1),L[4]);
  s.mac(get_val(2),L[3]);
  s.mac(get_val(3),L[2]);
  s.msub(get_val(4),L[1]);
  s.mac(get_val(5),L[0]);
  return s.convert(2,3); // convert sum of 6 products -> logsums=3
}


// class matrix_T2
inline bool matrix_T2::isnull(void) const {
  return v.iszero();
}
inline scalar matrix_T2::get_val(int row, int col) const {
  return v.getval(row,col);
}
inline void matrix_T2::get_col(int col, base_T2& b) const {
  b.set_base(v.getval(0,col),
	     v.getval(1,col),
	     v.getval(2,col)
	     );
}
inline void matrix_T2::get_row(int row, base_T2& b) const {
  b.set_base(v.getval(row,0),
	     v.getval(row,1),
	     v.getval(row,2)
	     );
}
inline void matrix_T2::set_null(void) {
  v.zero();
}
inline void matrix_T2::set_unit(void) {
  v.zero();
  for(int i=0; i<3; i++)  v.setval(i,i, scalar::maxval);
}
inline void matrix_T2::set_cols(const base_T2 &w,
				const base_T2 &x,
				const base_T2 &y) {
  v.colref(0)= w.v;
  v.colref(1)= x.v;
  v.colref(2)= y.v;
}
inline void matrix_T2::set_rows(const base_T2 &w,
				const base_T2 &x,
				const base_T2 &y) {
  set_cols(w,x,y);
  transpose();
}

inline void matrix_T2::set_neg(const matrix_T2 &A) {
  v.setneg(A.v);
}
inline void matrix_T2::negate(void) {
  v.negate();
}

inline void matrix_T2::set_transpose(const matrix_T2 &A) {
  v.set_transpose(A.v);
}
inline void matrix_T2::transpose(void) {
  v.transpose();
}

inline void matrix_T2::set_product(const matrix_T2 &A, const matrix_T2 &B) {
  // *this= matrix multiply:  T_ij =  A_ik B_kj 
  // (there may be a positive scaling factor applied to the result)
  //
  // notation, um, notes:
  // - when applied to column vector c,
  //   (T c)_i = T_ij c_j = A_ik B_kj c_j  ->  T c = A (B c) = (A B) c
  // - when applied to row vector r,
  //   (r T)_j = r_i T_ij = r_i A_ik B_kj  ->  r T = (r A) B = r (A B)
  v.prod_matrix<3,2>(A.v,B.v);
}
inline void matrix_T2::xform_row(const base_T2& r_in, base_T2& r_out) const
  { r_out.v.xform_row<3,2>(r_in.v, v); }
inline void matrix_T2::xform_col(const base_T2& c_in, base_T2& c_out) const
  { c_out.v.xform_col<3,2>(v, c_in.v); }
inline scalar matrix_T2::det(void) const {
  base_T2 a,b,c;

  get_col(0,a);
  get_col(1,b);
  get_col(2,c);
  return a.det(b,c);
}


// class matrix_T3
inline bool matrix_T3::isnull(void) const {
  return v.iszero();
}
inline void matrix_T3::get_col(int col, base_T3& b) const {
  b.set_base(v.getval(0,col),
	     v.getval(1,col),
	     v.getval(2,col),
	     v.getval(3,col)
	     );
}
inline void matrix_T3::get_row(int row, base_T3& b) const {
  b.set_base(v.getval(row,0),
	     v.getval(row,1),
	     v.getval(row,2),
	     v.getval(row,3)
	     );
}
inline scalar matrix_T3::get_val(int row, int col) const {
  return v.getval(row,col);
}

inline void matrix_T3::set_null(void) {
  v.zero();
}
inline void matrix_T3::set_unit(void) {
  v.zero();
  for(int i=0; i<4; i++)  v.setval(i,i, scalar::maxval);
}
inline void matrix_T3::set_cols(const base_T3 &w, const base_T3 &x,
				const base_T3 &y, const base_T3 &z) {
  v.colref(0)= w.v;
  v.colref(1)= x.v;
  v.colref(2)= y.v;
  v.colref(3)= z.v;
}
inline void matrix_T3::set_rows(const base_T3 &w, const base_T3 &x,
				const base_T3 &y, const base_T3 &z) {
  // cheesy expedient... when optimizing, write it out
  set_cols(w,x,y,z);
  v.transpose();
}

inline void matrix_T3::set_neg(const matrix_T3 &A) {
  v.setneg(A.v);
}
inline void matrix_T3::negate(void) {
  v.negate();
}

inline void matrix_T3::set_product(const matrix_T3 &A, const matrix_T3 &B) {
  // *this= matrix multiply:  T_ij= A_ik B_kj
  // (there may be a positive scaling factor applied to the result)
  //
  // in conventional matrix notation:
  // - when applied to column vector c,
  //   T c = T_ij c_j = A_ik B_kj c_j  ->  T c = A (B c) = (A B) c
  // - when applied to row vector r,
  //   r T = r_i T_ij = r_i A_ik B_kj  ->  r T = (r A) B = r (A B)
  v.prod_matrix<4,2>(A.v,B.v);
}
inline void matrix_T3::xform_row(const base_T3& r_in, base_T3& r_out) const
  { r_out.v.xform_row<4,2>(r_in.v, v); }
inline void matrix_T3::xform_col(const base_T3& c_in, base_T3& c_out) const
  { c_out.v.xform_col<4,2>(v, c_in.v); }

inline scalar matrix_T3::det(void) const {
  base_T3 a,b,c,d;

  get_col(0,a);
  get_col(1,b);
  get_col(2,c);
  get_col(3,d);

  return a.det(b,c,d);
}


// class matrix_T5:
inline bool matrix_T5::isnull(void) const {
  return v.iszero();
}
inline void matrix_T5::set_null(void) {
  v.zero();
}
inline void matrix_T5::set_neg(const matrix_T5 &A) {
  v.setneg(A.v);
}
inline void matrix_T5::negate(void) {
  v.negate();
}
inline void matrix_T5::set_transpose(const matrix_T5 &A) {
  v.set_transpose(A.v);
}
inline void matrix_T5::transpose(void) {
  v.transpose();
}
inline void matrix_T5::xform_row(const base_T5& r_in, base_T5& r_out) const {
  r_out.v.xform_row<6,3>(r_in.v, v);
}
inline void matrix_T5::xform_col(const base_T5& c_in, base_T5& c_out) const {
  c_out.v.xform_col<6,3>(v, c_in.v);
}


// class gmap_T2:
inline void gmap_T2::set_null(void) {
  ptrans.set_null();
  dtrans.set_null();
  pos= false;
  neg= false;
}
inline void gmap_T2::set_unit(void) {
  ptrans.set_unit();
  dtrans.set_unit();
  pos= true;
  neg= false;
}
inline void gmap_T2::set_gmap(const gmap_T2& M) {
  *this= M;
}
inline void gmap_T2::set_p2p(const matrix_T2& A) {
  ptrans= A;
  scalar det= dtrans.set_cofactor(A);
  pos= det.ispos();
  neg= det.isneg();
}
inline void gmap_T2::set_inv(const gmap_T2& M) {
  assert(this != &M);
  pos= M.pos;
  neg= M.neg;
  ptrans.set_transpose(M.dtrans);
  dtrans.set_transpose(M.ptrans);
  if(pos) {
  }
  else if(neg) {
    ptrans.negate();
    dtrans.negate();
  } else {
    // there may be a better thing to do in this case...
  }
}
inline void gmap_T2::set_inv_p2p(const matrix_T2& A) {
  scalar det= ptrans.set_cofactor(A);
  ptrans.transpose();
  dtrans= A;
  dtrans.transpose();
  if(det.ispos()) {
    pos= true;
    neg= false;
  } else if(det.isneg()) {
    pos= false;
    neg= true;
    ptrans.negate();
    dtrans.negate();
  } else {
    pos= false;
    neg= false;
    // there may be a better thing to do in this case...
  }
}
inline void
gmap_T2::gmap_point(const point_T2& p_in, point_T2& p_out) const {
  // p_out= M(p_in) = p_in ptrans   [M === (*this)]
  ptrans.xform_row(p_in, p_out);
}
inline void
gmap_T2::gmap_inv_point(const point_T2& p_in, point_T2& p_out) const {
  // p_out= M^-1(p_in) = p_in (sgn()*dtrans^t)   [M === (*this)]
  base_T2 tmp;
  dtrans.xform_col(p_in, tmp);
  if(isneg()) {
    p_out.base_T2::set_anti(tmp);
  } else {
    p_out.set_base(tmp);
  }
}
inline void
gmap_T2::gmap_line(const line_T2& d_in, line_T2& d_out) const {
  // d_out= M(d_in) = d_in dtrans   [M === *this]
  dtrans.xform_row(d_in, d_out);
}
inline void
gmap_T2::gmap_inv_line(const line_T2& d_in, line_T2& d_out) const {
  // d_out= M^-1(d_in) = d_in (sgn()*ptrans^t)   [M === *this]
  base_T2 tmp;
  ptrans.xform_col(d_in, tmp);
  if(isneg()) {
    d_out.base_T2::set_anti(tmp);
  } else {
    d_out.set_base(tmp);
  }
}

// class gmap_T3
inline void gmap_T3::set_null() {
  ptrans.set_null();
  dtrans.set_null();
  pos= false;
  neg= false;
}
inline void gmap_T3::set_unit() {
  ptrans.set_unit();
  dtrans.set_unit();
  pos= true;
  neg= false;
}
inline void gmap_T3::set_gmap(const gmap_T3& M) {
  *this= M;
}
inline void gmap_T3::set_inv(const gmap_T3& M) {
  assert(this!= &M);
  pos= M.pos;
  neg= M.neg;
  ptrans.set_transpose(M.dtrans);
  dtrans.set_transpose(M.ptrans);
  ltrans.set_transpose(M.ltrans);
  if(pos) {}
  else if(neg) {
    ptrans.negate();
    dtrans.negate();
    ltrans.negate();
  } else {  // near-singular case
    // there may be a better thing to do here...
  }
}
inline void gmap_T3::set_p2p(const matrix_T3& A) {
  ptrans= A;
  ltrans.set_expand(ptrans);  // this makes set_cofactor somewhat redundant...
  scalar det= dtrans.set_cofactor(ptrans);
  pos= det.ispos();
  neg= det.isneg();
}
inline void gmap_T3::set_inv_p2p(const matrix_T3& A) {
  gmap_T3 tmp;
  tmp.set_p2p(A);
  set_inv(tmp);
}
inline void
gmap_T3::gmap_point(const point_T3& p_in, point_T3& p_out) const {
  // p_out= M(p_in) = p_in ptrans   [M === (*this)]
  ptrans.xform_row(p_in, p_out);
}
inline void
gmap_T3::gmap_inv_point(const point_T3& p_in, point_T3& p_out) const {
  // p_out= M^-1(p_in) = p_in (sgn()*dtrans^t)   [M === (*this)]
  base_T3 tmp;
  dtrans.xform_col(p_in, tmp);
  if(isneg()) {
    p_out.base_T3::set_anti(tmp);
  } else {
    p_out.base_T3::set_base(tmp);
  }
}
inline void
gmap_T3::gmap_plane(const plane_T3& d_in, plane_T3& d_out) const {
  // d_out= M(d_in) = d_in dtrans   [M === (*this)]
  dtrans.xform_row(d_in, d_out);
}
inline void
gmap_T3::gmap_inv_plane(const plane_T3& d_in, plane_T3& d_out) const {
  // d_out= M^-1(d_in) = d_in (sgn()*ptrans^t)   [M === (*this)]
  base_T3 tmp;
  ptrans.xform_col(d_in, d_out);
  if(isneg()) {
    d_out.base_T3::set_anti(tmp);
  } else {
    d_out.base_T3::set_base(tmp);
  }
}
inline void
gmap_T3::gmap_line(const line_T3& L_in, line_T3& L_out) const {
  // L_out= M(L_in) = L_in ltrans   [M === (*this)]
  ltrans.xform_row(L_in, L_out);
}
inline void
gmap_T3::gmap_inv_line(const line_T3& L_in, line_T3& L_out) const {
}



// class dmap_T2
inline void dmap_T2::set_null() {
  ptrans.set_null();
  dtrans.set_null();
}
inline void dmap_T2::set_unit() {
  ptrans.set_unit();
  dtrans.set_unit();
}
inline void dmap_T2::set_dmap(const dmap_T2& D) {
  *this= D;
}
inline void dmap_T2::set_inv(const dmap_T2& D) {
  assert(this != &D);
  ptrans.set_transpose(D.ptrans);
  dtrans.set_transpose(D.dtrans);
}
inline void dmap_T2::set_p2d(const matrix_T2& A) {
  // rep(D(point)) === rep(point) A
  ptrans= A;
  scalar det= dtrans.set_cofactor(ptrans);
  if(det.isneg()) {
    dtrans.negate();
  }
}
inline void dmap_T2::set_d2p(const matrix_T2& B) {
  // rep(D(dual)) === rep(dual) B
  dtrans= B;
  scalar det= ptrans.set_cofactor(ptrans);
  if(det.isneg()) {
    ptrans.negate();
  }
}
inline void
dmap_T2::dmap_line(const line_T2& d_in, point_T2& p_out) const {
  // p_out= D(d_in) = d_in dtrans   [D === *this]
  dtrans.xform_row(d_in, p_out);
}
inline void
dmap_T2::dmap_inv_line(const line_T2& d_in, point_T2& p_out) const {
  // p_out= D^-1(d_in) = d_in (dtrans^t)   [D === *this]
  dtrans.xform_col(d_in, p_out);
}
inline void
dmap_T2::dmap_point(const point_T2& p_in, line_T2& d_out) const {
  ptrans.xform_row(p_in, d_out);
}
inline void
dmap_T2::dmap_inv_point(const point_T2& p_in, line_T2& d_out) const {
  ptrans.xform_col(p_in, d_out);
}


// class dmap_T3
inline void dmap_T3::set_null() {
  ptrans.set_null();
  dtrans.set_null();
}
inline void dmap_T3::set_unit() {
  ptrans.set_unit();
  dtrans.set_unit();
  dtrans.negate();
}
inline void dmap_T3::set_dmap(const dmap_T3& D) {
  *this= D;
}
inline void dmap_T3::set_p2d(const matrix_T3& A) {
  // rep(D(point)) === rep(point) A
  ptrans= A;
  scalar det= dtrans.set_cofactor(ptrans);
  if(det.ispos()) {
    dtrans.negate();
  }
}
inline void dmap_T3::set_inv(const dmap_T3& D) {
  assert(this != &D);
  ptrans.set_transpose(D.ptrans);
  ptrans.negate();
  dtrans.set_transpose(D.dtrans);
  ptrans.negate();
}
inline void dmap_T3::set_d2p(const matrix_T3& B) {
  // rep(D(dual)) === rep(dual) B
  dtrans= B;
  scalar det= ptrans.set_cofactor(ptrans);
  if(det.ispos()) {
    ptrans.negate();
  }
}




#endif // GEOMETRY_HH