Initial commit

.pre-commit-config.yaml

@@ -0,0 +1,10 @@
+repos:
+  - repo: https://github.com/pre-commit/mirrors-clang-format
+    rev: v17.0.6
+    hooks:
+      - id: clang-format
+  - repo: https://github.com/cheshirekow/cmake-format-precommit
+    rev: v0.6.13
+    hooks:
+      - id: cmake-format
+      - id: cmake-lint

.vscode/settings.json

@@ -0,0 +1,77 @@
+{
+  "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

@@ -0,0 +1,13 @@
+cmake_minimum_required(VERSION 3.16)
+project(trigdx LANGUAGES CXX)
+
+set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+option(USE_MKL "Enable Intel MKL backend" OFF)
+
+include_directories(${PROJECT_SOURCE_DIR}/include)
+
+add_subdirectory(src)
+add_subdirectory(tests)
+add_subdirectory(benchmarks)

benchmarks/CMakeLists.txt

@@ -0,0 +1,10 @@
+add_executable(benchmark_reference benchmark_reference.cpp)
+target_link_libraries(benchmark_reference PRIVATE trigdx)
+
+add_executable(benchmark_lookup benchmark_lookup.cpp)
+target_link_libraries(benchmark_lookup PRIVATE trigdx)
+
+if(USE_MKL)
+  add_executable(benchmark_mkl benchmark_mkl.cpp)
+  target_link_libraries(benchmark_mkl PRIVATE trigdx)
+endif()

benchmarks/benchmark_lookup.cpp

@@ -0,0 +1,13 @@
+#include <trigdx/lookup.hpp>
+
+#include "benchmark_utils.hpp"
+
+int main() {
+  benchmark_sinf<LookupBackend<16384>>();
+  benchmark_cosf<LookupBackend<16384>>();
+  benchmark_sincosf<LookupBackend<16384>>();
+
+  benchmark_sinf<LookupBackend<32768>>();
+  benchmark_cosf<LookupBackend<32768>>();
+  benchmark_sincosf<LookupBackend<32768>>();
+}

benchmarks/benchmark_mkl.cpp

@@ -0,0 +1,9 @@
+#include <trigdx/mkl.hpp>
+
+#include "benchmark_utils.hpp"
+
+int main() {
+  benchmark_sinf<MKLBackend>();
+  benchmark_cosf<MKLBackend>();
+  benchmark_sincosf<MKLBackend>();
+}

benchmarks/benchmark_reference.cpp

@@ -0,0 +1,9 @@
+#include <trigdx/reference.hpp>
+
+#include "benchmark_utils.hpp"
+
+int main() {
+  benchmark_sinf<ReferenceBackend>();
+  benchmark_cosf<ReferenceBackend>();
+  benchmark_sincosf<ReferenceBackend>();
+}

benchmarks/benchmark_utils.hpp

@@ -0,0 +1,76 @@
+#pragma once
+
+#include <chrono>
+#include <iomanip>
+#include <iostream>
+#include <vector>
+
+const size_t N = 1e7;
+
+inline void report(const std::string &name, double sec, double throughput) {
+  std::ios state(nullptr);
+  state.copyfmt(std::cout);
+  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() {
+  std::vector<float> x(N), s(N);
+
+  for (size_t i = 0; i < N; ++i)
+    x[i] = (i % 360) * 0.0174533f; // degrees to radians
+
+  Backend backend;
+  backend.init();
+
+  auto start = std::chrono::high_resolution_clock::now();
+  backend.compute_sinf(N, x.data(), s.data());
+  auto end = std::chrono::high_resolution_clock::now();
+
+  double sec = std::chrono::duration<double>(end - start).count();
+  double throughput = N / sec / 1e6;
+
+  report("sinf", sec, throughput);
+}
+
+template <typename Backend> inline void benchmark_cosf() {
+  std::vector<float> x(N), c(N);
+
+  for (size_t i = 0; i < N; ++i)
+    x[i] = (i % 360) * 0.0174533f; // degrees to radians
+
+  Backend backend;
+  backend.init();
+
+  auto start = std::chrono::high_resolution_clock::now();
+  backend.compute_cosf(N, x.data(), c.data());
+  auto end = std::chrono::high_resolution_clock::now();
+
+  double sec = std::chrono::duration<double>(end - start).count();
+  double throughput = N / sec / 1e6;
+
+  report("cosf", sec, throughput);
+}
+
+template <typename Backend> inline void benchmark_sincosf() {
+  std::vector<float> x(N), s(N), c(N);
+
+  for (size_t i = 0; i < N; ++i)
+    x[i] = (i % 360) * 0.0174533f; // degrees to radians
+
+  Backend backend;
+  backend.init();
+
+  auto start = std::chrono::high_resolution_clock::now();
+  backend.compute_sincosf(N, x.data(), s.data(), c.data());
+  auto end = std::chrono::high_resolution_clock::now();
+
+  double sec = std::chrono::duration<double>(end - start).count();
+  double throughput = N / sec / 1e6;
+
+  report("sincosf", sec, throughput);
+}

include/trigdx/interface.hpp

@@ -0,0 +1,22 @@
+#pragma once
+
+#include <cstddef>
+
+// Base interface for all math backends
+class Backend {
+public:
+  virtual ~Backend() = default;
+
+  // Optional initialization
+  virtual void init() {}
+
+  // Compute sine for n elements
+  virtual void compute_sinf(size_t n, const float *x, float *s) const = 0;
+
+  // Compute cosine for n elements
+  virtual void compute_cosf(size_t n, const float *x, float *c) const = 0;
+
+  // Compute sine and cosine for n elements
+  virtual void compute_sincosf(size_t n, const float *x, float *s,
+                               float *c) const = 0;
+};

include/trigdx/lookup.hpp

@@ -0,0 +1,22 @@
+#pragma once
+
+#include <cmath>
+#include <vector>
+
+#include "interface.hpp"
+
+template <size_t NR_SAMPLES> class LookupBackend : public Backend {
+public:
+  void init() override;
+  void compute_sinf(size_t n, const float *x, float *s) const override;
+  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;
+
+private:
+  std::vector<float> lookup;
+  static constexpr size_t MASK = NR_SAMPLES - 1;
+  static constexpr float SCALE = NR_SAMPLES / (2.0f * float(M_PI));
+};
+
+#include "lookup.tpp" // include implementation

include/trigdx/lookup.tpp

@@ -0,0 +1,44 @@
+#pragma once
+
+#include <cmath>
+
+template <size_t NR_SAMPLES>
+void LookupBackend<NR_SAMPLES>::init() {
+    lookup.resize(NR_SAMPLES);
+    for (size_t i = 0; i < NR_SAMPLES; ++i)
+        lookup[i] = std::sinf(i * (2.0f * float(M_PI) / NR_SAMPLES));
+}
+
+template <size_t NR_SAMPLES>
+void LookupBackend<NR_SAMPLES>::compute_sinf(size_t n,
+                                               const float* x,
+                                               float* s) const {
+    for (size_t i = 0; i < n; ++i) {
+        size_t idx = static_cast<size_t>(x[i] * SCALE) & MASK;
+        s[i] = lookup[idx];
+    }
+}
+
+template <size_t NR_SAMPLES>
+void LookupBackend<NR_SAMPLES>::compute_cosf(size_t n,
+                                               const float* x,
+                                               float* c) const {
+    for (size_t i = 0; i < n; ++i) {
+        size_t idx = static_cast<size_t>(x[i] * SCALE) & MASK;
+        size_t idx_cos = (idx + NR_SAMPLES / 4) & MASK;
+        c[i] = lookup[idx_cos];
+    }
+}
+
+template <size_t NR_SAMPLES>
+void LookupBackend<NR_SAMPLES>::compute_sincosf(size_t n,
+                                                  const float* x,
+                                                  float* s,
+                                                  float* c) const {
+    for (size_t i = 0; i < n; ++i) {
+        size_t idx = static_cast<size_t>(x[i] * SCALE) & MASK;
+        size_t idx_cos = (idx + NR_SAMPLES / 4) & MASK;
+        s[i] = lookup[idx];
+        c[i] = lookup[idx_cos];
+    }
+}

include/trigdx/mkl.hpp

@@ -0,0 +1,13 @@
+#pragma once
+
+#include "interface.hpp"
+
+class MKLBackend : public Backend {
+public:
+  void compute_sinf(size_t n, const float *x, float *s) const override;
+
+  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;
+};

include/trigdx/reference.hpp

@@ -0,0 +1,13 @@
+#pragma once
+
+#include "interface.hpp"
+
+class ReferenceBackend : public Backend {
+public:
+  void compute_sinf(size_t n, const float *x, float *s) const override;
+
+  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;
+};

src/CMakeLists.txt

@@ -0,0 +1,9 @@
+add_library(trigdx reference.cpp)
+
+target_include_directories(trigdx PUBLIC ${PROJECT_SOURCE_DIR}/include)
+
+if(USE_MKL)
+  find_package(MKL REQUIRED)
+  target_sources(trigdx PRIVATE mkl.cpp)
+  target_link_libraries(trigdx PRIVATE MKL::MKL)
+endif()

src/mkl.cpp

@@ -0,0 +1,16 @@
+#include <mkl_vml.h>
+
+#include "trigdx/mkl.hpp"
+
+void MKLBackend::compute_sinf(size_t n, const float *x, float *s) const {
+  vmsSin(static_cast<MKL_INT>(n), x, s, VML_HA);
+}
+
+void MKLBackend::compute_cosf(size_t n, const float *x, float *c) const {
+  vmsCos(static_cast<MKL_INT>(n), x, c, VML_HA);
+}
+
+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);
+}

