sinhatanh2

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <fenv.h>
#include <typeinfo>
#include <mpreal.h>
#include <iomanip>
#include <immintrin.h>

#define PREC 50
#define ULP_THRESH 1

#define debug(msg, x) std::cout << std::scientific << std::setprecision(PREC) << msg << " = " << x << std::endl;

template<class R>
R sinhatanh(R x)
{
    if (typeid(R) == typeid(float))
        return sinhf(atanhf(x));
    else if (typeid(R) == typeid(double))
        return sinh(atanh(x));
    else if (typeid(R) == typeid(long double))
        return sinhl(atanhl(x));

    __builtin_unreachable();
}

template<class R>
R sinhatanh2_dijkstra(R x)
{
    if (typeid(R) == typeid(float))
    {
        float u = 1 - fabsf(x);
        return x/sqrtf(2*u - u*u);
    }
    else if (typeid(R) == typeid(double))
    {
        double u = 1 - fabs(x);
        return x/sqrt(2*u - u*u);
    }
    else if (typeid(R) == typeid(long double))
    {
        long double u = 1 - fabsl(x);
        return x/sqrtl(2*u - u*u);
    }
    __builtin_unreachable();
}

template<class R>
R sinhatanh2_fma(R x)
{
    if (typeid(R) == typeid(float))
    {
        __m128 mx  = _mm_set_ps(x  , x  ,  x ,  x );
        __m128 one = _mm_set_ps(1.f, 1.f, 1.f, 1.f);
        __m128 sq  = _mm_fnmadd_ps(mx, mx, one);
        float r = * (float*) &sq;
        return x/sqrtf(r);
    }
    else if (typeid(R) == typeid(double))
    {
        __m128d mx  = _mm_set_pd(x , x);
        __m128d one = _mm_set_pd(1., 1.);
        __m128d sq  = _mm_fnmadd_pd(mx, mx, one);
        double r = * (double*) &sq;
        return x/sqrt(r);
    }
    else if (typeid(R) == typeid(long double))
        return x/sqrtl(1 - x*x);

    __builtin_unreachable();
}

template<class R>
R coshatanh(R x)
{
    if (typeid(R) == typeid(float))
        return coshf(atanhf(x));
    else if (typeid(R) == typeid(double))
        return cosh(atanh(x));
    else if (typeid(R) == typeid(long double))
        return coshl(atanhl(x));

    __builtin_unreachable();
}

template<class R>
R coshatanh2_dijkstra(R x)
{
    R ret;
    if (typeid(R) == typeid(float))
    {
        double u = 1 - fabs(x);
        double sqrtr = sqrtf(2*u - u*u);
        ret =  1/sqrtr;
    }
    else if (typeid(R) == typeid(double))
    {
        double u = 1 - fabs(x);
        ret =  1/sqrt(2*u - u*u);
    }
    else if (typeid(R) == typeid(long double))
    {
        long double u = 1 - fabsl(x);
        ret = 1/sqrtl(2*u - u*u);
    }
    return ret;
}

template<class R>
R coshatanh2_fma(R x)
{
    if (typeid(R) == typeid(float))
    {
        __m128 mx  = _mm_set_ps(x  , x  ,  x ,  x );
        __m128 one = _mm_set_ps(1.f, 1.f, 1.f, 1.f);
        __m128 sq  = _mm_fnmadd_ps(mx, mx, one);
        float r = * (float*) &sq;
        return 1/sqrtf(r);
    }
    else if (typeid(R) == typeid(double))
    {
        __m128d mx  = _mm_set_pd(x , x);
        __m128d one = _mm_set_pd(1., 1.);
        __m128d sq  = _mm_fnmadd_pd(mx, mx, one);
        double r = * (double*) &sq;
        return 1/sqrt(r);
    }
    else if (typeid(R) == typeid(long double))
        return 1/sqrtl(1 - x*x);

    __builtin_unreachable();
}


template<class R>
R nextafter_(R x, R y)
{
    if (typeid(R) == typeid(float))
        return nextafterf(x, y);
    else if (typeid(R) == typeid(double))
        return nextafter(x, y);
    else if (typeid(R) == typeid(long double))
        return nextafterl(x, y);

    __builtin_unreachable();
}

//not efficient, but OK for our needs.
template<class R>
int ulp_diff(R x, R y)
{
    int ulp = 0;

    while (x != y)
    {
        x = nextafter_(x, y);
        ulp++;
    }
    return ulp;
}

template<class R>
mpfr::mpreal mpfr_sinhatanh(R x)
{
    mpfr::mpreal mpfr_x = x;
    return sinh(atanh(mpfr_x));
}

template<class R>
mpfr::mpreal mpfr_coshatanh(R x)
{
    mpfr::mpreal mpfr_x = x;
    return cosh(atanh(mpfr_x));
}

#define err(x) fabs(x - ms)

template<class R>
void sinh_test()
{
    R xi = 1, s1, s2, sp, mone=-1;
    mpfr::mpreal ms;

    do
    {
        xi = nextafter_(xi, mone);
        s1 = sinhatanh(xi);
        s2 = sinhatanh2_dijkstra(xi);
        ms = mpfr_sinhatanh(xi);
        sp = (R) ms;
    } while(isnan(s1-s2) || isinf(s2-s1) || ulp_diff(s1, sp) + ULP_THRESH >= ulp_diff(s2, sp));

    debug("with input  : ", xi);
    debug("sinh: before: ", s1);
    debug("sinh: after : ", s2);
    debug("sinh: mpfr  : ", ms);
    debug("ulp err befr: ", ulp_diff(s1, sp));
    debug("ulp err aftr: ", ulp_diff(s2, sp));
}
#undef err


#define err(x) fabs(x - mc)
template<class R>
void cosh_test()
{
    R xi = 1, c1, c2, cp, mone=-1;
    mpfr::mpreal mc;

    do
    {
        xi = nextafter_(xi, mone);
        c1 = coshatanh(xi);
        c2 = coshatanh2_dijkstra(xi);
        mc = mpfr_coshatanh(xi);
        cp = (R) mc;

    } while (isnan(c1-c2) || isinf(c2-c1) || ulp_diff(c1, cp) + ULP_THRESH >= ulp_diff(c2, cp));

    debug("with input  : ", xi);
    debug("cosh: before: ", c1);
    debug("cosh: after : ", c2);
    debug("cosh: mpfr  : ", mc);
    debug("ulp err befr: ", ulp_diff(c1, cp));
    debug("ulp err aftr: ", ulp_diff(c2, cp));
}
#undef err

int main()
{
    mpfr_set_default_prec(512);

    std::cout << "With floats: " << std::endl;
    sinh_test<float>();
//    cosh_test<float>();
    std::cout << std::endl;

    std::cout << "With doubles: " << std::endl;
    sinh_test<double>();
//    cosh_test<double>();
    std::cout << std::endl;

    std::cout << "With long doubles: " << std::endl;
    sinh_test<long double>();
//    cosh_test<long double>();
    std::cout << std::endl;

    return 0;
}