Add implementations of sinh, cosh, and tan

README.md

@@ -4,10 +4,10 @@ A C++ library for math approximations.
 
 Currently supported:
 
-- sin/cos
+- sin/cos/tan
 - exp/exp2/exp10
 - log/log2/log10
-- tanh
+- sinh/cosh/tanh
 - sigmoid
 - Wright-Omega function
 - Dilogarithm function

include/math_approx/math_approx.hpp

@@ -6,10 +6,11 @@ namespace math_approx
 
 #include "src/basic_math.hpp"
 
-#include "src/tanh_approx.hpp"
-#include "src/sigmoid_approx.hpp"
 #include "src/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/sigmoid_approx.hpp"
 #include "src/wright_omega_approx.hpp"
 #include "src/polylogarithm_approx.hpp"

include/math_approx/src/pow_approx.hpp

@@ -10,7 +10,7 @@ namespace pow_detail
     template <typename T, int order, bool C1_continuous>
     constexpr T pow2_approx (T x)
     {
-        static_assert (order >= 3 && order <= 6);
+        static_assert (order >= 3 && order <= 7);
         using S = scalar_of_t<T>;
 
         const auto x_sq = x * x;
@@ -46,6 +46,17 @@ namespace pow_detail
                 const auto x_1_2_3_4_5_6 = x_1_2 + x_3_4_5_6 * x_sq;
                 return (S) 1 + x_1_2_3_4_5_6 * x;
             }
+            else if constexpr (order == 7)
+            {
+                // doesn't seem to help at single-precision
+                const auto x_6_7 = (S) 0.000133154170702612 + (S) 0.0000245778949916153 * x;
+                const auto x_4_5 = (S) 0.00960612128901630 + (S) 0.00135551454943593 * x;
+                const auto x_2_3 = (S) 0.240226202240181 + (S) 0.0555072492957270 * x;
+                const auto x_0_1 = (S) 1 + (S) 0.693147180559945 * x;
+                const auto x_4_5_6_7 = x_4_5 + x_6_7 * x_sq;
+                const auto x_0_1_2_3 = x_0_1 + x_2_3 * x_sq;
+                return x_0_1_2_3 + x_4_5_6_7 * x_sq * x_sq;
+            }
             else
             {
                 return {};
@@ -83,6 +94,17 @@ namespace pow_detail
                 const auto x_1_2_3_4_5_6 = x_1_2 + x_3_4_5_6 * x_sq;
                 return (S) 1 + x_1_2_3_4_5_6 * x;
             }
