Trig approx improvements

.github/workflows/run_tests.yml

@@ -42,4 +42,4 @@ jobs:
       - name: CMake Test
         run: |
           ctest --test-dir build -C RelWithDebInfo --show-only
-          ctest --test-dir build -C RelWithDebInfo -j4 --output-on-failure
+          ctest --test-dir build -C RelWithDebInfo -j2 --output-on-failure

CMakeLists.txt

@@ -15,6 +15,9 @@ endif()
 
 add_library(math_approx INTERFACE)
 target_include_directories(math_approx INTERFACE include)
+if(MSVC)
+    target_compile_definitions(math_approx INTERFACE _USE_MATH_DEFINES=1)
+endif()
 if (TARGET xsimd)
     message(STATUS "math_approx -- Linking with XSIMD...")
     target_link_libraries(math_approx INTERFACE xsimd)

include/math_approx/math_approx.hpp

@@ -8,4 +8,4 @@ namespace math_approx
 
 #include "src/tanh_approx.hpp"
 #include "src/sigmoid_approx.hpp"
-#include "src/sin_approx.hpp"
+#include "src/trig_approx.hpp"

include/math_approx/src/sin_approx.hpp → include/math_approx/src/trig_approx.hpp

@@ -99,11 +99,16 @@ T cos_mpi_pi (T x)
 
     using S = scalar_of_t<T>;
     static constexpr auto pi = static_cast<S> (M_PI);
-    static constexpr auto pi_o_2 = pi * (S) 0.5;;
+    static constexpr auto pi_sq = pi * pi;
+    static constexpr auto pi_o_2 = pi * (S) 0.5;
+
+    using std::abs;
+#if defined(XSIMD_HPP)
+    using xsimd::abs;
+#endif
+    x = abs (x);
 
-    const auto hpmx = (x > (S) 0 ? (S) 1 : (S) -1) * pi_o_2 - x;
-    const auto thpmx = (x > (S) 0 ? (S) 3 : (S) -3) * pi_o_2 - x;
-    const auto nhpmx = (x > (S) 0 ? (S) -1 : (S) 1) * pi_o_2 - x;
+    const auto hpmx = pi_o_2 - x;
     const auto hpmx_sq = hpmx * hpmx;
 
     T x_poly {};
@@ -114,7 +119,7 @@ T cos_mpi_pi (T x)
     else if constexpr (order == 5)
         x_poly = sin_detail::sin_poly_5 (hpmx, hpmx_sq);
 
-    return thpmx * nhpmx * (x > (S) 0 ? (S) -1 : (S) 1) * x_poly;
+    return (pi_sq - hpmx_sq) * x_poly;
 }
 
 template <int order, typename T>

test/CMakeLists.txt

@@ -30,5 +30,4 @@ endfunction(setup_catch_test)
 
 setup_catch_test(tanh_approx_test)
 setup_catch_test(sigmoid_approx_test)
-setup_catch_test(sin_approx_test)
-setup_catch_test(cos_approx_test)
+setup_catch_test(trig_approx_test)

test/src/sin_approx_test.cpp → test/src/trig_approx_test.cpp

@@ -44,3 +44,44 @@ TEST_CASE ("Sine Approx Test")
                      7.5e-4f);
     }
 }
+
+TEST_CASE ("Cosine Approx Test")
+{
+#if ! defined(WIN32)
+    const auto all_floats = test_helpers::all_32_bit_floats (-10.0f, 10.0f, 1.0e-3f);
+#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 (all_floats, [] (auto x)
+                                                    { return std::cos (x); });
+
+    const auto test_approx = [&all_floats, &y_exact] (auto&& f_approx, float err_bound)
+    {
+        const auto y_approx = test_helpers::compute_all (all_floats, f_approx);
+
+        const auto error = test_helpers::compute_error (y_exact, y_approx);
+        const auto max_error = test_helpers::abs_max (error);
+
+        // std::cout << max_error << std::endl;
+        REQUIRE (std::abs (max_error) < err_bound);
+    };
+
+    SECTION ("9th-Order")
+    {
+        test_approx ([] (auto x)
+                     { return math_approx::cos<9> (x); },
+                     7.5e-7f);
+    }
+    SECTION ("7th-Order")
+    {
+        test_approx ([] (auto x)
+                     { return math_approx::cos<7> (x); },
+                     1.8e-5f);
+    }
+    SECTION ("5th-Order")
+    {
+        test_approx ([] (auto x)
+                     { return math_approx::cos<5> (x); },
+                     7.5e-4f);
+    }
+}

