/* class Vec2, Vec3, Vec4, and Vec<N>
 *
 * These classes represent a DISPLACEMENT in R^2, R^3, R^4, and R^N
 * space, respectively. For example, the origin SHOULD NOT be described
 * as a Vec3(0,0,0). It is a point, and as such, should be represented
 * by a point object. There is NO way to cast a Vec into a Pt. Needing
 * to do so means you're thinking about it wrong. The way to do it would
 * be (PtN::origin + vecN) if you really have to.
 *
 * Things to note:
 *  -static member e[] is the set of canonical basis vectors. i.e.
 *   pt2.e[0]*Vec2::e[0] + pt2.e[1]*Vec2::e[1] =~= pt2 - Pt2::origin
 *  -The components of a vector can be extracted by the [] operator.
 *  -There are no overloaded operators for vector by vector
 *   multiplication. Vec3Dot(A,B) returns the dot product, while
 *   Vec3Cross(C, A, B) sets C = A x B. In 2D, the cross product is just
 *   the signed magnitude of C.
 *
 * Affine geometry notes:
 *   Point +/- Vector = Point
 *   Vector + Vector  = Vector
 *   Point - Point    = Vector
 *   Point + Point    = Point (so that we can do affine combinations of points)
 *
 * Originally based off Matthew Fisher's Vec3 class.
 *
 * Last modified: 2008-02-10
 */

#ifndef _VECTOR_H_
#define _VECTOR_H_

#include <vDEC/include/vdec_types.h>
#include <cstdlib>
#include <cmath>

struct Vec2{
	/// The vector component data.
	real x, y;

	/// Returns the zero vector.
	inline static const Vec2 Zero();
	/// Returns the which-th canonical basis vector (the vector of all zeroes and a single 1 in the which-th slot).
	inline static const Vec2 Basis(unsigned int which);
	
	// Default constructor initializes all components to zero.
	inline Vec2();
	inline Vec2(const Vec2 &V);
	inline Vec2(const real &_x, const real &_y);

	inline Vec2& operator=(const Vec2 &V);
	/// Sets the current vector to be a point randomly chosen on the unit circle.
	inline static Vec2 RandomNormal();

	inline Vec2& operator *= (const real &C);
	inline Vec2& operator /= (const real &C);
	inline Vec2& operator += (const Vec2 &Right);
	inline Vec2& operator -= (const Vec2 &Right);

	/// Indexes into the components.
	inline real& operator [] (unsigned int axis);
	/// Indexes into the components.
	inline real operator [] (unsigned int axis) const;

	/// Normalizes the vector. If zero vector, produces an infinite vector.
	inline void Normalize();

	/// Returns the L^2 norm of the vector.
	inline real Length() const;
	/// Returns the square of the L^2 norm of the vector.
	inline real LengthSq() const;

	/// Performs the dot (inner) product on two vectors (Euclidean R^2 metric).
	inline static real Dot(const Vec2 &Left, const Vec2 &Right);
	/// Performs the scalar product (length of cross product).
	inline static real Cross(const Vec2 &Left, const Vec2 &Right);
};

//Vec2 overloaded operators
inline Vec2 operator * (const Vec2 &Left, real Right);
inline Vec2 operator * (real Left, const Vec2 &Right);
inline Vec2 operator / (const Vec2 &Left, real Right);
inline Vec2 operator + (const Vec2 &Left, const Vec2 &Right);
inline Vec2 operator - (const Vec2 &Left, const Vec2 &Right);
inline Vec2 operator - (const Vec2 &V);

struct Vec3 {
	/// The vector component data.
	real x, y, z;

	/// Returns the zero vector.
	inline static const Vec3 Zero();
	/// Returns the which-th canonical basis vector (the vector of all zeroes and a single 1 in the which-th slot).
	inline static const Vec3 Basis(unsigned int which);
	
	// Default constructor initializes all components to zero.
	inline Vec3();
	inline Vec3(const Vec3 &V);
	inline Vec3(const real &_x, const real &_y, const real &_z);

	inline Vec3& operator=(const Vec3 &V);
	/// Sets the current vector to be a point randomly chosen on the unit circle.
	inline static Vec3 RandomNormal();

	inline Vec3& operator *= (const real &C);
	inline Vec3& operator /= (const real &C);
	inline Vec3& operator += (const Vec3 &Right);
	inline Vec3& operator -= (const Vec3 &Right);

	/// Indexes into the components.
	inline real& operator [] (unsigned int axis);
	/// Indexes into the components.
	inline real operator [] (unsigned int axis) const;

	/// Normalizes the vector. If zero vector, produces an infinite vector.
	inline void Normalize();