+            else if constexpr (order == 7)
+            {
+                // doesn't seem to help at single-precision
+                const auto x_6_7 = (S) 0.000136898688977877 + (S) 0.0000234440812713967 * x;
+                const auto x_4_5 = (S) 0.00960825566419915 + (S) 0.00135107295099880 * x;
+                const auto x_2_3 = (S) 0.240226092549669 + (S) 0.0555070350342468 * x;
+                const auto x_0_1 = (S) 1 + (S) 0.693147201030637 * x;
+                const auto x_4_5_6_7 = x_4_5 + x_6_7 * x_sq;
+                const auto x_0_1_2_3 = x_0_1 + x_2_3 * x_sq;
+                return x_0_1_2_3 + x_4_5_6_7 * x_sq * x_sq;
+            }
             else
             {
                 return {};

include/math_approx/src/sinh_cosh_approx.hpp

@@ -0,0 +1,53 @@
+#pragma once
+
+#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>
+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>
+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>
+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 math_approx

include/math_approx/src/trig_approx.hpp

@@ -4,7 +4,7 @@
 
 namespace math_approx
 {
-namespace sin_detail
+namespace trig_detail
 {
     template <typename T>
     T truncate (T x)
@@ -34,6 +34,20 @@ namespace sin_detail
         return select (x >= (T) 0, mod, mod + two_pi) - pi;
     }
 
+    /** Fast method to wrap a value into the range [-pi/2, pi/2] */
+    template <typename T>
+    T fast_mod_mhalfpi_halfpi (T x)
+    {
+        using S = scalar_of_t<T>;
+        static constexpr auto half_pi = static_cast<S> (M_PI) * (S) 0.5;
+        static constexpr auto pi = static_cast<S> (M_PI);
+        static constexpr auto recip_pi = (S) 1 / pi;
+
+        x += half_pi;
+        const auto mod = x - pi * truncate (x * recip_pi);
+        return select (x >= (T) 0, mod, mod + pi) - half_pi;
+    }
+
     template <typename T>
     T sin_poly_9 (T x, T x_sq)
     {
@@ -68,7 +82,7 @@ namespace sin_detail
 template <int order, typename T>
 T sin_mpi_pi (T x)
 {
-    static_assert (order % 2 == 1 && order <= 9 && order >= 5, "Order must e an odd number within [5, 9]");
+    static_assert (order % 2 == 1 && order <= 9 && order >= 5, "Order must be an odd number within [5, 9]");
 
     using S = scalar_of_t<T>;
     static constexpr auto pi = static_cast<S> (M_PI);
@@ -77,25 +91,25 @@ T sin_mpi_pi (T x)
 
     T x_poly {};
     if constexpr (order == 9)
-        x_poly = sin_detail::sin_poly_9 (x, x_sq);
+        x_poly = trig_detail::sin_poly_9 (x, x_sq);
     else if constexpr (order == 7)
-        x_poly = sin_detail::sin_poly_7 (x, x_sq);
+        x_poly = trig_detail::sin_poly_7 (x, x_sq);
     else if constexpr (order == 5)
-        x_poly = sin_detail::sin_poly_5 (x, x_sq);
+        x_poly = trig_detail::sin_poly_5 (x, x_sq);
 
     return (pi_sq - x_sq) * x_poly;
 }
 
 template <int order, typename T>
 T sin (T x)
 {
-    return sin_mpi_pi<order, T> (sin_detail::fast_mod_mpi_pi (x));
+    return sin_mpi_pi<order, T> (trig_detail::fast_mod_mpi_pi (x));
 }
 
 template <int order, typename T>
 T cos_mpi_pi (T x)
 {
-    static_assert (order % 2 == 1 && order <= 9 && order >= 5, "Order must e an odd number within [5, 9]");
+    static_assert (order % 2 == 1 && order <= 9 && order >= 5, "Order must be an odd number within [5, 9]");
 
     using S = scalar_of_t<T>;
     static constexpr auto pi = static_cast<S> (M_PI);
@@ -113,18 +127,116 @@ T cos_mpi_pi (T x)
 
     T x_poly {};
     if constexpr (order == 9)
-        x_poly = sin_detail::sin_poly_9 (hpmx, hpmx_sq);
+        x_poly = trig_detail::sin_poly_9 (hpmx, hpmx_sq);
     else if constexpr (order == 7)
-        x_poly = sin_detail::sin_poly_7 (hpmx, hpmx_sq);
+        x_poly = trig_detail::sin_poly_7 (hpmx, hpmx_sq);
     else if constexpr (order == 5)
-        x_poly = sin_detail::sin_poly_5 (hpmx, hpmx_sq);
+        x_poly = trig_detail::sin_poly_5 (hpmx, hpmx_sq);
 
     return (pi_sq - hpmx_sq) * x_poly;
 }
 
 template <int order, typename T>
 T cos (T x)
 {
-    return cos_mpi_pi<order, T> (sin_detail::fast_mod_mpi_pi (x));
+    return cos_mpi_pi<order, T> (trig_detail::fast_mod_mpi_pi (x));
+}
+
+/** Approximation of tan(x) on the range [-pi/4, pi/4] */
+template <int order, typename T>
+T tan_mquarterpi_quarterpi (T x)
+{
+    static_assert (order % 2 == 1 && order >= 3 && order <= 15, "Order must be an odd number within [3, 15]");
+
+    using S = scalar_of_t<T>;
+    const auto x_sq = x * x;
+    if constexpr (order == 3)
+    {
+        const auto x_1_3 = (S) 1 + (S) 0.442959265447 * x_sq;
+        return x * x_1_3;
+    }
+    else if constexpr (order == 5)
+    {
+        const auto x_3_5 = (S) 0.317574684334 + (S) 0.203265826702 * x_sq;
+        const auto x_1_3_5 = (S) 1 + x_3_5 * x_sq;
+        return x * x_1_3_5;
+    }
+    else if constexpr (order == 7)
+    {
+        const auto x_5_7 = (S) 0.116406244996 + (S) 0.0944480566104 * x_sq;
+        const auto x_1_3 = (S) 1 + (S) 0.335216153138 * x_sq;
+        const auto x_1_3_5_7 = x_1_3 + x_5_7 * x_sq * x_sq;
+        return x * x_1_3_5_7;
+    }
+    else if constexpr (order == 9)
+    {
+        const auto x_7_9 = (S) 0.0405232529373 + (S) 0.0439292071029 * x_sq;
+        const auto x_3_5 = (S) 0.333131667276 + (S) 0.136333765649 * x_sq;
+        const auto x_3_5_7_9 = x_3_5 + x_7_9 * x_sq * x_sq;
+        return x * ((S) 1 + x_3_5_7_9 * x_sq);
+    }
+    else if constexpr (order == 11)
+    {
+        const auto x_q = x_sq * x_sq;
+        const auto x_9_11 = (S) 0.0126603694551 + (S) 0.0203633469693 * x_sq;
+        const auto x_5_7 = (S) 0.132897195017 + (S) 0.0570525279731 * x_sq;
+        const auto x_1_3 = (S) 1 + (S) 0.333353019629 * x_sq;
+        const auto x_5_7_9_11 = x_5_7 + x_9_11 * x_q;
+        const auto x_1_3_5_7_9_11 = x_1_3 + x_5_7_9_11 * x_q;
+        return x * x_1_3_5_7_9_11;
+    }
+    else if constexpr (order == 13)
+    {
+        const auto x_q = x_sq * x_sq;
+        const auto x_6 = x_q * x_sq;
+        const auto x_11_13 = (S) 0.00343732283737 + (S) 0.00921082294855 * x_sq;
+        const auto x_7_9 = (S) 0.0534743904687 + (S) 0.0242183751709 * x_sq;
+        const auto x_3_5 = (S) 0.333331890901 + (S) 0.133379954680 * x_sq;
+        const auto x_7_9_11_13 = x_7_9 + x_11_13 * x_q;
+        const auto x_1_3_5 = (S) 1 + x_3_5 * x_sq;
+        return x * (x_1_3_5 + x_7_9_11_13 * x_6);
+    }
+    else if constexpr (order == 15)
+    {
+        // doesn't seem to help much at single-precision, but here it is:
+        const auto x_q = x_sq * x_sq;
+        const auto x_8 = x_q * x_q;
+        const auto x_13_15 = (S) 0.000292958045126 + (S) 0.00427933470414 * x_sq;
+        const auto x_9_11 = (S) 0.0213477960960 + (S) 0.0106702896251 * x_sq;
+        const auto x_5_7 = (S) 0.133327796402 + (S) 0.0540469276103* x_sq;
+        const auto x_1_3 = (S) 1 + (S) 0.333333463757 * x_sq;
+        const auto x_9_11_13_15 = x_9_11 + x_13_15 * x_q;
+        const auto x_1_3_5_7 = x_1_3 + x_5_7 * x_q;
+        const auto x_1_3_5_7_9_11_13_15 = x_1_3_5_7 + x_9_11_13_15 * x_8;
+        return x * x_1_3_5_7_9_11_13_15;
+    }
+    else
+    {
+        return {};
+    }
+}
+
+/**
+ * Approximation of tan(x) on the range [-pi/2, pi/2]
+ *
+ * Accuracy may suffer as x approaches ±pi/2.
+ */
+template <int order, typename T>
+T tan_mhalfpi_halfpi (T x)
+{
+    using S = scalar_of_t<T>;
+    const auto h_x = tan_mquarterpi_quarterpi<order> ((S) 0.5 * x);
+    return (S) 2 * h_x / ((S) 1 - h_x * h_x);
+}
+
+/**
+ * Full-range approximation of tan(x)
+ *
+ * Accuracy may suffer as x approaches values for which tan(x) approaches ±Inf.
+ */
+template <int order, typename T>
+T tan (T x)
+{
+    return tan_mhalfpi_halfpi<order> (trig_detail::fast_mod_mhalfpi_halfpi (x));
 }
 } // namespace math_approx

test/CMakeLists.txt

@@ -35,3 +35,4 @@ setup_catch_test(pow_approx_test)
 setup_catch_test(log_approx_test)
 setup_catch_test(wright_omega_approx_test)
 setup_catch_test(polylog_approx_test)
+setup_catch_test(sinh_cosh_approx_test)