/*
    Hyperbolic Tessellations
    Copyright (C) 2007 Dmitry Brant
    http://dmitrybrant.com

    Based largely on code from
    http://www.plunk.org/~hatch/HyperbolicApplet/
*/
#ifndef _myvector_h
#define _myvector_h


template <class T>
class Vector {
public:
    T* data;
    int length;

    Vector():length(0), data(NULL) { }
    ~Vector(){ if(data){ delete [] data; data = NULL; } }

    Vector(int _length):length(_length) {
        data = new T[_length];
    }
    Vector(T* _data, int _length){
        length = _length;
        data = new T[length];
        for(int i=0; i<length; i++)
            data[i] = _data[i];
    }
    Vector(const Vector& a){
        length = a.length;
        data = new T[length];
        for(int i=0; i<length; i++)
            data[i] = a.data[i];
    }

    inline
    Vector<T>& operator=(const Vector<T>& a){
        if(a.length != length){
            if(data){ delete [] data; data = NULL; }
            length = a.length;
            data = new T[length];
        }
        for(int i=0; i<length; i++)
            data[i] = a.data[i];
        return *this;
    }


    inline
    void Resize(int newLength){
        if(newLength == length) return;
        if(data){ delete [] data; data = NULL; }
        data = new T[newLength];
        length = newLength;
    }

    inline
    T& operator[](int index){ return data[index]; }

    inline
    const T& operator[](int index) const{ return data[index]; }

    inline
    Vector& operator*=(const T& t){
        for(int i=length-1; i>=0; i--) data[i] *= t;
        return *this;
    }

    inline
    Vector& operator/=(const T& t){
        for(int i=length-1; i>=0; i--) data[i] *= t;
        return *this;
    }

    inline
    Vector& operator+=(const Vector<T>& a){
        for(int i=length-1; i>=0; i--) data[i] += a.data[i];
        return *this;
    }

    inline
    Vector& operator-=(const Vector<T>& a){
        for(int i=length-1; i>=0; i--) data[i] -= a.data[i];
        return *this;
    }

    /*
    String toString(){
        String s;
        for(int i=length-1; i>=0; i--) s += String(data[i]) + " ";
        return s;
    }
    */

};



//template vector functions

//dot product of two vectors (any length, as long as they're equal)
inline
template <class T>
T dot(const Vector<T>& v0, const Vector<T>& v1){
    T result=0;
    for(int i=v0.length-1; i>=0; i--)
        result += v0[i]*v1[i];
    return result;
}

// perp dot
// only makes sense for length 2
inline
template <class T>
Vector<T> xv2(const Vector<T>& v){
    return complex(-v[1],v[0]);
}

// determinant of matrix having the two vectors as rows
// or cross product in two dimensions
// only makes sense for length 2
inline
template <class T>
T vxv2(const Vector<T>& v, const Vector<T>& w){
    return v[0]*w[1] - v[1]*w[0];
}

//cross product in three dimensions
//only makes sense for length 3
inline
template <class T>
Vector<T> vxv3(const Vector<T>& v1, const Vector<T>& v2){
    Vector<T> result(3);
    result[0] = v1[1]*v2[2] - v1[2]*v2[1];
    result[1] = v1[2]*v2[0] - v1[0]*v2[2];
    result[2] = v1[0]*v2[1] - v1[1]*v2[0];
    return result;
}


// determinant of matrix having the three vectors as rows

// only makes sense for length 3

inline
template <class T>
T vxvxv3(const Vector<T>& v0, const Vector<T>& v1, const Vector<T>& v2){
    return v0[0] * (v1[1]*v2[2] - v1[2]*v2[1])
         + v0[1] * (v1[2]*v2[0] - v1[0]*v2[2])
         + v0[2] * (v1[0]*v2[1] - v1[1]*v2[0]);
}



// ||b-a||^2 ||c-a||^2 - (b-a)dot(c-a)
inline
template <class T>
T sqrdTwiceTriangleArea(const Vector<T>& a, const Vector<T>& b, const Vector<T>& c){
    T ab_dot_ab = 0.;
    T ab_dot_ac = 0.;
    T ac_dot_ac = 0.;
    for(int i=a.length-1; i>=0; i--){
        T ai = a[i];
        T ab = b[i] - ai;
        T ac = c[i] - ai;
        ab_dot_ab += ab*ab;
        ab_dot_ac += ab*ac;
        ac_dot_ac += ac*ac;
    }
    return ab_dot_ab*ac_dot_ac - (ab_dot_ac*ab_dot_ac);
} // sqrdTwiceTriangleArea


inline
template <class T>
T normsqrd(const Vector<T>& v){
    return dot(v,v);
}

inline
template <class T>
T norm(const Vector<T>& v){
    return sqrt(dot(v,v));
}

inline
template <class T>
void normalize(Vector<T>& result, const Vector<T>& v){
    T len = norm(v);
    result = v;
    result *= (len==0. ? 1. : 1./len);
}


inline

template <class T>

T distsqrd(const Vector<T>& v0, const Vector<T>& v1){

    T result = 0.;
    for(int i=v0.length-1; i>=0; i--){
        T diff = v1[i]-v0[i];
        result += diff*diff;
    }
    return result;
}


inline

template <class T>

T dist(const Vector<T>& v0, const Vector<T>& v1){

    return sqrt(distsqrd(v0, v1));
}


inline
template <class T>
T angleBetweenUnitVectors(const Vector<T>& u0, const Vector<T>& u1){
    // XXX can we cleverly avoid the square root?
    if (dot(u0, u1) > 0)
        return 2*asin(.5*dist(u0,u1));
    else{
        Vector<T> minusU1 = u1;
        minusU1 *= -1.;
        return PI - 2*asin(.5*dist(u0,minusU1));
    }
}





#define DoubleVector        Vector<double>
#define IntVector           Vector<int>


#endif