From f40c44d5dd5fcb399908ac9bddc89683bff06501 Mon Sep 17 00:00:00 2001
From: Bram Veenboer <bram.veenboer@gmail.com>
Date: Tue, 12 Aug 2025 14:55:28 +0200
Subject: [PATCH 1/4] Add Python interface

---
 CMakeLists.txt        | 10 +++++-
 python/CMakeLists.txt |  8 +++++
 python/bindings.cpp   | 80 +++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 97 insertions(+), 1 deletion(-)
 create mode 100644 python/CMakeLists.txt
 create mode 100644 python/bindings.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index e807f5e..ef1df5e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -3,18 +3,22 @@ project(trigdx LANGUAGES CXX)
 
 set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 
 option(TRIGDX_USE_MKL "Enable Intel MKL backend" OFF)
 option(TRIGDX_USE_GPU "Enable GPU 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)
 
 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)
+if(TRIGDX_BUILD_TESTS
+   OR TRIGDX_BUILD_BENCHMARKS
+   OR TRIGDX_BUILD_PYTHON)
   include(FetchContent)
 endif()
 
@@ -31,3 +35,7 @@ endif()
 if(TRIGDX_BUILD_BENCHMARKS)
   add_subdirectory(benchmarks)
 endif()
+
+if(TRIGDX_BUILD_PYTHON)
+  add_subdirectory(python)
+endif()
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
new file mode 100644
index 0000000..81e1494
--- /dev/null
+++ b/python/CMakeLists.txt
@@ -0,0 +1,8 @@
+FetchContent_Declare(
+  pybind11
+  GIT_REPOSITORY https://github.com/pybind/pybind11.git
+  GIT_TAG v3.0.0)
+FetchContent_MakeAvailable(pybind11)
+
+pybind11_add_module(pytrigdx bindings.cpp)
+target_link_libraries(pytrigdx PRIVATE trigdx)
diff --git a/python/bindings.cpp b/python/bindings.cpp
new file mode 100644
index 0000000..5c57756
--- /dev/null
+++ b/python/bindings.cpp
@@ -0,0 +1,80 @@
+#include <memory>
+
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <trigdx/trigdx.hpp>
+
+namespace py = pybind11;
+
+py::array_t<float> compute_sinf_py(
+    const Backend &backend,
+    py::array_t<float, py::array::c_style | py::array::forcecast> x) {
+  ssize_t n = x.shape(0);
+  auto x_ptr = x.data();
+
+  py::array_t<float> s(n);
+  auto s_ptr = s.mutable_data();
+
+  backend.compute_sinf(static_cast<size_t>(n), x_ptr, s_ptr);
+
+  return s;
+}
+
+py::array_t<float> compute_cosf_py(
+    const Backend &backend,
+    py::array_t<float, py::array::c_style | py::array::forcecast> x) {
+  ssize_t n = x.shape(0);
+  auto x_ptr = x.data();
+
+  py::array_t<float> c(n);
+  auto c_ptr = c.mutable_data();
+
+  backend.compute_cosf(static_cast<size_t>(n), x_ptr, c_ptr);
+
+  return c;
+}
+
+std::tuple<py::array_t<float>, py::array_t<float>> compute_sincosf_py(
+    const Backend &backend,
+    py::array_t<float, py::array::c_style | py::array::forcecast> x) {
+  ssize_t n = x.shape(0);
+  auto x_ptr = x.data();
+
+  py::array_t<float> s(n);
+  py::array_t<float> c(n);
+
+  backend.compute_sincosf(static_cast<size_t>(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_sinf_py)
+      .def("compute_cosf", &compute_cosf_py)
+      .def("compute_sincosf", &compute_sincosf_py);
+}
+
+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
+}
\ No newline at end of file

From 97692cface1911070a0a8ee94fcedf5d433b7f98 Mon Sep 17 00:00:00 2001
From: Bram Veenboer <bram.veenboer@gmail.com>
Date: Thu, 14 Aug 2025 11:02:30 +0200
Subject: [PATCH 2/4] Use Type template for helper functions

---
 python/bindings.cpp | 43 +++++++++++++++++++++++--------------------
 1 file changed, 23 insertions(+), 20 deletions(-)

diff --git a/python/bindings.cpp b/python/bindings.cpp
index 5c57756..1988e98 100644
--- a/python/bindings.cpp
+++ b/python/bindings.cpp
@@ -6,42 +6,45 @@
 
 namespace py = pybind11;
 
