
 

// Emacs: -*- C++ -*-

#ifndef IA_NO_RCS_STATIC_DATA
static char rcsid[] = "$Id: Point_c.m4,v 1.25.1.2 1994/12/08 14:31:41 thoth Exp $";
#endif

/*
        Copyright 1993, Center for Computer Vision and Visualization,
        University of Florida.  All rights reserved.

// $Log: Point_c.m4,v $
// Revision 1.25.1.2  1994/12/08  14:31:41  thoth
// All throws are now of class objects.
//
// Revision 1.25.1.1  1994/11/30  16:16:38  thoth
// Added cross product and dot product.
//
// Revision 1.25  1994/11/10  18:49:44  ljr
// Added IA_NO_RCS_STATIC_DATA guard.
//
// Revision 1.24  1994/10/17  19:09:42  fjsoria
// *** empty log message ***
//
// Revision 1.23  1994/09/13  21:02:51  fjsoria
// correction to derange for no_loop_unrolling
//
// Revision 1.22  1994/08/25  19:28:53  fjsoria
// abs is not defined in c-front redefined abs, fixed bug in derange
//
// Revision 1.21  1994/08/25  18:07:25  fjsoria
// removed abs for scalars since it is supported by UNIX, DOS, and Mac
// platforms
//
// Revision 1.20  1994/08/25  03:24:11  gmt
// Added Mac specific pragma (template_access) to tell the compiler that the
// generated files are actually template instances via template instantiation.
//
// Revision 1.19  1994/08/22  14:39:02  thoth
// fixed conflict with system's abs.
//
// Revision 1.18  1994/08/18  17:51:28  fjsoria
// corrected DoublePoint.h .c to DblPoint.h .c for DOS compat.
//
// Revision 1.17  1994/08/16  21:03:13  fjsoria
// added derange and enorm as non-member functions, keeping member versions
// added inorm, mnorm, abs, fabs as non-member functions
//
// Revision 1.16  1994/07/25  16:42:26  thoth
// FloatPoint is now either Point<double> or DoublePoint.
//
// Revision 1.15  1994/03/14  18:44:44  jnw
// Modified sum and product
//
// Revision 1.14  1994/03/14  18:33:16  jnw
// Added point reduce functions 'sum' and 'product'
//
// Revision 1.13  1993/11/16  13:27:10  thoth
// Pointcmp is now pointcmp
//
// Revision 1.12  1993/09/15  12:25:58  thoth
// Point is now a template with specializations only.
//
// Revision 1.11  93/05/28  15:43:31  thoth
// fix comment errors.
// 
// Revision 1.10  1993/05/26  16:59:57  rjj
// Added copyright notice
//
// Revision 1.9  93/04/06  22:28:17  rjj
// Changed name of file from "Point.C.m4" to "Point.c.m4"
// Fixed lattice bugs with != by defining functions to invoke ==
// Changed "const {int|double}" to "{int|double}"
// Fixed a leftover FloatPoint name ref in trunc()
// Changed bitwise unsigned parameters to int
// 
//Revision 1.8  93/04/05  23:18:00  rjj
//Many changes to support changes necessary to mutate to the
//  interface defined by newpoint.  Please see the list for
//  Point.h.m4 revison 2.10
//  Error in operator ~ () corrected
//  Lattice operator != is probably still wrong, need to check !
//  m4 defines *_*_RELATION changed to LATTICE_*_*
//  Postfix ++ and -- code added
//  file rearanged to suit my tastes, your milage my vary ;-)
//
//Revision 1.7  92/12/11  01:10:35  rjj
//Modified to support the IA_ prefix naming convention.
//Also put #ifdef IA_PRII around the #pragma stmt
//
//Revision 1.6  92/09/02  14:51:27  rjj
//Corrected code in FloatPoint( IntPoint ) to allocated array of type int
//Changed enorm() cast from (float) to (double)
//
//Revision 1.5  92/09/01  21:08:17  rjj
//Changed max(IA_Point<int>,IA_Point<int>) to be supremum(IA_Point<int>,IA_Point<int>)
//Changed min(IA_Point<int>,IA_Point<int>) to be infimum(IA_Point<int>,IA_Point<int>)
//
//Revision 1.4  92/09/01  20:46:24  rjj
//Added support to ensure float op IntPoint and IntPoint op float yield FloatPoint.
//
//Revision 1.3  92/08/27  12:11:55  jnw
//Added C style comments around Log Messages to fixx gcc-cpp bug
//
//Revision 1.2  92/08/23  13:33:19  thoth
//Cfrontified
//didn't hint to inline for loop stuff any more (braindead compiler).
//hack-around for switch statements that CFront's flow analysis
//can't fully understand.
//throw changed to ia_throw.
//removed duplicate default values for formal parameters (ARM spec).
//
// Revision 1.1  92/08/20  22:24:54  rjj
// Initial revision
// 
// 
// 

*/

#include <math.h>

#ifdef IA_PRII
#pragma implementation "IntPoint.h"

#endif
#include "IntPoint.h"

#include "point_errors.h"

#ifdef __SC__
#pragma template_access public
#endif

////////////////////////////////////////////////////////////
// C O N S T R U C T O R S
////////////////////////////////////////////////////////////

IA_Point<int> extend_to_point (int initial_value, unsigned length) {
    IA_Point<int> pt;

    pt.coordinates = new int[pt.dimen = length];

#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < length; i++)
	pt.coordinates[i] = initial_value;
#else
    switch (pt.dim()) {
    case 5: pt.coordinates[4] = initial_value;
    case 4: pt.coordinates[3] = initial_value;
    case 3: pt.coordinates[2] = initial_value;
    case 2: pt.coordinates[1] = initial_value;
    case 1: pt.coordinates[0] = initial_value;
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < length; i++)
		    pt.coordinates[i] = initial_value;
	    }
    }
#endif    
    return pt;
}

IA_Point<int>::IA_Point (const IA_Point<int> &ip)
{

// new int[0] creates a pointer to one int; we never use it.
// This is to keep from having to do "if (dim() == 0) ..." everywhere.
// The destructor will delete it, though.
    
    this->coordinates = new int[this->dimen = ip.dim()];

#ifdef NO_LOOP_UNROLLING    
    for (unsigned i=0; i<this->dim(); i++)
	this->coordinates[i] = ip.coordinates[i];
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] = ip.coordinates[4];
    case 4: this->coordinates[3] = ip.coordinates[3];
    case 3: this->coordinates[2] = ip.coordinates[2];
    case 2: this->coordinates[1] = ip.coordinates[1];
    case 1: this->coordinates[0] = ip.coordinates[0];
    case 0: break;
    default:
	    {
		for (unsigned i=0; i<this->dim(); i++)
		    this->coordinates[i] = ip.coordinates[i];
	    }
    }
#endif
}

IA_Point<int>::IA_Point (int *vec, unsigned length, int giveaway)
{
    this->dimen = length;
    if (giveaway) {
	this->coordinates = vec;
	return;
    }
    
    this->coordinates = new int[this->dim()];
    
#ifdef NO_LOOP_UNROLLING    
    for (unsigned length=0; length<this->dim(); length++)
	this->coordinates[length] = vec[length];
#else	
    switch (this->dim()) {
    case 5: this->coordinates[4] = vec[4];
    case 4: this->coordinates[3] = vec[3];
    case 3: this->coordinates[2] = vec[2];
    case 2: this->coordinates[1] = vec[1];
    case 1: this->coordinates[0] = vec[0];
    case 0: break;
    default:
	    {
		for (unsigned length=0; length<this->dim(); length++)
		    this->coordinates[length] = vec[length];
	    }
    }
#endif
}

// IA_Point<int> = concat (IA_Point<int>, IA_Point<int>)
//
// Given a combination of Points and numbers, create a new Point that is
// the concatenation of the arguments.  Lets the user construct larger
// points than the usual 5 arg constructor.

// TO DO -- These could be modified to not check for zero-length
// points (since we now know that "new foo[0]" returns a deletable
// pointer to one foo that we will never access), and also we could
// do some loop unrolling here as well.

