//
// ======================================================
//
//       Interface for Class
//
//          Vec4
//
// ======================================================
//

#ifndef JL_VEC4_H
#define JL_VEC4_H

#include "assert.h"
#include <iostream.h>
#include <math.h>
#include "jllib_include.h"
#include "matrix.h"

#define SQRT_TABLE_SIZE 1000

typedef Vec4 Plane;
typedef Vec4 Point3d;

class Vec4 {

	friend class Quaternion;

	private:
		real _data[4];

		static float _float_array[3];
		static real _sqrt_table[SQRT_TABLE_SIZE+2];

	public:

// *** Constructors/Destructor: ***
		inline Vec4() { _data[0] = _data[1] = _data[2] = _data[3] = 0; }

		inline Vec4(real x, real y, real z, real w=0) {
			_data[0]=x; _data[1]=y; _data[2]=z; _data[3]=w; }

		Vec4(unsigned long);

		inline Vec4(const Vec4 &v) {
				_data[0] = v._data[0]; _data[1] = v._data[1];
				_data[2] = v._data[2]; _data[3] = v._data[3];
			}

		Vec4(const Vec4 &v1, const Vec4 &v2) {
			_data[0] = v1._data[0] - v2._data[0];
			_data[1] = v1._data[1] - v2._data[1];
			_data[2] = v1._data[2] - v2._data[2];
			_data[3] = v1._data[3] - v2._data[3];
		}

		inline ~Vec4() { }

// *** Operators: ***
		inline Vec4& operator=(const Vec4& v) {
				_data[0] = v._data[0]; _data[1] = v._data[1];
				_data[2] = v._data[2]; _data[3] = v._data[3];
				return *this;
			}
		Vec4* Copy() const { return new Vec4(*this); }
		friend int operator==(const Vec4&, const Vec4&);
		friend int operator!=(const Vec4&, const Vec4&);
		friend ostream& operator<<(ostream&, const Vec4&);
		friend Vec4 operator+(const Vec4&, const Vec4&);
		friend Vec4 operator-(const Vec4&, const Vec4&);
		friend Vec4 operator*(Vec4, const real&);
		friend Vec4 operator*(const real&, Vec4);
//		friend Vec4 operator/(Vec4, const real&);
		Vec4& operator+=(const Vec4 &v) {
				_data[0] += v._data[0];
				_data[1] += v._data[1];
				_data[2] += v._data[2];
				_data[3] += v._data[3];
				return *this;
			}
	Vec4& operator-=(const Vec4&);
		Vec4& operator*=(const JLmatrix&);
	inline void Negate() {
			_data[0] = -_data[0];
			_data[1] = -_data[1];
			_data[2] = -_data[2];
			// *** I don't think I want to negate w... ***
		}


	inline void TotalNegate() {
			_data[0] = -_data[0];
			_data[1] = -_data[1];
			_data[2] = -_data[2];
                        _data[3] = -_data[3]; //I do. jpa
		}

		inline void operator*=(real d) {
				_data[0] *= d, _data[1] *= d, _data[2] *= d; _data[3] *= d;
			}

	inline void operator/=(real d) {
		_data[0] /= d, _data[1] /= d, _data[2] /= d; _data[3] /= d;
		}

	inline Vec4& SetWto1() { _data[3] = 1; return *this; }

	inline friend void
		Vec4FastAdd(Vec4 &sum, const Vec4 &v1, const Vec4 &v2) {
			sum._data[0] = v1._data[0] + v2._data[0];
			sum._data[1] = v1._data[1] + v2._data[1];
			sum._data[2] = v1._data[2] + v2._data[2];
			sum._data[3] = 1;
		}
	inline friend void
		Vec4FastSub(Vec4 &diff, const Vec4 &v1, const Vec4 &v2) {
				diff._data[0] = v1._data[0] - v2._data[0];
				diff._data[1] = v1._data[1] - v2._data[1];
				diff._data[2] = v1._data[2] - v2._data[2];
				diff._data[3] = 1;
			}

