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

    Based largely on code from
    http://www.plunk.org/~hatch/HyperbolicApplet/
*/
//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop

#include "Isometry.h"



//actually implement some stuff


Complex Isometry::apply(Complex z){
    if (R < 0)
        z = conj(z);
    return (T*z + P) / (conj(P)*T*z + 1.0);
}
//------------------------------------------------------------


// Note order: (F1*F0)(x,y) = F1(F0(x,y))
// Attempted to optimize somewhat, but ... temporaries are created
// and destroyed.
Isometry Isometry::mul(Isometry F1, Isometry F0)
{
    if (F1.R < 0){
        F0.T = conj(F0.T);
        F0.P = conj(F0.P);
    }
    Complex denom = F1.T*F0.P*conj(F1.P) + 1.0;
    Isometry result;
    result.R = F1.R * F0.R;
    result.P = (F1.T*F0.P + F1.P) / denom;
    result.T = (F0.T*F1.T + F0.T*conj(F0.P)*F1.P) / denom;
    // renormalize, as recommended in the paper.
    // use the fact that sqrt(1+eps) is approximately 1+eps/2 for small eps.
    /* double invLenT = 1./hypot(result.T.x, result.T.y); */
    double invLenT = 1.0 / ((1+result.T.real()*result.T.real()+result.T.imag()*result.T.imag())*0.5);
    result.T = Complex(result.T.real()*invLenT, result.T.imag()*invLenT);
    return result;
}
//------------------------------------------------------------


Isometry Isometry::pow(Isometry F, int e){
    if (e > 0){
        Isometry temp = pow(F, e/2);
        if ((e % 2) == 0)
            return (temp*temp);
        else
            return ((F*temp)*temp);
    }
    else if (e < 0){
        return pow(F.inverse(), -e);
    }
    else{
        return Isometry(Complex(1.0,0.0),Complex(0.0,0.0),1);
    }
}
//------------------------------------------------------------

// not optimized, shouldn't be used in inner loops
Isometry Isometry::inverse(){
    Isometry F0(1,-P), F1(conj(T),0), F2(1,0,R);
    return F2 * F1 * F0;
}
//------------------------------------------------------------

Isometry Isometry::pureRotation(double angle){
    return Isometry(Complex(cos(angle),sin(angle)), Complex(0.0,0.0), 1);

}
//------------------------------------------------------------

// the simple isometry that takes 0,0 to x1,y1, and takes -x1,-y1 to 0,0
Isometry Isometry::pureTranslation(double x1, double y1){
    return Isometry(Complex(1.0,0.0), Complex(x1,y1), 1);
}
//------------------------------------------------------------

//
// Find the pure translation (i.e. moves the origin straight in some
// direction) that takes p0 to p1.
// This is a result of calculation on paper, and it can probably
// be simplified.
// It's ambiguous if both points are on the circumference of the circle.
// (Note, this is completely non-optimized and shouldn't be used
// in heavy computation.  Mostly it's good for interaction.)
//
Isometry Isometry::pureTranslation(Complex p0, Complex p1){
    Complex A = p1-p0;
    Complex B = p1*p0;
    double denom = 1.0 - (B.real()*B.real() + B.imag()*B.imag());
    return Isometry(Complex(1.0,0.0), Complex((A.real()*(1+B.real())+A.imag()*B.imag())/denom,
                                               (A.imag()*(1-B.real())+A.real()*B.imag())/denom),
                                                1);
}
//------------------------------------------------------------

// utility... not optimized XXX probably a much more direct way
Isometry Isometry::reflectionAcrossLine(Complex l0, Complex l1){
    Isometry takel0ToOrigin = Isometry::pureTranslation(-l0.real(),-l0.imag());
    Isometry takeOriginTol0 = Isometry::pureTranslation(l0.real(),l0.imag());
    Complex l1_ = takel0ToOrigin.apply(l1);
    double angle = atan2(l1_.imag(), l1_.real());
    Isometry rotatel1_ToXAxis = Isometry::pureRotation(-angle);
    Isometry rotateXAxisTol1_ = Isometry::pureRotation(angle);
    Isometry reflectAboutXAxis(Complex(1.0, 0.0), Complex(0.0, 0.0), -1);
    Isometry result = takel0ToOrigin;
    result = rotatel1_ToXAxis * result;
    result = reflectAboutXAxis * result;
    result = rotateXAxisTol1_ * result;
    result = takeOriginTol0 * result;
    return result;
}
//------------------------------------------------------------


// This isn't foolproof, but is a quick and dirty way
// to hash a tuple of doubles in the range [-1..1] (approximately).
// The algorithm is:
//    for each element
//        1. Add a small uncommon irrational constant > 2, to make sure
//           the result is > 1 and to minimize the chance that truncating
//           the low bits will be unstable.
//           We use the constant PI*E/4 =~ 2.1349336...
//        2. The result is a small double > 1.
//           Take the high 16 bits of the fractional part,
//           and throw the result into the hash.
//           (If we use too many bits, we raise the chance of
//           instability which can make us violate
//           the contract that objects for which equal() is true
//           should have the same hash code; I don't know what will
//           happen then.  But if we use too few
//           bits, we will get hash collisions, which aren't fatal
//           but which degrade performance).
int Isometry::hashCode(){
    int hash = 0;
    int i = 0;
    hash = _accumulateHash(hash, _hashCode1(T.real()), i++);
    hash = _accumulateHash(hash, _hashCode1(T.imag()), i++);
    hash = _accumulateHash(hash, _hashCode1(P.real()), i++);
    hash = _accumulateHash(hash, _hashCode1(P.imag()), i++);
    hash = _accumulateHash(hash, R, i++);
    return hash;
}
//-----------------------------------------------------------








Isometry operator*(Isometry F1, Isometry F0)

{
    if (F1.R < 0){
        F0.T = conj(F0.T);
        F0.P = conj(F0.P);
    }
    Complex denom = F1.T*F0.P*conj(F1.P) + 1.0;
    Isometry result;
    result.R = F1.R * F0.R;
    result.P = (F1.T*F0.P + F1.P) / denom;
    result.T = (F0.T*F1.T + F0.T*conj(F0.P)*F1.P) / denom;
    // renormalize, as recommended in the paper.
    // use the fact that sqrt(1+eps) is approximately 1+eps/2 for small eps.
    /* double invLenT = 1./hypot(result.T.x, result.T.y); */
    double invLenT = 1.0 / ((1+result.T.real()*result.T.real()+result.T.imag()*result.T.imag())*0.5);
    result.T = Complex(result.T.real()*invLenT, result.T.imag()*invLenT);
    return result;
}

//---------------------------------------------------------


bool operator==(Isometry& first, Isometry& other){

    // XXX clean up-- this is error-prone!
    double eps=1e-6;
    return (first.R == other.R
        && first.T.real()-other.T.real() <= eps
        && other.T.real()-first.T.real() <= eps
        && first.T.imag()-other.T.imag() <= eps
        && other.T.imag()-first.T.imag() <= eps
        && first.P.real()-other.P.real() <= eps
        && other.P.real()-first.P.real() <= eps
        && first.P.imag()-other.P.imag() <= eps
        && other.P.imag()-first.P.imag() <= eps);
}
//----------------------------------------------------------













//---------------------------------------------------------------------------
#pragma package(smart_init)