#ifndef __VECTOR4_H__
#define __VECTOR4_H__

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

#include "public.h"


namespace xrt {



// The Vector4 class is for 4-component spatial values.
//
// We GUARANTEE that the memory layout and size of a Vector4 is
// identical to a generic float[4], and thus it's presumed relatively
// safe to cast between a (float *) and a (Vector4 *).
//
// Semantics are as expected: construction from one or four floats or
// a float pointer (presumed to point to 4 floats), [] for accessing
// individual components, the usual arithmetic (+, -, *, /).  Note
// that math ops between two vectors are always component-by-component.
// 
// There are a handful of useful friends -- routines for normalizing,
// taking distance between two points, dot and cross products, etc.
//
// It's up to the user to keep track of whether they represent
// points, vectors, or normals, and to only use the operators where
// their semantics make sense.
//
// See also the Matrix4 class for other handy friends.

class PUBLIC Vector4 {

public:

    // Constructors
    Vector4 (void) { }
    explicit Vector4 (float x) {
#ifdef VECMAT_LOOPS
	for (int i = 0;  i < 4;  ++i)
            v[i] = x;
#else
        v[0] = x; v[1] = x; v[2] = x; v[3] = x;
#endif
    }
    Vector4 (float x, float y, float z, float w) {
	v[0] = x; v[1] = y; v[2] = z; v[3] = w;
    }
    explicit Vector4 (const float *xyzw) {
	assert (xyzw != NULL);
#ifdef VECMAT_LOOPS
	for (int i = 0;  i < 4;  ++i)
            v[i] = xyzw[i];
#else
	v[0] = xyzw[0]; v[1] = xyzw[1]; v[2] = xyzw[2]; v[3] = xyzw[3];
#endif
    }

    // Direct access to the data -- use with caution
    const float *data(void) const { return v; }

    // Simple assignment
    const Vector4& operator= (float x) {
#ifdef VECMAT_LOOPS
	for (int i = 0;  i < 4;  ++i)
            v[i] = x;
#else
	v[0] = x; v[1] = x; v[2] = x; v[3] = x;
#endif
	return *this;
    }
    // Use the default assignment operator (bitwise copy)


    // Component access
    float operator[] (int i) const {
	assert(i>=0 && i < 4);  // range check -- only in DEBUG mode
	return v[i];
    }
    float& operator[] (int i) {
	assert(i>=0 && i < 4);  // range check -- only in DEBUG mode
	return v[i];
    }


    // Comparisons between points, and with floats
    friend bool operator== (const Vector4 &a, const Vector4 &b) {
	return a[0] == b[0] && a[1] == b[1] && a[2] == b[2] && a[3] == b[3];
    }
    friend bool operator!= (const Vector4 &a, const Vector4 &b) {
	return a[0] != b[0] || a[1] != b[1] || a[2] != b[2] || a[3] != b[3];
    }
    friend bool operator== (const Vector4 &a, float b) {
	return a[0] == b && a[1] == b && a[2] == b && a[3] == b;
    }
    friend bool operator!= (const Vector4 &a, float b) {
	return a[0] != b || a[1] != b || a[2] != b || a[3] != b;
    }
    friend bool operator== (float b, const Vector4 &a) { return (a==b); }
    friend bool operator!= (float b, const Vector4 &a) { return (a!=b); }

    // less-than operator: this isn't mathematically meaningful,
    // but lets you use certain STL classes.
    friend bool operator< (const Vector4 &a, const Vector4 &b) {
	if (a[0] < b[0]) return true;
	else if (a[0] > b[0]) return false;
	if (a[1] < b[1]) return true;
	else if (a[1] > b[1]) return false;
	if (a[2] < b[2]) return true;
	else if (a[2] > b[2]) return false;
	if (a[3] < b[3]) return true;
	else return false;
    }

    // Stream output -- very handy for debugging
    friend std::ostream& operator<< (std::ostream& out, const Vector4& a) {
	return (out << '[' << a[0] << ' ' << a[1] << ' ' 
                    << a[2] << ' ' << a[3] << ']');
    }


