#ifndef VECTORS_H
#define VECTORS_H

// K-3D
// Copyright (c) 1995-2004, Timothy M. Shead
//
// Contact: tshead@k-3d.com
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation; either
// version 2 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

/** \file
		\brief Vectors and algebra routines
		\author Tim Shead (tshead@k-3d.com)
*/

/****************************************************************
*
* C++ Vector and Matrix Algebra routines
* Author: Jean-Francois DOUE
* Version 3.1 --- October 1993
*
****************************************************************/

//
//	From "Graphics Gems IV / Edited by Paul S. Heckbert
//	Academic Press, 1994, ISBN 0-12-336156-9
//	"You are free to use and modify this code in any way
//	you like." (p. xv)
//
//	Modified by J. Nagle, March 1997
//	-	All functions are inline.
//	-	All functions are const-correct.
//	-	All checking is via the standard "assert" macro.
//

// Modified by Tim Shead for use with K-3D, January 1998

#include <cassert>
#include <cmath>
#include <iostream>
#include <vector>

namespace k3d
{

// Forward declarations ...
class vector2;
class vector3;
class vector4;

/// Axes array indices
enum {VX, VY, VZ, VW};
/// Plane array indices
enum {PA, PB, PC, PD};

/// Returns the maximum of two values
template<class Type> Type sdpMax(const Type& A, const Type& B)
{ return A >= B ? A : B; }

/// Returns the minimum of two values
template<class Type> Type sdpMin(const Type& A, const Type& B)
{ return A <= B ? A : B; }

/////////////////////////////////////////////////////////////////////////////
// vector2

/// A two-dimensional vector
class vector2
{
public:
	/// Stores the vector values
	double n[2];

	// Constructors
	vector2();
	vector2(const double x, const double y);
	vector2(const double d);
	vector2(const double[2]);
	/// Copy constructor
	vector2(const vector2& v);
	/// Casts an vector3 to an vector2, dropping the third dimension
	vector2(const vector3& v);
	/// Casts an vector3 to an vector2, dropping the given dimension
	vector2(const vector3& v, int dropAxis);
	/// Assigns an vector2
	vector2& operator	= ( const vector2& v );
	/// Assigns a C/C++ style array
	vector2& operator  = ( const double[2]);
	/// Adds an vector2
	vector2& operator += ( const vector2& v );
	/// Subtracts an vector2
	vector2& operator -= ( const vector2& v );
	/// Multiplication by a constant
	vector2& operator *= ( const double d );
	/// Division by a constant
	vector2& operator /= ( const double d );
	/// Returns an indexed dimension by reference
	double& operator [] ( int i);
	/// Returns an indexed dimension by value
	double vector2::operator [] ( int i) const;
	/// Returns the vector length
	double Length() const;
	/// Returns the square of the vector length
	double Length2() const;
	/// Normalizes the vector in place
	vector2& Normalize();
	/// Returns the perpendicular to this vector
	vector2 Perpendicular() const;
	/// Returns the quadrant containing this vector
	long Quadrant();
	/// Returns the angle of this vector with respect to the positive X axis in radians [-pi, pi]
	double Angle() const;
	/// Casting to a C/C++ style array pointer
	operator double*() { return &n[0]; }
	/// Copies the vector into a C/C++ style array
	void CopyArray(float f[2]) const;
	/// Copies the vector into a C/C++ style array
	void CopyArray(double d[2]) const;

	friend std::ostream& operator<<(std::ostream& Stream, const vector2& RHS)
	{
		Stream << RHS.n[0] << " " << RHS.n[1];
		return Stream;
	}