IA_Point<int> concat(const IA_Point<int>& lhs, const IA_Point<int>& rhs)
{
    IA_Point<int> result_point;
    
    if (lhs.dim() == 0) {
	if (rhs.dim() == 0)
	    return result_point;
	else 
	    return rhs;
    } else if (rhs.dim() == 0)
	return lhs;

    unsigned int total_dimen = lhs.dim() + rhs.dim();

    result_point.coordinates = new int[result_point.dimen = total_dimen];
    
    for (unsigned int i = 0; i < lhs.dim(); i++) 
	result_point.coordinates[i] = lhs.coordinates[i];

    for (unsigned int j = 0; j < rhs.dim(); j++) 
	result_point.coordinates[i+j] = rhs.coordinates[j];

    return result_point;
}

IA_Point<int> concat(const IA_Point<int>& lhs, int rhs)
{
    IA_Point<int> result_point;

    if (lhs.dim() == 0) 
	return IA_Point<int> (rhs);
    
    unsigned int total_dimen = lhs.dim() + 1;

    result_point.coordinates = new int[total_dimen];
    result_point.dimen = total_dimen;
    
    for (unsigned int i = 0; i < lhs.dim(); i++) 
	result_point.coordinates[i] = lhs.coordinates[i];
    
    result_point.coordinates[i] = rhs;
    
    return result_point;
}

IA_Point<int> concat(int lhs, const IA_Point<int>& rhs)
{
    IA_Point<int> result_point;
    
    if (rhs.dim() == 0)
	return IA_Point<int> (lhs);

    unsigned int total_dimen = rhs.dim() + 1;

    result_point.coordinates = new int[total_dimen];
    result_point.dimen = total_dimen;

    result_point.coordinates[0] = lhs;

    for (unsigned int i = 0; i < rhs.dim(); i++) 
	result_point.coordinates[i+1] = rhs.coordinates[i];

    return result_point;
}

//  cdr(IA_Point<int>, skip=1):
//
//  Given a point of dimension n, return a point made of the n-skip last
//  coordinates (or if skip is negative, the n+skip first coordinates),
//  where skip defaults to 1.

IA_Point<int> cdr(const IA_Point<int>& pt, int skip) {

    unsigned  abs_skip;

    abs_skip = (skip < 0) ? -skip : skip;
    
    if (pt.dim() < abs_skip) {
	ia_throw(Point_BadIndex_Exception(__FILE__, __LINE__));
	return IA_Point<int> ();
    } else if (pt.dim() == skip)
	return IA_Point<int> ();

    if (skip < 0)
	return IA_Point<int> (pt.coordinates, pt.dim()+skip);
    else
	return IA_Point<int> (pt.coordinates+skip, pt.dim()-skip);
}

////////////////////////////////////////////////////////////
// M I S C    O P E R A T I O N S
////////////////////////////////////////////////////////////

IA_Point<int> &IA_Point<int>::operator = (const IA_Point<int> &point) {
    // return directly when self-assignment
    //  (not sure the if is worth it in general)

    if (this == &point)
	return *this;

    // We need to deallocate/reallocate if differing dimens.
    
    if (point.dim() != this->dim()) {
        delete[] this->coordinates;
        this->coordinates = new int[this->dimen = point.dim()];
    }

#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)
	this->coordinates[i] = point.coordinates[i];
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] = point.coordinates[4];
    case 4: this->coordinates[3] = point.coordinates[3];
    case 3: this->coordinates[2] = point.coordinates[2];
    case 2: this->coordinates[1] = point.coordinates[1];
    case 1: this->coordinates[0] = point.coordinates[0];
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)
		    this->coordinates[i] = point.coordinates[i];
	    }
    }
#endif
    return *this;
}

// inf returns an empty IA_Point<int> if passed two empty points.

IA_Point<int> inf (const IA_Point<int>& lhs, const IA_Point<int>& rhs) 
{ 
    if (lhs.dim() != rhs.dim()) {  
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return IA_Point<int> ();  
    }  

    IA_Point<int> p = rhs;

#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < p.dim(); i++) 
	p.coordinates[i]
	    = (lhs.coordinates[i] < rhs.coordinates[i])
		?  lhs.coordinates[i] : rhs.coordinates[i];
#else
    switch (p.dim()) {
    case 5: p.coordinates[4] =	
	(lhs.coordinates[4] < rhs.coordinates[4])
	    ? lhs.coordinates[4] : rhs.coordinates[4];
    case 4: p.coordinates[3] = 
	(lhs.coordinates[3] < rhs.coordinates[3])
	    ? lhs.coordinates[3] : rhs.coordinates[3];
    case 3: p.coordinates[2] = 
	(lhs.coordinates[2] < rhs.coordinates[2])
	    ? lhs.coordinates[2] : rhs.coordinates[2];
    case 2: p.coordinates[1] = 
	(lhs.coordinates[1] < rhs.coordinates[1])
	    ? lhs.coordinates[1] : rhs.coordinates[1];
    case 1: p.coordinates[0] = 
	(lhs.coordinates[0] < rhs.coordinates[0])
	    ? lhs.coordinates[0] : rhs.coordinates[0];
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < p.dim(); i++) 
		    p.coordinates[i]
			= (lhs.coordinates[i] < rhs.coordinates[i])
			    ?  lhs.coordinates[i] : rhs.coordinates[i];
	    }
    }
#endif    
    return p; 
}

// sup returns an empty IA_Point<int> if passed two empty points.

IA_Point<int> sup (const IA_Point<int>& lhs, const IA_Point<int>& rhs) 
{ 
    if (lhs.dim() != rhs.dim()) {  
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return IA_Point<int> ();  
    }  

    IA_Point<int> p = rhs;

#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < p.dim(); i++) 
	p.coordinates[i]
	    = (lhs.coordinates[i] > rhs.coordinates[i])
		?  lhs.coordinates[i] : rhs.coordinates[i];
#else
    switch (p.dim()) {
    case 5: p.coordinates[4] =	
	(lhs.coordinates[4] > rhs.coordinates[4])
	    ? lhs.coordinates[4] : rhs.coordinates[4];
    case 4: p.coordinates[3] = 
	(lhs.coordinates[3] > rhs.coordinates[3])
	    ? lhs.coordinates[3] : rhs.coordinates[3];
    case 3: p.coordinates[2] = 
	(lhs.coordinates[2] > rhs.coordinates[2])
	    ? lhs.coordinates[2] : rhs.coordinates[2];
    case 2: p.coordinates[1] = 
	(lhs.coordinates[1] > rhs.coordinates[1])
	    ? lhs.coordinates[1] : rhs.coordinates[1];
    case 1: p.coordinates[0] = 
	(lhs.coordinates[0] > rhs.coordinates[0])
	    ? lhs.coordinates[0] : rhs.coordinates[0];
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < p.dim(); i++) 
		    p.coordinates[i]
			= (lhs.coordinates[i] > rhs.coordinates[i])
			    ?  lhs.coordinates[i] : rhs.coordinates[i];
	    }
    }
#endif    
    return p; 
}

int min (const IA_Point<int>& point)
{

    if (point.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return 0;
    }

    int result = point.coordinates[0];

#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 1; i < point.dim(); i++) 
	result =  (point.coordinates[i] < result)
			? point.coordinates[i] : result;
#else
    switch (point.dim()) {
    case 5: result = (point.coordinates[4] < result)
			? point.coordinates[4] : result;
    case 4: result = (point.coordinates[3] < result)
			? point.coordinates[3] : result;
    case 3: result = (point.coordinates[2] < result)
			? point.coordinates[2] : result;
    case 2: result = (point.coordinates[1] < result)
			? point.coordinates[1] : result;
    case 1: break;
    default:
	    {
		for (unsigned i = 1; i < point.dim(); i++) 
		    result =  (point.coordinates[i] < result)
				? point.coordinates[i] : result;
	    }
    }
#endif
    return result; 
}

int max (const IA_Point<int>& point)
{

    if (point.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return 0;
    }

    int result = point.coordinates[0];

#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 1; i < point.dim(); i++) 
	result =  (point.coordinates[i] > result)
			? point.coordinates[i] : result;
#else
    switch (point.dim()) {
    case 5: result = (point.coordinates[4] > result)
			? point.coordinates[4] : result;
    case 4: result = (point.coordinates[3] > result)
			? point.coordinates[3] : result;
    case 3: result = (point.coordinates[2] > result)
			? point.coordinates[2] : result;
    case 2: result = (point.coordinates[1] > result)
			? point.coordinates[1] : result;
    case 1: break;
    default:
	    {
		for (unsigned i = 1; i < point.dim(); i++) 
		    result =  (point.coordinates[i] > result)
				? point.coordinates[i] : result;
	    }
    }
#endif
    return result; 
}

