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

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

#include "mycomplex.h"

//linear interpolation
double LERP(double a, double b, double t){
    return ((a) + (t)*((b)-(a)));
}
//-----------------------------------------------------

int MIN(int a, int b){
    return ((a)<=(b)?(a):(b));
}
//-----------------------------------------------------

double MIN(double a, double b){
    return ((a)<=(b)?(a):(b));
}
//-----------------------------------------------------

int MAX(int a, int b){
    return ((a)>=(b)?(a):(b));
}
//-----------------------------------------------------

double MAX(double a, double b){
    return ((a)>=(b)?(a):(b));
}
//-----------------------------------------------------

double MAX4(double a, double b, double c, double d){
    if(a>=b){
        if(a>=c){
            if(a>=d) return a; else return d;
        }else{
            if(c>=d) return c; else return d;
        }
    }else{
        if(b>=c){
            if(b>=d) return b; else return d;
        }else{
            if(c>=d) return c; else return d;
        }
    }
}
//-----------------------------------------------------



double expm1(double x){
    double u =exp(x);
    if(u == 1.) return x;
    if(u-1. == -1.) return -1;
    return (u-1.)*x/log(u);
}


double log1p(double x){
    double u = 1.+x;
    return log(u) - ((u-1.)-x)/u; // cancels errors with IEEE arithmetic
}


double __sinh(double x){
    // sinh(x) = (e^x - e^-x) / 2
    //         = (e^x - 1)(e^x + 1) / e^x / 2;
    //         = expm1(x) * (expm1(x)+2) / (expm1(x)+1) / 2
    double u = expm1(x);
    return .5 * u / (u+1) * (u+2); // ordered to avoid overflow when big
}

double __cosh(double x){
    // cosh(x) = (e^x + e^-x) / 2
    // I don't think there are any cancellation issues
    // (though probably coshm1, below, is more useful).
    double e_x = exp(x);
    return (e_x + 1./e_x) * .5;
}

// cosh(x)-1, accurate even when x is small.
double coshm1(double x){
    // cosh(x) - 1 = (e^x + e^-x) / 2 - 1
    //             = (e^x - 2 + e^-x) / 2
    //             = (e^2x - 2*e^x + 1) / e^x / 2
    //             = (e^x - 1)^2 / e^x / 2
    //             = expm1(x)^2 / (expm1(x)+1) / 2
    double u = expm1(x);
    return .5 * u / (u+1) * u; // ordered to avoid overflow when big
}

double __tanh(double x){
    // tanh(x) = sinh(x) / cosh(x)
    //         = (e^x - e^-x) / (e^x + e^-x)
    //         = (e^2x - 1) / (e^2x + 1)
    //         = expm1(2*x) / (expm1(2*x) + 2)
    // That works great but overflows prematurely, so do it
    // this way instead:
    // tanh(x) = (e^2x - 1) / (e^2x + 1)
    //         = (e^x - 1)(e^x + 1) / ((e^x - 1)(e^x + 1) + 2)
    //         = expm1(x)(expm1(x)+1) / (expm1(x)(expm1(x)+2) + 2)
    double u = expm1(x);
    return u / (u*(u+2.)+2.) * (u+2.); // ordered to avoid overflow when big
    // XXX oops, that doesn't avoid overflow
    // XXX since u*(u+2) can overflow... can we find another formulation?
}


double asinh(double x){
    // asinh(x) = log(x + sqrt(x^2 + 1))
    //          = log1p(x + sqrt(x^2 + 1) - 1)
    //          = log1p(x + (sqrt(x^2+1)-1)*(sqrt(x^2+1)+1)/(sqrt(x^2+1)+1))
    //          = log1p(x + (x^2+1 - 1)/(sqrt(x^2+1)+1))
    //          = log1p(x + x^2 / (sqrt(x^2+1)+1))
    //          = log1p(x * (1 + x / (sqrt(x^2+1)+1) ))
    return log1p(x * (1. + x / (sqrt(x*x+1.)+1.)));
}

double acosh(double x){
    // Only defined for x >= 1.

    // acosh(x) = log(x + sqrt(x^2 - 1))
    // Use the formula given by Kahan in
    // "Branch Cuts for Complex Elementary Functions",
    // as quoted by Cleve Moler in:
    // http://www.mathworks.com/company/newsletter/clevescorner/sum98cleve.shtml
    return 2 * log(sqrt((x+1)*.5) + sqrt((x-1)*.5));
}