	/// Returns the L^2 norm of the vector.
	inline real Length() const;
	/// Returns the square of the L^2 norm of the vector.
	inline real LengthSq() const;

	/// Returns the cross product of two vectors Left x Right.
	inline static Vec3 Cross(const Vec3 &Left, const Vec3 &Right);
	/// Returns the dot (inner) product of the two vectors (Euclidean R^3 metric).
	inline static real Dot(const Vec3 &Left, const Vec3 &Right);
	/// Returns the angle between two vectors. Always in the range [0, Pi].
	inline static real AngleBetween(const Vec3 &Left, const Vec3 &Right);
	
	/// Returns an orthonormal basis, given v0, such that: v0 x v1 = v2.
	inline static void MakeBasis(const Vec3 &v0, Vec3 &v1, Vec3 &v2);
};

//Vec3 overloaded operators
inline Vec3 operator * (const Vec3 &Left, real Right);
inline Vec3 operator * (real Left, const Vec3 &Right);
inline Vec3 operator / (const Vec3 &Left, real Right);
inline Vec3 operator + (const Vec3 &Left, const Vec3 &Right);
inline Vec3 operator - (const Vec3 &Left, const Vec3 &Right);
inline Vec3 operator - (const Vec3 &V);


struct Vec4 {
	/// The vector component data.
	real x, y, z, w;
	
	// Default constructor initializes all components to zero.
	inline Vec4();
	inline Vec4(const Vec4 &V);
	inline Vec4(const real &_x, const real &_y, const real &_z, const real &_w);

	inline Vec4& operator = (const Vec4 &V);

	inline Vec4& operator *= (const real &C);
	inline Vec4& operator /= (const real &C);
	inline Vec4& operator += (const Vec4 &Right);
	inline Vec4& operator -= (const Vec4 &Right);

	/// Indexes into the components.
	inline real& operator [] (unsigned int axis);
	/// Indexes into the components.
	inline real operator [] (unsigned int axis) const;

	/// Normalizes the vector. If zero vector, produces an infinite vector.
	inline void Normalize();

	/// Returns the L^2 norm of the vector.
	inline real Length() const;
	/// Returns the square of the L^2 norm of the vector.
	inline real LengthSq() const;

	/// Returns the dot (inner) product of the two vectors (Euclidean R^4 metric).
	inline static real Dot(const Vec4 &Left, const Vec4 &Right);
};


//Vec4 overloaded operators
inline Vec4 operator * (const Vec4 &Left, real Right);
inline Vec4 operator * (real Left, const Vec4 &Right);
inline Vec4 operator / (const Vec4 &Left, real Right);
inline Vec4 operator + (const Vec4 &Left, const Vec4 &Right);
inline Vec4 operator - (const Vec4 &Left, const Vec4 &Right);
inline Vec4 operator - (const Vec4 &V);

/// Generic multidimensional vector class
/// This class has no constructors in order to make things like this work:
///   Vec<3> u = {1,0,0};
template <unsigned int N>
struct Vec{
	real e[N];
	
	/*
	Vec(){ for(unsigned int i = 0; i < N; ++i){ e[i] = 0; } }
	Vec(const real V[N]){ for(unsigned int i = 0; i < N; ++i){ e[i] = V[i]; } }
	Vec(const Vec<N> &V){ for(unsigned int i = 0; i < N; ++i){ e[i] = V.e[i]; } }
	*/

	Vec<N>& operator = (const Vec<N> &V){ for(unsigned int i = 0; i < N; ++i){ e[i] = V.e[i]; } return *this; }

	Vec<N>& operator *= (const real &C){ for(unsigned int i = 0; i < N; ++i){ e[i] *= C; } return *this; }
	Vec<N>& operator /= (const real &C){ real f = (real)1.0/C; for(unsigned int i = 0; i < N; ++i){ e[i] *= f; } return *this; }
	Vec<N>& operator += (const Vec<N> &Right){ for(unsigned int i = 0; i < N; ++i){ e[i] += Right.e[i]; } return *this; }
	Vec<N>& operator -= (const Vec<N> &Right){ for(unsigned int i = 0; i < N; ++i){ e[i] -= Right.e[i]; } return *this; }

	/// Indexes into the components.
	real& operator [] (unsigned int axis){ return e[axis]; }
	/// Indexes into the components.
	real operator [] (unsigned int axis) const{ return e[axis]; }

	/// Normalizes the vector. If zero vector, produces an infinite vector.
	void Normalize(){ real l = Length(); for(unsigned int i = 0; i < N; ++i){ e[i] /= l; } }

