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

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

#ifndef _complex_h
#define _complex_h

#include <math.h>
#include "myvector.h"



//------- SLOWER (BUT MORE POLITICALLY CORRECT) COMPLEX CLASS ----------
//
// (could easily transition to >2 dimensions)
//
// uncomment the following stuff and comment out the rest


/*
template <class T>
class Complex : public Vector<T> {
public:

    // Constructors...
    Complex(): Vector<T>(2) {
     // data[0]=0.; data[1]=0.;
     }

    Complex(const T& x, const T& y): Vector<T>(2) {
        data[0] = x;
        data[1] = y;
    }

    inline
    Complex<T>& operator =(Complex<T>& a){
        data[0] = a.data[0];
        data[1] = a.data[1];
        return *this;
    }

    inline
    Complex<T>& operator =(const Complex<T>& a){
        data[0] = a.data[0];
        data[1] = a.data[1];
        return *this;
    }


    inline
    Complex& operator*=(const Complex& b){
        T temp = data[0]*b.data[0] - data[1]*b.data[1];
        data[1] = data[0]*b.data[1] + data[1]*b.data[0];
        data[0] = temp;
        return *this;
    }

    inline
    Complex& operator*=(const T& b){
        data[0] *= b;
        data[1] *= b;
        return *this;
    }

};



inline
template <class T>
Complex<T> operator+(const Complex<T>& a, const Complex<T>& b){
    return Complex<T>(a.data[0]+b.data[0], a.data[1]+b.data[1]);
}

inline
template <class T>
Complex<T> operator+(const Complex<T>& a, const T& b){
    return Complex<T>(a.data[0]+b, a.data[1]);
}

inline
template <class T>
Complex<T> operator+(const T& a, const Complex<T>& b){
    return Complex<T>(b.data[0]+a, b.data[1]);
}

inline
template <class T>
Complex<T> operator-(const Complex<T>& a, const Complex<T>& b){
    return Complex<T>(a.data[0]-b.data[0], a.data[1]-b.data[1]);
}

inline
template <class T>
Complex<T> operator-(Complex<T>& a){
    return Complex<T>(-a.data[0], -a.data[1]);
}

//conjugate operator!
inline
template <class T>
Complex<T> operator~(const Complex<T>& a){
    return Complex<T>(a.data[0], -a.data[1]);
}

inline
template <class T>
Complex<T> operator*(const Complex<T> a, const Complex<T> b){
    Complex<T> result;
    result.data[0] = a.data[0]*b.data[0] - a.data[1]*b.data[1];
    result.data[1] = a.data[0]*b.data[1] + a.data[1]*b.data[0];
    return result;
}

inline
template <class T>
Complex<T> operator*(const Complex<T>& a, const T& b){
    return Complex<T>(a.data[0]*b, a.data[1]*b);
}

inline
template <class T>
Complex<T> operator*(const T& b, const Complex<T>& a){
    return Complex<T>(a.data[0]*b, a.data[1]*b);
}

inline
template <class T>
Complex<T> operator/(const Complex<T>& a, const Complex<T>& b){
    T invLenSqrd = (b.data[0]*b.data[0] + b.data[1]*b.data[1]);
    return Complex<T> ((a.data[0]*b.data[0] + a.data[1]*b.data[1]) / invLenSqrd,
                (a.data[1]*b.data[0] - a.data[0]*b.data[1]) / invLenSqrd);
}

// Change from Klein model to Poincare model, in place.
inline
template <class T>
void k2p(Vector<T>& p){
    T temp=0;
    for(int i=p.length-1; i>=0; i--)
        temp += p.data[i]*p.data[i];

    //slow?
    if(temp > 1.) return;

    T scale = 1./(1. + sqrt(1. - temp));
    try{
        for(int i=p.length-1; i>=0; i--)
            p.data[i] *= scale;
    }catch(...){
    // could be NaN exception
    }
}

// Change from Poincare disk model to Klein disk model, in place.
inline
template <class T>
void p2k(Vector<T>& p){
    T temp=0;
    for(int i=p.length-1; i>=0; i--)
        temp += p.data[i]*p.data[i];
    T scale = 2. / (1. + temp);
    for(int i=p.length-1; i>=0; i--)
        p.data[i] *= scale;
}

*/





// ------- ALTERNATE (FASTER) COMPLEX CLASS ------------
//
// (but stuck in 2 dimentions)
//
//



template <class T>
class Complex {
public:

    T x, y;

    // Constructors...
    Complex(){
     // data[0]=0.; data[1]=0.;
     }

    Complex(const T& _x, const T& _y){
        x = _x;
        y = _y;
    }