ostream &operator << (ostream &stream, const IA_Point<int>& point)
{
    int last_index = point.dim() - 1;
    stream << "<";
    for  ( int i = 0; i < last_index; ++i ) {
	stream << point.coordinates[i] << " ";
    }
    if ( point.dim() > 0 ) {
	stream << point.coordinates[last_index];
    };
    stream << ">";
    return stream;
}

////////////////////////////////////////////////////////////
// I N D E X I N G    O P E R A T I O N S
////////////////////////////////////////////////////////////
// Index derangement: "IA_Point<int> (point_index, point_index..)"
//
// Construct new points from old points, possibly with different
// dimensions and indices rearanged.
// For instance, let  a=Point(10, 20, 30).
// Then b=a.derange(2,0) is equivalent to b=Point(30,10).

IA_Point<int> derange (const IA_Point<int> &pt, unsigned i0)
{
      
    if (i0 >= pt.dim()) {
	ia_throw(Point_BadIndex_Exception(__FILE__, __LINE__));
	return IA_Point<int> ();
    }
    
    int *p = new int [1];

    p[0] = pt(i0);
    return IA_Point<int> ( p, 1, 1);
}

IA_Point<int> derange (const IA_Point<int> &pt, unsigned i0, unsigned i1) 
{
    

    if (i0 >= pt.dim() || i1 >= pt.dim()) {
	ia_throw(Point_BadIndex_Exception(__FILE__, __LINE__));
	return IA_Point<int> ();
    }

    int *p = new int [2];
    
    p[0] = pt(i0);
    p[1] = pt(i1);
    return IA_Point<int> ( p, 2, 1);
}

IA_Point<int> derange (const IA_Point<int> &pt, unsigned i0, unsigned i1,unsigned i2)
{
    

    if (i0 >= pt.dim() || i1 >= pt.dim() || i2 >= pt.dim()) {
	ia_throw(Point_BadIndex_Exception(__FILE__, __LINE__));
	return IA_Point<int> ();
    }

    int *p = new int [3];
    p[0] = pt(i0);
    p[1] = pt(i1);
    p[2] = pt(i2);
    return IA_Point<int> ( p, 3, 1);
}

IA_Point<int> derange (const IA_Point<int> &pt, unsigned i0, unsigned i1, unsigned i2, unsigned i3)
{
    
    if (i0 >= pt.dim() || i1 >= pt.dim() || i2 >= pt.dim() || i3 >= pt.dim()) {
	ia_throw(Point_BadIndex_Exception(__FILE__, __LINE__));
	return IA_Point<int> ();
    }

    int *p = new int [4];
    p[0]= pt(i0);
    p[1] = pt(i1);
    p[2] = pt(i2);
    p[3] = pt(i3);
    return IA_Point<int> ( p, 4, 1);
}

IA_Point<int> derange (const IA_Point<int> &pt, unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4)
{
    
    if (i0 >= pt.dim() || i1 >= pt.dim() || i2 >= pt.dim() ||
	i3 >= pt.dim() || i4 >= pt.dim()) {
	ia_throw(Point_BadIndex_Exception(__FILE__, __LINE__));
	return IA_Point<int> ();
    }

    int *p = new int [5];
    p[0] = pt(i0);
    p[1] = pt(i1);
    p[2] = pt(i2);
    p[3] = pt(i3);
    p[4] = pt(i4);
    return IA_Point<int> ( p, 5, 1);
}

    
IA_Point<int> derange (const IA_Point<int> &pt, unsigned *vec, const unsigned vecLen)
{
    // We need to deallocate/reallocate if differing dimens.
    unsigned i;
    for (i=0; i<vecLen; i++) {
	if (vec[i] >= pt.dim()) {
	    ia_throw(Point_BadIndex_Exception(__FILE__, __LINE__));
	    return IA_Point<int> ();  // An empty point
	}
    }

    int *p = new int [vecLen];

#ifdef NO_LOOP_UNROLLING
    for (i = 0; i < vecLen; i++)
	p[i] = pt(vec[i]);
#else
    switch (vecLen) {
    case 5: p[4] = pt(vec[4]);
    case 4: p[3] = pt(vec[3]);
    case 3: p[2] = pt(vec[2]);
    case 2: p[1] = pt(vec[1]);
    case 1: p[0] = pt(vec[0]);
    case 0: break;
    default:
	    {
		for (i = 0; i < vecLen; i++)
		    p[i] = pt(vec[i]);
	    }
    }
#endif
    return IA_Point<int> ( p, vecLen, 1);
}

  inline static int IA_abs ( int t)
{
  t = (t < 0) ? -t : t;
  return t;
}

// abs will work on double and int points, is a friend of point

IA_Point<int> abs (const IA_Point<int> &p)
{
    
  if (p.dim() == 0) {
    ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;
    return IA_Point<int> (); //if invalid return a point of 0 dim
  }

  int *result = new int[p.dim()];
  
#ifdef NO_LOOP_UNROLLING
  for (unsigned i = 0; i < p.dim(); i++)	
    result[i] = IA_abs(p(i));
#else
  switch (p.dim()) {
  case 5: result[4] =IA_abs(p(4));  
  case 4: result[3] =IA_abs(p(3));  
  case 3: result[2] =IA_abs(p(2));  
  case 2: result[1] =IA_abs(p(1));
  case 1: result[0] =IA_abs(p(0));  
  case 0: break;
  default:
    {
      for (unsigned i = 0; i < p.dim(); i++)
	result[i] = IA_abs(p(i));
    }
  
  }
#endif 
  return IA_Point<int> (result, p.dim(), 1);
}

// fabs to follow stdlib.h will only work on double

// inorm and enorm uses abs() for both points of double and ints
// enorm: Euclidean distance from origin.

double enorm (const IA_Point<int> &p)
{
    double result_number = 0;

    if (p.dim() == 0)
      {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;
	return result_number;
      }

#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < p.dim(); i++)
      result_number = result_number +
	(double) (p(i) * p(i));
#else
    switch (p.dim()) {
    case 5: result_number += (double) (p(4) * p(4));
    case 4: result_number += (double) (p(3) * p(3));
    case 3: result_number += (double) (p(2) * p(2));
    case 2: result_number += (double) (p(1) * p(1));
    case 1: result_number += (double) (p(0) * p(0));
    case 0: break;
    default:
      {
	for (unsigned i = 0; i < p.dim(); i++)
	  result_number += (double) (p(i) * p(i));
      }
    }
#endif
    return sqrt(result_number);
}

// inorm: infinity norm , checkerboard distance

int inorm(const IA_Point<int> &p)
{
  int result = 0;
  
  if (p.dim() == 0) {
       ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;
       return result;
   }
  result =  max(abs(p));
  return result;
}

// mnorm: manhattan distance

int mnorm(const IA_Point<int> &p)
{
  int result = 0;
  
  if (p.dim() == 0) {
       ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;
       return result;
    }
  result = sum(abs(p));
  return result;
}