	friend std::istream& operator>>(std::istream& Stream, vector2& RHS)
	{
		Stream >> RHS.n[0];
		RHS.n[1] = RHS.n[0];
		Stream >> RHS.n[1];

		return Stream;
	}

};

/// Negation
vector2 operator - (const vector2& v);
/// Addition
vector2 operator + (const vector2& a, const vector2& b);
/// Subtraction
vector2 operator - (const vector2& a, const vector2& b);
/// Multiplication by a constant
vector2 operator * (const vector2& a, const double d);
/// Multiplication by a constant
vector2 operator * (const double d, const vector2& a);

/// Returns the dot product of two vectors
double operator * (const vector2& a, const vector2& b);
/// Division by a constant
vector2 operator / (const vector2& a, const double d);
/// Returns the cross product of two vectors
vector3 operator ^ (const vector2& a, const vector2& b);
/// Tests for equality
bool operator == (const vector2& a, const vector2& b);
/// Tests for non-equality
bool operator != (const vector2& a, const vector2& b);
/// Returns the term-by-term minimum of two vectors
vector2 vectorMin(const vector2& a, const vector2& b);
/// Returns the term-by-term maximum of two vectors
vector2 vectorMax(const vector2& a, const vector2& b);
/// Returns the term-by-term product of two vectors
vector2 Product(const vector2& a, const vector2& b);

/// An array of two-dimensional vectors
typedef std::vector<vector2> vector2Array;

/// Tests to see if the point is "within" the area defined by the array of points
bool IsInside(const vector2 Point, vector2Array Array);
/// Returns the array point closest to the given coordinate
vector2 ClosestPoint(vector2& Point, vector2Array& Array);

/////////////////////////////////////////////////////////////////////////////
// vector3

// We have an unfortunate clash with <windows.h> ... go figure
#ifdef WIN32
#ifdef RGB
#undef RGB
#endif // RGB
#endif // WIN32

/// A three-dimensional vector
class vector3
{
public:
	/// Stores the vector values
	double n[3];

	// Constructors
	vector3();
	vector3(const double x, const double y, const double z);
	vector3(const double d);
	vector3(const double d[3]);
	/// Copy constructor
	vector3(const vector3& v);
	/// Casts an vector2 to an vector3 (with zero third dimension)
	vector3(const vector2& v);
	/// Casts an vector2 to an vector3 with the given third dimension
	vector3(const vector2& v, double d);
	/// Casts an vector4 to an vector3, dropping the fourth dimension
	vector3(const vector4& v);
	/// Casts an vector4 to an vector3, dropping the given dimension
	vector3(const vector4& v, int dropAxis);
	/// Assignment of an vector3
	vector3& operator	= ( const vector3& v );
	/// Assignment of a C/C++ style array
	vector3& operator  = ( const double d[3]);
	/// Addition
	vector3& operator += ( const vector3& v );
	/// Subtraction
	vector3& operator -= ( const vector3& v );
	/// Multiplication by a constant
	vector3& operator *= ( const double d );
	/// Division by a constant
	vector3& operator /= ( const double d );
	/// Returns the given dimension by reference
	double& operator [] ( int i);
	/// Returns the given dimension by value
	double operator[] (int i) const;
	/// Returns the vector length
	double Length() const;
	/// Returns the squared vector length
	double Length2() const;
	/// Normalizes the vector in place
	vector3& Normalize();
	/// Converts Cartesian coordinates to spherical coordinates
	vector3 Spherical() const;
	/// Sets each dimension to its absolute value in place
	vector3& Abs();
	/// Casting to a C/C++ style array pointer
	/** \deprecated */
	operator double*() { return &n[0]; }
	/// A better way to integrate with legacy array-oriented APIs
	const double* data() const { return &n[0]; }
	/// A better way to integrate with legacy array-oriented APIs
	double* data() { return &n[0]; }
	/// Copies the vector into a C/C++ style array
	void CopyArray(float f[3]) const;
	/// Copies the vector into a C/C++ style array
	void CopyArray(double d[3]) const;

	friend std::ostream& operator<<(std::ostream& Stream, const vector3& RHS)
	{
		Stream << RHS.n[0] << " " << RHS.n[1] << " " << RHS.n[2];
		return Stream;
	}

