#ifndef POINTS_DECL
#define POINTS_DECL

#include <vector>
#include <cmath>
#include <iostream>
#include "Auxiliary.cpp"
using namespace std;



#pragma mark -
#pragma mark XY Point
#pragma mark 


class XYPoint{
public:
	double X,Y;
	static const unsigned int dimension=2;
	
	XYPoint(){ X=Y=0; };
	XYPoint(double _X, double _Y){ X=_X; Y=_Y; }
	explicit XYPoint(double value){ X=Y=value; } 
	explicit XYPoint(int value){ X=Y=value; }
	XYPoint &operator=(const double &value){ X=Y=value; return *this; } //assignment operator
	
	//vector space operations
	XYPoint& operator+=(const XYPoint& p){ 	X += p.X; 		Y += p.Y; 		return *this; }
	XYPoint& operator-=(const XYPoint& p){ 	X -= p.X; 		Y -= p.Y; 		return *this; }
	XYPoint& operator*=(const XYPoint& p){ 	X *= p.X; 		Y *= p.Y; 		return *this; }
	XYPoint& operator/=(const XYPoint& p){ 	X /= p.X; 		Y /= p.Y; 		return *this; }
	XYPoint& operator*=(double scalar){ 	X *= scalar; 	Y *= scalar; 	return *this; }
	XYPoint& operator/=(double scalar){ 	X /= scalar; 	Y /= scalar; 	return *this; }
	
	
	//imitating STL vectors
	double& operator[](unsigned int k){
		switch(k){
		case 0: return X;
		case 1: return Y;
		}
	}
	
	const double& operator[](unsigned int k) const{
		switch(k){
		case 0: return X;
		case 1: return Y;
		}
	}
	
	unsigned int size() const{ return XYPoint::dimension; }
	
	//text-descriptions of state and type
	string currentValueDescription() const{ return "X=" + makeString(X) + ", Y=" + makeString(Y); }
	static string listComponentNames(const string &separator){ return "X" + separator + "Y"; }
	static string componentName(unsigned int c){ return (c==0 ? "X" : "Y"); }
	
};



inline XYPoint operator+(XYPoint p1, const XYPoint &p2){
	p1.X += p2.X; 
	p1.Y += p2.Y; 
	return p1;
}

inline XYPoint operator-(XYPoint p1, const XYPoint &p2){
	p1.X -= p2.X;
	p1.Y -= p2.Y;
	return p1;
}

inline XYPoint operator*(double scalar, XYPoint p){
	p.X *= scalar;
	p.Y *= scalar;
	return p;
}

inline XYPoint operator*(XYPoint p, double scalar){
	p.X *= scalar;
	p.Y *= scalar;
	return p;
}

inline XYPoint operator/(XYPoint p, double scalar){
	p.X /= scalar;
	p.Y /= scalar;
	return p;
}

//component-wise multiplication
inline XYPoint operator*(XYPoint p1, const XYPoint &p2){
	p1.X *= p2.X;
	p1.Y *= p2.Y;
	return p1;
}

//component-wise division
inline XYPoint operator/(XYPoint p1, const XYPoint &p2){
	p1.X /= p2.X;
	p1.Y /= p2.Y;
	return p1;
}

inline XYPoint min(XYPoint p1, const XYPoint &p2){
	p1.X = std::min(p1.X, p2.X);
	p1.Y = std::min(p1.Y, p2.Y);
	return p1;
}

inline XYPoint max(XYPoint p1, const XYPoint &p2){
	p1.X = std::max(p1.X, p2.X);
	p1.Y = std::max(p1.Y, p2.Y);
	return p1;
}

inline XYPoint max(XYPoint p1, double b){
	p1.X = std::max(p1.X, b);
	p1.Y = std::max(p1.Y, b);
	return p1;
}

inline XYPoint sqrt(XYPoint p){
	p.X = std::sqrt(p.X);
	p.Y = std::sqrt(p.Y);
	return p;
}

inline XYPoint log(XYPoint p){
	p.X = std::log(p.X);
	p.Y = std::log(p.Y);
	return p;
}

inline XYPoint exp(XYPoint p){
	p.X = std::exp(p.X);
	p.Y = std::exp(p.Y);
	return p;
}