////////////////////////////////////////////////////////////
// A R I T H M E T I C   O P E R A T I O N S
////////////////////////////////////////////////////////////

  
IA_Point<int>& IA_Point<int>::operator += (const IA_Point<int>& rhs) {
    if (this->dim() != rhs.dim()) {  
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return *this;  
    }  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] += rhs.coordinates[i];  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] += rhs.coordinates[4];  
    case 4: this->coordinates[3] += rhs.coordinates[3];  
    case 3: this->coordinates[2] += rhs.coordinates[2];  
    case 2: this->coordinates[1] += rhs.coordinates[1];  
    case 1: this->coordinates[0] += rhs.coordinates[0];  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    coordinates[i] += rhs.coordinates[i];  
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator -= (const IA_Point<int>& rhs) {
    if (this->dim() != rhs.dim()) {  
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return *this;  
    }  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] -= rhs.coordinates[i];  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] -= rhs.coordinates[4];  
    case 4: this->coordinates[3] -= rhs.coordinates[3];  
    case 3: this->coordinates[2] -= rhs.coordinates[2];  
    case 2: this->coordinates[1] -= rhs.coordinates[1];  
    case 1: this->coordinates[0] -= rhs.coordinates[0];  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    coordinates[i] -= rhs.coordinates[i];  
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator *= (const IA_Point<int>& rhs) {
    if (this->dim() != rhs.dim()) {  
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return *this;  
    }  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] *= rhs.coordinates[i];  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] *= rhs.coordinates[4];  
    case 4: this->coordinates[3] *= rhs.coordinates[3];  
    case 3: this->coordinates[2] *= rhs.coordinates[2];  
    case 2: this->coordinates[1] *= rhs.coordinates[1];  
    case 1: this->coordinates[0] *= rhs.coordinates[0];  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    coordinates[i] *= rhs.coordinates[i];  
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator /= (const IA_Point<int>& rhs) {
    if (this->dim() != rhs.dim()) {  
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return *this;  
    }  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] /= rhs.coordinates[i];  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] /= rhs.coordinates[4];  
    case 4: this->coordinates[3] /= rhs.coordinates[3];  
    case 3: this->coordinates[2] /= rhs.coordinates[2];  
    case 2: this->coordinates[1] /= rhs.coordinates[1];  
    case 1: this->coordinates[0] /= rhs.coordinates[0];  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    coordinates[i] /= rhs.coordinates[i];  
	    }
    }
#endif
    return *this;  
}

  
IA_Point<int>& IA_Point<int>::operator += (int rhs) { 
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] += rhs;  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] += rhs;  
    case 4: this->coordinates[3] += rhs;  
    case 3: this->coordinates[2] += rhs;  
    case 2: this->coordinates[1] += rhs;  
    case 1: this->coordinates[0] += rhs;  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] += rhs;  
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator -= (int rhs) { 
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] -= rhs;  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] -= rhs;  
    case 4: this->coordinates[3] -= rhs;  
    case 3: this->coordinates[2] -= rhs;  
    case 2: this->coordinates[1] -= rhs;  
    case 1: this->coordinates[0] -= rhs;  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] -= rhs;  
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator *= (int rhs) { 
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] *= rhs;  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] *= rhs;  
    case 4: this->coordinates[3] *= rhs;  
    case 3: this->coordinates[2] *= rhs;  
    case 2: this->coordinates[1] *= rhs;  
    case 1: this->coordinates[0] *= rhs;  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] *= rhs;  
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator /= (int rhs) { 
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] /= rhs;  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] /= rhs;  
    case 4: this->coordinates[3] /= rhs;  
    case 3: this->coordinates[2] /= rhs;  
    case 2: this->coordinates[1] /= rhs;  
    case 1: this->coordinates[0] /= rhs;  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] /= rhs;  
	    }
    }
#endif
    return *this;  
}
    

IA_Point<int> IA_Point<int>::operator + (const IA_Point<int>& rhs) const  
{ 
   IA_Point<int> result_point = *this;    
   return result_point += rhs;  
}

IA_Point<int> IA_Point<int>::operator - (const IA_Point<int>& rhs) const  
{ 
   IA_Point<int> result_point = *this;    
   return result_point -= rhs;  
}

IA_Point<int> IA_Point<int>::operator * (const IA_Point<int>& rhs) const  
{ 
   IA_Point<int> result_point = *this;    
   return result_point *= rhs;  
}

IA_Point<int> IA_Point<int>::operator / (const IA_Point<int>& rhs) const  
{ 
   IA_Point<int> result_point = *this;    
   return result_point /= rhs;  
}
    

IA_Point<int> operator + (int lhs, const IA_Point<int>& rhs)  
{ 
  IA_Point<int> result_point = extend_to_point (lhs, rhs.dim()); 
  return result_point += rhs;  
}

IA_Point<int> operator - (int lhs, const IA_Point<int>& rhs)  
{ 
  IA_Point<int> result_point = extend_to_point (lhs, rhs.dim()); 
  return result_point -= rhs;  
}

IA_Point<int> operator * (int lhs, const IA_Point<int>& rhs)  
{ 
  IA_Point<int> result_point = extend_to_point (lhs, rhs.dim()); 
  return result_point *= rhs;  
}

IA_Point<int> operator / (int lhs, const IA_Point<int>& rhs)  
{ 
  IA_Point<int> result_point = extend_to_point (lhs, rhs.dim()); 
  return result_point /= rhs;  
}
    

IA_Point<int> IA_Point<int>::operator + (int rhs) const
{ 
    IA_Point<int> result_point = *this;
    return result_point += rhs; 
}

IA_Point<int> IA_Point<int>::operator - (int rhs) const
{ 
    IA_Point<int> result_point = *this;
    return result_point -= rhs; 
}

IA_Point<int> IA_Point<int>::operator * (int rhs) const
{ 
    IA_Point<int> result_point = *this;
    return result_point *= rhs; 
}

IA_Point<int> IA_Point<int>::operator / (int rhs) const
{ 
    IA_Point<int> result_point = *this;
    return result_point /= rhs; 
}
    

IA_Point<int> IA_Point<int>::operator - () const
{
    IA_Point<int> result_point;

    if (this->dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;
	return result_point;
    }

    result_point.coordinates = new int[result_point.dimen = this->dim()];

#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < result_point.dim(); i++)
	result_point.coordinates[i] =  -this->coordinates[i];
#else
    switch (result_point.dim()) {
    case 5: result_point.coordinates[4] = -this->coordinates[4];
    case 4: result_point.coordinates[3] = -this->coordinates[3];
    case 3: result_point.coordinates[2] = -this->coordinates[2];
    case 2: result_point.coordinates[1] = -this->coordinates[1];
    case 1: result_point.coordinates[0] = -this->coordinates[0];
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < result_point.dim(); i++)
		    result_point.coordinates[i] =  - this->coordinates[i];
	    }
    }
#endif
    return result_point;
}

// Prefix ++
IA_Point<int>& IA_Point<int>::operator ++ () {
    if (this->dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;
	return *this;
    }
    *this += 1;
    return  *this;
}

// Prefix --
IA_Point<int>& IA_Point<int>::operator -- () {
    if (this->dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;
	return *this;
    }
    *this -= 1;
    return  *this;
}

// Postfix ++
IA_Point<int>  IA_Point<int>::operator ++ (int) {
    if (this->dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;
	return *this;
    }
    IA_Point<int>  result = *this;
    *this += 1;
    return  result;
}

// Postfix --
IA_Point<int>  IA_Point<int>::operator -- (int) {
    if (this->dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;
	return *this;
    }
    IA_Point<int>  result = *this;
    *this -= 1;
    return  result;
}

// Point reduce operations

int sum(const IA_Point<int> &p)
{
    int result;
    
    if (p.dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;
	return result;
    }
    
    result = p[0];
    
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 1; i < p.dim(); i++)
	result += p[i];
#else
    switch (p.dim()) {
    case 5: result += p[4];
    case 4: result += p[3];
    case 3: result += p[2];
    case 2: result += p[1];
    case 1: break;
    default:
	{
	    for (unsigned i = 1; i < p.dim(); i++)
		result += p[i];
	}
    }
#endif
    return result;
}

int product(const IA_Point<int> &p)
{
    int result;
    
    if (p.dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;
	return result;
    }
    
    result = p[0];
    
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 1; i < p.dim(); i++)
	result *= p[i];
#else
    switch (p.dim()) {
    case 5: result *= p[4];
    case 4: result *= p[3];
    case 3: result *= p[2];
    case 2: result *= p[1];
    case 1: break;
    default:
	{
	    for (unsigned i = 1; i < p.dim(); i++)
		result *= p[i];
	}
    }
#endif
    return result;
}

////////////////////////////////////////////////////////////
// L A T T I C E   O P E R A T I O N S
////////////////////////////////////////////////////////////

// Relational operators