    inline
    Complex<T>& operator =(Complex<T>& a){
        x = a.x;
        y = a.y;
        return *this;
    }

    inline
    Complex<T>& operator =(const Complex<T>& a){
        x = a.x;
        y = a.y;
        return *this;
    }

    inline
    Complex& operator*=(const Complex& b){
        T temp = x*b.x - y*b.y;
        y = x*b.y + y*b.x;
        x = temp;
        return *this;
    }

    inline
    T& operator[](int b){
        if(b == 0) return x;
        else return y;
    }

    inline
    const T& operator[](int b) const{
        if(b == 0) return x;
        else return y;
    }

    inline
    Complex& operator+(const Complex& b){
        x += b.x;
        y += b.y;
        return *this;
    }

    inline
    Complex& operator-=(const Complex& b){
        x -= b.x;
        y -= b.y;
        return *this;
    }

    inline
    Complex& operator*=(const T& b){
        x *= b;
        y *= b;
        return *this;
    }

};

inline
template <class T>
Complex<T> operator+(const Complex<T>& a, const Complex<T>& b){
    return Complex<T>(a.x+b.x, a.y+b.y);
}

inline
template <class T>
Complex<T> operator+(const Complex<T>& a, const T& b){
    return Complex<T>(a.x+b, a.y);
}

inline
template <class T>
Complex<T> operator+(const T& a, const Complex<T>& b){
    return Complex<T>(b.x+a, b.y);
}

inline
template <class T>
Complex<T> operator-(const Complex<T>& a, const Complex<T>& b){
    return Complex<T>(a.x-b.x, a.y-b.y);
}

inline
template <class T>
Complex<T> operator-(Complex<T>& a){
    return Complex<T>(-a.x, -a.y);
}

//conjugate operator!
inline
template <class T>
Complex<T> operator~(const Complex<T>& a){
    return Complex<T>(a.x, -a.y);
}

inline
template <class T>
Complex<T> operator*(const Complex<T> a, const Complex<T> b){
    Complex<T> result;
    result.x = a.x*b.x - a.y*b.y;
    result.y = a.x*b.y + a.y*b.x;
    return result;
}

inline
template <class T>
Complex<T> operator*(const Complex<T>& a, const T& b){
    return Complex<T>(a.x*b, a.y*b);
}

inline
template <class T>
Complex<T> operator*(const T& b, const Complex<T>& a){
    return Complex<T>(a.x*b, a.y*b);
}

inline
template <class T>
Complex<T> operator/(const Complex<T>& a, const Complex<T>& b){
    T invLenSqrd = (b.x*b.x + b.y*b.y);
    return Complex<T> ((a.x*b.x + a.y*b.y) / invLenSqrd,
                (a.y*b.x - a.x*b.y) / invLenSqrd);
}

// Change from Klein model to Poincare model, in place.
inline
template <class T>
void k2p(Complex<T>& p){
    T temp=0;
    temp += p.x*p.x;
    temp += p.y*p.y;

    //slow?
    if(temp > 1.) return;

    T scale = 1./(1. + sqrt(1. - temp));
    try{
        p.x *= scale;
        p.y *= scale;
    }catch(...){
    // could be NaN exception
    }
}

// Change from Poincare disk model to Klein disk model, in place.
inline
template <class T>
void p2k(Complex<T>& p){
    T temp=0;
    temp += p.x*p.x;
    temp += p.y*p.y;
    T scale = 2. / (1. + temp);
    p.x *= scale;
    p.y *= scale;
}

inline
template <class T>

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

    T result = 0.;
    T diff = v1.x-v0.x;
    result += diff*diff;
    diff = v1.y-v0.y;
    result += diff*diff;
    return result;
}


void vxm(Complex<double>& result, Complex<double>& v, Vector<Complex<double> >& m){
    double re = 0.0, im = 0.0;
    re += v[0] * m[0][0];
    im += v[0] * m[0][1];
    re += v[1] * m[1][0];
    im += v[1] * m[1][1];
    result[0] = re;
    result[1] = im;
}
//-----------------------------------------------------

inline
template <class T>
T vxv2(const Complex<T>& v, const Complex<T>& w){
    return v.x*w.y - v.y*w.x;
}

inline
template <class T>

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

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


inline

template <class T>
T dot(const Complex<T>& v0, const Complex<T>& v1){
    T result=0;
    result += v0.x*v1.x;
    result += v0.y*v1.y;
    return result;
}

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

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

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





// ------- END ALTERNATE (FASTER) COMPLEX CLASS ------------


















#define complex     Complex<double>


#define PI          3.141592653589793238
#define EULER_E     2.718281828459045235


extern complex zero;
extern complex one;

#endif