Re-organizing and adding acosh approximation

README.md

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

include/math_approx/math_approx.hpp

@@ -10,9 +10,8 @@ namespace math_approx
 #include "src/inverse_trig_approx.hpp"
 #include "src/pow_approx.hpp"
 #include "src/log_approx.hpp"
-#include "src/tanh_approx.hpp"
-#include "src/sinh_cosh_approx.hpp"
-#include "src/asinh_approx.hpp"
+#include "src/hyperbolic_trig_approx.hpp"
+#include "src/inverse_hyperbolic_trig_approx.hpp"
 #include "src/sigmoid_approx.hpp"
 #include "src/wright_omega_approx.hpp"
 #include "src/polylogarithm_approx.hpp"

include/math_approx/src/tanh_approx.hpp → include/math_approx/src/hyperbolic_trig_approx.hpp

@@ -1,9 +1,56 @@
 #pragma once
 
-#include "basic_math.hpp"
+#include "pow_approx.hpp"
 
 namespace math_approx
 {
+// ref: https://en.wikipedia.org/wiki/Hyperbolic_functions#Definitions
+// sinh = (e^(2x) - 1) / (2e^x), cosh = (e^(2x) + 1) / (2e^x)
+// let B = e^x, then sinh = (B^2 - 1) / (2B), cosh = (B^2 + 1) / (2B)
+// simplifying, we get: sinh = 0.5 (B - 1/B), cosh = 0.5 (B + 1/B)
+
+/** Approximation of sinh(x), using exp(x) internally */
+template <int order, typename T>
+constexpr T sinh (T x)
+{
+    using S = scalar_of_t<T>;
+    auto B = exp<order> (x);
+    auto Br = (S) 0.5 / B;
+    B *= (S) 0.5;
+    return B - Br;
+}
+
+/** Approximation of cosh(x), using exp(x) internally */
+template <int order, typename T>
+constexpr T cosh (T x)
+{
+    using S = scalar_of_t<T>;
+    auto B = exp<order> (x);
+    auto Br = (S) 0.5 / B;
+    B *= (S) 0.5;
+    return B + Br;
+}
+
+/**
+ * Simultaneous pproximation of sinh(x) and cosh(x),
+ * using exp(x) internally.
+ *
+ * For more information see the comments above.
+ */
+template <int order, typename T>
+constexpr auto sinh_cosh (T x)
+{
+    using S = scalar_of_t<T>;
+    auto B = exp<order> (x);
+    auto Br = (S) 0.5 / B;
+    B *= (S) 0.5;
+
+    auto sinh = B - Br;
+    auto cosh = B + Br;
+
+    return std::make_pair (sinh, cosh);
+}
+
 namespace tanh_detail
 {
     // These polynomial fits were generated from: https://www.wolframcloud.com/obj/chowdsp/Published/tanh_approx.nb

include/math_approx/src/asinh_approx.hpp → include/math_approx/src/inverse_hyperbolic_trig_approx.hpp

@@ -1,5 +1,6 @@
 #pragma once
 
+#include "basic_math.hpp"
 #include "log_approx.hpp"
 
 namespace math_approx
@@ -55,11 +56,9 @@ template <int order, typename T>
 constexpr T asinh (T x)
 {
     using S = scalar_of_t<T>;
-    using std::abs;
-    using std::sqrt;
+    using std::abs, std::sqrt;
 #if defined(XSIMD_HPP)
-    using xsimd::abs;
-    using xsimd::sqrt;
+    using xsimd::abs, xsimd::sqrt;
 #endif
 
     const auto sign = select (x > (S) 0, (T) (S) 1, select (x < (S) 0, (T) (S) -1, (T) (S) 0));
@@ -76,4 +75,18 @@ constexpr T asinh (T x)
 
     return sign * y;
 }
-} // namespace math_approx
+
+template <int order, typename T>
+constexpr T acosh (T x)
+{
+    using S = scalar_of_t<T>;
+    using std::sqrt;
+#if defined(XSIMD_HPP)
+    using xsimd::sqrt;
+#endif
+
+    const auto z0 = x - (S) 1;
+    const auto z1 = z0 + sqrt (z0 + z0 + z0 * z0);
+    return log1p<order> (z1);
+}
+}