// IA_Point<int> < IA_Point<int>    

 
int IA_Point<int>::operator < (const IA_Point<int>& rhs) const  
{  
// this is suspicious -XXX-
    if (this->dim() == 0 || rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  
    if (this->dim() != rhs.dim())  
	return  0;  

#ifdef NO_LOOP_UNROLLING    
    for (unsigned i = 0; i < this->dim() ; i++)
	if (! (this->coordinates[i] < rhs.coordinates[i])) return 0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] < rhs.coordinates[4])) return 0;  
    case 4: if (! (this->coordinates[3] < rhs.coordinates[3])) return 0;  
    case 3: if (! (this->coordinates[2] < rhs.coordinates[2])) return 0;  
    case 2: if (! (this->coordinates[1] < rhs.coordinates[1])) return 0;  
    case 1: if (! (this->coordinates[0] < rhs.coordinates[0])) return 0;  
    case 0: break;
    default: 
	    {
		for (unsigned i = 0; i < this->dim() ; i++)   
		    if (! (this->coordinates[i] < rhs.coordinates[i]))
			return 0;  
	    }
    }
#endif
    return 1;  
}
 
int IA_Point<int>::operator <= (const IA_Point<int>& rhs) const  
{  
// this is suspicious -XXX-
    if (this->dim() == 0 || rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  
    if (this->dim() != rhs.dim())  
	return  0;  

#ifdef NO_LOOP_UNROLLING    
    for (unsigned i = 0; i < this->dim() ; i++)
	if (! (this->coordinates[i] <= rhs.coordinates[i])) return 0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] <= rhs.coordinates[4])) return 0;  
    case 4: if (! (this->coordinates[3] <= rhs.coordinates[3])) return 0;  
    case 3: if (! (this->coordinates[2] <= rhs.coordinates[2])) return 0;  
    case 2: if (! (this->coordinates[1] <= rhs.coordinates[1])) return 0;  
    case 1: if (! (this->coordinates[0] <= rhs.coordinates[0])) return 0;  
    case 0: break;
    default: 
	    {
		for (unsigned i = 0; i < this->dim() ; i++)   
		    if (! (this->coordinates[i] <= rhs.coordinates[i]))
			return 0;  
	    }
    }
#endif
    return 1;  
}
 
int IA_Point<int>::operator > (const IA_Point<int>& rhs) const  
{  
// this is suspicious -XXX-
    if (this->dim() == 0 || rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  
    if (this->dim() != rhs.dim())  
	return  0;  

#ifdef NO_LOOP_UNROLLING    
    for (unsigned i = 0; i < this->dim() ; i++)
	if (! (this->coordinates[i] > rhs.coordinates[i])) return 0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] > rhs.coordinates[4])) return 0;  
    case 4: if (! (this->coordinates[3] > rhs.coordinates[3])) return 0;  
    case 3: if (! (this->coordinates[2] > rhs.coordinates[2])) return 0;  
    case 2: if (! (this->coordinates[1] > rhs.coordinates[1])) return 0;  
    case 1: if (! (this->coordinates[0] > rhs.coordinates[0])) return 0;  
    case 0: break;
    default: 
	    {
		for (unsigned i = 0; i < this->dim() ; i++)   
		    if (! (this->coordinates[i] > rhs.coordinates[i]))
			return 0;  
	    }
    }
#endif
    return 1;  
}
 
int IA_Point<int>::operator >= (const IA_Point<int>& rhs) const  
{  
// this is suspicious -XXX-
    if (this->dim() == 0 || rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  
    if (this->dim() != rhs.dim())  
	return  0;  

#ifdef NO_LOOP_UNROLLING    
    for (unsigned i = 0; i < this->dim() ; i++)
	if (! (this->coordinates[i] >= rhs.coordinates[i])) return 0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] >= rhs.coordinates[4])) return 0;  
    case 4: if (! (this->coordinates[3] >= rhs.coordinates[3])) return 0;  
    case 3: if (! (this->coordinates[2] >= rhs.coordinates[2])) return 0;  
    case 2: if (! (this->coordinates[1] >= rhs.coordinates[1])) return 0;  
    case 1: if (! (this->coordinates[0] >= rhs.coordinates[0])) return 0;  
    case 0: break;
    default: 
	    {
		for (unsigned i = 0; i < this->dim() ; i++)   
		    if (! (this->coordinates[i] >= rhs.coordinates[i]))
			return 0;  
	    }
    }
#endif
    return 1;  
}
 
int IA_Point<int>::operator == (const IA_Point<int>& rhs) const  
{  
// this is suspicious -XXX-
    if (this->dim() == 0 || rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  
    if (this->dim() != rhs.dim())  
	return  0;  

#ifdef NO_LOOP_UNROLLING    
    for (unsigned i = 0; i < this->dim() ; i++)
	if (! (this->coordinates[i] == rhs.coordinates[i])) return 0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] == rhs.coordinates[4])) return 0;  
    case 4: if (! (this->coordinates[3] == rhs.coordinates[3])) return 0;  
    case 3: if (! (this->coordinates[2] == rhs.coordinates[2])) return 0;  
    case 2: if (! (this->coordinates[1] == rhs.coordinates[1])) return 0;  
    case 1: if (! (this->coordinates[0] == rhs.coordinates[0])) return 0;  
    case 0: break;
    default: 
	    {
		for (unsigned i = 0; i < this->dim() ; i++)   
		    if (! (this->coordinates[i] == rhs.coordinates[i]))
			return 0;  
	    }
    }
#endif
    return 1;  
}

int IA_Point<int>::operator != (const IA_Point<int>& rhs) const
{
    return ! (*this == rhs);
}

// IA_Point<int> < int

int  IA_Point<int>::operator < (int rhs) const
{  
    if (this->dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < this->dim(); i++)  
	if (! (this->coordinates[i] < rhs)) return  0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] < rhs)) return  0;  
    case 4: if (! (this->coordinates[3] < rhs)) return  0;  
    case 3: if (! (this->coordinates[2] < rhs)) return  0;  
    case 2: if (! (this->coordinates[1] < rhs)) return  0;  
    case 1: if (! (this->coordinates[0] < rhs)) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < this->dim(); i++)  
		    if (! (this->coordinates[i] < rhs)) return  0;  
	    }
    }
#endif  
    return  1;  
}

int  IA_Point<int>::operator <= (int rhs) const
{  
    if (this->dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < this->dim(); i++)  
	if (! (this->coordinates[i] <= rhs)) return  0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] <= rhs)) return  0;  
    case 4: if (! (this->coordinates[3] <= rhs)) return  0;  
    case 3: if (! (this->coordinates[2] <= rhs)) return  0;  
    case 2: if (! (this->coordinates[1] <= rhs)) return  0;  
    case 1: if (! (this->coordinates[0] <= rhs)) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < this->dim(); i++)  
		    if (! (this->coordinates[i] <= rhs)) return  0;  
	    }
    }
#endif  
    return  1;  
}

int  IA_Point<int>::operator > (int rhs) const
{  
    if (this->dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < this->dim(); i++)  
	if (! (this->coordinates[i] > rhs)) return  0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] > rhs)) return  0;  
    case 4: if (! (this->coordinates[3] > rhs)) return  0;  
    case 3: if (! (this->coordinates[2] > rhs)) return  0;  
    case 2: if (! (this->coordinates[1] > rhs)) return  0;  
    case 1: if (! (this->coordinates[0] > rhs)) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < this->dim(); i++)  
		    if (! (this->coordinates[i] > rhs)) return  0;  
	    }
    }
#endif  
    return  1;  
}

int  IA_Point<int>::operator >= (int rhs) const
{  
    if (this->dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < this->dim(); i++)  
	if (! (this->coordinates[i] >= rhs)) return  0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] >= rhs)) return  0;  
    case 4: if (! (this->coordinates[3] >= rhs)) return  0;  
    case 3: if (! (this->coordinates[2] >= rhs)) return  0;  
    case 2: if (! (this->coordinates[1] >= rhs)) return  0;  
    case 1: if (! (this->coordinates[0] >= rhs)) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < this->dim(); i++)  
		    if (! (this->coordinates[i] >= rhs)) return  0;  
	    }
    }
#endif  
    return  1;  
}

int  IA_Point<int>::operator == (int rhs) const
{  
    if (this->dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < this->dim(); i++)  
	if (! (this->coordinates[i] == rhs)) return  0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] == rhs)) return  0;  
    case 4: if (! (this->coordinates[3] == rhs)) return  0;  
    case 3: if (! (this->coordinates[2] == rhs)) return  0;  
    case 2: if (! (this->coordinates[1] == rhs)) return  0;  
    case 1: if (! (this->coordinates[0] == rhs)) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < this->dim(); i++)  
		    if (! (this->coordinates[i] == rhs)) return  0;  
	    }
    }