	friend std::istream& operator>>(std::istream& Stream, vector3& RHS)
	{
		Stream >> RHS.n[0];
		RHS.n[1] = RHS.n[2] = RHS.n[0];
		Stream >> RHS.n[1] >> RHS.n[2];

		return Stream;
	}

};

/// Negation
vector3 operator - (const vector3& v);
/// Addition
vector3 operator + (const vector3& a, const vector3& b);
/// Subtraction
vector3 operator - (const vector3& a, const vector3& b);
/// Multiplication by a constant
vector3 operator * (const vector3& a, const double d);
/// Multiplication by a constant
vector3 operator * (const double d, const vector3& a);
/// Returns the dot product of two vectors
double operator * (const vector3& a, const vector3& b);
/// Division by a constant
vector3 operator / (const vector3& a, const double d);
/// Returns the cross product of two vectors
vector3 operator ^ (const vector3& a, const vector3& b);
/// Equality
bool operator == (const vector3& a, const vector3& b);
/// Non-equality
bool operator != (const vector3& a, const vector3& b);
/// Returns the term-by-term minimum of two vectors
vector3 vectorMin(const vector3& a, const vector3& b);
/// Returns the term-by-term maximum of two vectors
vector3 vectorMax(const vector3& a, const vector3& b);
/// Returns the term-by-term product of two vectors
vector3 Product(const vector3& a, const vector3& b);

/// An array of three-dimensional vectors
typedef std::vector<vector3> vector3Array;

/////////////////////////////////////////////////////////////////////////////
// vector4

/// A four-dimensional vector
class vector4
{
public:
	/// Stores the vector values
	double n[4];

	// Constructors
	vector4();
	vector4(const double x, const double y, const double z, const double w);
	vector4(const double d);
	vector4(const double d[4]);
	/// Copy Constructor
	vector4(const vector4& v);
	/// Casts an vector3 to an vector4 with zero fourth dimension
	vector4(const vector3& v);
	/// Casts an vector3 to an vector4 with the given fourth dimension
	vector4(const vector3& v, const double d);
	/// Assignment of an vector4
	vector4& operator	= ( const vector4& v );
	/// Assignment of a C/C++ array
	vector4& operator  = ( const double d[4]);
	/// Addition
	vector4& operator += ( const vector4& v );
	/// Subtraction
	vector4& operator -= ( const vector4& v );
	/// Multiplication by a constant
	vector4& operator *= ( const double d );
	/// Division by a constant
	vector4& operator /= ( const double d );
	/// Returns a vector element by reference
	double& operator [] ( int i);
	/// Returns a vector element by value
	double operator[] (int i) const;
	/// Returns the length of a vector
	double Length() const;
	/// Returns the square of the vector length
	double Length2() const;
	/// Normalizes the vector in place
	vector4& Normalize();
	/// Casts the contents of the vector to a C/C++ array
	operator double*() { return &n[0]; }

	friend std::ostream& operator<<(std::ostream& Stream, const vector4& RHS)
	{
		Stream << RHS.n[0] << " " << RHS.n[1] << " " << RHS.n[2] << " " << RHS.n[3];
		return Stream;
	}