-py::array_t<float> compute_sinf_py(
-    const Backend &backend,
-    py::array_t<float, py::array::c_style | py::array::forcecast> x) {
+template <typename T>
+py::array_t<T>
+compute_sin(const Backend &backend,
+            py::array_t<T, py::array::c_style | py::array::forcecast> x) {
   ssize_t n = x.shape(0);
-  auto x_ptr = x.data();
+  const T *x_ptr = x.data();
 
   py::array_t<float> s(n);
-  auto s_ptr = s.mutable_data();
+  T *s_ptr = s.mutable_data();
 
   backend.compute_sinf(static_cast<size_t>(n), x_ptr, s_ptr);
 
   return s;
 }
 
-py::array_t<float> compute_cosf_py(
-    const Backend &backend,
-    py::array_t<float, py::array::c_style | py::array::forcecast> x) {
+template <typename T>
+py::array_t<T>
+compute_cos(const Backend &backend,
+            py::array_t<T, py::array::c_style | py::array::forcecast> x) {
   ssize_t n = x.shape(0);
-  auto x_ptr = x.data();
+  const T *x_ptr = x.data();
 
-  py::array_t<float> c(n);
-  auto c_ptr = c.mutable_data();
+  py::array_t<T> c(n);
+  T *c_ptr = c.mutable_data();
 
   backend.compute_cosf(static_cast<size_t>(n), x_ptr, c_ptr);
 
   return c;
 }
 
-std::tuple<py::array_t<float>, py::array_t<float>> compute_sincosf_py(
-    const Backend &backend,
-    py::array_t<float, py::array::c_style | py::array::forcecast> x) {
+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) {
   ssize_t n = x.shape(0);
-  auto x_ptr = x.data();
+  const T *x_ptr = x.data();
 
-  py::array_t<float> s(n);
-  py::array_t<float> c(n);
+  py::array_t<T> s(n);
+  py::array_t<T> c(n);
 
   backend.compute_sincosf(static_cast<size_t>(n), x_ptr, s.mutable_data(),
                           c.mutable_data());
@@ -53,9 +56,9 @@ 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_sinf_py)
-      .def("compute_cosf", &compute_cosf_py)
-      .def("compute_sincosf", &compute_sincosf_py);
+      .def("compute_sinf", &compute_sin<float>)
+      .def("compute_cosf", &compute_cos<float>)
+      .def("compute_sincosf", &compute_sincos<float>);
 }
 
 PYBIND11_MODULE(pytrigdx, m) {

From 79dc7b4285ea274649a7a6d0295d8b2183fc84ce Mon Sep 17 00:00:00 2001
From: Bram Veenboer <bram.veenboer@gmail.com>
Date: Thu, 14 Aug 2025 11:05:31 +0200
Subject: [PATCH 3/4] Add dimension checks

---
 python/bindings.cpp | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/python/bindings.cpp b/python/bindings.cpp
index 1988e98..ba63955 100644
--- a/python/bindings.cpp
+++ b/python/bindings.cpp
@@ -11,6 +11,10 @@ py::array_t<T>
 compute_sin(const Backend &backend,
             py::array_t<T, py::array::c_style | py::array::forcecast> x) {
   ssize_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);
@@ -26,6 +30,10 @@ py::array_t<T>
 compute_cos(const Backend &backend,
             py::array_t<T, py::array::c_style | py::array::forcecast> x) {
   ssize_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);
@@ -41,6 +49,10 @@ 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) {
   ssize_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);

From b3467840f965a9792465c2a456a98387ff94dc27 Mon Sep 17 00:00:00 2001
From: Bram Veenboer <bram.veenboer@gmail.com>
Date: Thu, 14 Aug 2025 11:07:45 +0200
Subject: [PATCH 4/4] Use const size_t for n

---
 python/bindings.cpp | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/python/bindings.cpp b/python/bindings.cpp
index ba63955..c8153df 100644
--- a/python/bindings.cpp
+++ b/python/bindings.cpp
@@ -10,7 +10,7 @@ template <typename T>
 py::array_t<T>
 compute_sin(const Backend &backend,
             py::array_t<T, py::array::c_style | py::array::forcecast> x) {
-  ssize_t n = x.shape(0);
+  const size_t n = x.shape(0);
   if (x.ndim() != 1) {
     throw py::value_error("Input array must be 1-dimensional");
   }
@@ -20,7 +20,7 @@ compute_sin(const Backend &backend,
   py::array_t<float> s(n);
   T *s_ptr = s.mutable_data();
 
-  backend.compute_sinf(static_cast<size_t>(n), x_ptr, s_ptr);
+  backend.compute_sinf(n, x_ptr, s_ptr);
 
   return s;
 }
@@ -29,7 +29,7 @@ template <typename T>
 py::array_t<T>
 compute_cos(const Backend &backend,
             py::array_t<T, py::array::c_style | py::array::forcecast> x) {
-  ssize_t n = x.shape(0);
+  const size_t n = x.shape(0);
   if (x.ndim() != 1) {
     throw py::value_error("Input array must be 1-dimensional");
   }
@@ -39,7 +39,7 @@ compute_cos(const Backend &backend,
   py::array_t<T> c(n);
   T *c_ptr = c.mutable_data();
 
-  backend.compute_cosf(static_cast<size_t>(n), x_ptr, c_ptr);
+  backend.compute_cosf(n, x_ptr, c_ptr);
 
   return c;
 }
@@ -48,7 +48,7 @@ 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) {
-  ssize_t n = x.shape(0);
+  const size_t n = x.shape(0);
   if (x.ndim() != 1) {
     throw py::value_error("Input array must be 1-dimensional");
   }
@@ -58,8 +58,7 @@ compute_sincos(const Backend &backend,
   py::array_t<T> s(n);
   py::array_t<T> c(n);
 
-  backend.compute_sincosf(static_cast<size_t>(n), x_ptr, s.mutable_data(),
-                          c.mutable_data());
+  backend.compute_sincosf(n, x_ptr, s.mutable_data(), c.mutable_data());
 
   return std::make_tuple(s, c);
 }