inline XYPoint pow(XYPoint p, double power){
	p.X = std::pow(p.X, power);
	p.Y = std::pow(p.Y, power);
	return p;
}

inline ostream& operator << (ostream &sout, const XYPoint &point){
	return sout << point.X << " " << point.Y;
}

inline istream& operator >> (istream &sIn, XYPoint &p){
	return sIn >> p.X >> p.Y;
}

			

#pragma mark -
#pragma mark Generic N-dimensional points
#pragma mark



template<unsigned int DIM>
class NPoint{
public:
	double X[DIM];
	static const unsigned int dimension=DIM;

	NPoint(){ for(unsigned int i=0; i<DIM; ++i){ X[i] = 0; } }
	explicit NPoint(const double _X[]){ for(unsigned int i=0; i<DIM; ++i){ X[i] = _X[i]; } }
	explicit NPoint(double value){ for(unsigned int i=0; i<DIM; ++i){ X[i] = value; } } 
	explicit NPoint(int value){ for(unsigned int i=0; i<DIM; ++i){ X[i] = value; } }
	NPoint &operator=(const double &value){ for(unsigned int i=0; i<DIM; ++i){ X[i] = value; } return *this; } //assignment operator
	
	//vector space operations
	NPoint<DIM>& operator+=(const NPoint<DIM>& p){ for(unsigned int i=0; i<DIM; ++i){ X[i] += p.X[i]; } return *this; }
	NPoint<DIM>& operator-=(const NPoint<DIM>& p){ for(unsigned int i=0; i<DIM; ++i){ X[i] -= p.X[i]; } return *this; }
	NPoint<DIM>& operator*=(const NPoint<DIM>& p){ for(unsigned int i=0; i<DIM; ++i){ X[i] *= p.X[i]; } return *this; }
	NPoint<DIM>& operator/=(const NPoint<DIM>& p){ for(unsigned int i=0; i<DIM; ++i){ X[i] /= p.X[i]; } return *this; }
	NPoint<DIM>& operator*=(double scalar){ for(unsigned int i=0; i<DIM; ++i){ X[i] *= scalar; } return *this; }
	NPoint<DIM>& operator/=(double scalar){ for(unsigned int i=0; i<DIM; ++i){ X[i] /= scalar; } return *this; }
	
	//imitating STL vectors
	double& operator[](unsigned int k){
		return X[k];
	}
	
	const double& operator[](unsigned int k) const{
		return X[k];
	}
	
	unsigned int size() const{ return NPoint::dimension; }
	
	//text-descriptions of state and type
	string currentValueDescription() const{ ostringstream s; for(unsigned int i=0; i<DIM; ++i){ s << (i==0 ? "" : ", ") << "X" << (i+1) << " = " << X[i]; } return s.str(); }
	static string listComponentNames(const string &separator){ string s; for(unsigned int i=0; i<DIM; ++i){ s+= (i>0 ? separator : "") + "X" + makeString(i+1); } return s; }
	static string componentName(unsigned int c){ return "X" + makeString(c+1); }
	
};


template<unsigned int DIM>
inline NPoint<DIM> operator+(NPoint<DIM> p1, const NPoint<DIM> &p2){
	for(unsigned int i=0; i<DIM; ++i){
		p1.X[i] += p2.X[i];
	}
	return p1;
}

template<unsigned int DIM>
inline NPoint<DIM> operator-(NPoint<DIM> p1, const NPoint<DIM> &p2){
	for(unsigned int i=0; i<DIM; ++i){
		p1.X[i] -= p2.X[i];
	}
	return p1;
}

template<unsigned int DIM>
inline NPoint<DIM> operator*(double scalar, NPoint<DIM> p){
	for(unsigned int i=0; i<DIM; ++i){
		p.X[i] *= scalar;
	}
	return p;
}

template<unsigned int DIM>
inline NPoint<DIM> operator*(NPoint<DIM> p, double scalar){
	for(unsigned int i=0; i<DIM; ++i){
		p.X[i] *= scalar;
	}
	return p;
}

template<unsigned int DIM>
inline NPoint<DIM> operator/(NPoint<DIM> p, double scalar){
	for(unsigned int i=0; i<DIM; ++i){
		p.X[i] /= scalar;
	}
	return p;
}