src/reference.cpp

@@ -0,0 +1,23 @@
+#include <cmath>
+
+#include "trigdx/reference.hpp"
+
+void ReferenceBackend::compute_sinf(size_t n, const float *x, float *s) const {
+  for (size_t i = 0; i < n; ++i) {
+    s[i] = std::sinf(x[i]);
+  }
+}
+
+void ReferenceBackend::compute_cosf(size_t n, const float *x, float *c) const {
+  for (size_t i = 0; i < n; ++i) {
+    c[i] = std::cosf(x[i]);
+  }
+}
+
+void ReferenceBackend::compute_sincosf(size_t n, const float *x, float *s,
+                                       float *c) const {
+  for (size_t i = 0; i < n; ++i) {
+    s[i] = std::sinf(x[i]);
+    c[i] = std::cosf(x[i]);
+  }
+}

tests/CMakeLists.txt

@@ -0,0 +1,24 @@
+include(FetchContent)
+
+FetchContent_Declare(
+  catch2
+  GIT_REPOSITORY https://github.com/catchorg/Catch2.git
+  GIT_TAG v3.9.0)
+FetchContent_MakeAvailable(catch2)
+
+# Lookup backend test
+add_executable(test_lookup test_lookup.cpp)
+target_link_libraries(test_lookup 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)
+  add_test(NAME test_mkl COMMAND test_mkl)
+endif()