	inline friend void
		Vec4FastAdd4(Vec4 &sum, const Vec4 &v1, const Vec4 &v2) {
			sum._data[0] = v1._data[0] + v2._data[0];
			sum._data[1] = v1._data[1] + v2._data[1];
			sum._data[2] = v1._data[2] + v2._data[2];
			sum._data[3] = v1._data[3] + v2._data[3];
		}
	inline friend void
		Vec4FastSub4(Vec4 &diff, const Vec4 &v1, const Vec4 &v2) {
				diff._data[0] = v1._data[0] - v2._data[0];
				diff._data[1] = v1._data[1] - v2._data[1];
				diff._data[2] = v1._data[2] - v2._data[2];
				diff._data[3] = v1._data[3] - v2._data[3];
			}

	inline void Vec4FastScale(real x, real y, real z) {
			_data[0] *= x;
			_data[1] *= y;
			_data[2] *= z;
		}

	inline void operator *=(const Vec4 &v2) {
			_data[0] *= v2._data[0];
			_data[1] *= v2._data[1];
			_data[2] *= v2._data[2];
		}

        inline Vec4 operator*(const Vec4 &v2) {
                return Vec4 (_data[0] * v2._data[0],
                _data[1] * v2._data[1],
                _data[2] * v2._data[2], 1.);
		}

	inline friend void Vec4FastCopyScale(Vec4 &v1, const Vec4 &v2, real d) {
			v1._data[0] = v2._data[0] * d;
			v1._data[1] = v2._data[1] * d;
			v1._data[2] = v2._data[2] * d;
		}

	inline friend void Vec4FastAddScale(
		Vec4 &sum, const Vec4 &v1, const Vec4 &v2, real d) {
			sum._data[0] = v1._data[0] + v2._data[0] * d;
			sum._data[1] = v1._data[1] + v2._data[1] * d;
			sum._data[2] = v1._data[2] + v2._data[2] * d;
			sum._data[3] = 1;
		}

	inline friend void Vec3_FastWeightedSum(
		Vec4 &sum, const Vec4 &v1, real w1, const Vec4 &v2, real w2) {
			sum._data[0] = v1._data[0]*w1 + v2._data[0]*w2;
			sum._data[1] = v1._data[1]*w1 + v2._data[1]*w2;
			sum._data[2] = v1._data[2]*w1 + v2._data[2]*w2;
			sum._data[3] = 1;
		}

	inline friend void Vec4_FastWeightedSum(
		Vec4 &sum, const Vec4 &v1, real w1, const Vec4 &v2, real w2) {
			sum._data[0] = v1._data[0]*w1 + v2._data[0]*w2;
			sum._data[1] = v1._data[1]*w1 + v2._data[1]*w2;
			sum._data[2] = v1._data[2]*w1 + v2._data[2]*w2;
			sum._data[3] = v1._data[3]*w1 + v2._data[3]*w2;
		}

	inline void FastJLmatrixMultiply(const JLmatrix &m) {
			real v0 = 0, v1 = 0, v2 = 0, v3 = 0;
			v0 = _data[0] * m._data[0] +
			 _data[1] * m._data[4] +
			 _data[2] * m._data[8] +
			 m._data[12];
					v1 = _data[0] * m._data[1] +
						 _data[1] * m._data[5] +
						 _data[2] * m._data[9] +
						 m._data[13];
					v2 = _data[0] * m._data[2] +
						 _data[1] * m._data[6] +
						 _data[2] * m._data[10] +
						 m._data[14];
					v3 = _data[0] * m._data[3] +
						 _data[1] * m._data[7] +
						 _data[2] * m._data[11] +
						 m._data[15];
			_data[0] = v0; _data[1] = v1; _data[2] = v2; _data[3] = v3;
		}

	// v = *this * m, then divide by w
	inline void
		FastJLmatrixMultiplyAndDivideByW(Vec4 &v, const JLmatrix &m) const {
				real w = _data[0] * m._data[3] +
						 _data[1] * m._data[7] +
						 _data[2] * m._data[11] +
						 m._data[15];
				if (w) {
					v._data[0] = (_data[0] * m._data[0] +
					 			  _data[1] * m._data[4] +
					 			  _data[2] * m._data[8] +
					 			  m._data[12]) / w;
					v._data[1] = (_data[0] * m._data[1] +
						 		  _data[1] * m._data[5] +
					 			  _data[2] * m._data[9] +
					 			  m._data[13]) / w;
					v._data[2] = (_data[0] * m._data[2] +
					 			  _data[1] * m._data[6] +
					 			  _data[2] * m._data[10] +
					 			  m._data[14]) / w;
				}
				else v._data[0] = v._data[1] = v._data[2] = 0;
				v._data[3] = 1;
			}