//component-wise multiplication
template<unsigned int DIM>
inline NPoint<DIM> operator*(NPoint<DIM> p1, const NPoint<DIM> &p2){
	for(unsigned int i=0; i<DIM; ++i){
		p1.X[i] *= p2.X[i];
	}
	return p1;
}

//component-wise division
template<unsigned int DIM>
inline NPoint<DIM> operator/(NPoint<DIM> p1, const NPoint<DIM> &p2){
	for(unsigned int i=0; i<DIM; ++i){
		p1.X[i] /= p2.X[i];
	}
	return p1;
}

template<unsigned int DIM>
inline NPoint<DIM> min(NPoint<DIM> p1, const NPoint<DIM> &p2){
	for(unsigned int i=0; i<DIM; ++i){
		p1.X[i] = std::min(p1.X[i], p2.X[i]);
	}
	return p1;
}

template<unsigned int DIM>
inline NPoint<DIM> max(NPoint<DIM> p1, const NPoint<DIM> &p2){
	for(unsigned int i=0; i<DIM; ++i){
		p1.X[i] = std::max(p1.X[i], p2.X[i]);
	}
	return p1;
}

template<unsigned int DIM>
inline NPoint<DIM> max(NPoint<DIM> p1, double b){
	for(unsigned int i=0; i<DIM; ++i){
		p1.X[i] = std::max(p1.X[i], b);
	}
	return p1;
}

template<unsigned int DIM>
inline NPoint<DIM> sqrt(NPoint<DIM> p){
	for(unsigned int i=0; i<DIM; ++i){
		p.X[i] = std::sqrt(p.X[i]);
	}
	return p;
}

template<unsigned int DIM>
inline NPoint<DIM> log(NPoint<DIM> p){
	for(unsigned int i=0; i<DIM; ++i){
		p.X[i] = std::log(p.X[i]);
	}
	return p;
}

template<unsigned int DIM>
inline NPoint<DIM> exp(NPoint<DIM> p){
	for(unsigned int i=0; i<DIM; ++i){
		p.X[i] = std::exp(p.X[i]);
	}
	return p;
}

template<unsigned int DIM>
inline NPoint<DIM> pow(NPoint<DIM> p, double power){
	for(unsigned int i=0; i<DIM; ++i){
		p.X[i] = std::pow(p.X[i],power);
	}
	return p;
}

template<unsigned int DIM>
inline ostream& operator << (ostream &sout, const NPoint<DIM> &point){
	for(unsigned int i=0; i<DIM; ++i){
		sout << (i==0 ? "" : " ") << point.X[i];
	}
	return sout;
}

template<unsigned int DIM>
istream& operator >> (istream &sIn, NPoint<DIM> &p){
	for(unsigned int i=0; i<DIM; ++i){
		sIn >> p.X[i];
	}
	return sIn;
}






#pragma mark -
#pragma mark Generic 1-dimensional point


class XPoint{
public:
	double X;
	static const unsigned int dimension = 1;
	
	XPoint(){ X = 0; }
	explicit XPoint(double _X[]){ X = _X[0]; }
	explicit XPoint(double value){ X = value; } 
	explicit XPoint(int value){ X = value; }
	XPoint &operator=(const double &value){ X = value; return *this; } //assignment operator
	
	//vector space operations
	XPoint& operator+=(const XPoint& p){ X += p.X; return *this; }
	XPoint& operator-=(const XPoint& p){ X -= p.X; return *this; }
	XPoint& operator*=(const XPoint& p){ X *= p.X; return *this; }
	XPoint& operator/=(const XPoint& p){ X /= p.X; return *this; }
	XPoint& operator*=(double scalar){ X *= scalar; return *this; }
	XPoint& operator/=(double scalar){ X /= scalar;  return *this; }
	
	//imitating STL vectors
	double& operator[](unsigned int k){
		return X;
	}
	
	const double& operator[](unsigned int k) const{
		return X;
	}
	
	unsigned int size() const{ return XPoint::dimension; }
	
	//text-descriptions of state and type
	string currentValueDescription() const{ return "X = " + makeString(X); }
	static string listComponentNames(const string &separator){ return "X"; }
	static string componentName(unsigned int c){ return (c==0 ? "X" : ""); }
	
};