double atanh(double x){
    // Only defined for x < 1.

    // atanh(x) = log((1+x)/(1-x))/2
    //          = log((1 - x + 2x) / (1-x)) / 2
    //          = log(1 + 2x/(1-x)) / 2
    //          = log1p(2x/(1-x)) / 2

    if(x>=1.) x=0.9999999;

    return .5 * log1p(2.*x/(1.-x));
}

// 1-cos(x), doesn't lose accuracy for small x
double cosf1(double x) {
    double sinHalfX = sin(.5*x);
    return 2.*sinHalfX*sinHalfX;
}

// acos(1-x), doesn't lose accuracy for small x
double acos1m(double x){
    return 2.*asin(sqrt(.5*x));
}

// sin(x)/x, but stable when small
double sin_over_x(double x){
    //
    // It's 1 - x^2/3! + x^4/5! - ...
    // so if |x| is so small that 1-x^2/6 is indistinguishable from 1,
    // then the result will be indistinguishable from 1 too.
    //
    if (1. - x*x*(1/6.) == 1.){
        //System.out.println("Ha! sin(x)/x("+x+") returning 1. instead of "+Math.sin(x)+"/"+x);
        return 1.;
    }
    else
        return sin(x)/x;
}

// asin(x)/x, but stable when small
double asin_over_x(double x){
    // x + 1/6*x^3 + 3/40*x^5 + 5/112*x^7 + 35/1152*x^9 + 63/2816*x^11 + 231/13312*x^13 + ...
    if (1. + x*x*(1/3.) == 1.)
        return 1.;
    else
        return asin(x)/x;
}

// (1-cos(x))/x, but stable when small
double cosf1_over_x(double x){
    //
    // It's x/2! - x^3/4! + x^5/6! - x^7/8! + ...
    // So if x is so small that x/2 - x^3/24
    // is indistinguishable from x/2,
    // then the result will be indistinguishable from .5*x.
    // This is pretty much the same as .5 - x^2/24
    // being indistinguishable from .5,
    // but maybe off by as much as 1 in the least-significant bit
    // (if x is just on the wrong side of a power of 2, I think)
    // so be conservative and test whether
    // .5 - x^2/12 is indistinguishable from .5 instead.
    // This is exactly the same as whether 1 - x^2/6
    // is indistingishable from 1.
    //
    if (1. - x*x*(1/6.) == 1.){
        //System.out.println("Ha! (1-cos(x))/x("+x+") returning "+(.5*x)+"instead of "+cosf1(x)+"/"+x+" = "+(cosf1(x) / x)+"");
        return .5*x;
    }
    else
        return cosf1(x) / x;
}

double myhypot(double x, double y){
    x = abs(x);
    y = abs(y);
    double min, max;
    if (x < y){
        min = x;
        max = y;
    }else{
        min = y;
        max = x;
    }
    if(min == 0.) return max;
    double min_max = min/max;
    return max * sqrt(1. + min_max * min_max);
}







//multiply a vector of complex numbers by a vector of real numbers
//(the result is just one complex number)
void vxm(complex& result, Vector<double>& v, Vector<complex>& m){
    if(!v.length) return;
    double re = 0.0, im = 0.0;
    for(int j = v.length-1; j >= 0; j--){
        re += v[j] * m[j][0];
        im += v[j] * m[j][1];
    }
    result[0] = re;
    result[1] = im;
}
//-----------------------------------------------------



void mxv(double result[4], double m[4][4], double v[4]){
    for (int i = 0; i < 4; ++i){
        double *m_i = m[i];
        double sum = 0.0;
        for (int j = 0; j < 4; ++j){
            sum += m_i[j] * v[j];
        }
        result[i] = sum;
    }
}
//-----------------------------------------------------

void mxm(double result[4][4], double m0[4][4], double m1[4][4]){
    for (int i = 0; i < 4; ++i){
        double *m0_i = m0[i];
        double *result_i = result[i];
        for (int j = 0; j < 4; ++j){
            double sum = 0.0;
            for (int k = 0; k < 4; ++k){
                sum += m0_i[k]*m1[k][j];
            }
            result_i[j] = sum;
        }
    }
}
//-----------------------------------------------------


//---------------------------------------------------------------------------
#endif