    //////////
    // The obvious vector operations:
    // 

    // Addition
    friend Vector4 operator+ (const Vector4& a, const Vector4& b) {
#ifdef VECMAT_LOOPS
        Vector4 r;
        for (int i = 0;  i < 4;  ++i)
            r[i] = a[i] + b[i];
        return r;
#else
	return Vector4(a.v[0]+b.v[0], a.v[1]+b.v[1],
                       a.v[2]+b.v[2], a.v[3]+b.v[3]);
#endif
    }
    const Vector4& operator+= (const Vector4& b) {
#ifdef VECMAT_LOOPS
	for (int i = 0;  i < 4;  ++i)
            v[i] += b.v[i];
#else
	v[0] += b.v[0];  v[1] += b.v[1];  v[2] += b.v[2];  v[3] += b.v[3];
#endif
	return *this;
    }

    // Subtraction
    friend Vector4 operator- (const Vector4& a, const Vector4& b) {
#ifdef VECMAT_LOOPS
        Vector4 r;
        for (int i = 0;  i < 4;  ++i)
            r[i] = a[i] - b[i];
        return r;
#else
	return Vector4(a.v[0]-b.v[0], a.v[1]-b.v[1],
                       a.v[2]-b.v[2], a.v[3]-b.v[3]);
#endif
    }
    const Vector4& operator-= (const Vector4& b) {
#ifdef VECMAT_LOOPS
	for (int i = 0;  i < 4;  ++i)
            v[i] -= b.v[i];
#else
	v[0] -= b.v[0];  v[1] -= b.v[1];  v[2] -= b.v[2];  v[3] -= b.v[3];
#endif
	return *this;
    }

    // Negation
    friend Vector4 operator- (const Vector4& b) {
	return Vector4(-b.v[0], -b.v[1], -b.v[2], -b.v[3]);
    }

    // Scalar multiplication
    friend Vector4 operator* (const Vector4& a, float b) {
	return Vector4(a.v[0]*b, a.v[1]*b, a.v[2]*b, a.v[3]*b);
    }
    friend Vector4 operator* (float b, const Vector4& a) {
	return Vector4(a.v[0]*b, a.v[1]*b, a.v[2]*b, a.v[3]*b);
    }
    const Vector4& operator*= (float b) {
#ifdef VECMAT_LOOPS
	for (int i = 0;  i < 4;  ++i)
            v[i] *= b;
#else
	v[0] *= b;  v[1] *= b;  v[2] *= b;  v[3] *= b;
#endif
	return *this;
    }

    // Component-by-component multiplication
    friend Vector4 operator* (const Vector4& a, const Vector4& b) {
#ifdef VECMAT_LOOPS
        Vector4 r;
        for (int i = 0;  i < 4;  ++i)
            r[i] = a[i] * b[i];
        return r;
#else
	return Vector4(a.v[0]*b.v[0], a.v[1]*b.v[1],
                       a.v[2]*b.v[2], a.v[3]*b.v[3]);
#endif
    }
    const Vector4& operator*= (const Vector4& b) {
#ifdef VECMAT_LOOPS
	for (int i = 0;  i < 4;  ++i)
            v[i] *= b.v[i];
#else
	v[0] *= b.v[0];  v[1] *= b.v[1];  v[2] *= b.v[2];  v[3] *= b.v[3];
#endif
	return *this;
    }

    // Scalar division -- watch out for divide-by-zero
    friend Vector4 operator/ (const Vector4& a, float b) {
	if (b == 0.0f)
	    return Vector4(0.0f);
	else {
            float binv = 1.0f/b;
	    return Vector4(a.v[0]*binv, a.v[1]*binv, a.v[2]*binv, a.v[3]*binv);
        }
    }
    friend Vector4 operator/ (float a, const Vector4& b) {
        Vector4 dst;
#ifdef VECMAT_LOOPS
	for (int i = 0;  i < 4;  ++i)
            if (b[i] == 0.0f) dst[i] = 0.0f;
            else dst[i] = a / b[i];
#else
        if (b[0] == 0.0f) dst[0] = 0;
        else dst[0] = a / b[0];
        if (b[1] == 0.0f) dst[1] = 0;
        else dst[1] = a / b[1];
        if (b[2] == 0.0f) dst[2] = 0;
        else dst[2] = a / b[2];
        if (b[3] == 0.0f) dst[3] = 0;
        else dst[3] = a / b[3];
#endif
        return dst;
    }
    const Vector4& operator/= (float b) {
	if (b == 0.0f)
	    *this = 0.0f;
	else {
            float binv = 1.0f/b;
#ifdef VECMAT_LOOPS
            for (int i = 0;  i < 4;  ++i)
                v[i] *= binv;
#else
	    v[0] *= binv;  v[1] *= binv;  v[2] *= binv;  v[3] *= binv;
#endif
	}
	return *this;
    }