inline XPoint operator+(XPoint p1, const XPoint &p2){
	p1.X += p2.X;
	return p1;
}


inline XPoint operator-(XPoint p1, const XPoint &p2){
	p1.X -= p2.X;
	return p1;
}


inline XPoint operator*(double scalar, XPoint p){
	p.X *= scalar;
	return p;
}


inline XPoint operator*(XPoint p, double scalar){
	p.X *= scalar;
	return p;
}


inline XPoint operator/(XPoint p, double scalar){
	p.X /= scalar;
	return p;
}

//component-wise multiplication

inline XPoint operator*(XPoint p1, const XPoint &p2){
	p1.X *= p2.X;
	return p1;
}

//component-wise division

inline XPoint operator/(XPoint p1, const XPoint &p2){
	p1.X /= p2.X;
	return p1;
}


inline XPoint min(XPoint p1, const XPoint &p2){
	p1.X = std::min(p1.X, p2.X);
	return p1;
}


inline XPoint max(XPoint p1, const XPoint &p2){
	p1.X = std::max(p1.X, p2.X);
	return p1;
}


inline XPoint max(XPoint p1, double b){
	p1.X = std::max(p1.X, b);
	return p1;
}


inline XPoint sqrt(XPoint p){
	p.X = std::sqrt(p.X);
	return p;
}


inline XPoint log(XPoint p){
	p.X = std::log(p.X);
	return p;
}

inline XPoint exp(XPoint p){
	p.X = std::exp(p.X);
	return p;
}


inline XPoint pow(XPoint p, double power){
	p.X = std::pow(p.X,power);
	return p;
}


inline ostream& operator << (ostream &sout, const XPoint &point){
	sout << point.X;
	return sout;
}

inline istream& operator >> (istream &sIn, XPoint &p){
	return sIn >> p.X;
}




#pragma mark -
#pragma mark Tensor products between points and vectors
#pragma mark 

// tensor product between N-point and a vector of size DIM2
// resulting vector respresents a DIM1*DIM2 dimensional matrix stored in row-major order
// That is, matrix[r*DIM2+c] = p1[r]*p2[c]
template<unsigned int DIM1>
inline vector<double> TensorProduct(const NPoint<DIM1> &p1, const vector<double> &p2){
	unsigned int r, c, DIM2=p2.size();
	vector<double> p(DIM1*DIM2);
	for(r=0; r<DIM1; ++r){
		for(c=0; c<DIM2; ++c){
			p[r*DIM2 + c] = p1.X[r] * p2[c];
		}
	}
	return p;
}


// tensor product between XYpoint and a vector of size DIM2
// resulting vector respresents a 2*DIM2 dimensional matrix stored in row-major order
// That is, matrix[r*DIM2+c] = p1[r]*p2[c]
inline vector<double> TensorProduct(const XYPoint &p1, const vector<double> &p2){
	unsigned int c, DIM2=p2.size();
	vector<double> p(2*DIM2);
	for(c=0; c<DIM2; ++c){
		p[0*DIM2 + c] = p1.X * p2[c];
		p[1*DIM2 + c] = p1.Y * p2[c];
	}
	return p;
}



inline vector<double> TensorProduct(const XYPoint &p1, const XYPoint &p2){
	vector<double> p(2*2);
	p[0] = p1.X * p2.X;
	p[1] = p1.X * p2.Y;
	p[2] = p1.Y * p2.X;
	p[3] = p1.Y * p2.Y;
	return p;
}


// tensor product between Xpoint and a vector of size DIM2
// resulting vector respresents a 1*DIM2 dimensional matrix stored in row-major order
// That is, matrix[r*DIM2+c] = p1[r]*p2[c]
inline vector<double> TensorProduct(const XPoint &p1, const vector<double> &p2){
	unsigned int c, DIM2=p2.size();
	vector<double> p(1*DIM2);
	for(c=0; c<DIM2; ++c){
		p[c] = p1.X * p2[c];
	}
	return p;
}

inline vector<double> TensorProduct(const XPoint &p1, const XPoint &p2){
	vector<double> p(1*1);
	p[0] = p1.X * p2.X;
	return p;
}




			

#endif