# Computing max t such that tanh(pi*sinh(t)/2) <1

I'm looking for a way to quickly determine that maximum value of t such that tanh(pi/2 sinh(t)) is strictly less than 1, for a given precision.

I have reproduced the code I am currently using to do this below:

#include <cmath>
#include <iostream>
#include <random>
#include <boost/math/special_functions/next.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <benchmark/benchmark.h>

using std::tanh;
using std::sinh;
using std::asinh;
using std::atanh;
using std::isnormal;
using boost::math::float_prior;
using boost::math::constants::half_pi;
using boost::math::constants::two_div_pi;
using boost::multiprecision::float128;
using boost::multiprecision::cpp_bin_float_50;
using boost::multiprecision::cpp_bin_float_100;

template<class Real>
Real g(Real t)
{
Real tmp = half_pi<Real>()*sinh(t);
return tanh(tmp);
}

template<class Real>
Real g_inv(Real x)
{
return asinh(two_div_pi<Real>()*atanh(x));
}

template<class Real>
Real find_max_t()
{
Real x = float_prior((Real) 1);
Real t_max = g_inv(x);
while(!isnormal(t_max))
{
x = float_prior(x);
t_max = g_inv(x);
}
// This occurs once on 100 digit arithmetic:
while(!(g(t_max) < (Real) 1))
{
x = float_prior(x);
t_max = g_inv(x);
}
return t_max;
}

template<class Real>
static void BM_max_t(benchmark::State& state)
{
Real t_max;
while(state.KeepRunning())
{
benchmark::DoNotOptimize(t_max = find_max_t<Real>());
}
std::cout << "t_max = " << t_max << std::endl;
}

BENCHMARK_TEMPLATE(BM_max_t, float);
BENCHMARK_TEMPLATE(BM_max_t, double);
BENCHMARK_TEMPLATE(BM_max_t, long double);
BENCHMARK_TEMPLATE(BM_max_t, float128);
BENCHMARK_TEMPLATE(BM_max_t, cpp_bin_float_50);
BENCHMARK_TEMPLATE(BM_max_t, cpp_bin_float_100);

BENCHMARK_MAIN();


The build sequence is:

INC = -I/path/to/boost
all: perf.x

perf.x: perf.o
g++ -o $@$< -lbenchmark -pthread -lquadmath

perf.o: perf.cpp
g++ --std=c++14 -fext-numeric-literals -Wfatal-errors -g -O3 $(INC) -c$< -o \$@


This computes the correct value of t_max in 153 ns (float), 214 ns (double), 246 ns (long double), 5280 ns (float128), 53us (cpp_bin_float_50), and 455 us (cpp_bin_float_100).

How can this code be made faster? For reference, at double precision, I have found 30% of the time is spent calling __expm1, 24% calling __ieee754_log_avx, and 15% calling __asinh.

• You do not handle domain_error and overflow_error (at least in this code). Can you assume they won't happen (or an undefined result is acceptable in those cases)? In general I'd not suggest this but you may well re-write float_prev() removing all special cases (you already have a call to isnormal()...) and error handling you do not want to have. It's a tiny overhead but if you're really squeezing the ns... – Adriano Repetti Oct 6 '17 at 8:09

I don't have any suggestions for making the code faster (other than pre-computing at compile-time, which I assume isn't an option for you), but I do have some general observations on the code style.

Instead of bringing definitions of the std and boost names at global scope, they can be isolated in the functions that use them (still with using, so ADL will continue to work the same as in the original code):

#include <cmath>
#include <boost/math/constants/constants.hpp>

template<class Real>
Real g(Real t)
{
using std::tanh;
using std::sinh;
using boost::math::constants::half_pi;
Real tmp = half_pi<Real>()*sinh(t);
return tanh(tmp);
}

template<class Real>
Real g_inv(Real x)
{
using std::asinh;
using std::atanh;
using boost::math::constants::two_div_pi;
return asinh(two_div_pi<Real>()*atanh(x));
}


Standard C++ provides nexttoward(), which can be used in place of boost::math::float_prior:

template<class Real>
Real find_max_t()
{
using std::isnormal;
using std::nexttoward;

static const Real ONE{1};

Real x = nexttoward(ONE, 0);
Real t_max = g_inv(x);
while (!isnormal(t_max))
{
x = float_prior(x);
t_max = g_inv(x);
}
// This occurs once on 100 digit arithmetic:
while (!(g(t_max) < ONE))
{
x = nexttoward(x, 0);
t_max = g_inv(x);
}
return t_max;
}


Unfortunately, the Boost numeric types don't provide any nexttoward(), so you may be stuck using boost::math::float_prior after all. :-(

The benchmark code can now be separable:

#include <iostream>
#include <benchmark/benchmark.h>

template<class Real>
static void BM_max_t(benchmark::State& state)
{
Real t_max;
while(state.KeepRunning())
{
benchmark::DoNotOptimize(t_max = find_max_t<Real>());
}
std::cout << "t_max = " << t_max << std::endl;
}

#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
BENCHMARK_TEMPLATE(BM_max_t, float);
BENCHMARK_TEMPLATE(BM_max_t, double);
BENCHMARK_TEMPLATE(BM_max_t, long double);
BENCHMARK_TEMPLATE(BM_max_t, boost::multiprecision::float128);
BENCHMARK_TEMPLATE(BM_max_t, boost::multiprecision::cpp_bin_float_50);
BENCHMARK_TEMPLATE(BM_max_t, boost::multiprecision::cpp_bin_float_100);

BENCHMARK_MAIN();