	inline friend void
		Vec4FastAverage(Vec4 &a, const Vec4 &v1, const Vec4 &v2) {
			a._data[0] = ((v1._data[0] + v2._data[0]) * 0.5);
			a._data[1] = ((v1._data[1] + v2._data[1]) * 0.5);
			a._data[2] = ((v1._data[2] + v2._data[2]) * 0.5);
			a._data[3] = 1;
		}

// *** Convert from homogeneous coordinates to 3D: ***
	inline void DivideByW() {
			if (_data[3]) {
			_data[0] /= _data[3];
			_data[1] /= _data[3];
			_data[2] /= _data[3];
			}
			else
			_data[0] = _data[1] = _data[2] = 0;
			_data[3] = 1;
		}

	inline void Mult_XYZ_By_W() {
			_data[0] *= _data[3];
			_data[1] *= _data[3];
			_data[2] *= _data[3];
		}

// *** Some nice standard vector functions: ***

	inline void MakeUnit() {
			real l = fsqrt(_data[0]*_data[0] +
					_data[1]*_data[1] +
					_data[2]*_data[2]);
			if (l) {
			_data[0] /= l;
			_data[1] /= l;
			_data[2] /= l;
			}
			//else
			//_data[0] = _data[1] = _data[2] = 0;
		}

	inline void FastMakeUnit3() { // SJT
			real l = 1.0F / (fsqrt(_data[0]*_data[0] +
					_data[1]*_data[1] +
					_data[2]*_data[2]));
			_data[0] *= l;
			_data[1] *= l;
			_data[2] *= l;
		}

//divide everything by a factor s.t. the first three are normalized.
inline void PlaneNormalize (void) {
			real l = 1.0F / (fsqrt(_data[0]*_data[0] +
					_data[1]*_data[1] +
					_data[2]*_data[2]));
			_data[0] *= l;
			_data[1] *= l;
			_data[2] *= l;
			_data[3] *= l;
		}


	inline void FastMakeUnit4() { // SJT/JPA
			real l = 1.0F / (fsqrt(_data[0]*_data[0] +
					_data[1]*_data[1] +
					_data[2]*_data[2] + 
					_data[3]*_data[3]));
			_data[0] *= l;
			_data[1] *= l;
			_data[2] *= l;
			_data[3] *= l;
		}

	inline Vec4& MakeAndGetUnit() { MakeUnit(); return *this; }

	static void MakeSqrtTable();

// # define FastMakeUnit MakeUnit
/*
	void FastMakeUnit() {
		real approx_l = ABS(_data[0]) + ABS(_data[1]) + ABS(_data[2]);
		real w = ((_data[0] * _data[0]) +
				  (_data[1] * _data[1]) +
				  (_data[2] * _data[2])) /
					(approx_l * approx_l) *
					SQRT_TABLE_SIZE;
		int i = (int) w;
		w = w-i;
		real l = ((Vec4::_sqrt_table[i+1] * w) +
				  (Vec4::_sqrt_table[i] * (1-w))) * approx_l;
		if (l) {
			_data[0] /= l;
			_data[1] /= l;
			_data[2] /= l;
		}
		//else _data[0] = _data[1] = _data[2] = 0;
	}
*/

		inline real Length() const {
			return fsqrt(_data[0] * _data[0] +
				_data[1] * _data[1] +
				_data[2] * _data[2]);
		}

		inline real DistanceTo(const Vec4 &v) const {
			Vec4 tmp(*this, v);
			return tmp.Length();
		}
			

		inline real Dot3(const Vec4 &v) const {
			return _data[0] * v._data[0] +
			   _data[1] * v._data[1] +
			   _data[2] * v._data[2];
		}