	/// Returns the L^2 norm of the vector.
	real Length() const{ return (real)VDEC_SQRT(LengthSq()); }
	/// Returns the square L^2 norm of the vector.
	real LengthSq() const{ real s = 0; for(unsigned int i = 0; i < N; ++i){ s += e[i]*e[i]; } return s; }

	/// Returns the dot (inner) product of the two vectors (Euclidean R^N metric).
	static real Dot(const Vec<N> &Left, const Vec<N> &Right){
		real s = 0;
		for(unsigned int i = 0; i < N; ++i){ s += Left.e[i] * Right.e[i]; }
		return s;
	}
};

template <unsigned int N>
Vec<N> operator * (const Vec<N> &V, const real &C){
	Vec<N> Result(V);
	return (Result *= C);
}

template <unsigned int N>
Vec<N> operator * (const real &C, const Vec<N> &V){
	Vec<N> Result(V);
	return (Result *= C);
}

template <unsigned int N>
Vec<N> operator / (const Vec<N> &V, const real &C){
	Vec<N> Result(V);
	return (Result /= C);
}

template <unsigned int N>
Vec<N> operator + (const Vec<N> &Left, const Vec<N> &Right){
	Vec<N> Result(Left);
	return (Result += Right);
}

template <unsigned int N>
Vec<N> operator - (const Vec<N> &Left, const Vec<N> &Right){
	Vec<N> Result(Left);
	return (Result -= Right);
}

template <unsigned int N>
Vec<N> operator - (const Vec<N> &V){
	Vec<N> Result(V);
	return (Result *= -1);
}

/// Same as Vec2::Cross
inline real operator^(const Vec<2> &Left, const Vec<2> &Right);
/// Same as Vec3::Cross
inline Vec<3> operator^(const Vec<3> &Left, const Vec<3> &Right);

/********************** Implementation of Vec2 ************************/


const Vec2 Vec2::Zero(){
	static const Vec2 z(0,0);
	return z;
}

const Vec2 Vec2::Basis(unsigned int which){
	static const Vec2 bases[2] = {Vec2(1,0), Vec2(0,1)};
	return bases[which];
}

Vec2::Vec2():x(0),y(0){}
Vec2::Vec2(const Vec2 &V):x(V.x),y(V.y){}
Vec2::Vec2(const real &_x, const real &_y):x(_x),y(_y){}

Vec2& Vec2::operator=(const Vec2 &V){
	x = V.x;
	y = V.y;
	return *this;
}

Vec2 Vec2::RandomNormal(){
	real theta = (real)(real(rand()) * 2.0 * M_PI / (real)RAND_MAX);
	return Vec2(VDEC_COS(theta), VDEC_SIN(theta));
}

Vec2& Vec2::operator *= (const real &C){ x *= C; y *= C; return *this; }
Vec2& Vec2::operator /= (const real &C){ x /= C; y /= C; return *this; }
Vec2& Vec2::operator += (const Vec2 &Right){ x += Right.x; y += Right.y; return *this; }
Vec2& Vec2::operator -= (const Vec2 &Right){ x -= Right.x; y -= Right.y; return *this; }

real& Vec2::operator [] (unsigned int axis)      { return ((real*)this)[axis]; }
real  Vec2::operator [] (unsigned int axis) const{ return ((real*)this)[axis]; }

void Vec2::Normalize(){
	real s = (real)1.0/Length();
	x *= s; y *= s;
}

real Vec2::Length() const{ return (real)VDEC_SQRT(x*x+y*y); }
real Vec2::LengthSq() const{ return (real)(x*x+y*y); }

real Vec2::Dot(const Vec2 &Left, const Vec2 &Right){ return (Left.x * Right.x + Left.y * Right.y); }
real Vec2::Cross(const Vec2 &Left, const Vec2 &Right){ return Left.x * Right.y - Left.y * Right.x; }

Vec2 operator * (const Vec2 &Left, real Right){
	Vec2 Return;
	Return.x = Left.x * Right;
	Return.y = Left.y * Right;
	return Return;
}

Vec2 operator * (real Right, const Vec2 &Left){
	Vec2 Return;
	Return.x = Left.x * Right;
	Return.y = Left.y * Right;
	return Return;
}

Vec2 operator / (const Vec2 &Left, real Right){
	Vec2 Return;
	Return.x = Left.x / Right;
	Return.y = Left.y / Right;
	return Return;
}

Vec2 operator + (const Vec2 &Left, const Vec2 &Right){
	Vec2 Return;
	Return.x = Left.x + Right.x;
	Return.y = Left.y + Right.y;
	return Return;
}