#endif  
    return  1;  
}

int  IA_Point<int>::operator != (int rhs) const
{
    return ! (*this == rhs);
}

// int < IA_Point<int> 

 
int operator < (int lhs, const IA_Point<int>& rhs)  
{  
    if (rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < rhs.dim(); i++)  
	if (! (lhs < rhs.coordinates[i])) return  0;  
#else
    switch (rhs.dim()) {
    case 5: if (! (lhs < rhs.coordinates[4])) return  0;  
    case 4: if (! (lhs < rhs.coordinates[3])) return  0;  
    case 3: if (! (lhs < rhs.coordinates[2])) return  0;  
    case 2: if (! (lhs < rhs.coordinates[1])) return  0;  
    case 1: if (! (lhs < rhs.coordinates[0])) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < rhs.dim(); i++)  
		    if (! (lhs < rhs.coordinates[i])) return  0;  
	    }
    }
#endif  

    return  1;  
}
 
int operator <= (int lhs, const IA_Point<int>& rhs)  
{  
    if (rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < rhs.dim(); i++)  
	if (! (lhs <= rhs.coordinates[i])) return  0;  
#else
    switch (rhs.dim()) {
    case 5: if (! (lhs <= rhs.coordinates[4])) return  0;  
    case 4: if (! (lhs <= rhs.coordinates[3])) return  0;  
    case 3: if (! (lhs <= rhs.coordinates[2])) return  0;  
    case 2: if (! (lhs <= rhs.coordinates[1])) return  0;  
    case 1: if (! (lhs <= rhs.coordinates[0])) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < rhs.dim(); i++)  
		    if (! (lhs <= rhs.coordinates[i])) return  0;  
	    }
    }
#endif  

    return  1;  
}
 
int operator > (int lhs, const IA_Point<int>& rhs)  
{  
    if (rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < rhs.dim(); i++)  
	if (! (lhs > rhs.coordinates[i])) return  0;  
#else
    switch (rhs.dim()) {
    case 5: if (! (lhs > rhs.coordinates[4])) return  0;  
    case 4: if (! (lhs > rhs.coordinates[3])) return  0;  
    case 3: if (! (lhs > rhs.coordinates[2])) return  0;  
    case 2: if (! (lhs > rhs.coordinates[1])) return  0;  
    case 1: if (! (lhs > rhs.coordinates[0])) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < rhs.dim(); i++)  
		    if (! (lhs > rhs.coordinates[i])) return  0;  
	    }
    }
#endif  

    return  1;  
}
 
int operator >= (int lhs, const IA_Point<int>& rhs)  
{  
    if (rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < rhs.dim(); i++)  
	if (! (lhs >= rhs.coordinates[i])) return  0;  
#else
    switch (rhs.dim()) {
    case 5: if (! (lhs >= rhs.coordinates[4])) return  0;  
    case 4: if (! (lhs >= rhs.coordinates[3])) return  0;  
    case 3: if (! (lhs >= rhs.coordinates[2])) return  0;  
    case 2: if (! (lhs >= rhs.coordinates[1])) return  0;  
    case 1: if (! (lhs >= rhs.coordinates[0])) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < rhs.dim(); i++)  
		    if (! (lhs >= rhs.coordinates[i])) return  0;  
	    }
    }
#endif  

    return  1;  
}
 
int operator == (int lhs, const IA_Point<int>& rhs)  
{  
    if (rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < rhs.dim(); i++)  
	if (! (lhs == rhs.coordinates[i])) return  0;  
#else
    switch (rhs.dim()) {
    case 5: if (! (lhs == rhs.coordinates[4])) return  0;  
    case 4: if (! (lhs == rhs.coordinates[3])) return  0;  
    case 3: if (! (lhs == rhs.coordinates[2])) return  0;  
    case 2: if (! (lhs == rhs.coordinates[1])) return  0;  
    case 1: if (! (lhs == rhs.coordinates[0])) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < rhs.dim(); i++)  
		    if (! (lhs == rhs.coordinates[i])) return  0;  
	    }
    }
#endif  

    return  1;  
}

int operator != (int lhs, const IA_Point<int>& rhs)  
{
    return !(rhs == lhs);
}

// Intpoint op double
 
int IA_Point<int>::operator <  (double rhs) const
{  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < this->dim(); i++)  
	if (! (this->coordinates[i] < rhs)) return  0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] < rhs)) return  0;  
    case 4: if (! (this->coordinates[3] < rhs)) return  0;  
    case 3: if (! (this->coordinates[2] < rhs)) return  0;  
    case 2: if (! (this->coordinates[1] < rhs)) return  0;  
    case 1: if (! (this->coordinates[0] < rhs)) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < this->dim(); i++)  
		    if (! (this->coordinates[i] < rhs)) return  0;  
	    }
    }
#endif
    return  1;  
}
 
int IA_Point<int>::operator <=  (double rhs) const
{  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < this->dim(); i++)  
	if (! (this->coordinates[i] <= rhs)) return  0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] <= rhs)) return  0;  
    case 4: if (! (this->coordinates[3] <= rhs)) return  0;  
    case 3: if (! (this->coordinates[2] <= rhs)) return  0;  
    case 2: if (! (this->coordinates[1] <= rhs)) return  0;  
    case 1: if (! (this->coordinates[0] <= rhs)) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < this->dim(); i++)  
		    if (! (this->coordinates[i] <= rhs)) return  0;  
	    }
    }
#endif
    return  1;  
}
 
int IA_Point<int>::operator >  (double rhs) const
{  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < this->dim(); i++)  
	if (! (this->coordinates[i] > rhs)) return  0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] > rhs)) return  0;  
    case 4: if (! (this->coordinates[3] > rhs)) return  0;  
    case 3: if (! (this->coordinates[2] > rhs)) return  0;  
    case 2: if (! (this->coordinates[1] > rhs)) return  0;  
    case 1: if (! (this->coordinates[0] > rhs)) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < this->dim(); i++)  
		    if (! (this->coordinates[i] > rhs)) return  0;  
	    }
    }
#endif
    return  1;  
}
 
int IA_Point<int>::operator >=  (double rhs) const
{  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < this->dim(); i++)  
	if (! (this->coordinates[i] >= rhs)) return  0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] >= rhs)) return  0;  
    case 4: if (! (this->coordinates[3] >= rhs)) return  0;  
    case 3: if (! (this->coordinates[2] >= rhs)) return  0;  
    case 2: if (! (this->coordinates[1] >= rhs)) return  0;  
    case 1: if (! (this->coordinates[0] >= rhs)) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < this->dim(); i++)  
		    if (! (this->coordinates[i] >= rhs)) return  0;  
	    }
    }
#endif
    return  1;  
}
 
int IA_Point<int>::operator ==  (double rhs) const
{  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < this->dim(); i++)  
	if (! (this->coordinates[i] == rhs)) return  0;  
#else
    switch (this->dim()) {
    case 5: if (! (this->coordinates[4] == rhs)) return  0;  
    case 4: if (! (this->coordinates[3] == rhs)) return  0;  
    case 3: if (! (this->coordinates[2] == rhs)) return  0;  
    case 2: if (! (this->coordinates[1] == rhs)) return  0;  
    case 1: if (! (this->coordinates[0] == rhs)) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < this->dim(); i++)  
		    if (! (this->coordinates[i] == rhs)) return  0;  
	    }
    }
#endif
    return  1;  
}

int IA_Point<int>::operator !=  (double rhs) const
{
    return ! (*this == rhs);
}

// double < 

int operator < (double lhs, const IA_Point<int>& rhs)  
{  
    if (rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < rhs.dim(); i++)  
	if (! (lhs < rhs.coordinates[i]))  
	    return  0;  

#else
    switch (rhs.dim()) {
    case 5: if (! (lhs < rhs.coordinates[4])) return  0;  
    case 4: if (! (lhs < rhs.coordinates[3])) return  0;  
    case 3: if (! (lhs < rhs.coordinates[2])) return  0;  
    case 2: if (! (lhs < rhs.coordinates[1])) return  0;  
    case 1: if (! (lhs < rhs.coordinates[0])) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < rhs.dim(); i++)  
		    if (! (lhs < rhs.coordinates[i])) return  0;  
	    }
    }