tools/bench/CMakeLists.txt

@@ -10,3 +10,6 @@ target_link_libraries(tanh_approx_bench PRIVATE benchmark::benchmark math_approx
 
 add_executable(sigmoid_approx_bench sigmoid_bench.cpp)
 target_link_libraries(sigmoid_approx_bench PRIVATE benchmark::benchmark math_approx)
+
+add_executable(trig_approx_bench trig_bench.cpp)
+target_link_libraries(trig_approx_bench PRIVATE benchmark::benchmark math_approx)

tools/bench/trig_bench.cpp

@@ -0,0 +1,63 @@
+#include <math_approx/math_approx.hpp>
+#include <benchmark/benchmark.h>
+
+static constexpr size_t N = 2000;
+const auto data = []
+{
+    std::vector<float> x;
+    x.resize (N, 0.0f);
+    for (size_t i = 0; i < N; ++i)
+        x[i] = -10.0f + 20.0f * (float) i / (float) N;
+    return x;
+}();
+
+#define TRIG_BENCH(name, func) \
+void name (benchmark::State& state) \
+{ \
+for (auto _ : state) \
+{ \
+for (auto& x : data) \
+{ \
+auto y = func (x); \
+benchmark::DoNotOptimize (y); \
+} \
+} \
+} \
+BENCHMARK (name);
+
+TRIG_BENCH (cos_std, std::cos)
+TRIG_BENCH (cos_approx9, math_approx::cos<9>)
+TRIG_BENCH (cos_approx7, math_approx::cos<7>)
+TRIG_BENCH (cos_approx5, math_approx::cos<5>)
+
+TRIG_BENCH (sin_std, std::sin)
+TRIG_BENCH (sin_approx9, math_approx::sin<9>)
+TRIG_BENCH (sin_approx7, math_approx::sin<7>)
+TRIG_BENCH (sin_approx5, math_approx::sin<5>)
+
+#define TRIG_SIMD_BENCH(name, func) \
+void name (benchmark::State& state) \
+{ \
+for (auto _ : state) \
+{ \
+for (auto& x : data) \
+{ \
+auto y = func (xsimd::broadcast (x)); \
+static_assert (std::is_same_v<xsimd::batch<float>, decltype(y)>); \
+benchmark::DoNotOptimize (y); \
+} \
+} \
+} \
+BENCHMARK (name);
+
+TRIG_SIMD_BENCH (sin_xsimd, xsimd::sin)
+TRIG_SIMD_BENCH (sin_simd_approx9, math_approx::sin<9>)
+TRIG_SIMD_BENCH (sin_simd_approx7, math_approx::sin<7>)
+TRIG_SIMD_BENCH (sin_simd_approx5, math_approx::sin<5>)
+
+TRIG_SIMD_BENCH (cos_xsimd, xsimd::cos)
+TRIG_SIMD_BENCH (cos_simd_approx9, math_approx::cos<9>)
+TRIG_SIMD_BENCH (cos_simd_approx7, math_approx::cos<7>)
+TRIG_SIMD_BENCH (cos_simd_approx5, math_approx::cos<5>)
+
+BENCHMARK_MAIN();

tools/plotter/plotter.cpp

@@ -66,8 +66,8 @@ int main()
 
     // // plot_error (all_floats, y_exact, [] (float x) { return math_approx::sin<5> (x); }, "Sin-5");
     // // plot_error (all_floats, y_exact, [] (float x) { return math_approx::sin<7> (x); }, "Sin-7");
-    plot_ulp_error (all_floats, y_exact, [] (float x) { return math_approx::cos_mpi_pi<9> (x); }, "Cos-9");
-    // plot_function (all_floats, [] (float x) { return math_approx::cos_mpi_pi<9> (x); }, "Cos-9");
+    // plot_ulp_error (all_floats, y_exact, [] (float x) { return math_approx::cos_mpi_pi<9> (x); }, "Cos-9");
+    plot_function (all_floats, [] (float x) { return math_approx::cos_mpi_pi<9> (x); }, "Cos-9");
 
     plt::legend ({ { "loc", "upper right" } });
     plt::xlim (range.first, range.second);