	friend std::istream& operator>>(std::istream& Stream, vector4& RHS)
	{
		Stream >> RHS.n[0];
		RHS.n[1] = RHS.n[2] = RHS.n[3] = RHS.n[0];
		Stream >> RHS.n[1] >> RHS.n[2] >> RHS.n[3];

		return Stream;
	}

};

/// Negation
vector4 operator - (const vector4& v);
/// Addition
vector4 operator + (const vector4& a, const vector4& b);
/// Subtraction
vector4 operator - (const vector4& a, const vector4& b);
/// Multiplication by a constant
vector4 operator * (const vector4& a, const double d);
/// Multiplication by a constant
vector4 operator * (const double d, const vector4& a);
/// Returns the dot product of two vectors
double operator * (const vector4& a, const vector4& b);
/// Division by a constant
vector4 operator / (const vector4& a, const double d);
/// Equality
bool operator == (const vector4& a, const vector4& b);
/// Non-equality
bool operator != (const vector4& a, const vector4& b);
/// Returns the element-by-element minimum of two vectors
vector4 vectorMin(const vector4& a, const vector4& b);
/// Returns the element-by-element maximum of two vectors
vector4 vectorMax(const vector4& a, const vector4& b);
/// Returns the element-by-element product of two vectors
vector4 Product(const vector4& a, const vector4& b);


/// Returns the normalized form of a vector ...
template<typename T>
const T normalize(const T& Vector)
{
	return Vector / Vector.Length();
}

/////////////////////////////////////////////////////////////////////////////
// vector2 implementation

inline vector2::vector2()
{ n[VX] = n[VY] = 0.0; }

inline vector2::vector2(const double x, const double y)
{ n[VX] = x; n[VY] = y; }

inline vector2::vector2(const double d)
{ n[VX] = n[VY] = d; }

inline vector2::vector2(const double d[2])
{ n[VX] = d[0]; n[VY] = d[1]; }

inline vector2::vector2(const vector2& v)
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; }

inline vector2::vector2(const vector3& v) // it is up to caller to avoid divide-by-zero
{ n[VX] = v.n[VX]/v.n[VZ]; n[VY] = v.n[VY]/v.n[VZ]; }

inline vector2::vector2(const vector3& v, int dropAxis) {
 switch (dropAxis) {
	case VX: n[VX] = v.n[VY]; n[VY] = v.n[VZ]; break;
	case VY: n[VX] = v.n[VX]; n[VY] = v.n[VZ]; break;
	default: n[VX] = v.n[VX]; n[VY] = v.n[VY]; break;
 }
}

inline vector2& vector2::operator = (const vector2& v)
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; return *this; }

inline vector2& vector2::operator = (const double d[2])
{ n[VX] = d[0]; n[VY] = d[1]; return *this; }

inline vector2& vector2::operator += ( const vector2& v )
{ n[VX] += v.n[VX]; n[VY] += v.n[VY]; return *this; }

inline vector2& vector2::operator -= ( const vector2& v )
{ n[VX] -= v.n[VX]; n[VY] -= v.n[VY]; return *this; }

inline vector2& vector2::operator *= ( const double d )
{ n[VX] *= d; n[VY] *= d; return *this; }

inline vector2& vector2::operator /= ( const double d )
{ double d_inv = 1./d; n[VX] *= d_inv; n[VY] *= d_inv; return *this; }

inline double& vector2::operator [] ( int i)
{
 assert(!(i < VX || i > VY));
 return n[i];
}

inline double vector2::operator [] ( int i) const
{
 assert(!(i < VX || i > VY));
 return n[i];
}

inline double vector2::Length() const
{ return sqrt(Length2()); }

inline double vector2::Length2() const
{ return n[VX]*n[VX] + n[VY]*n[VY]; }

inline vector2& vector2::Normalize()
{
	if(const double length = Length())
		*this /= length;

	return *this;
}

inline vector2 vector2::Perpendicular() const
{
	return vector2(-n[1], n[0]);
}

inline long vector2::Quadrant()
{
	if(n[VX] < 0.0)
		return n[VY] < 0.0 ? 2 : 1;

	return n[VY] < 0.0 ? 3 : 0;
}

inline double vector2::Angle() const
{
	return atan2(n[VY], n[VX]);
}

inline void vector2::CopyArray(float f[2]) const
{	f[0] = float(n[0]); f[1] = float(n[1]); }

inline void vector2::CopyArray(double d[2]) const
{	d[0] = n[0]; d[1] = n[1]; }