test/CMakeLists.txt

@@ -28,13 +28,12 @@ function(setup_catch_test target)
     catch_discover_tests(${target} TEST_PREFIX ${target}_)
 endfunction(setup_catch_test)
 
-setup_catch_test(tanh_approx_test)
-setup_catch_test(sigmoid_approx_test)
 setup_catch_test(trig_approx_test)
 setup_catch_test(inverse_trig_approx_test)
 setup_catch_test(pow_approx_test)
 setup_catch_test(log_approx_test)
+setup_catch_test(hyperbolic_trig_approx_test)
+setup_catch_test(inverse_hyperbolic_trig_approx_test)
+setup_catch_test(sigmoid_approx_test)
 setup_catch_test(wright_omega_approx_test)
 setup_catch_test(polylog_approx_test)
-setup_catch_test(sinh_cosh_approx_test)
-setup_catch_test(asinh_approx_test)

test/src/sinh_cosh_approx_test.cpp → test/src/hyperbolic_trig_approx_test.cpp

@@ -145,3 +145,74 @@ TEST_CASE ("Cosh Approx Test")
                      0);
     }
 }
+
+TEST_CASE ("Tanh Approx Test")
+{
+#if ! defined(WIN32)
+    const auto all_floats = test_helpers::all_32_bit_floats (-10.0f, 10.0f, 1.0e-2f);
+#else
+    const auto all_floats = test_helpers::all_32_bit_floats (-10.0f, 10.0f, 1.0e-1f);
+#endif
+    const auto y_exact = test_helpers::compute_all<float> (all_floats, [] (auto x)
+                                                    { return std::tanh (x); });
+
+    const auto test_approx = [&all_floats, &y_exact] (auto&& f_approx, float err_bound, float rel_err_bound, uint32_t ulp_bound)
+    {
+        const auto y_approx = test_helpers::compute_all<float> (all_floats, f_approx);
+
+        const auto error = test_helpers::compute_error<float> (y_exact, y_approx);
+        const auto rel_error = test_helpers::compute_rel_error<float> (y_exact, y_approx);
+        const auto ulp_error = test_helpers::compute_ulp_error (y_exact, y_approx);
+
+        const auto max_error = test_helpers::abs_max<float> (error);
+        const auto max_rel_error = test_helpers::abs_max<float> (rel_error);
+        const auto max_ulp_error = *std::max_element (ulp_error.begin(), ulp_error.end());
+
+        std::cout << max_error << ", " << max_rel_error << ", " << max_ulp_error << std::endl;
+        REQUIRE (std::abs (max_error) < err_bound);
+        REQUIRE (std::abs (max_rel_error) < rel_err_bound);
+        if (ulp_bound > 0)
+            REQUIRE (max_ulp_error < ulp_bound);
+    };
+
+    SECTION ("11th-Order")
+    {
+        test_approx ([] (auto x)
+                     { return math_approx::tanh<11> (x); },
+                     2.5e-7f,
+                     4.0e-7f,
+                     7);
+    }
+    SECTION ("9th-Order")
+    {
+        test_approx ([] (auto x)
+                     { return math_approx::tanh<9> (x); },
+                     1.5e-6f,
+                     1.5e-6f,
+                     20);
+    }
+    SECTION ("7th-Order")
+    {
+        test_approx ([] (auto x)
+                     { return math_approx::tanh<7> (x); },
+                     1.5e-5f,
+                     1.5e-5f,
+                     230);
+    }
+    SECTION ("5th-Order")
+    {
+        test_approx ([] (auto x)
+                     { return math_approx::tanh<5> (x); },
+                     2.5e-4f,
+                     2.5e-4f,
+                     0);
+    }
+    SECTION ("3th-Order")
+    {
+        test_approx ([] (auto x)
+                     { return math_approx::tanh<3> (x); },
+                     4.0e-3f,
+                     4.0e-3f,
+                     0);
+    }
+}

test/src/asinh_approx_test.cpp → test/src/inverse_hyperbolic_trig_approx_test.cpp

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