tests/test_lookup.cpp

@@ -0,0 +1,19 @@
+#include <catch2/catch_test_macros.hpp>
+#include <trigdx/lookup.hpp>
+
+#include "test_utils.hpp"
+
+TEST_CASE("sinf") {
+  test_sinf<LookupBackend<16384>>(1e-2f);
+  test_sinf<LookupBackend<32768>>(1e-2f);
+}
+
+TEST_CASE("cosf") {
+  test_cosf<LookupBackend<16384>>(1e-2f);
+  test_cosf<LookupBackend<32768>>(1e-2f);
+}
+
+TEST_CASE("sincosf") {
+  test_sincosf<LookupBackend<16384>>(1e-2f);
+  test_sincosf<LookupBackend<32768>>(1e-2f);
+}

tests/test_mkl.cpp

@@ -0,0 +1,10 @@
+#include <catch2/catch_test_macros.hpp>
+#include <trigdx/mkl.hpp>
+
+#include "test_utils.hpp"
+
+TEST_CASE("sinf") { test_sinf<MKLBackend>(1e-6f); }
+
+TEST_CASE("cosf") { test_cosf<MKLBackend>(1e-6f); }
+
+TEST_CASE("sincosf") { test_sincosf<MKLBackend>(1e-6f); }

tests/test_utils.hpp

@@ -0,0 +1,68 @@
+#pragma once
+
+#include <algorithm>
+#include <cmath>
+#include <vector>
+
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <trigdx/reference.hpp>
+
+const size_t N = 1e7;
+
+template <typename Backend> inline void test_sinf(float tol) {
+  std::vector<float> x(N), s_ref(N), s(N);
+
+  for (size_t i = 0; i < N; ++i) {
+    x[i] = float(i) * 0.01f;
+  }
+
+  ReferenceBackend ref;
+  Backend backend;
+  backend.init();
+
+  ref.compute_sinf(N, x.data(), s_ref.data());
+  backend.compute_sinf(N, x.data(), s.data());
+
+  for (size_t i = 0; i < N; ++i) {
+    REQUIRE_THAT(s[i], Catch::Matchers::WithinAbs(s_ref[i], tol));
+  }
+}
+
+template <typename Backend> inline void test_cosf(float tol) {
+  std::vector<float> x(N), c_ref(N), c(N);
+
+  for (size_t i = 0; i < N; ++i) {
+    x[i] = float(i) * 0.01f;
+  }
+
+  ReferenceBackend ref;
+  Backend backend;
+  backend.init();
+
+  ref.compute_cosf(N, x.data(), c_ref.data());
+  backend.compute_cosf(N, x.data(), c.data());
+
+  for (size_t i = 0; i < N; ++i) {
+    REQUIRE_THAT(c[i], Catch::Matchers::WithinAbs(c_ref[i], 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);
+
+  for (size_t i = 0; i < N; ++i) {
+    x[i] = float(i) * 0.01f;
+  }
+
+  ReferenceBackend ref;
+  Backend backend;
+  backend.init();
+
+  ref.compute_sincosf(N, x.data(), s_ref.data(), c_ref.data());
+  backend.compute_sincosf(N, x.data(), s.data(), c.data());
+
+  for (size_t i = 0; i < N; ++i) {
+    REQUIRE_THAT(s[i], Catch::Matchers::WithinAbs(s_ref[i], tol));
+    REQUIRE_THAT(c[i], Catch::Matchers::WithinAbs(c_ref[i], tol));
+  }
+}