inline bool IsInside(const vector2 Point, vector2Array Array)
{
	// Sanity check ...
	if(!Array.size())
		{
			std::cerr << __PRETTY_FUNCTION__ << " : empty array." << std::endl;
			return false;
		}

	long winding = 0;

	// Close the set defined by this list ...
	Array.push_back(Array[0]);

	// Get the first point ...
	vector2 testpoint(Array[0] - Point);
	long lastquadrant = testpoint.Quadrant();

	// For each remaining point ...
	for(unsigned long i = 1; i < Array.size(); ++i)
	{
		testpoint = vector2(Array[i] - Point);
		long quadrant = testpoint.Quadrant();
		if(quadrant == lastquadrant) continue;

		if(((lastquadrant + 1) & 3) == quadrant)
		{
			++winding;
		}
		else if(((quadrant + 1) & 3) == lastquadrant)
		{
			--winding;
		}
		else
			{
				double right = Array[i-1][0] - Array[i][0];
				double left = (Array[i-1][1] - Array[i][1]) * ((Point)[0] - Array[i-1][0]) + (Array[i-1][1] * right);
				winding += (left >= right * (Point)[1]) ? 2 : -2;
			}

		lastquadrant = quadrant;
	}

	return winding ? true : false;
}

inline vector2 ClosestPoint(vector2& Point, vector2Array& Array)
{
	vector2 result;
	double distance2;
	double distance;

	assert(Array.size());

	for(unsigned long i = 0; i < Array.size(); ++i)
	{
		distance2 = (Array[i] - Point).Length2();
		if(i == 0 || distance2 <= distance)
		{
			distance = distance2;
			result = Array[i];
		}
	}

	return result;
}

inline vector2 operator - (const vector2& a)
{ return vector2(-a.n[VX],-a.n[VY]); }

inline vector2 operator + (const vector2& a, const vector2& b)
{ return vector2(a.n[VX]+ b.n[VX], a.n[VY] + b.n[VY]); }

inline vector2 operator - (const vector2& a, const vector2& b)
{ return vector2(a.n[VX]-b.n[VX], a.n[VY]-b.n[VY]); }

inline vector2 operator * (const vector2& a, const double d)
{ return vector2(d*a.n[VX], d*a.n[VY]); }

inline vector2 operator * (const double d, const vector2& a)
{ return a*d; }

inline double operator * (const vector2& a, const vector2& b)
{ return (a.n[VX]*b.n[VX] + a.n[VY]*b.n[VY]); }

inline vector2 operator / (const vector2& a, const double d)
{ double d_inv = 1./d; return vector2(a.n[VX]*d_inv, a.n[VY]*d_inv); }

inline vector3 operator ^ (const vector2& a, const vector2& b)
{ return vector3(0.0, 0.0, a.n[VX] * b.n[VY] - b.n[VX] * a.n[VY]); }

inline bool operator == (const vector2& a, const vector2& b)
{ return (a.n[VX] == b.n[VX]) && (a.n[VY] == b.n[VY]); }

inline bool operator != (const vector2& a, const vector2& b)
{ return !(a == b); }

inline vector2 vectorMin(const vector2& a, const vector2& b)
{ return vector2(sdpMin(a.n[VX], b.n[VX]), sdpMin(a.n[VY], b.n[VY])); }

inline vector2 vectorMax(const vector2& a, const vector2& b)
{ return vector2(sdpMax(a.n[VX], b.n[VX]), sdpMax(a.n[VY], b.n[VY])); }

inline vector2 Product(const vector2& a, const vector2& b)
{ return vector2(a.n[VX] * b.n[VX], a.n[VY] * b.n[VY]); }

/////////////////////////////////////////////////////////////////////////////
// vector3 implementation

inline vector3::vector3()
{ n[VX] = n[VY] = n[VZ] = 0.0; }

inline vector3::vector3(const double x, const double y, const double z)
{ n[VX] = x; n[VY] = y; n[VZ] = z; }