		inline real Dot4(const Vec4 &v) const {
			return _data[0] * v._data[0] +
			   _data[1] * v._data[1] +
			   _data[2] * v._data[2] +
			   _data[3] * v._data[3];
		}

        // interpret this as 4-coefficient plane equation
        // interpret pt as 3D point x, y, z; 
        // compute signed distance of point above plane
		inline real PointPlaneDist(const Vec4 &pt) const {
			return _data[0] * pt._data[0] +
			   _data[1] * pt._data[1] +
			   _data[2] * pt._data[2] +
			   _data[3]; // assumes pt._data[3] = 1.0
		}

	Vec4 Cross3(const Vec4&);
	friend inline void Vec4FastCross3(Vec4 &c, const Vec4 &v1, const Vec4 &v2) {
			c.Set(  v1._data[1] * v2._data[2] - v1._data[2] * v2._data[1],
				  - v1._data[0] * v2._data[2] + v1._data[2] * v2._data[0],
					v1._data[0] * v2._data[1] - v1._data[1] * v2._data[0]);
		}


// *** Access functions: ***
	Vec4& Set(const Vec4 &v) { return (*this = v); }
	Vec4& Set(const float v[3]) {
			_data[0] = v[0]; _data[1] = v[1]; _data[2] = v[2];
			return *this;
		}
	inline Vec4& Set(real x, real y, real z, real w = 1) {
		_data[0] = x; _data[1] = y; _data[2] = z; _data[3] = w;
		return *this;
		}
	inline Vec4& Inc(real x, real y, real z, real w = 0) {
		_data[0] += x; _data[1] += y; _data[2] += z; _data[3] += w;
		return *this;
		}
	inline void Vec4FastGet3(real &x, real &y, real &z) const {
		x = _data[0]; y = _data[1]; z = _data[2]; }

        inline real operator [](int i) const { assert ((i >= 0) && (i <= 3)); return _data[i];}
        inline real &operator [](int i) { assert ((i >= 0) && (i <= 3)); return _data[i];}


	inline real x() const { return _data[0]; }
	inline real y() const { return _data[1]; }
	inline real z() const { return _data[2]; }
	inline real w() const { return _data[3]; }
	inline real r() const { return _data[0]; }
	inline real g() const { return _data[1]; }
	inline real b() const { return _data[2]; }
	inline real v(int i) const { return _data[i]; }
	inline void set_x(float x) { _data[0]=x; }
	inline void set_y(float y) { _data[1]=y; }
	inline void set_z(float z) { _data[2]=z; }
	inline void set_w(float w) { _data[3]=w; }
	inline void set_r(float r) { _data[0]=r; }
	inline void set_g(float g) { _data[1]=g; }
	inline void set_b(float b) { _data[2]=b; }

	inline void AugmentBound(real min_max[2][3]) {
			if		(_data[0] < min_max[0][0]) min_max[0][0] = _data[0];
			else if (_data[0] > min_max[1][0]) min_max[1][0] = _data[0];
			if		(_data[1] < min_max[0][1]) min_max[0][1] = _data[1];
			else if (_data[1] > min_max[1][1]) min_max[1][1] = _data[1];
			if		(_data[2] < min_max[0][2]) min_max[0][2] = _data[2];
			else if (_data[2] > min_max[1][2]) min_max[1][2] = _data[2];
		}

// *** If you're using a Vec4 to store a color: ***
	inline unsigned long RGBColor() {
		return int(_data[0] * 255) +
			(0x100 * int(_data[1] * 255)) +
			(0x10000 * int(_data[2] * 255));
		}

// *** For use with gl "v" function:
	inline float* ArrayForGL() {
			_float_array[0] = _data[0];
			_float_array[1] = _data[1];
			_float_array[2] = _data[2];
			return _float_array;
		}

	inline void Clamp(float min, float max) {
			if (_data[0] < min) _data[0] = min;
			if (_data[1] < min) _data[1] = min;
			if (_data[2] < min) _data[2] = min;
			if (_data[0] > max) _data[0] = max;
			if (_data[1] > max) _data[1] = max;
			if (_data[2] > max) _data[2] = max;
		}
};

#endif