Vec2 operator - (const Vec2 &Left, const Vec2 &Right){
	Vec2 Return;
	Return.x = Left.x - Right.x;
	Return.y = Left.y - Right.y;
	return Return;
}

Vec2 operator - (const Vec2 &V){
	Vec2 Result;
	Result.x = -V.x;
	Result.y = -V.y;
	return Result;
}


/********************** Implementation of Vec3 ************************/

const Vec3 Vec3::Zero(){
	static const Vec3 z(0,0,0);
	return z;
}

const Vec3 Vec3::Basis(unsigned int which){
	static const Vec3 bases[3] = {Vec3(1,0,0), Vec3(0,1,0), Vec3(0,0,1)};
	return bases[which];
}

Vec3::Vec3():x(0),y(0),z(0){}
Vec3::Vec3(const Vec3 &V):x(V.x),y(V.y),z(V.z){}
Vec3::Vec3(const real &_x, const real &_y, const real &_z):x(_x),y(_y),z(_z){}

Vec3& Vec3::operator=(const Vec3 &V){
	x = V.x;
	y = V.y;
	z = V.z;
	return *this;
}

Vec3 Vec3::RandomNormal(){
	Vec3 ret(0,0, real(rand()) / RAND_MAX);
	real theta = (real)(real(rand()) * 2.0 * M_PI / (real)RAND_MAX);
	ret.x = VDEC_SQRT(1-ret.z*ret.z);
	ret.y = ret.x * VDEC_SIN(theta);
	ret.x *= VDEC_COS(theta);
	return ret;
}

//Overloaded operators
Vec3& Vec3::operator *= (const real &C){ x *= C; y *= C; z *= C; return *this; }
Vec3& Vec3::operator /= (const real &C){ x /= C; y /= C; z /= C; return *this; }
Vec3& Vec3::operator += (const Vec3 &Right){ x += Right.x; y += Right.y; z += Right.z; return *this; }
Vec3& Vec3::operator -= (const Vec3 &Right){ x -= Right.x; y -= Right.y; z -= Right.z; return *this; }

real& Vec3::operator [] (unsigned int axis){ return ((real*)this)[axis]; }
real  Vec3::operator [] (unsigned int axis) const{ return ((real*)this)[axis]; }

void Vec3::Normalize(){
	real s = (real)1.0/Length();
	x *= s;
	y *= s;
	z *= s;
}

real Vec3::Length() const{ return (real)VDEC_SQRT(x*x+y*y+z*z); }
real Vec3::LengthSq() const{ return (real)(x*x+y*y+z*z); }

Vec3 Vec3::Cross(const Vec3 &Left, const Vec3 &Right){
	return Vec3(
		Left.y * Right.z - Left.z * Right.y,
		Left.z * Right.x - Left.x * Right.z,
		Left.x * Right.y - Left.y * Right.x);
}

real Vec3::Dot(const Vec3 &Left, const Vec3 &Right){ return (Left.x * Right.x + Left.y * Right.y + Left.z * Right.z); }

real Vec3::AngleBetween(const Vec3 &Left, const Vec3 &Right){
	Vec3 V1(Left), V2(Right);
	V1.Normalize();
	V2.Normalize();
	real Dot = Vec3::Dot(V1, V2);
	if(VDEC_FABS(Dot) >= 1.0f){ Dot = 1; }
	return (real)VDEC_ACOS(Dot);
}

void Vec3::MakeBasis(const Vec3 &v0, Vec3 &v1, Vec3 &v2){
	Vec3 Given = v0;
	Given.Normalize();
	if(VDEC_FABS(Given.x) > VDEC_FABS(Given.y)){
		real invLen = 1 / VDEC_SQRT(Given.x*Given.x + Given.z*Given.z);
		v1.x = -Given.z * invLen;
		v1.y = 0;
		v1.z = Given.x * invLen;
	}else{
		real invLen = 1 / VDEC_SQRT(Given.y*Given.y + Given.z*Given.z);
		v1.x = 0;
		v1.y = Given.z * invLen;
		v1.z = -Given.y * invLen;
	}
	v2 = Vec3::Cross(Given, v1);
}

Vec3 operator * (const Vec3 &Left, real Right){
	Vec3 Return;
	Return.x = Left.x * Right;
	Return.y = Left.y * Right;
	Return.z = Left.z * Right;
	return Return;
}

Vec3 operator * (real Right, const Vec3 &Left){
	Vec3 Return;
	Return.x = Left.x * Right;
	Return.y = Left.y * Right;
	Return.z = Left.z * Right;
	return Return;
}

Vec3 operator / (const Vec3 &Left, real Right){
	Vec3 Return;
	Return.x = Left.x / Right;
	Return.y = Left.y / Right;
	Return.z = Left.z / Right;
	return Return;
}