inline vector3::vector3(const double d)
{ n[VX] = n[VY] = n[VZ] = d; }

inline vector3::vector3(const double d[3])
{ n[VX] = d[0]; n[VY] = d[1]; n[VZ] = d[2]; }

inline vector3::vector3(const vector3& v)
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; }

inline vector3::vector3(const vector2& v)
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = 1.0; }

inline vector3::vector3(const vector2& v, double d)
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = d; }

inline vector3::vector3(const vector4& v) // it is up to caller to avoid divide-by-zero
{ n[VX] = v.n[VX] / v.n[VW]; n[VY] = v.n[VY] / v.n[VW];
 n[VZ] = v.n[VZ] / v.n[VW]; }

inline vector3::vector3(const vector4& v, int dropAxis) {
 switch (dropAxis) {
	case VX: n[VX] = v.n[VY]; n[VY] = v.n[VZ]; n[VZ] = v.n[VW]; break;
	case VY: n[VX] = v.n[VX]; n[VY] = v.n[VZ]; n[VZ] = v.n[VW]; break;
	case VZ: n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VW]; break;
	default: n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; break;
 }
}

inline vector3& vector3::operator = (const vector3& v)
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; return *this; }

inline vector3& vector3::operator = (const double d[3])
{ n[VX] = d[0]; n[VY] = d[1]; n[VZ] = d[2]; return *this; }

inline vector3& vector3::operator += ( const vector3& v )
{ n[VX] += v.n[VX]; n[VY] += v.n[VY]; n[VZ] += v.n[VZ]; return *this; }

inline vector3& vector3::operator -= ( const vector3& v )
{ n[VX] -= v.n[VX]; n[VY] -= v.n[VY]; n[VZ] -= v.n[VZ]; return *this; }

inline vector3& vector3::operator *= ( const double d )
{ n[VX] *= d; n[VY] *= d; n[VZ] *= d; return *this; }

inline vector3& vector3::operator /= ( const double d )
{ double d_inv = 1./d; n[VX] *= d_inv; n[VY] *= d_inv; n[VZ] *= d_inv;
 return *this; }

inline double& vector3::operator [] ( int i) {
 assert(! (i < VX || i > VZ));
 return n[i];
}

inline double vector3::operator [] ( int i) const {
 assert(! (i < VX || i > VZ));
 return n[i];
}

inline double vector3::Length() const
{ return sqrt(Length2()); }

inline double vector3::Length2() const
{ return n[VX]*n[VX] + n[VY]*n[VY] + n[VZ]*n[VZ]; }

inline vector3& vector3::Normalize()
{
	if(const double length = Length())
		*this /= length;

	return *this;
}

inline void vector3::CopyArray(float f[3]) const
{	f[0] = (float)n[0];	f[1] = (float)n[1];	f[2] = (float)n[2]; }

inline void vector3::CopyArray(double d[3]) const
{	d[0] = n[0]; d[1] = n[1]; d[2] = n[2]; }

inline vector3 vector3::Spherical() const
{
	return vector3(Length(), atan2(n[VX], n[VZ]), atan2(n[VY], sqrt(n[VX] * n[VX] + n[VZ] * n[VZ])));
/*
	// Handle the singularity at the poles
	if(0.0 == n[VX] && 0.0 == n[VZ])
		return vector3(n[VY] > 0.0 ? sdpPiOver2 : -sdpPiOver2, 0.0, Length());


	return vector3(atan2(n[VY], sqrt(n[VX] * n[VX] + n[VZ] * n[VZ])), atan2(n[VX], n[VZ]), Length());
*/
}

inline vector3& vector3::Abs()
{	n[0] = fabs(n[0]); n[1] = fabs(n[1]); n[2] = fabs(n[2]); return *this; }

inline vector3 operator - (const vector3& a)
{ return vector3(-a.n[VX],-a.n[VY],-a.n[VZ]); }

