Merge pull request #9 from barrett-ruth/sieve

sieve of eratosthenes

include/bmath/bmath.hh

@@ -2,5 +2,6 @@
 #define BMATH_BMATH_HH
 
 #include "bmath/mint.hh"
+#include "bmath/sieve.hh"
 
 #endif

include/bmath/mint.hh

@@ -1,14 +1,14 @@
 #ifndef BMATH_MINT_HH
 #define BMATH_MINT_HH
 
-namespace bmath {
-
 #include <format>
 #include <iterator>
 #include <ostream>
 #include <stdexcept>
 #include <string>
 
+namespace bmath {
+
 inline constexpr uint64_t DEFAULT_MOD = 1'000'000'007;
 
 template <std::integral T, T M = static_cast<T>(DEFAULT_MOD)>

include/bmath/sieve.hh

@@ -0,0 +1,62 @@
+#ifndef BMATH_SIEVE_HH
+#define BMATH_SIEVE_HH
+
+#include <bitset>
+#include <cstdint>
+#include <format>
+#include <stdexcept>
+
+namespace bmath {
+
+template <typename Derived, size_t Limit>
+  requires(Limit > 0)
+class Sieve {
+ public:
+  constexpr explicit Sieve() noexcept { static_cast<Derived *>(this)->build(); }
+
+ protected:
+  std::bitset<Limit + 1> sieve;
+};
+
+template <size_t Limit>
+  requires(Limit > 0)
+class Eratosthenes : public Sieve<Eratosthenes<Limit>, Limit> {
+ public:
+  constexpr void build() noexcept {
+    this->sieve.set();
+    this->sieve[0] = false;
+    this->sieve[1] = false;
+
+    for (size_t x = 2; x * x <= Limit; ++x) {
+      if (!this->sieve[x]) {
+        continue;
+      }
+
+      for (size_t current = x + x; current <= Limit; current += x) {
+        this->sieve[current] = false;
+      }
+    }
+  }
+
+  [[nodiscard]] constexpr bool operator[](size_t const number) const {
+    if consteval {
+      if (number > Limit) {
+        throw std::out_of_range(
+            std::format("cannot determine primality of {} > size of sieve {}",
+                        number, Limit));
+      }
+    } else {
+      if (number > Limit) {
+        throw std::out_of_range(
+            std::format("cannot determine primality of {} > size of sieve {}",
+                        number, Limit));
+      }
+    }
+
+    return this->sieve[number];
+  };
+};
+
+}  // namespace bmath
+
+#endif

readme.md

@@ -4,6 +4,7 @@ header-only c++ 23 math library
 
 ## todo
 
-[] sieves
-[] factorization
-[x] `std::unsigned_integral` as modulus type
+- sieves (pritchard)
+- factorization (pollard-rho, rabin-miller)
+- profiling
+- saturating arithmetic operations

tests/CMakeLists.txt

@@ -2,16 +2,19 @@ include(FetchContent)
 FetchContent_Declare(
   googletest
   GIT_REPOSITORY https://github.com/google/googletest.git
-  GIT_TAG        v1.14.0
-)
+  GIT_TAG v1.14.0)
 set(gtest_force_shared_crt
     ON
     CACHE BOOL "" FORCE)
 FetchContent_MakeAvailable(googletest)
 
-add_executable(test_bmath test_bmath.cc)
+file(GLOB TEST_SOURCES "*.cc")
 
-target_link_libraries(test_bmath PRIVATE bmath::bmath gtest_main)
+foreach(TEST_SOURCE ${TEST_SOURCES})
+  get_filename_component(TEST_NAME ${TEST_SOURCE} NAME_WE)
+  add_executable(${TEST_NAME} ${TEST_SOURCE})
+  target_link_libraries(${TEST_NAME} PRIVATE bmath::bmath gtest_main)
 
-include(GoogleTest)
-gtest_discover_tests(test_bmath)
+  include(GoogleTest)
+  gtest_discover_tests(${TEST_NAME})
+endforeach()

tests/test_sieve.cc

@@ -0,0 +1,35 @@
+#include <gtest/gtest.h>
+
+#include <cstdint>
+#include <format>
+
+#include "../include/bmath/bmath.hh"
+
+using namespace bmath;
+
+bool naive_prime(size_t n) {
+  if (n < 2) {
+    return false;
+  }
+  if (n == 2) {
+    return true;
+  }
+  if (n % 2 == 0) {
+    return false;
+  }
+  for (size_t i = 3; i * i <= n; i += 2) {
+    if (n % i == 0) return false;
+  }
+  return true;
+}
+
+TEST(SieveTest, BasicConstruction) { Eratosthenes<100> sieve; }
+
+TEST(SieveTest, ValidateAgainstNaive) {
+  Eratosthenes<1000> sieve;
+
+  for (size_t i = 0; i <= 1000; ++i) {
+    EXPECT_EQ(sieve[i], naive_prime(i))
+        << "Sieve disagrees with naive checker for " << i;
+  }
+}