Vec3 operator + (const Vec3 &Left, const Vec3 &Right){
	Vec3 Return;
	Return.x = Left.x + Right.x;
	Return.y = Left.y + Right.y;
	Return.z = Left.z + Right.z;
	return Return;
}

Vec3 operator - (const Vec3 &Left, const Vec3 &Right){
	Vec3 Return;
	Return.x = Left.x - Right.x;
	Return.y = Left.y - Right.y;
	Return.z = Left.z - Right.z;
	return Return;
}

Vec3 operator - (const Vec3 &V){
	Vec3 Result;
	Result.x = -V.x;
	Result.y = -V.y;
	Result.z = -V.z;
	return Result;
}

/********************** Implementation of Vec4 ************************/

Vec4::Vec4():x(0),y(0),z(0),w(0){}
Vec4::Vec4(const Vec4 &V):x(V.x),y(V.y),z(V.z),w(V.w){}
Vec4::Vec4(const real &_x, const real &_y, const real &_z, const real &_w):x(_x),y(_y),z(_z),w(_w){}

Vec4& Vec4::operator = (const Vec4 &V){ x = V.x; y = V.y; z = V.z; w = V.w; return *this; }

Vec4& Vec4::operator *= (const real &C){ x *= C; y *= C; z *= C; w *= C; return *this; }
Vec4& Vec4::operator /= (const real &C){ x /= C; y /= C; z /= C; w /= C; return *this; }
Vec4& Vec4::operator += (const Vec4 &Right){ x += Right.x; y += Right.y; z += Right.z; w += Right.w; return *this; }
Vec4& Vec4::operator -= (const Vec4 &Right){ x -= Right.x; y -= Right.y; z -= Right.z; w -= Right.w; return *this; }

real& Vec4::operator [] (unsigned int axis){ return ((real*)this)[axis]; }
real  Vec4::operator [] (unsigned int axis) const{ return ((real*)this)[axis]; }

void Vec4::Normalize(){ real s = (real)1.0/Length(); x *= s; y *= s; z *= s; w *= s; }

real Vec4::Length() const{ return (real)VDEC_SQRT(x*x+y*y+z*z+w*w); }
real Vec4::LengthSq() const{ return (real)(x*x+y*y+z*z+w*w); }

real Vec4::Dot(const Vec4 &Left, const Vec4 &Right){ return (Left.x * Right.x + Left.y * Right.y + Left.z * Right.z + Left.w * Right.w); }

Vec4 operator * (const Vec4 &Left, real Right){
	Vec4 Return;
	Return.x = Left.x * Right;
	Return.y = Left.y * Right;
	Return.z = Left.z * Right;
	Return.w = Left.w * Right;
	return Return;
}

Vec4 operator * (real Right, const Vec4 &Left){
	Vec4 Return;
	Return.x = Left.x * Right;
	Return.y = Left.y * Right;
	Return.z = Left.z * Right;
	Return.w = Left.w * Right;
	return Return;
}

Vec4 operator / (const Vec4 &Left, real Right){
	Vec4 Return;
	Return.x = Left.x / Right;
	Return.y = Left.y / Right;
	Return.z = Left.z / Right;
	Return.w = Left.w / Right;
	return Return;
}

Vec4 operator + (const Vec4 &Left, const Vec4 &Right){
	Vec4 Return;
	Return.x = Left.x + Right.x;
	Return.y = Left.y + Right.y;
	Return.z = Left.z + Right.z;
	Return.w = Left.w + Right.w;
	return Return;
}

Vec4 operator - (const Vec4 &Left, const Vec4 &Right){
	Vec4 Return;
	Return.x = Left.x - Right.x;
	Return.y = Left.y - Right.y;
	Return.z = Left.z - Right.z;
	Return.w = Left.w - Right.w;
	return Return;
}

Vec4 operator - (const Vec4 &V){
	Vec4 Result;
	Result.x = -V.x;
	Result.y = -V.y;
	Result.z = -V.z;
	Result.w = -V.w;
	return Result;
}

/********************* Implementation of Vec<N> ***********************/

real operator^(const Vec<2> &Left, const Vec<2> &Right){
	return (Left[0] * Right[1] - Left[1] * Right[0]);
}

Vec<3> operator^(const Vec<3> &Left, const Vec<3> &Right){
	Vec<3> ret;
	ret[0] = Left[1] * Right[2] - Left[2] * Right[1];
	ret[1] = Left[2] * Right[0] - Left[0] * Right[2];
	ret[2] = Left[0] * Right[1] - Left[1] * Right[0];
	return ret;
}
#endif