TODO: first changes to add expf

* Fixed broken pybind11 target check

* Update python/CMakeLists.txt

Co-authored-by: Bram Veenboer <bram.veenboer@gmail.com>

---------

Co-authored-by: Bram Veenboer <bram.veenboer@gmail.com>

.github/workflows/linting.yml

@@ -0,0 +1,15 @@
+name: linting
+on:
+  push:
+jobs:
+  format:
+    name: Linting
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: '3.12'
+      - run: sudo apt-get update
+      - run: sudo apt-get install -y clang-format-14 cmake-format pre-commit
+      - run: pre-commit run -a

.github/workflows/test.yml

@@ -0,0 +1,89 @@
+name: test
+on:
+  push:
+
+jobs:
+  build:
+    runs-on: [slurm]
+    strategy:
+      matrix: &matrix
+        include:
+          - name: GNU, NVIDIA A4000
+            partition: defq
+            gres: gpu:A4000
+            cmake_flags: "-DTRIGDX_USE_MKL=1 -DTRIGDX_USE_GPU=1 -DTRIGDX_USE_MKL=1 -DTRIGDX_USE_XSIMD=1 -DCMAKE_CUDA_ARCHITECTURES=86"
+            environment_modules: "spack/20250403 intel-oneapi-mkl cuda python"
+          - name: Intel, NVIDIA A4000
+            partition: defq
+            gres: gpu:A4000
+            cmake_flags: "-DTRIGDX_USE_MKL=1 -DTRIGDX_USE_GPU=1 -DTRIGDX_USE_MKL=1 -DTRIGDX_USE_XSIMD=1 -DCMAKE_CUDA_ARCHITECTURES=86"
+            environment_modules: "spack/20250403 intel-oneapi-compilers intel-oneapi-mkl cuda python"
+          - name: NVIDIA GH200
+            partition: ghq
+            gres: gpu:GH200
+            cmake_flags: "-DTRIGDX_USE_GPU=1 -DTRIGDX_USE_XSIMD=1 -DCMAKE_CUDA_ARCHITECTURES=90a -DCMAKE_CUDA_COMPILER=/usr/local/cuda/bin/nvcc"
+    env:
+      DEVICE_NAME: ${{ matrix.name }}
+      PARTITION_NAME: ${{ matrix.partition }}
+      CMAKE_FLAGS: ${{ matrix.cmake_flags }}
+      ENVIRONMENT_MODULES: ${{ matrix.environment_modules }}
+    steps:
+      - &cleanup
+        name: "Remove old workspace"
+        run: |
+          rm -rf $GITHUB_WORKSPACE/*
+          mkdir -p $GITHUB_WORKSPACE
+      - uses: actions/checkout@v4
+      - uses: astron-rd/slurm-action@v1.2
+        with:
+          partition: ${{ matrix.partition }}
+          gres: ${{ matrix.gres }}
+          commands: |
+            module load ${ENVIRONMENT_MODULES}
+            cmake -S . -B build ${CMAKE_FLAGS}
+            make -C build -j
+      - name: Upload build
+        uses: pyTooling/upload-artifact@v4
+        with:
+          name: build-${{ matrix.name }}
+          path: build
+
+  test:
+    runs-on: [slurm]
+    needs: build
+    strategy:
+      matrix: *matrix
+    env:
+      DEVICE_NAME: ${{ matrix.name }}
+      PARTITION_NAME: ${{ matrix.partition }}
+    steps:
+      - *cleanup
+      - uses: pyTooling/download-artifact@v4
+        with:
+          name: build-${{ matrix.name }}
+      - uses: astron-rd/slurm-action@v1.2
+        with:
+          partition: ${{ matrix.partition }}
+          gres: ${{ matrix.gres }}
+          commands: |
+            find build/tests -type f -executable -exec {} \;
+
+  benchmark:
+    runs-on: [slurm]
+    needs: build
+    strategy:
+      matrix: *matrix
+    env:
+      DEVICE_NAME: ${{ matrix.name }}
+      PARTITION_NAME: ${{ matrix.partition }}
+    steps:
+      - *cleanup
+      - uses: pyTooling/download-artifact@v4
+        with:
+          name: build-${{ matrix.name }}
+      - uses: astron-rd/slurm-action@v1.2
+        with:
+          partition: ${{ matrix.partition }}
+          gres: ${{ matrix.gres }}
+          commands: |
+            find build/benchmarks -type f -executable -exec {} \;

.vscode/settings.json

@@ -1,77 +0,0 @@
-{
-  "files.associations": {
-    "*.cc": "cpp",
-    "*.tpp": "cpp",
-    "array": "cpp",
-    "atomic": "cpp",
-    "bit": "cpp",
-    "bitset": "cpp",
-    "cctype": "cpp",
-    "charconv": "cpp",
-    "chrono": "cpp",
-    "clocale": "cpp",
-    "cmath": "cpp",
-    "compare": "cpp",
-    "concepts": "cpp",
-    "condition_variable": "cpp",
-    "csignal": "cpp",
-    "cstdarg": "cpp",
-    "cstddef": "cpp",
-    "cstdint": "cpp",
-    "cstdio": "cpp",
-    "cstdlib": "cpp",
-    "cstring": "cpp",
-    "ctime": "cpp",
-    "cwchar": "cpp",
-    "cwctype": "cpp",
-    "deque": "cpp",
-    "list": "cpp",
-    "map": "cpp",
-    "set": "cpp",
-    "string": "cpp",
-    "unordered_map": "cpp",
-    "vector": "cpp",
-    "exception": "cpp",
-    "algorithm": "cpp",
-    "functional": "cpp",
-    "iterator": "cpp",
-    "memory": "cpp",
-    "memory_resource": "cpp",
-    "numeric": "cpp",
-    "optional": "cpp",
-    "random": "cpp",
-    "ratio": "cpp",
-    "regex": "cpp",
-    "string_view": "cpp",
-    "system_error": "cpp",
-    "tuple": "cpp",
-    "type_traits": "cpp",
-    "utility": "cpp",
-    "format": "cpp",
-    "fstream": "cpp",
-    "future": "cpp",
-    "initializer_list": "cpp",
-    "iomanip": "cpp",
-    "iosfwd": "cpp",
-    "iostream": "cpp",
-    "istream": "cpp",
-    "limits": "cpp",
-    "mutex": "cpp",
-    "new": "cpp",
-    "numbers": "cpp",
-    "ostream": "cpp",
-    "ranges": "cpp",
-    "semaphore": "cpp",
-    "span": "cpp",
-    "sstream": "cpp",
-    "stdexcept": "cpp",
-    "stop_token": "cpp",
-    "streambuf": "cpp",
-    "text_encoding": "cpp",
-    "thread": "cpp",
-    "typeinfo": "cpp",
-    "variant": "cpp",
-    "queue": "cpp",
-    "stack": "cpp"
-  }
-}

CMakeLists.txt

@@ -3,13 +3,40 @@ project(trigdx LANGUAGES CXX)
 
 set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 
-option(USE_MKL "Enable Intel MKL backend" OFF)
+option(TRIGDX_USE_MKL "Enable Intel MKL backend" OFF)
-option(USE_GPU "Enable GPU backend" OFF)
+option(TRIGDX_USE_GPU "Enable GPU backend" OFF)
-option(USE_XSIMD "Enable XSIMD backend" OFF)
+option(TRIGDX_USE_XSIMD "Enable XSIMD backend" OFF)
+option(TRIGDX_BUILD_TESTS "Build tests" ON)
+option(TRIGDX_BUILD_BENCHMARKS "Build tests" ON)
+option(TRIGDX_BUILD_PYTHON "Build Python interface" ON)
+
+list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake")
+configure_file(
+  ${CMAKE_CURRENT_SOURCE_DIR}/cmake/trigdx_config.hpp.in
+  ${CMAKE_CURRENT_BINARY_DIR}/include/trigdx/trigdx_config.hpp @ONLY)
+
+if(TRIGDX_BUILD_TESTS
+   OR TRIGDX_BUILD_BENCHMARKS
+   OR TRIGDX_BUILD_PYTHON)
+  include(FetchContent)
+endif()
 
 include_directories(${PROJECT_SOURCE_DIR}/include)
 
 add_subdirectory(src)
-add_subdirectory(tests)
+
-add_subdirectory(benchmarks)
+if(TRIGDX_BUILD_TESTS)
+  include(CTest)
+  enable_testing()
+  add_subdirectory(tests)
+endif()
+
+if(TRIGDX_BUILD_BENCHMARKS)
+  add_subdirectory(benchmarks)
+endif()
+
+if(TRIGDX_BUILD_PYTHON)
+  add_subdirectory(python)
+endif()

benchmarks/CMakeLists.txt

@@ -1,23 +1,36 @@
+FetchContent_Declare(
+  benchmark
+  GIT_REPOSITORY https://github.com/google/benchmark.git
+  GIT_TAG v1.9.4)
+set(BENCHMARK_ENABLE_TESTING
+    OFF
+    CACHE BOOL "" FORCE)
+FetchContent_MakeAvailable(benchmark)
+
 add_executable(benchmark_reference benchmark_reference.cpp)
-target_link_libraries(benchmark_reference PRIVATE trigdx)
+target_link_libraries(benchmark_reference PRIVATE trigdx benchmark::benchmark)
 
 add_executable(benchmark_lookup benchmark_lookup.cpp)
-target_link_libraries(benchmark_lookup PRIVATE trigdx)
+target_link_libraries(benchmark_lookup PRIVATE trigdx benchmark::benchmark)
 
-add_executable(benchmark_lookup_avx benchmark_lookup_avx.cpp)
+if(HAVE_AVX)
-target_link_libraries(benchmark_lookup_avx PRIVATE trigdx)
+  add_executable(benchmark_lookup_avx benchmark_lookup_avx.cpp)
+  target_link_libraries(benchmark_lookup_avx PRIVATE trigdx
+                                                     benchmark::benchmark)
+endif()
 
-if(USE_MKL)
+if(TRIGDX_USE_MKL)
   add_executable(benchmark_mkl benchmark_mkl.cpp)
-  target_link_libraries(benchmark_mkl PRIVATE trigdx)
+  target_link_libraries(benchmark_mkl PRIVATE trigdx benchmark::benchmark)
 endif()
 
-if(USE_GPU)
+if(TRIGDX_USE_GPU)
   add_executable(benchmark_gpu benchmark_gpu.cpp)
-  target_link_libraries(benchmark_gpu PRIVATE trigdx gpu)
+  target_link_libraries(benchmark_gpu PRIVATE trigdx gpu benchmark::benchmark)
 endif()
 
-if(USE_XSIMD)
+if(TRIGDX_USE_XSIMD)
   add_executable(benchmark_lookup_xsimd benchmark_lookup_xsimd.cpp)
-  target_link_libraries(benchmark_lookup_xsimd PRIVATE trigdx)
+  target_link_libraries(benchmark_lookup_xsimd PRIVATE trigdx
+                                                       benchmark::benchmark)
 endif()

benchmarks/benchmark_gpu.cpp

@@ -1,9 +1,21 @@
-#include <trigdx/gpu.hpp>
+#include <trigdx/trigdx.hpp>
 
 #include "benchmark_utils.hpp"
 
-int main() {
+BENCHMARK_TEMPLATE(benchmark_sinf, GPUBackend)
-  benchmark_sinf<GPUBackend>();
+    ->Unit(benchmark::kMillisecond)
-  benchmark_cosf<GPUBackend>();
+    ->Arg(1e5)
-  benchmark_sincosf<GPUBackend>();
+    ->Arg(1e6)
-}
+    ->Arg(1e7);
+BENCHMARK_TEMPLATE(benchmark_cosf, GPUBackend)
+    ->Unit(benchmark::kMillisecond)
+    ->Arg(1e5)
+    ->Arg(1e6)
+    ->Arg(1e7);
+BENCHMARK_TEMPLATE(benchmark_sincosf, GPUBackend)
+    ->Unit(benchmark::kMillisecond)
+    ->Arg(1e5)
+    ->Arg(1e6)
+    ->Arg(1e7);
+
+BENCHMARK_MAIN();

benchmarks/benchmark_lookup.cpp

@@ -1,13 +1,30 @@
-#include <trigdx/lookup.hpp>
+#include <trigdx/trigdx.hpp>
 
 #include "benchmark_utils.hpp"
 
-int main() {
+template <typename Backend> void register_benchmarks() {
-  benchmark_sinf<LookupBackend<16384>>();
+  BENCHMARK_TEMPLATE(benchmark_sinf, Backend)
-  benchmark_cosf<LookupBackend<16384>>();
+      ->Unit(benchmark::kMillisecond)
-  benchmark_sincosf<LookupBackend<16384>>();
+      ->Arg(1e5)
-
+      ->Arg(1e6)
-  benchmark_sinf<LookupBackend<32768>>();
+      ->Arg(1e7);
-  benchmark_cosf<LookupBackend<32768>>();
+  BENCHMARK_TEMPLATE(benchmark_cosf, Backend)
-  benchmark_sincosf<LookupBackend<32768>>();
+      ->Unit(benchmark::kMillisecond)
+      ->Arg(1e5)
+      ->Arg(1e6)
+      ->Arg(1e7);
+  BENCHMARK_TEMPLATE(benchmark_sincosf, Backend)
+      ->Unit(benchmark::kMillisecond)
+      ->Arg(1e5)
+      ->Arg(1e6)
+      ->Arg(1e7);
 }
+
+int main(int argc, char **argv) {
+  ::benchmark::Initialize(&argc, argv);
+
+  register_benchmarks<LookupBackend<16384>>();
+  register_benchmarks<LookupBackend<32768>>();
+
+  return ::benchmark::RunSpecifiedBenchmarks();
+}

benchmarks/benchmark_lookup_avx.cpp

@@ -1,13 +1,30 @@
-#include <trigdx/lookup_avx.hpp>
+#include <trigdx/trigdx.hpp>
 
 #include "benchmark_utils.hpp"
 
-int main() {
+template <typename Backend> void register_benchmarks() {
-  benchmark_sinf<LookupAVXBackend<16384>>();
+  BENCHMARK_TEMPLATE(benchmark_sinf, Backend)
-  benchmark_cosf<LookupAVXBackend<16384>>();
+      ->Unit(benchmark::kMillisecond)
-  benchmark_sincosf<LookupAVXBackend<16384>>();
+      ->Arg(1e5)
-
+      ->Arg(1e6)
-  benchmark_sinf<LookupAVXBackend<32768>>();
+      ->Arg(1e7);
-  benchmark_cosf<LookupAVXBackend<32768>>();
+  BENCHMARK_TEMPLATE(benchmark_cosf, Backend)
-  benchmark_sincosf<LookupAVXBackend<32768>>();
+      ->Unit(benchmark::kMillisecond)
+      ->Arg(1e5)
+      ->Arg(1e6)
+      ->Arg(1e7);
+  BENCHMARK_TEMPLATE(benchmark_sincosf, Backend)
+      ->Unit(benchmark::kMillisecond)
+      ->Arg(1e5)
+      ->Arg(1e6)
+      ->Arg(1e7);
 }
+
+int main(int argc, char **argv) {
+  ::benchmark::Initialize(&argc, argv);
+
+  register_benchmarks<LookupAVXBackend<16384>>();
+  register_benchmarks<LookupAVXBackend<32768>>();
+
+  return ::benchmark::RunSpecifiedBenchmarks();
+}

benchmarks/benchmark_lookup_xsimd.cpp

@@ -1,13 +1,30 @@
-#include <trigdx/lookup_xsimd.hpp>
+#include <trigdx/trigdx.hpp>
 
 #include "benchmark_utils.hpp"
 
-int main() {
+template <typename Backend> void register_benchmarks() {
-  benchmark_sinf<LookupXSIMDBackend<16384>>();
+  BENCHMARK_TEMPLATE(benchmark_sinf, Backend)
-  benchmark_cosf<LookupXSIMDBackend<16384>>();
+      ->Unit(benchmark::kMillisecond)
-  benchmark_sincosf<LookupXSIMDBackend<16384>>();
+      ->Arg(1e5)
-
+      ->Arg(1e6)
-  benchmark_sinf<LookupXSIMDBackend<32768>>();
+      ->Arg(1e7);
-  benchmark_cosf<LookupXSIMDBackend<32768>>();
+  BENCHMARK_TEMPLATE(benchmark_cosf, Backend)
-  benchmark_sincosf<LookupXSIMDBackend<32768>>();
+      ->Unit(benchmark::kMillisecond)
+      ->Arg(1e5)
+      ->Arg(1e6)
+      ->Arg(1e7);
+  BENCHMARK_TEMPLATE(benchmark_sincosf, Backend)
+      ->Unit(benchmark::kMillisecond)
+      ->Arg(1e5)
+      ->Arg(1e6)
+      ->Arg(1e7);
 }
+
+int main(int argc, char **argv) {
+  ::benchmark::Initialize(&argc, argv);
+
+  register_benchmarks<LookupXSIMDBackend<16384>>();
+  register_benchmarks<LookupXSIMDBackend<32768>>();
+
+  return ::benchmark::RunSpecifiedBenchmarks();
+}

benchmarks/benchmark_mkl.cpp

@@ -1,9 +1,21 @@
-#include <trigdx/mkl.hpp>
+#include <trigdx/trigdx.hpp>
 
 #include "benchmark_utils.hpp"
 
-int main() {
+BENCHMARK_TEMPLATE(benchmark_sinf, MKLBackend)
-  benchmark_sinf<MKLBackend>();
+    ->Unit(benchmark::kMillisecond)
-  benchmark_cosf<MKLBackend>();
+    ->Arg(1e5)
-  benchmark_sincosf<MKLBackend>();
+    ->Arg(1e6)
-}
+    ->Arg(1e7);
+BENCHMARK_TEMPLATE(benchmark_cosf, MKLBackend)
+    ->Unit(benchmark::kMillisecond)
+    ->Arg(1e5)
+    ->Arg(1e6)
+    ->Arg(1e7);
+BENCHMARK_TEMPLATE(benchmark_sincosf, MKLBackend)
+    ->Unit(benchmark::kMillisecond)
+    ->Arg(1e5)
+    ->Arg(1e6)
+    ->Arg(1e7);
+
+BENCHMARK_MAIN();

benchmarks/benchmark_reference.cpp

@@ -1,9 +1,21 @@
-#include <trigdx/reference.hpp>
+#include <trigdx/trigdx.hpp>
 
 #include "benchmark_utils.hpp"
 
-int main() {
+BENCHMARK_TEMPLATE(benchmark_sinf, ReferenceBackend)
-  benchmark_sinf<ReferenceBackend>();
+    ->Unit(benchmark::kMillisecond)
-  benchmark_cosf<ReferenceBackend>();
+    ->Arg(1e5)
-  benchmark_sincosf<ReferenceBackend>();
+    ->Arg(1e6)
-}
+    ->Arg(1e7);
+BENCHMARK_TEMPLATE(benchmark_cosf, ReferenceBackend)
+    ->Unit(benchmark::kMillisecond)
+    ->Arg(1e5)
+    ->Arg(1e6)
+    ->Arg(1e7);
+BENCHMARK_TEMPLATE(benchmark_sincosf, ReferenceBackend)
+    ->Unit(benchmark::kMillisecond)
+    ->Arg(1e5)
+    ->Arg(1e6)
+    ->Arg(1e7);
+
+BENCHMARK_MAIN();

benchmarks/benchmark_utils.hpp

@@ -1,76 +1,90 @@
 #pragma once
 
 #include <chrono>
-#include <iomanip>
+#include <cmath>
-#include <iostream>
+#include <string>
 #include <vector>
 
-const size_t N = 1e7;
+#include <benchmark/benchmark.h>
 
-inline void report(const std::string &name, double sec, double throughput) {
+void init_x(std::vector<float> &x) {
-  std::ios state(nullptr);
+  for (size_t i = 0; i < x.size(); ++i) {
-  state.copyfmt(std::cout);
+    x[i] = (i % 360) * 0.0174533f; // degrees to radians
-  std::cout << std::setw(7) << name << " -> ";
+  }
-  std::cout << "time: ";
-  std::cout << std::fixed << std::setprecision(3) << std::setfill('0');
-  std::cout << sec << " s, ";
-  std::cout << "throughput: " << throughput << " M elems/sec\n";
-  std::cout.copyfmt(state);
 }
 
-template <typename Backend> inline void benchmark_sinf() {
+template <typename Backend>
+static void benchmark_sinf(benchmark::State &state) {
+  const size_t N = static_cast<size_t>(state.range(0));
   std::vector<float> x(N), s(N);
-
+  init_x(x);
-  for (size_t i = 0; i < N; ++i)
-    x[i] = (i % 360) * 0.0174533f; // degrees to radians
 
   Backend backend;
-  backend.init(N);
 
   auto start = std::chrono::high_resolution_clock::now();
-  backend.compute_sinf(N, x.data(), s.data());
+  backend.init(N);
   auto end = std::chrono::high_resolution_clock::now();
+  state.counters["init_ms"] =
+      std::chrono::duration_cast<std::chrono::microseconds>(end - start)
+          .count() /
+      1.e3;
 
-  double sec = std::chrono::duration<double>(end - start).count();
+  for (auto _ : state) {
-  double throughput = N / sec / 1e6;
+    backend.compute_sinf(N, x.data(), s.data());
+    benchmark::DoNotOptimize(s);
+  }
 
-  report("sinf", sec, throughput);
+  state.SetItemsProcessed(static_cast<int64_t>(state.iterations()) *
+                          static_cast<int64_t>(N));
 }
 
-template <typename Backend> inline void benchmark_cosf() {
+template <typename Backend>
+static void benchmark_cosf(benchmark::State &state) {
+  const size_t N = static_cast<size_t>(state.range(0));
   std::vector<float> x(N), c(N);
-
+  init_x(x);
-  for (size_t i = 0; i < N; ++i)
-    x[i] = (i % 360) * 0.0174533f; // degrees to radians
 
   Backend backend;
-  backend.init(N);
 
   auto start = std::chrono::high_resolution_clock::now();
-  backend.compute_cosf(N, x.data(), c.data());
+  backend.init(N);
   auto end = std::chrono::high_resolution_clock::now();
+  state.counters["init_ms"] =
+      std::chrono::duration_cast<std::chrono::microseconds>(end - start)
+          .count() /
+      1.e3;
 
-  double sec = std::chrono::duration<double>(end - start).count();
+  for (auto _ : state) {
-  double throughput = N / sec / 1e6;
+    backend.compute_cosf(N, x.data(), c.data());
+    benchmark::DoNotOptimize(c);
+  }
 
-  report("cosf", sec, throughput);
+  state.SetItemsProcessed(static_cast<int64_t>(state.iterations()) *
+                          static_cast<int64_t>(N));
 }
 
-template <typename Backend> inline void benchmark_sincosf() {
+template <typename Backend>
+static void benchmark_sincosf(benchmark::State &state) {
+  const size_t N = static_cast<size_t>(state.range(0));
   std::vector<float> x(N), s(N), c(N);
-
+  init_x(x);
-  for (size_t i = 0; i < N; ++i)
-    x[i] = (i % 360) * 0.0174533f; // degrees to radians
 
   Backend backend;
-  backend.init(N);
 
   auto start = std::chrono::high_resolution_clock::now();
-  backend.compute_sincosf(N, x.data(), s.data(), c.data());
+  backend.init(N);
   auto end = std::chrono::high_resolution_clock::now();
+  state.counters["init_ms"] =
+      std::chrono::duration_cast<std::chrono::microseconds>(end - start)
+          .count() /
+      1.e3;
 
-  double sec = std::chrono::duration<double>(end - start).count();
+  for (auto _ : state) {
-  double throughput = N / sec / 1e6;
+    backend.compute_sincosf(N, x.data(), s.data(), c.data());
+    benchmark::DoNotOptimize(s);
+    benchmark::DoNotOptimize(c);
+  }
 
-  report("sincosf", sec, throughput);
+  state.SetItemsProcessed(static_cast<int64_t>(state.iterations()) *
+                          static_cast<int64_t>(N));
 }

cmake/FindAVX.cmake

@@ -0,0 +1,42 @@
+include(CheckCXXSourceRuns)
+
+set(SUPPORTED_COMPILERS Clang;GNU;Intel;IntelLLVM)
+
+if(CMAKE_CXX_COMPILER_ID IN_LIST SUPPORTED_COMPILERS)
+  if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
+    set(CMAKE_REQUIRED_FLAGS "-xHost") # ICC
+  elseif(CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM")
+    set(CMAKE_REQUIRED_FLAGS "-march=native") # ICX
+  else()
+    set(CMAKE_REQUIRED_FLAGS "-march=native") # GCC/Clang
+  endif()
+else()
+  message(FATAL_ERROR "Unsupported compiler: ${CMAKE_CXX_COMPILER_ID}.")
+endif()
+
+# AVX check
+check_cxx_source_runs(
+  "
+#include <immintrin.h>
+int main() {
+    __m256 a = _mm256_setzero_ps(); // AVX
+    (void) a;
+    return 0;
+}
+"
+  HAVE_AVX)
+
+if(HAVE_AVX)
+  # AVX2 check
+  check_cxx_source_runs(
+    "
+#include <immintrin.h>
+int main() {
+    __m256i a = _mm256_set1_epi32(-1);
+    __m256i b = _mm256_abs_epi32(a); // AVX2
+    (void) b;
+    return 0;
+}
+"
+    HAVE_AVX2)
+endif()

cmake/trigdx_config.hpp.in

@@ -0,0 +1,5 @@
+#pragma once
+
+#cmakedefine TRIGDX_USE_MKL
+#cmakedefine TRIGDX_USE_GPU
+#cmakedefine TRIGDX_USE_XSIMD

include/trigdx/gpu.hpp

@@ -16,6 +16,7 @@ public:
   void compute_cosf(size_t n, const float *x, float *c) const override;
   void compute_sincosf(size_t n, const float *x, float *s,
                        float *c) const override;
+  void compute_expf(size_t n, const float *x, float *e) const override;
 
 private:
   struct Impl;

include/trigdx/interface.hpp

@@ -19,4 +19,7 @@ public:
   // Compute sine and cosine for n elements
   virtual void compute_sincosf(size_t n, const float *x, float *s,
                                float *c) const = 0;
+
+  // Compute exponent for n elements
+  virtual void compute_expf(size_t n, const float *x, float *e) const = 0;
 };

include/trigdx/mkl.hpp

@@ -10,4 +10,6 @@ public:
 
   void compute_sincosf(size_t n, const float *x, float *s,
                        float *c) const override;
+
+  void compute_expf(size_t n, const float *x, float *e) const override;
 };

include/trigdx/reference.hpp

@@ -10,4 +10,6 @@ public:
 
   void compute_sincosf(size_t n, const float *x, float *s,
                        float *c) const override;
+
+  void compute_expf(size_t n, const float *x, float *e) const override;
 };

include/trigdx/trigdx.hpp

@@ -0,0 +1,19 @@
+#pragma once
+
+#include <trigdx/trigdx_config.hpp>
+
+#include <trigdx/lookup.hpp>
+#include <trigdx/lookup_avx.hpp>
+#include <trigdx/reference.hpp>
+
+#if defined(TRIGDX_USE_MKL)
+#include <trigdx/mkl.hpp>
+#endif
+
+#if defined(TRIGDX_USE_GPU)
+#include <trigdx/gpu.hpp>
+#endif
+
+#if defined(TRIGDX_USE_XSIMD)
+#include <trigdx/lookup_xsimd.hpp>
+#endif

python/CMakeLists.txt

@@ -0,0 +1,12 @@
+find_package(pybind11 CONFIG QUIET)
+
+if(NOT pybind11_FOUND)
+  FetchContent_Declare(
+    pybind11
+    GIT_REPOSITORY https://github.com/pybind/pybind11.git
+    GIT_TAG v3.0.0)
+  FetchContent_MakeAvailable(pybind11)
+endif()
+
+pybind11_add_module(pytrigdx bindings.cpp)
+target_link_libraries(pytrigdx PRIVATE trigdx)

python/bindings.cpp

@@ -0,0 +1,94 @@
+#include <memory>
+
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <trigdx/trigdx.hpp>
+
+namespace py = pybind11;
+
+template <typename T>
+py::array_t<T>
+compute_sin(const Backend &backend,
+            py::array_t<T, py::array::c_style | py::array::forcecast> x) {
+  const size_t n = x.shape(0);
+  if (x.ndim() != 1) {
+    throw py::value_error("Input array must be 1-dimensional");
+  }
+
+  const T *x_ptr = x.data();
+
+  py::array_t<float> s(n);
+  T *s_ptr = s.mutable_data();
+
+  backend.compute_sinf(n, x_ptr, s_ptr);
+
+  return s;
+}
+
+template <typename T>
+py::array_t<T>
+compute_cos(const Backend &backend,
+            py::array_t<T, py::array::c_style | py::array::forcecast> x) {
+  const size_t n = x.shape(0);
+  if (x.ndim() != 1) {
+    throw py::value_error("Input array must be 1-dimensional");
+  }
+
+  const T *x_ptr = x.data();
+
+  py::array_t<T> c(n);
+  T *c_ptr = c.mutable_data();
+
+  backend.compute_cosf(n, x_ptr, c_ptr);
+
+  return c;
+}
+
+template <typename T>
+std::tuple<py::array_t<T>, py::array_t<T>>
+compute_sincos(const Backend &backend,
+               py::array_t<T, py::array::c_style | py::array::forcecast> x) {
+  const size_t n = x.shape(0);
+  if (x.ndim() != 1) {
+    throw py::value_error("Input array must be 1-dimensional");
+  }
+
+  const T *x_ptr = x.data();
+
+  py::array_t<T> s(n);
+  py::array_t<T> c(n);
+
+  backend.compute_sincosf(n, x_ptr, s.mutable_data(), c.mutable_data());
+
+  return std::make_tuple(s, c);
+}
+
+template <typename BackendType>
+void bind_backend(py::module &m, const char *name) {
+  py::class_<BackendType, Backend, std::shared_ptr<BackendType>>(m, name)
+      .def(py::init<>())
+      .def("compute_sinf", &compute_sin<float>)
+      .def("compute_cosf", &compute_cos<float>)
+      .def("compute_sincosf", &compute_sincos<float>);
+}
+
+PYBIND11_MODULE(pytrigdx, m) {
+  py::class_<Backend, std::shared_ptr<Backend>>(m, "Backend")
+      .def("init", &Backend::init);
+
+  bind_backend<ReferenceBackend>(m, "Reference");
+  bind_backend<LookupBackend<16384>>(m, "Lookup16K");
+  bind_backend<LookupBackend<32768>>(m, "Lookup32K");
+  bind_backend<LookupAVXBackend<16384>>(m, "LookupAVX16K");
+  bind_backend<LookupAVXBackend<32768>>(m, "LookupAVX32K");
+#if defined(TRIGDX_USE_MKL)
+  bind_backend<MKLBackend>(m, "MKL");
+#endif
+#if defined(TRIGDX_USE_GPU)
+  bind_backend<GPUBackend>(m, "GPU");
+#endif
+#if defined(TRIGDX_USE_XSIMD)
+  bind_backend<LookupXSIMDBackend<16384>>(m, "LookupXSIMD16K");
+  bind_backend<LookupXSIMDBackend<32768>>(m, "LookupXSIMD32K");
+#endif
+}

src/CMakeLists.txt

@@ -1,16 +1,20 @@
-add_library(trigdx reference.cpp lookup.cpp lookup_avx.cpp)
+include(FetchContent)
+include(FindAVX)
+add_library(trigdx reference.cpp lookup.cpp)
 
 target_include_directories(trigdx PUBLIC ${PROJECT_SOURCE_DIR}/include)
 
-target_compile_options(trigdx PRIVATE -O3 -march=native)
+if(HAVE_AVX)
+  target_sources(trigdx PRIVATE lookup_avx.cpp)
+endif()
 
-if(USE_MKL)
+if(TRIGDX_USE_MKL)
   find_package(MKL REQUIRED)
   target_sources(trigdx PRIVATE mkl.cpp)
   target_link_libraries(trigdx PRIVATE MKL::MKL)
 endif()
 
-if(USE_GPU)
+if(TRIGDX_USE_GPU)
   enable_language(CUDA)
   find_package(CUDAToolkit REQUIRED)
   add_library(gpu SHARED gpu/gpu.cu)
@@ -19,8 +23,20 @@ if(USE_GPU)
   target_link_libraries(trigdx PRIVATE gpu)
 endif()
 
-if(USE_XSIMD)
+if(TRIGDX_USE_XSIMD)
-  find_package(xsimd REQUIRED)
+  # Requires XSIMD > 13 for architecture independent dispatching
+  find_package(xsimd 13 QUIET)
+  if(NOT TARGET xsimd)
+    FetchContent_Declare(
+      xsimd
+      GIT_REPOSITORY https://github.com/xtensor-stack/xsimd.git
+      GIT_TAG 13.2.0)
+    FetchContent_MakeAvailable(xsimd)
+  endif()
   target_sources(trigdx PRIVATE lookup_xsimd.cpp)
   target_link_libraries(trigdx PRIVATE xsimd)
 endif()
+
+target_include_directories(
+  trigdx INTERFACE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include>
+                   $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>)

src/gpu.cpp

@@ -20,6 +20,9 @@ struct GPUBackend::Impl {
     if (h_c) {
       cudaFreeHost(h_c);
     }
+    if (h_e) {
+      cudaFreeHost(h_e);
+    }
     if (d_x) {
       cudaFree(d_x);
     }
@@ -29,6 +32,9 @@ struct GPUBackend::Impl {
     if (d_c) {
       cudaFree(d_c);
     }
+    if (d_e) {
+      cudaFree(d_e);
+    }
   }
 
   void init(size_t n) {
@@ -36,42 +42,59 @@ struct GPUBackend::Impl {
     cudaMallocHost(&h_x, bytes);
     cudaMallocHost(&h_s, bytes);
     cudaMallocHost(&h_c, bytes);
+    cudaMallocHost(&h_e, bytes);
     cudaMalloc(&d_x, bytes);
     cudaMalloc(&d_s, bytes);
     cudaMalloc(&d_c, bytes);
+    cudaMalloc(&d_e, bytes);
   }
 
   void compute_sinf(size_t n, const float *x, float *s) const {
     const size_t bytes = n * sizeof(float);
     std::memcpy(h_x, x, bytes);
     cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice);
-    launch_sincosf_kernel(d_x, d_s, d_c, n);
+    launch_sinf_kernel(d_x, d_s, n);
     cudaMemcpy(h_s, d_s, bytes, cudaMemcpyDeviceToHost);
     std::memcpy(s, h_s, bytes);
   }
 
   void compute_cosf(size_t n, const float *x, float *c) const {
     const size_t bytes = n * sizeof(float);
+    std::memcpy(h_x, x, bytes);
     cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice);
-    launch_sincosf_kernel(d_x, d_s, d_c, n);
+    launch_cosf_kernel(d_x, d_c, n);
     cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);
     std::memcpy(c, h_c, bytes);
   }
 
   void compute_sincosf(size_t n, const float *x, float *s, float *c) const {
     const size_t bytes = n * sizeof(float);
+    std::memcpy(h_x, x, bytes);
     cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice);
     launch_sincosf_kernel(d_x, d_s, d_c, n);
     cudaMemcpy(h_s, d_s, bytes, cudaMemcpyDeviceToHost);
     cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);
+    std::memcpy(s, h_s, bytes);
+    std::memcpy(c, h_c, bytes);
+  }
+
+  void compute_expf(size_t n, const float *x, float *e) const {
+    const size_t bytes = n * sizeof(float);
+    std::memcpy(h_x, x, bytes);
+    cudaMemcpy(d_x, h_x, bytes, cudaMemcpyHostToDevice);
+    launch_expf_kernel(d_x, d_e, n);
+    cudaMemcpy(h_e, d_e, bytes, cudaMemcpyDeviceToHost);
+    std::memcpy(e, h_e, bytes);
   }
 
   float *h_x = nullptr;
   float *h_s = nullptr;
   float *h_c = nullptr;
+  float *h_e = nullptr;
   float *d_x = nullptr;
   float *d_s = nullptr;
   float *d_c = nullptr;
+  float *d_e = nullptr;
 };
 
 GPUBackend::GPUBackend() : impl(std::make_unique<Impl>()) {}
@@ -91,4 +114,8 @@ void GPUBackend::compute_cosf(size_t n, const float *x, float *c) const {
 void GPUBackend::compute_sincosf(size_t n, const float *x, float *s,
                                  float *c) const {
   impl->compute_sincosf(n, x, s, c);
-}
+}
+
+void GPUBackend::compute_expf(size_t n, const float *x, float *e) const {
+  impl->compute_expf(n, x, e);
+}

src/gpu/gpu.cu

@@ -31,6 +31,15 @@ __global__ void kernel_sincosf(const float *__restrict__ x,
   }
 }
 
+__global__ void kernel_expf(const float *__restrict__ x, float *__restrict__ e,
+                            size_t n) {
+  size_t idx = blockIdx.x * blockDim.x + threadIdx.x;
+  if (idx < n) {
+    // e[idx] = __expf(x[idx]);
+    e[idx] = expf(x[idx]);
+  }
+}
+
 namespace {
 inline dim3 make_grid(size_t n, size_t threadsPerBlock = 256) {
   return dim3((n + threadsPerBlock - 1) / threadsPerBlock);
@@ -54,3 +63,9 @@ void launch_sincosf_kernel(const float *d_x, float *d_s, float *d_c, size_t n) {
   dim3 grid = make_grid(n, blocks.x);
   kernel_sincosf<<<grid, blocks>>>(d_x, d_s, d_c, n);
 }
+
+void launch_expf_kernel(const float *d_x, float *d_e, size_t n) {
+  dim3 blocks(256);
+  dim3 grid = make_grid(n, blocks.x);
+  kernel_expf<<<grid, blocks>>>(d_x, d_e, n);
+}

src/gpu/gpu.cuh

@@ -6,3 +6,4 @@ void launch_sinf_kernel(const float *d_x, float *d_s, size_t n);
 void launch_cosf_kernel(const float *d_x, float *d_c, size_t n);
 void launch_sincosf_kernel(const float *d_x, float *d_s, float *d_c,
                            std::size_t n);
+void launch_expf_kernel(const float *d_x, float *d_e, size_t n);

src/lookup_avx.cpp

@@ -160,7 +160,6 @@ template <std::size_t NR_SAMPLES> struct LookupAVXBackend<NR_SAMPLES>::Impl {
     for (std::size_t i = 0; i < n; ++i) {
       std::size_t idx = static_cast<std::size_t>(x[i] * SCALE) & MASK;
       std::size_t idx_cos = (idx + NR_SAMPLES / 4) & MASK;
-      s[i] = lookup[idx];
       c[i] = lookup[idx_cos];
     }
 #endif

src/lookup_xsimd.cpp

@@ -8,12 +8,20 @@
 template <std::size_t NR_SAMPLES> struct lookup_table {
   static constexpr std::size_t MASK = NR_SAMPLES - 1;
   static constexpr float SCALE = NR_SAMPLES / (2.0f * float(M_PI));
-  lookup_table() : values{} {
+  static constexpr float PI_FRAC = 2.0f * M_PIf32 / NR_SAMPLES;
+  static constexpr float TERM1 = 1.0f;         // 1
+  static constexpr float TERM2 = 0.5f;         // 1/2!
+  static constexpr float TERM3 = 1.0f / 6.0f;  // 1/3!
+  static constexpr float TERM4 = 1.0f / 24.0f; // 1/4!
+
+  constexpr lookup_table() : sin_values{}, cos_values{} {
     for (uint_fast32_t i = 0; i < NR_SAMPLES; i++) {
-      values[i] = sinf(i * (2.0f * float(M_PI) / NR_SAMPLES));
+      sin_values[i] = sinf(i * PI_FRAC);
+      cos_values[i] = cosf(i * PI_FRAC);
     }
   }
-  std::array<float, NR_SAMPLES> values;
+  std::array<float, NR_SAMPLES> cos_values;
+  std::array<float, NR_SAMPLES> sin_values;
 };
 
 template <std::size_t NR_SAMPLES> struct cosf_dispatcher {
@@ -27,26 +35,55 @@ template <std::size_t NR_SAMPLES> struct cosf_dispatcher {
     const uint_fast32_t VS = n - n % VL;
     const uint_fast32_t Q_PI = NR_SAMPLES / 4U;
     const b_type scale = b_type::broadcast(lookup_table_.SCALE);
+    const b_type pi_frac = b_type::broadcast(lookup_table_.PI_FRAC);
     const m_type mask = m_type::broadcast(lookup_table_.MASK);
 
+    const b_type term1 = b_type::broadcast(lookup_table_.TERM1); // 1
+    const b_type term2 = b_type::broadcast(lookup_table_.TERM2); // 1/2!
+    const b_type term3 = b_type::broadcast(lookup_table_.TERM3); // 1/3!
+    const b_type term4 = b_type::broadcast(lookup_table_.TERM4); // 1/4!
     const m_type quarter_pi = m_type::broadcast(Q_PI);
     uint_fast32_t i;
     for (i = 0; i < VS; i += VL) {
       const b_type vx = b_type::load(a + i, Tag());
       const b_type scaled = xsimd::mul(vx, scale);
       m_type idx = xsimd::to_int(scaled);
-      m_type idx_cos = xsimd::add(idx, quarter_pi);
+      const b_type f_idx = xsimd::to_float(idx);
-      idx_cos = xsimd::bitwise_and(idx_cos, mask);
+      idx = xsimd::bitwise_and(idx, mask);
-      const b_type cosv = b_type::gather(lookup_table_.values.data(), idx_cos);
+
+      b_type cosv = b_type::gather(lookup_table_.cos_values.data(), idx);
+      b_type sinv = b_type::gather(lookup_table_.sin_values.data(), idx);
+
+      const b_type dx = xsimd::sub(vx, xsimd::mul(f_idx, pi_frac));
+      const b_type dx2 = xsimd::mul(dx, dx);
+      const b_type dx3 = xsimd::mul(dx2, dx);
+      const b_type dx4 = xsimd::mul(dx2, dx);
+      const b_type t2 = xsimd::mul(dx2, term2);
+      const b_type t3 = xsimd::mul(dx3, term3);
+      const b_type t4 = xsimd::mul(dx4, term3);
+
+      const b_type cosdx = xsimd::add(xsimd::sub(term1, t2), t4);
+
+      const b_type sindx = xsimd::sub(dx, t3);
+
+      cosv = xsimd::sub(xsimd::mul(cosv, cosdx), xsimd::mul(sinv, sindx));
 
       cosv.store(c + i, Tag());
     }
     for (; i < n; i++) {
-      std::size_t idx = static_cast<std::size_t>(a[i] * lookup_table_.SCALE) &
+      std::size_t idx = static_cast<std::size_t>(a[i] * lookup_table_.SCALE);
-                        lookup_table_.MASK;
-      std::size_t idx_cos = (idx + Q_PI) & lookup_table_.MASK;
 
-      c[i] = lookup_table_.values[idx_cos];
+      std::size_t masked = idx & lookup_table_.MASK;
+      const float cosv = lookup_table_.cos_values[masked];
+      const float sinv = lookup_table_.sin_values[masked];
+      const float dx = a[i] - idx * lookup_table_.PI_FRAC;
+      const float dx2 = dx * dx;
+      const float dx3 = dx2 * dx;
+      const float dx4 = dx3 * dx;
+      const float cosdx =
+          1.0f - lookup_table_.TERM2 * dx2 + lookup_table_.TERM4 * dx4;
+      const float sindx = dx - lookup_table_.TERM3 * dx3;
+      c[i] = cosv * cosdx - sinv * sindx;
     }
   }
   lookup_table<NR_SAMPLES> lookup_table_;
@@ -63,23 +100,52 @@ template <std::size_t NR_SAMPLES> struct sinf_dispatcher {
     const uint_fast32_t VS = n - n % VL;
     const uint_fast32_t Q_PI = NR_SAMPLES / 4U;
     const b_type scale = b_type::broadcast(lookup_table_.SCALE);
+    const b_type pi_frac = b_type::broadcast(lookup_table_.PI_FRAC);
     const m_type mask = m_type::broadcast(lookup_table_.MASK);
 
+    const b_type term1 = b_type::broadcast(lookup_table_.TERM1); // 1
+    const b_type term2 = b_type::broadcast(lookup_table_.TERM2); // 1/2!
+    const b_type term3 = b_type::broadcast(lookup_table_.TERM3); // 1/3!
+    const b_type term4 = b_type::broadcast(lookup_table_.TERM4); // 1/4!
     const m_type quarter_pi = m_type::broadcast(Q_PI);
     uint_fast32_t i;
     for (i = 0; i < VS; i += VL) {
       const b_type vx = b_type::load(a + i, Tag());
       const b_type scaled = xsimd::mul(vx, scale);
       m_type idx = xsimd::to_int(scaled);
-      idx = xsimd::bitwise_and(idx, mask);
+      b_type f_idx = xsimd::to_float(idx);
-      const b_type sinv = b_type::gather(lookup_table_.values.data(), idx);
+      const b_type dx = xsimd::sub(vx, xsimd::mul(f_idx, pi_frac));
+      const b_type dx2 = xsimd::mul(dx, dx);
+      const b_type dx3 = xsimd::mul(dx2, dx);
+      const b_type dx4 = xsimd::mul(dx2, dx);
+      const b_type t2 = xsimd::mul(dx2, term2);
+      const b_type t3 = xsimd::mul(dx3, term3);
+      const b_type t4 = xsimd::mul(dx4, term3);
 
+      const b_type cosdx = xsimd::add(xsimd::sub(term1, t2), t4);
+      const b_type sindx = xsimd::sub(dx, t3);
+
+      idx = xsimd::bitwise_and(idx, mask);
+      b_type sinv = b_type::gather(lookup_table_.sin_values.data(), idx);
+      const b_type cosv = b_type::gather(lookup_table_.cos_values.data(), idx);
+
+      sinv = xsimd::add(xsimd::mul(cosv, sindx), xsimd::mul(sinv, cosdx));
       sinv.store(s + i, Tag());
     }
     for (; i < n; i++) {
-      std::size_t idx = static_cast<std::size_t>(a[i] * lookup_table_.SCALE) &
+      std::size_t idx = static_cast<std::size_t>(a[i] * lookup_table_.SCALE);
-                        lookup_table_.MASK;
+      std::size_t masked = idx & lookup_table_.MASK;
-      s[i] = lookup_table_.values[idx];
+      const float cosv = lookup_table_.cos_values[masked];
+      const float sinv = lookup_table_.sin_values[masked];
+      const float dx = a[i] - idx * lookup_table_.PI_FRAC;
+      const float dx2 = dx * dx;
+      const float dx3 = dx2 * dx;
+      const float dx4 = dx3 * dx;
+      const float cosdx =
+          1.0f - lookup_table_.TERM2 * dx2 + lookup_table_.TERM4 * dx4;
+      const float sindx = dx - lookup_table_.TERM3 * dx3;
+
+      s[i] = sinv * cosdx + cosv * sindx;
     }
   }
   lookup_table<NR_SAMPLES> lookup_table_;
@@ -97,6 +163,12 @@ template <std::size_t NR_SAMPLES> struct sin_cosf_dispatcher {
     const uint_fast32_t Q_PI = NR_SAMPLES / 4U;
     const b_type scale = b_type::broadcast(lookup_table_.SCALE);
     const m_type mask = m_type::broadcast(lookup_table_.MASK);
+    const b_type pi_frac = b_type::broadcast(lookup_table_.PI_FRAC);
+
+    const b_type term1 = b_type::broadcast(lookup_table_.TERM1); // 1
+    const b_type term2 = b_type::broadcast(lookup_table_.TERM2); // 1/2!
+    const b_type term3 = b_type::broadcast(lookup_table_.TERM3); // 1/3!
+    const b_type term4 = b_type::broadcast(lookup_table_.TERM4); // 1/4!
 
     const m_type quarter_pi = m_type::broadcast(Q_PI);
     uint_fast32_t i;
@@ -104,21 +176,42 @@ template <std::size_t NR_SAMPLES> struct sin_cosf_dispatcher {
       const b_type vx = b_type::load(a + i, Tag());
       const b_type scaled = xsimd::mul(vx, scale);
       m_type idx = xsimd::to_int(scaled);
-      m_type idx_cos = xsimd::add(idx, quarter_pi);
+      b_type f_idx = xsimd::to_float(idx);
+      const b_type dx = xsimd::sub(vx, xsimd::mul(f_idx, pi_frac));
+      const b_type dx2 = xsimd::mul(dx, dx);
+      const b_type dx3 = xsimd::mul(dx2, dx);
+      const b_type dx4 = xsimd::mul(dx2, dx);
+      const b_type t2 = xsimd::mul(dx2, term2);
+      const b_type t3 = xsimd::mul(dx3, term3);
+      const b_type t4 = xsimd::mul(dx4, term3);
+
       idx = xsimd::bitwise_and(idx, mask);
-      idx_cos = xsimd::bitwise_and(idx_cos, mask);
+      b_type sinv = b_type::gather(lookup_table_.sin_values.data(), idx);
-      const b_type sinv = b_type::gather(lookup_table_.values.data(), idx);
+      b_type cosv = b_type::gather(lookup_table_.cos_values.data(), idx);
-      const b_type cosv = b_type::gather(lookup_table_.values.data(), idx_cos);
+
+      const b_type cosdx = xsimd::add(xsimd::sub(term1, t2), t4);
+      const b_type sindx = xsimd::sub(dx, t3);
+
+      sinv = xsimd::add(xsimd::mul(cosv, sindx), xsimd::mul(sinv, cosdx));
+      cosv = xsimd::sub(xsimd::mul(cosv, cosdx), xsimd::mul(sinv, sindx));
 
       sinv.store(s + i, Tag());
       cosv.store(c + i, Tag());
     }
     for (; i < n; i++) {
-      std::size_t idx = static_cast<std::size_t>(a[i] * lookup_table_.SCALE) &
+      std::size_t idx = static_cast<std::size_t>(a[i] * lookup_table_.SCALE);
-                        lookup_table_.MASK;
+      std::size_t masked = idx & lookup_table_.MASK;
-      std::size_t idx_cos = (idx + Q_PI) & lookup_table_.MASK;
+      const float cosv = lookup_table_.cos_values[masked];
-      s[i] = lookup_table_.values[idx];
+      const float sinv = lookup_table_.sin_values[masked];
-      c[i] = lookup_table_.values[idx_cos];
+      const float dx = a[i] - idx * lookup_table_.PI_FRAC;
+      const float dx2 = dx * dx;
+      const float dx3 = dx2 * dx;
+      const float dx4 = dx3 * dx;
+      const float cosdx =
+          1.0f - lookup_table_.TERM2 * dx2 + lookup_table_.TERM4 * dx4;
+      const float sindx = dx - lookup_table_.TERM3 * dx3;
+      s[i] = sinv * cosdx + cosv * sindx;
+      c[i] = cosv * cosdx - sinv * sindx;
     }
   }
   lookup_table<NR_SAMPLES> lookup_table_;

src/mkl.cpp

@@ -14,3 +14,7 @@ void MKLBackend::compute_sincosf(size_t n, const float *x, float *s,
                                  float *c) const {
   vmsSinCos(static_cast<MKL_INT>(n), x, s, c, VML_HA);
 }
+
+void MKLBackend::compute_expf(size_t n, const float *x, float *e) const {
+  vmsExp(static_cast<MKL_INT>(n), x, e, VML_HA);
+}

src/reference.cpp

@@ -21,3 +21,9 @@ void ReferenceBackend::compute_sincosf(size_t n, const float *x, float *s,
     c[i] = cosf(x[i]);
   }
 }
+
+void ReferenceBackend::compute_expf(size_t n, const float *x, float *e) const {
+  for (size_t i = 0; i < n; ++i) {
+    e[i] = expf(x[i]);
+  }
+}

tests/CMakeLists.txt

@@ -1,5 +1,3 @@
-include(FetchContent)
-
 FetchContent_Declare(
   catch2
   GIT_REPOSITORY https://github.com/catchorg/Catch2.git
@@ -9,31 +7,31 @@ FetchContent_MakeAvailable(catch2)
 # Lookup backend test
 add_executable(test_lookup test_lookup.cpp)
 target_link_libraries(test_lookup PRIVATE trigdx Catch2::Catch2WithMain)
-
-# LookupAVX backend test
-add_executable(test_lookup_avx test_lookup_avx.cpp)
-target_link_libraries(test_lookup_avx PRIVATE trigdx Catch2::Catch2WithMain)
-
-# MKL backend test
-if(USE_MKL)
-  add_executable(test_mkl test_mkl.cpp)
-  target_link_libraries(test_mkl PRIVATE trigdx Catch2::Catch2WithMain)
-endif()
-
-include(CTest)
 add_test(NAME test_lookup COMMAND test_lookup)
 
-if(USE_MKL)
+# LookupAVX backend test
+if(HAVE_AVX)
+  add_executable(test_lookup_avx test_lookup_avx.cpp)
+  target_link_libraries(test_lookup_avx PRIVATE trigdx Catch2::Catch2WithMain)
+  add_test(NAME test_lookup_avx COMMAND test_lookup_avx)
+endif()
+
+# MKL backend test
+if(TRIGDX_USE_MKL)
+  add_executable(test_mkl test_mkl.cpp)
+  target_link_libraries(test_mkl PRIVATE trigdx Catch2::Catch2WithMain)
   add_test(NAME test_mkl COMMAND test_mkl)
 endif()
 
-if(USE_GPU)
+# GPU backend test
+if(TRIGDX_USE_GPU)
   add_executable(test_gpu test_gpu.cpp)
   target_link_libraries(test_gpu PRIVATE trigdx Catch2::Catch2WithMain)
   add_test(NAME test_gpu COMMAND test_gpu)
 endif()
 
-if(USE_XSIMD)
+# XSIMD backend test
+if(TRIGDX_USE_XSIMD)
   add_executable(test_lookup_xsimd test_lookup_xsimd.cpp)
   target_link_libraries(test_lookup_xsimd PRIVATE trigdx Catch2::Catch2WithMain)
   add_test(NAME test_lookup_xsimd COMMAND test_lookup_xsimd)

tests/test_gpu.cpp

@@ -3,8 +3,8 @@
 
 #include "test_utils.hpp"
 
-TEST_CASE("sinf") { test_sinf<GPUBackend>(1e-1f); }
+TEST_CASE("sinf") { test_sinf<GPUBackend>(1e-6f); }
 
-TEST_CASE("cosf") { test_cosf<GPUBackend>(1e-1f); }
+TEST_CASE("cosf") { test_cosf<GPUBackend>(1e-6f); }
 
-TEST_CASE("sincosf") { test_sincosf<GPUBackend>(1e-1f); }
+TEST_CASE("sincosf") { test_sincosf<GPUBackend>(1e-6f); }

tests/test_lookup_xsimd.cpp

@@ -4,16 +4,16 @@
 #include "test_utils.hpp"
 
 TEST_CASE("sincosf") {
-  test_sincosf<LookupXSIMDBackend<16384>>(1e-2f);
+  test_sincosf<LookupXSIMDBackend<16384>>(1e-6f);
-  test_sincosf<LookupXSIMDBackend<32768>>(1e-2f);
+  test_sincosf<LookupXSIMDBackend<32768>>(1e-6f);
 }
 
 TEST_CASE("sinf") {
-  test_sinf<LookupXSIMDBackend<16384>>(1e-2f);
+  test_sinf<LookupXSIMDBackend<16384>>(1e-6f);
-  test_sinf<LookupXSIMDBackend<32768>>(1e-2f);
+  test_sinf<LookupXSIMDBackend<32768>>(1e-6f);
 }
 
 TEST_CASE("cosf") {
-  test_cosf<LookupXSIMDBackend<16384>>(1e-2f);
+  test_cosf<LookupXSIMDBackend<16384>>(1e-6f);
-  test_cosf<LookupXSIMDBackend<32768>>(1e-2f);
+  test_cosf<LookupXSIMDBackend<32768>>(1e-6f);
 }

tests/test_mkl.cpp

@@ -8,3 +8,5 @@ TEST_CASE("sinf") { test_sinf<MKLBackend>(1e-6f); }
 TEST_CASE("cosf") { test_cosf<MKLBackend>(1e-6f); }
 
 TEST_CASE("sincosf") { test_sincosf<MKLBackend>(1e-6f); }
+
+TEST_CASE("expf") { test_expf<MKLBackend>(1e-6f); }

tests/test_utils.hpp

@@ -7,14 +7,17 @@
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
 #include <trigdx/reference.hpp>
 
-const size_t N = 1e7;
+const size_t N = 1234;
+
+void init_x(std::vector<float> &x) {
+  for (size_t i = 0; i < x.size(); ++i) {
+    x[i] = (i % 360) * 0.0174533f; // degrees to radians
+  }
+}
 
 template <typename Backend> inline void test_sinf(float tol) {
   std::vector<float> x(N), s_ref(N), s(N);
-
+  init_x(x);
-  for (size_t i = 0; i < N; ++i) {
-    x[i] = float(i) * 0.01f;
-  }
 
   ReferenceBackend ref;
   Backend backend;
@@ -30,10 +33,7 @@ template <typename Backend> inline void test_sinf(float tol) {
 
 template <typename Backend> inline void test_cosf(float tol) {
   std::vector<float> x(N), c_ref(N), c(N);
-
+  init_x(x);
-  for (size_t i = 0; i < N; ++i) {
-    x[i] = float(i) * 0.01f;
-  }
 
   ReferenceBackend ref;
   Backend backend;
@@ -49,10 +49,7 @@ template <typename Backend> inline void test_cosf(float tol) {
 
 template <typename Backend> inline void test_sincosf(float tol) {
   std::vector<float> x(N), s_ref(N), c_ref(N), s(N), c(N);
-
+  init_x(x);
-  for (size_t i = 0; i < N; ++i) {
-    x[i] = float(i) * 0.01f;
-  }
 
   ReferenceBackend ref;
   Backend backend;
@@ -66,3 +63,19 @@ template <typename Backend> inline void test_sincosf(float tol) {
     REQUIRE_THAT(c[i], Catch::Matchers::WithinAbs(c_ref[i], tol));
   }
 }
+
+template <typename Backend> inline void test_expf(float tol) {
+  std::vector<float> x(N), e_ref(N), e(N);
+  init_x(x);
+
+  ReferenceBackend ref;
+  Backend backend;
+  backend.init(N);
+
+  ref.compute_expf(N, x.data(), e_ref.data());
+  backend.compute_expf(N, x.data(), e.data());
+
+  for (size_t i = 0; i < N; ++i) {
+    REQUIRE_THAT(e[i], Catch::Matchers::WithinAbs(e_ref[i], tol));
+  }
+}