    // Component-by-component division -- watch out for divide-by-zero
    friend Vector4 operator/ (const Vector4& a, const Vector4& b) {
	return Vector4((b.v[0] == 0.0f) ? 0.0f : a.v[0]/b.v[0],
                       (b.v[1] == 0.0f) ? 0.0f : a.v[1]/b.v[1],
                       (b.v[2] == 0.0f) ? 0.0f : a.v[2]/b.v[2],
                       (b.v[3] == 0.0f) ? 0.0f : a.v[3]/b.v[3]);
    }
    const Vector4& operator/= (const Vector4& b) {
#ifdef VECMAT_LOOPS
	for (int i = 0;  i < 4;  ++i)
            v[i] = (b.v[i] == 0.0f) ? 0.0f : v[i]/b.v[i];
#else
	v[0] = (b.v[0] == 0.0f) ? 0.0f : v[0]/b.v[0];
	v[1] = (b.v[1] == 0.0f) ? 0.0f : v[1]/b.v[1];
	v[2] = (b.v[2] == 0.0f) ? 0.0f : v[2]/b.v[2];
	v[3] = (b.v[3] == 0.0f) ? 0.0f : v[3]/b.v[3];
#endif
	return *this;
    }

    // Return the Euclidean length
    friend float length (const Vector4& a) {
	return sqrtf (a[0]*a[0] + a[1]*a[1] + a[2]*a[2] + a[3]*a[3]);
    }
    friend float squaredLength (const Vector4& a) {
	return (a[0]*a[0] + a[1]*a[1] + a[2]*a[2] + a[3]*a[3]);
    }

    // Member normalize alters the Vector4, returns original length
    float normalize (void) {
	float len = length(*this);
	if (len != 0.0f) {
            float s = 1.0f/len;
	    v[0] *= s; v[1] *= s; v[2] *= s;  v[3] *= s;
	}
	return len;
    }

    // Return the Euclidean distance between two points
    // NOTE: must be 'dist' since 'distance' is something else in STL
    friend float distSquared (const Vector4& a, const Vector4 &b) {
#ifdef VECMAT_LOOPS
        float x[4];
        for (int i = 0;  i < 4;  ++i) {
            x[i] = a[i] - b[i];
            x[i] *= x[i];
        }
        return x[0] + x[1] + x[2] + x[3];
#else
	float x = a[0] - b[0];
	float y = a[1] - b[1];
	float z = a[2] - b[2];
	float w = a[3] - b[3];
	return (x*x + y*y + z*z + w*w);
#endif
    }
    friend float dist (const Vector4& a, const Vector4 &b) {
	return sqrtf (distSquared (a, b));
    }

public:
    union {
		float v[4];
		struct { float x, y, z, w; };
	};

};


// Dot product
inline float dot (const Vector4& a, const Vector4& b) {
return a.v[0]*b.v[0] + a.v[1]*b.v[1] + a.v[2]*b.v[2] + a.v[3]*b.v[3];
}

// Functional normalize returns a normalized Vector4, leaves alone original
inline Vector4 normalize (const Vector4& P) {
float len = length(P);
return (len == 0.0f) ? Vector4(0.0f)
	                    : Vector4(P[0]/len, P[1]/len, P[2]/len, P[3]/len);
}


};  /* end namespace xrt */


#endif  /* !defined(__VECTOR4_H__) */