#endif    
    return  1;  
}

int operator <= (double lhs, const IA_Point<int>& rhs)  
{  
    if (rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < rhs.dim(); i++)  
	if (! (lhs <= rhs.coordinates[i]))  
	    return  0;  

#else
    switch (rhs.dim()) {
    case 5: if (! (lhs <= rhs.coordinates[4])) return  0;  
    case 4: if (! (lhs <= rhs.coordinates[3])) return  0;  
    case 3: if (! (lhs <= rhs.coordinates[2])) return  0;  
    case 2: if (! (lhs <= rhs.coordinates[1])) return  0;  
    case 1: if (! (lhs <= rhs.coordinates[0])) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < rhs.dim(); i++)  
		    if (! (lhs <= rhs.coordinates[i])) return  0;  
	    }
    }
#endif    
    return  1;  
}

int operator > (double lhs, const IA_Point<int>& rhs)  
{  
    if (rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < rhs.dim(); i++)  
	if (! (lhs > rhs.coordinates[i]))  
	    return  0;  

#else
    switch (rhs.dim()) {
    case 5: if (! (lhs > rhs.coordinates[4])) return  0;  
    case 4: if (! (lhs > rhs.coordinates[3])) return  0;  
    case 3: if (! (lhs > rhs.coordinates[2])) return  0;  
    case 2: if (! (lhs > rhs.coordinates[1])) return  0;  
    case 1: if (! (lhs > rhs.coordinates[0])) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < rhs.dim(); i++)  
		    if (! (lhs > rhs.coordinates[i])) return  0;  
	    }
    }
#endif    
    return  1;  
}

int operator >= (double lhs, const IA_Point<int>& rhs)  
{  
    if (rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < rhs.dim(); i++)  
	if (! (lhs >= rhs.coordinates[i]))  
	    return  0;  

#else
    switch (rhs.dim()) {
    case 5: if (! (lhs >= rhs.coordinates[4])) return  0;  
    case 4: if (! (lhs >= rhs.coordinates[3])) return  0;  
    case 3: if (! (lhs >= rhs.coordinates[2])) return  0;  
    case 2: if (! (lhs >= rhs.coordinates[1])) return  0;  
    case 1: if (! (lhs >= rhs.coordinates[0])) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < rhs.dim(); i++)  
		    if (! (lhs >= rhs.coordinates[i])) return  0;  
	    }
    }
#endif    
    return  1;  
}

int operator == (double lhs, const IA_Point<int>& rhs)  
{  
    if (rhs.dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return  0;  
    }	  

#ifdef NO_LOOP_UNROLLING
    for (unsigned int i = 0; i < rhs.dim(); i++)  
	if (! (lhs == rhs.coordinates[i]))  
	    return  0;  

#else
    switch (rhs.dim()) {
    case 5: if (! (lhs == rhs.coordinates[4])) return  0;  
    case 4: if (! (lhs == rhs.coordinates[3])) return  0;  
    case 3: if (! (lhs == rhs.coordinates[2])) return  0;  
    case 2: if (! (lhs == rhs.coordinates[1])) return  0;  
    case 1: if (! (lhs == rhs.coordinates[0])) return  0;
    case 0: break;
    default:
	    {
		for (unsigned int i = 0; i < rhs.dim(); i++)  
		    if (! (lhs == rhs.coordinates[i])) return  0;  
	    }
    }
#endif    
    return  1;  
}

int operator != (double lhs, const IA_Point<int>& rhs)  
{
    return ! (rhs == lhs);
}

//  Point comparison: pointcmp(IA_Point<int>, IA_Point<int>)
//
// A sort of lexicographic comparison of Points, except that
// points are zero extended for comparison, e.g., <1 0> is equal
// to <1>;  <> is equal <0 0 0>, etc. Returns 0 for equality, -1
// if the first is less than the second, 1 if the first is greater
// than the second,

int pointcmp(const IA_Point<int> &lhs, const IA_Point<int> &rhs)
{
    if (lhs.dim() > rhs.dim()) {
	for (unsigned i=0; i < rhs.dim(); i++) {
	    if (lhs[i] < rhs[i])
		return -1;
	    else if (lhs[i] > rhs[i])
		return 1;
	}
	for ( ; i<lhs.dim(); i++) {
	    if (lhs[i] < 0)
		return -1;
	    else if (lhs[i] > 0)
		return 1;
	}
    } else {
	for (unsigned i=0; i < lhs.dim(); i++) {
	    if (lhs[i] < rhs[i])
		return -1;
	    else if (lhs[i] > rhs[i])
		return 1;
	}
	for ( ; i < rhs.dim(); i++) {
	    if (0 < rhs[i])
		return  -1;
	    else if (0 > rhs[i])
		return  1;
	}
    }
    return 0;

    // If this were to be used for a lexicographic sorter, then we would
    // need the following additional code:
    //
    // At this point the zero-extended Coordinates do not differ in value,
    // so check for number of dimensions; shorter Points are considered
    // less than longer Points.
    // 
    //     if (lhs.dim() == rhs.dim())
    // 		return 0;
    //     else if (lhs.dim() > rhs.dim())
    // 		return 1;
    //     else
    // 		return -1;

}

////////////////////////////////////////////////////////////
// B I T W I S E   O P E R A T I O N S
////////////////////////////////////////////////////////////

IA_Point<int> IA_Point<int>::operator ~ () const
{
    IA_Point<int> result_point;
    
    if (this->dim() == 0) {
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;
	return *this;
    }

    result_point.coordinates = new int[result_point.dimen = this->dim()];

#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < result_point.dim(); i++)
	result_point.coordinates[i] =  ~(this->coordinates[i]);
#else
    switch (result_point.dim()) {
    case 5: result_point.coordinates[4] = ~(this->coordinates[4]);
    case 4: result_point.coordinates[3] = ~(this->coordinates[3]);
    case 3: result_point.coordinates[2] = ~(this->coordinates[2]);
    case 2: result_point.coordinates[1] = ~(this->coordinates[1]);
    case 1: result_point.coordinates[0] = ~(this->coordinates[0]);
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < result_point.dim(); i++)
		    result_point.coordinates[i] = ~(this->coordinates[i]);
	    }
    }
#endif    
    return result_point;
}

  
IA_Point<int>& IA_Point<int>::operator %= (const IA_Point<int>& rhs) {  
    if (this->dim() != rhs.dim()) {  
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return *this;  
    }  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] %= rhs.coordinates[i];  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] %= rhs.coordinates[4];  
    case 4: this->coordinates[3] %= rhs.coordinates[3];  
    case 3: this->coordinates[2] %= rhs.coordinates[2];  
    case 2: this->coordinates[1] %= rhs.coordinates[1];  
    case 1: this->coordinates[0] %= rhs.coordinates[0];  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] %= rhs.coordinates[i];  
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator <<= (const IA_Point<int>& rhs) {  
    if (this->dim() != rhs.dim()) {  
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return *this;  
    }  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] <<= rhs.coordinates[i];  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] <<= rhs.coordinates[4];  
    case 4: this->coordinates[3] <<= rhs.coordinates[3];  
    case 3: this->coordinates[2] <<= rhs.coordinates[2];  
    case 2: this->coordinates[1] <<= rhs.coordinates[1];  
    case 1: this->coordinates[0] <<= rhs.coordinates[0];  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] <<= rhs.coordinates[i];  
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator >>= (const IA_Point<int>& rhs) {  
    if (this->dim() != rhs.dim()) {  
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return *this;  
    }  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] >>= rhs.coordinates[i];  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] >>= rhs.coordinates[4];  
    case 4: this->coordinates[3] >>= rhs.coordinates[3];  
    case 3: this->coordinates[2] >>= rhs.coordinates[2];  
    case 2: this->coordinates[1] >>= rhs.coordinates[1];  
    case 1: this->coordinates[0] >>= rhs.coordinates[0];  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] >>= rhs.coordinates[i];  
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator &= (const IA_Point<int>& rhs) {  
    if (this->dim() != rhs.dim()) {  
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return *this;  
    }  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] &= rhs.coordinates[i];  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] &= rhs.coordinates[4];  
    case 4: this->coordinates[3] &= rhs.coordinates[3];  
    case 3: this->coordinates[2] &= rhs.coordinates[2];  
    case 2: this->coordinates[1] &= rhs.coordinates[1];  
    case 1: this->coordinates[0] &= rhs.coordinates[0];  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] &= rhs.coordinates[i];  
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator ^= (const IA_Point<int>& rhs) {  
    if (this->dim() != rhs.dim()) {  
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return *this;  
    }  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] ^= rhs.coordinates[i];  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] ^= rhs.coordinates[4];  
    case 4: this->coordinates[3] ^= rhs.coordinates[3];  
    case 3: this->coordinates[2] ^= rhs.coordinates[2];  
    case 2: this->coordinates[1] ^= rhs.coordinates[1];  
    case 1: this->coordinates[0] ^= rhs.coordinates[0];  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] ^= rhs.coordinates[i];  
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator |= (const IA_Point<int>& rhs) {  
    if (this->dim() != rhs.dim()) {  
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return *this;  
    }  
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] |= rhs.coordinates[i];  
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] |= rhs.coordinates[4];  
    case 4: this->coordinates[3] |= rhs.coordinates[3];  
    case 3: this->coordinates[2] |= rhs.coordinates[2];  
    case 2: this->coordinates[1] |= rhs.coordinates[1];  
    case 1: this->coordinates[0] |= rhs.coordinates[0];  
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] |= rhs.coordinates[i];  
	    }
    }
