Add arctanh implementation

README.md

@@ -9,7 +9,7 @@ Currently supported:
 - exp/exp2/exp10/expm1
 - log/log2/log10/log1p
 - sinh/cosh/tanh
-- arcsinh/arccosh/
+- arcsinh/arccosh/arctanh
 - sigmoid
 - Wright-Omega function
 - Dilogarithm function

include/math_approx/src/inverse_hyperbolic_trig_approx.hpp

@@ -89,4 +89,11 @@ constexpr T acosh (T x)
     const auto z1 = z0 + sqrt (z0 + z0 + z0 * z0);
     return log1p<order> (z1);
 }
+
+template <int order, typename T>
+constexpr T atanh (T x)
+{
+    using S = scalar_of_t<T>;
+    return (S) 0.5 * log<order> (((S) 1 + x) / ((S) 1 - x));
+}
 }

test/src/inverse_hyperbolic_trig_approx_test.cpp

@@ -84,3 +84,35 @@ TEMPLATE_TEST_CASE ("Acosh Approx Test", "", float, double)
                      6.5e-4f);
     }
 }
+
+TEMPLATE_TEST_CASE ("Atanh Approx Test", "", float, double)
+{
+    const auto all_floats = test_helpers::all_32_bit_floats<TestType> (-0.9999f, 0.9999f, 1.0e-2f);
+    const auto y_exact = test_helpers::compute_all<TestType> (all_floats, [] (auto x)
+                                                    { return std::atanh (x); });
+
+    SECTION ("6th-Order")
+    {
+        test_approx<TestType> (all_floats, y_exact, [] (auto x)
+                     { return math_approx::atanh<6> (x); },
+                     2.5e-6f);
+    }
+    SECTION ("5th-Order")
+    {
+        test_approx<TestType> (all_floats, y_exact, [] (auto x)
+                     { return math_approx::atanh<5> (x); },
+                     6.5e-6f);
+    }
+    SECTION ("4th-Order")
+    {
+        test_approx<TestType> (all_floats, y_exact, [] (auto x)
+                     { return math_approx::atanh<4> (x); },
+                     4.5e-5f);
+    }
+    SECTION ("3th-Order")
+    {
+        test_approx<TestType> (all_floats, y_exact, [] (auto x)
+                     { return math_approx::atanh<3> (x); },
+                     3.5e-4f);
+    }
+}

tools/bench/inverse_hyperbolic_trig_bench.cpp

@@ -20,6 +20,15 @@ const auto data_acosh = []
     return x;
 }();
 
+const auto data_atanh = []
+{
+    std::vector<float> x;
+    x.resize (N, 0.0f);
+    for (size_t i = 0; i < N; ++i)
+        x[i] = -1.0f + 2.0f * (float) i / (float) N;
+    return x;
+}();
+
 #define INV_HTRIG_BENCH(name, func, data) \
 void name (benchmark::State& state) \
 { \
@@ -47,6 +56,12 @@ INV_HTRIG_BENCH (acosh_approx5, math_approx::acosh<5>, data_acosh)
 INV_HTRIG_BENCH (acosh_approx4, math_approx::acosh<4>, data_acosh)
 INV_HTRIG_BENCH (acosh_approx3, math_approx::acosh<3>, data_acosh)
 
+INV_HTRIG_BENCH (atanh_std, std::atanh, data_atanh)
+INV_HTRIG_BENCH (atanh_approx6, math_approx::atanh<6>, data_atanh)
+INV_HTRIG_BENCH (atanh_approx5, math_approx::atanh<5>, data_atanh)
+INV_HTRIG_BENCH (atanh_approx4, math_approx::atanh<4>, data_atanh)
+INV_HTRIG_BENCH (atanh_approx3, math_approx::atanh<3>, data_atanh)
+
 #define INV_HTRIG_SIMD_BENCH(name, func, data) \
 void name (benchmark::State& state) \
 { \
@@ -75,4 +90,10 @@ INV_HTRIG_SIMD_BENCH (acosh_simd_approx5, math_approx::acosh<5>, data_acosh)
 INV_HTRIG_SIMD_BENCH (acosh_simd_approx4, math_approx::acosh<4>, data_acosh)
 INV_HTRIG_SIMD_BENCH (acosh_simd_approx3, math_approx::acosh<3>, data_acosh)
 
+INV_HTRIG_SIMD_BENCH (atanh_xsimd, xsimd::atanh, data_atanh)
+INV_HTRIG_SIMD_BENCH (atanh_simd_approx6, math_approx::atanh<6>, data_atanh)
+INV_HTRIG_SIMD_BENCH (atanh_simd_approx5, math_approx::atanh<5>, data_atanh)
+INV_HTRIG_SIMD_BENCH (atanh_simd_approx4, math_approx::atanh<4>, data_atanh)
+INV_HTRIG_SIMD_BENCH (atanh_simd_approx3, math_approx::atanh<3>, data_atanh)
+
 BENCHMARK_MAIN();

tools/plotter/plotter.cpp

@@ -61,12 +61,12 @@ void plot_function (std::span<const float> all_floats,
 int main()
 {
     plt::figure();
-    const auto range = std::make_pair (1.0f, 10.0f);
+    const auto range = std::make_pair (-0.99f, 0.99f);
     static constexpr auto tol = 1.0e-2f;
 
     const auto all_floats = test_helpers::all_32_bit_floats (range.first, range.second, tol);
-    const auto y_exact = test_helpers::compute_all<float> (all_floats, FLOAT_FUNC (std::acosh));
-    plot_ulp_error (all_floats, y_exact, FLOAT_FUNC ((math_approx::acosh<6>) ), "acosh-6");
+    const auto y_exact = test_helpers::compute_all<float> (all_floats, FLOAT_FUNC (std::atanh));
+    plot_error (all_floats, y_exact, FLOAT_FUNC ((math_approx::atanh<5>) ), "atanh_log-5");
 
     plt::legend ({ { "loc", "upper right" } });
     plt::xlim (range.first, range.second);