inline vector3 operator + (const vector3& a, const vector3& b)
{ return vector3(a.n[VX]+ b.n[VX], a.n[VY] + b.n[VY], a.n[VZ] + b.n[VZ]); }

inline vector3 operator - (const vector3& a, const vector3& b)
{ return vector3(a.n[VX]-b.n[VX], a.n[VY]-b.n[VY], a.n[VZ]-b.n[VZ]); }

inline vector3 operator * (const vector3& a, const double d)
{ return vector3(d*a.n[VX], d*a.n[VY], d*a.n[VZ]); }

inline vector3 operator * (const double d, const vector3& a)
{ return a*d; }

inline double operator * (const vector3& a, const vector3& b)
{ return (a.n[VX]*b.n[VX] + a.n[VY]*b.n[VY] + a.n[VZ]*b.n[VZ]); }

inline vector3 operator / (const vector3& a, const double d)
{ double d_inv = 1./d; return vector3(a.n[VX]*d_inv, a.n[VY]*d_inv,
 a.n[VZ]*d_inv); }

inline vector3 operator ^ (const vector3& a, const vector3& b) {
 return vector3(a.n[VY]*b.n[VZ] - a.n[VZ]*b.n[VY],
		a.n[VZ]*b.n[VX] - a.n[VX]*b.n[VZ],
		a.n[VX]*b.n[VY] - a.n[VY]*b.n[VX]);
}

inline bool operator == (const vector3& a, const vector3& b)
{ return (a.n[VX] == b.n[VX]) && (a.n[VY] == b.n[VY]) && (a.n[VZ] == b.n[VZ]);
}

inline bool operator != (const vector3& a, const vector3& b)
{ return !(a == b); }

inline vector3 vectorMin(const vector3& a, const vector3& b)
{ return vector3(sdpMin(a.n[VX], b.n[VX]), sdpMin(a.n[VY], b.n[VY]), sdpMin(a.n[VZ],
 b.n[VZ])); }

inline vector3 vectorMax(const vector3& a, const vector3& b)
{ return vector3(sdpMax(a.n[VX], b.n[VX]), sdpMax(a.n[VY], b.n[VY]), sdpMax(a.n[VZ],
 b.n[VZ])); }

inline vector3 Product(const vector3& a, const vector3& b)
{ return vector3(a.n[VX] * b.n[VX], a.n[VY] * b.n[VY], a.n[VZ] * b.n[VZ]); }

/////////////////////////////////////////////////////////////////////////////
// vector4 implementation

inline vector4::vector4()
{ n[VX] = n[VY] = n[VZ] = n[VW] = 0.0; }

inline vector4::vector4(const double x, const double y, const double z, const double w)
{ n[VX] = x; n[VY] = y; n[VZ] = z; n[VW] = w; }

inline vector4::vector4(const double d)
{ n[VX] = n[VY] = n[VZ] = n[VW] = d; }

inline vector4::vector4(const double d[4])
{ n[VX] = d[0]; n[VY] = d[1]; n[VZ] = d[2]; n[VW] = d[3]; }

inline vector4::vector4(const vector4& v)
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; n[VW] = v.n[VW]; }

inline vector4::vector4(const vector3& v)
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; n[VW] = 1.0; }

inline vector4::vector4(const vector3& v, const double d)
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; n[VW] = d; }

inline vector4& vector4::operator = (const vector4& v)
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; n[VW] = v.n[VW];
return *this; }

inline vector4& vector4::operator = (const double d[4])
{ n[VX] = d[0]; n[VY] = d[1]; n[VZ] = d[2]; n[VW] = d[3]; return *this; }

inline vector4& vector4::operator += ( const vector4& v )
{ n[VX] += v.n[VX]; n[VY] += v.n[VY]; n[VZ] += v.n[VZ]; n[VW] += v.n[VW];
return *this; }