#endif
    return *this;  
}

  
IA_Point<int>& IA_Point<int>::operator %= (int rhs) {
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] %= rhs;
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] %= rhs;
    case 4: this->coordinates[3] %= rhs;
    case 3: this->coordinates[2] %= rhs;
    case 2: this->coordinates[1] %= rhs;
    case 1: this->coordinates[0] %= rhs;
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] %= rhs;
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator <<= (int rhs) {
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] <<= rhs;
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] <<= rhs;
    case 4: this->coordinates[3] <<= rhs;
    case 3: this->coordinates[2] <<= rhs;
    case 2: this->coordinates[1] <<= rhs;
    case 1: this->coordinates[0] <<= rhs;
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] <<= rhs;
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator >>= (int rhs) {
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] >>= rhs;
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] >>= rhs;
    case 4: this->coordinates[3] >>= rhs;
    case 3: this->coordinates[2] >>= rhs;
    case 2: this->coordinates[1] >>= rhs;
    case 1: this->coordinates[0] >>= rhs;
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] >>= rhs;
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator &= (int rhs) {
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] &= rhs;
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] &= rhs;
    case 4: this->coordinates[3] &= rhs;
    case 3: this->coordinates[2] &= rhs;
    case 2: this->coordinates[1] &= rhs;
    case 1: this->coordinates[0] &= rhs;
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] &= rhs;
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator ^= (int rhs) {
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] ^= rhs;
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] ^= rhs;
    case 4: this->coordinates[3] ^= rhs;
    case 3: this->coordinates[2] ^= rhs;
    case 2: this->coordinates[1] ^= rhs;
    case 1: this->coordinates[0] ^= rhs;
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] ^= rhs;
	    }
    }
#endif
    return *this;  
}
  
IA_Point<int>& IA_Point<int>::operator |= (int rhs) {
    if (this->dim() == 0) {  
	ia_throw(Point_Uninitialized_Exception(__FILE__, __LINE__));;  
	return *this;  
    }  
#ifdef NO_LOOP_UNROLLING
    for (unsigned i = 0; i < this->dim(); i++)  
	this->coordinates[i] |= rhs;
#else
    switch (this->dim()) {
    case 5: this->coordinates[4] |= rhs;
    case 4: this->coordinates[3] |= rhs;
    case 3: this->coordinates[2] |= rhs;
    case 2: this->coordinates[1] |= rhs;
    case 1: this->coordinates[0] |= rhs;
    case 0: break;
    default:
	    {
		for (unsigned i = 0; i < this->dim(); i++)  
		    this->coordinates[i] |= rhs;
	    }
    }
#endif
    return *this;  
}

IA_Point<int> IA_Point<int>::operator % (const IA_Point<int>& rhs) const  
{ 
   IA_Point<int> result_point = *this;    
   return result_point %= rhs;  
}

IA_Point<int> IA_Point<int>::operator << (const IA_Point<int>& rhs) const  
{ 
   IA_Point<int> result_point = *this;    
   return result_point <<= rhs;  
}

IA_Point<int> IA_Point<int>::operator >> (const IA_Point<int>& rhs) const  
{ 
   IA_Point<int> result_point = *this;    
   return result_point >>= rhs;  
}

IA_Point<int> IA_Point<int>::operator & (const IA_Point<int>& rhs) const  
{ 
   IA_Point<int> result_point = *this;    
   return result_point &= rhs;  
}

IA_Point<int> IA_Point<int>::operator ^ (const IA_Point<int>& rhs) const  
{ 
   IA_Point<int> result_point = *this;    
   return result_point ^= rhs;  
}

IA_Point<int> IA_Point<int>::operator | (const IA_Point<int>& rhs) const  
{ 
   IA_Point<int> result_point = *this;    
   return result_point |= rhs;  
}

IA_Point<int> operator % (int lhs, const IA_Point<int>& rhs) {
    IA_Point<int> result_point = extend_to_point (lhs, rhs.dim());  
    return result_point %= rhs;  
}

IA_Point<int> operator << (int lhs, const IA_Point<int>& rhs) {
    IA_Point<int> result_point = extend_to_point (lhs, rhs.dim());  
    return result_point <<= rhs;  
}

IA_Point<int> operator >> (int lhs, const IA_Point<int>& rhs) {
    IA_Point<int> result_point = extend_to_point (lhs, rhs.dim());  
    return result_point >>= rhs;  
}

IA_Point<int> operator & (int lhs, const IA_Point<int>& rhs) {
    IA_Point<int> result_point = extend_to_point (lhs, rhs.dim());  
    return result_point &= rhs;  
}

IA_Point<int> operator ^ (int lhs, const IA_Point<int>& rhs) {
    IA_Point<int> result_point = extend_to_point (lhs, rhs.dim());  
    return result_point ^= rhs;  
}

IA_Point<int> operator | (int lhs, const IA_Point<int>& rhs) {
    IA_Point<int> result_point = extend_to_point (lhs, rhs.dim());  
    return result_point |= rhs;  
}

IA_Point<int> IA_Point<int>::operator % (int rhs) const
{ 
   IA_Point<int> result_point = *this;    
   return result_point %= rhs;  
}

IA_Point<int> IA_Point<int>::operator << (int rhs) const
{ 
   IA_Point<int> result_point = *this;    
   return result_point <<= rhs;  
}

IA_Point<int> IA_Point<int>::operator >> (int rhs) const
{ 
   IA_Point<int> result_point = *this;    
   return result_point >>= rhs;  
}

IA_Point<int> IA_Point<int>::operator & (int rhs) const
{ 
   IA_Point<int> result_point = *this;    
   return result_point &= rhs;  
}

IA_Point<int> IA_Point<int>::operator ^ (int rhs) const
{ 
   IA_Point<int> result_point = *this;    
   return result_point ^= rhs;  
}

IA_Point<int> IA_Point<int>::operator | (int rhs) const
{ 
   IA_Point<int> result_point = *this;    
   return result_point |= rhs;  
}

IA_Point<int> cross(const IA_Point<int> &a, const IA_Point<int> &b)
{
    if (a.dim()!=3 || b.dim()!=3) {
	ia_throw(Point_CrossProductRequireDimension3_Exception(__FILE__,__LINE__));
	return IA_Point<int> ();
    }

    return IA_Point<int> (a[1]*b[2]-a[2]*b[1],
			 a[2]*b[0]-a[0]*b[2],
			 a[0]*b[1]-a[1]*b[0]);
}

int dot(const IA_Point<int> &a, const IA_Point<int> &b)
{
    if (a.dim() != b.dim()) {
	ia_throw(Point_DimensionMismatch_Exception(__FILE__, __LINE__));  
	return int ();  
    }

    int	rval=0;
    for (int i=0; i<a.dim(); i++) {
	rval += a[i]*b[i];
    }
    return rval;
}