inline vector4& vector4::operator -= ( const vector4& v )
{ n[VX] -= v.n[VX]; n[VY] -= v.n[VY]; n[VZ] -= v.n[VZ]; n[VW] -= v.n[VW];
return *this; }

inline vector4& vector4::operator *= ( const double d )
{ n[VX] *= d; n[VY] *= d; n[VZ] *= d; n[VW] *= d; return *this; }

inline vector4& vector4::operator /= ( const double d )
{ double d_inv = 1./d; n[VX] *= d_inv; n[VY] *= d_inv; n[VZ] *= d_inv;
 n[VW] *= d_inv; return *this; }

inline double& vector4::operator [] ( int i) {
 assert(! (i < VX || i > VW));
 return n[i];
}

inline double vector4::operator [] ( int i) const {
 assert(! (i < VX || i > VW));
 return n[i];
}

inline double vector4::Length() const
{ return sqrt(Length2()); }

inline double vector4::Length2() const
{ return n[VX]*n[VX] + n[VY]*n[VY] + n[VZ]*n[VZ] + n[VW]*n[VW]; }

inline vector4& vector4::Normalize()
{
	if(const double length = Length())
		*this /= length;

	return *this;
}

inline vector4 operator - (const vector4& a)
{ return vector4(-a.n[VX],-a.n[VY],-a.n[VZ],-a.n[VW]); }

inline vector4 operator + (const vector4& a, const vector4& b)
{ return vector4(a.n[VX] + b.n[VX], a.n[VY] + b.n[VY], a.n[VZ] + b.n[VZ],
 a.n[VW] + b.n[VW]); }

inline vector4 operator - (const vector4& a, const vector4& b)
{ return vector4(a.n[VX] - b.n[VX], a.n[VY] - b.n[VY], a.n[VZ] - b.n[VZ],
 a.n[VW] - b.n[VW]); }

inline vector4 operator * (const vector4& a, const double d)
{ return vector4(d*a.n[VX], d*a.n[VY], d*a.n[VZ], d*a.n[VW] ); }

inline vector4 operator * (const double d, const vector4& a)
{ return a*d; }

inline double operator * (const vector4& a, const vector4& b)
{ return (a.n[VX]*b.n[VX] + a.n[VY]*b.n[VY] + a.n[VZ]*b.n[VZ] +
 a.n[VW]*b.n[VW]); }

inline vector4 operator / (const vector4& a, const double d)
{ double d_inv = 1./d; return vector4(a.n[VX]*d_inv, a.n[VY]*d_inv, a.n[VZ]*d_inv,
 a.n[VW]*d_inv); }

inline bool operator == (const vector4& a, const vector4& b)
{ return (a.n[VX] == b.n[VX]) && (a.n[VY] == b.n[VY]) && (a.n[VZ] == b.n[VZ])
 && (a.n[VW] == b.n[VW]); }

inline bool operator != (const vector4& a, const vector4& b)
{ return !(a == b); }

inline vector4 vectorMin(const vector4& a, const vector4& b)
{ return vector4(sdpMin(a.n[VX], b.n[VX]), sdpMin(a.n[VY], b.n[VY]), sdpMin(a.n[VZ],
 b.n[VZ]), sdpMin(a.n[VW], b.n[VW])); }

inline vector4 vectorMax(const vector4& a, const vector4& b)
{ return vector4(sdpMax(a.n[VX], b.n[VX]), sdpMax(a.n[VY], b.n[VY]), sdpMax(a.n[VZ],
 b.n[VZ]), sdpMax(a.n[VW], b.n[VW])); }

inline vector4 Product(const vector4& a, const vector4& b)
{ return vector4(a.n[VX] * b.n[VX], a.n[VY] * b.n[VY], a.n[VZ] * b.n[VZ],
 a.n[VW] * b.n[VW]); }

} // namespace k3d

#endif // VECTORS_H

