Initial commitHEADmaster

Sieve of Eratosthenes
- simple implementation and optimized versions
- Segmented sieve with pthreaded version

7 files changed, 322 insertions, 0 deletions

diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..af3901d
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,20 @@
+
+CXXFLAGS     += -Wall -Wextra -Wpedantic -g -Og --std=c++20
+TESTCXXFLAGS := $(shell pkg-config --cflags gtest gtest_main)
+TESTLDFLAGS  := $(shell pkg-config --libs gtest gtest_main) -lpthread
+
+SRCS := sieve_utils.cc sieve.cc segmented_sieve.cc threaded_sieve.cc
+OBJS := $(SRCS:.cc=.o)
+
+test: $(OBJS) test.o
+	$(CXX) $(TESTLDFLAGS) $^ -o $@
+
+clean:
+	rm -f *.o test
+
+%.o: %.cc
+	$(CXX) $(CXXFLAGS) -c $^ -o $@
+
+test.o: test.cc
+	$(CXX) $(CXXFLAGS) -c $^ -o $@
+
diff --git a/segmented_sieve.cc b/segmented_sieve.cc
new file mode 100644
index 0000000..b5974bb
--- /dev/null
+++ b/segmented_sieve.cc
@@ -0,0 +1,41 @@
+#include "sieve.hh"
+#include <cmath>
+#include <cstring>
+
+[[nodiscard]] std::vector<int> segmentedSieve(std::size_t size) {
+  const size_t s = sqrt(size);
+  const auto primes = sieve3(s);
+  std::vector<int> r{primes};
+  std::vector<int> r2;
+
+  size_t l = s;
+  size_t h = l + s;
+  unsigned char *n = (unsigned char *)malloc(s);
+
+  while (l < size) {
+
+    memset(n, false, s);
+
+    for (const auto prime : primes) {
+      // find the first number divisible by prime
+      std::size_t f = (l / prime) * prime;
+      if (f < l)
+        f += prime;
+
+      for (std::size_t i = f; i < h; i += prime)
+        n[i - l] = true;
+    }
+    for (std::size_t i = 0; i < s; ++i)
+      if (n[i] == false)
+        r2.emplace_back(i + l);
+
+    l += s;
+    h += s;
+  }
+
+  free(n);
+
+  r.reserve(r.size() + r2.size());
+  r.insert(r.end(), r2.begin(), r2.end());
+  return r;
+}
diff --git a/sieve.cc b/sieve.cc
new file mode 100644
index 0000000..4f7c84f
--- /dev/null
+++ b/sieve.cc
@@ -0,0 +1,58 @@
+#include "sieve.hh"
+#include <cmath>
+#include <cstring>
+
+[[nodiscard]] std::vector<int> sieve1(std::size_t size) {
+  std::vector<bool> n(size, false);
+  std::vector<int> r;
+
+  for (std::size_t i = 2; i < n.size(); ++i) {
+    if (n[i] == false) {
+      r.emplace_back(i);
+      for (std::size_t j = i + i; j < n.size(); j += i)
+        n[j] = true;
+    }
+  }
+
+  return r;
+}
+
+[[nodiscard]] std::vector<int> sieve2(std::size_t size) {
+  std::vector<bool> n(size, false);
+  std::vector<int> r{2};
+
+  for (std::size_t i = 3; i < n.size(); i += 2) {
+    if (n[i] == false) {
+      r.emplace_back(i);
+      for (std::size_t j = i + i; j < n.size(); j += i)
+        n[j] = true;
+    }
+  }
+
+  return r;
+}
+
+[[nodiscard]] std::vector<int> sieve3(std::size_t size) {
+  std::vector<bool> n(size, false);
+  std::vector<int> r{2};
+
+  std::size_t s = sqrt(n.size());
+
+  for (std::size_t i = 4; i <= s; i += 2)
+    n[i] = true;
+
+  for (std::size_t i = 3; i <= s; i += 2) {
+    if (n[i] == false) {
+      r.emplace_back(i);
+      for (std::size_t j = i * i; j < n.size(); j += i)
+        n[j] = true;
+    }
+  }
+
+  s += s % 2 ? 2 : 1;
+  for (std::size_t i = s; i < n.size(); i += 2)
+    if (n[i] == false)
+      r.emplace_back(i);
+
+  return r;
+}
diff --git a/sieve.hh b/sieve.hh
new file mode 100644
index 0000000..c3e6c14
--- /dev/null
+++ b/sieve.hh
@@ -0,0 +1,11 @@
+#pragma once
+
+#include <vector>
+
+[[nodiscard]] bool isPrime(std::size_t x);
+
+[[nodiscard]] std::vector<int> sieve1(std::size_t size);
+[[nodiscard]] std::vector<int> sieve2(std::size_t size);
+[[nodiscard]] std::vector<int> sieve3(std::size_t size);
+[[nodiscard]] std::vector<int> segmentedSieve(std::size_t size);
+[[nodiscard]] std::vector<int> pthreadSieve(std::size_t size);
diff --git a/sieve_utils.cc b/sieve_utils.cc
new file mode 100644
index 0000000..51bdcf1
--- /dev/null
+++ b/sieve_utils.cc
@@ -0,0 +1,10 @@
+#include "sieve.hh"
+#include <cmath>
+
+[[nodiscard]] bool isPrime(std::size_t x) {
+  for (size_t i = 2; i <= sqrt(x); ++i)
+    if (x % i == 0) {
+      return false;
+    }
+  return true;
+}
diff --git a/test.cc b/test.cc
new file mode 100644
index 0000000..ab31ba6
--- /dev/null
+++ b/test.cc
@@ -0,0 +1,97 @@
+#include "sieve.hh"
+#include <chrono>
+#include <functional>
+#include <gtest/gtest.h>
+#include <iostream>
+
+using ::testing::Combine;
+using ::testing::Values;
+
+TEST(Eratosthenes, isPrime) {
+  EXPECT_FALSE(isPrime(4));
+  EXPECT_FALSE(isPrime(9409));
+  EXPECT_TRUE(isPrime(9973));
+  EXPECT_FALSE(isPrime(53033));
+}
+
+constexpr auto sz = 1024 * 1024;
+static const auto n = sieve1(sz);
+
+TEST(Eratosthenes, validate) {
+  EXPECT_EQ(n[0], 2);
+
+  for (size_t i = 0; i < n.size(); ++i)
+    EXPECT_TRUE(isPrime(n[i])) << n[i];
+}
+
+typedef std::function<std::vector<int>(std::size_t)> SieveFunction;
+
+class SieveFunctionFixture
+    : public ::testing::TestWithParam<std::tuple<SieveFunction, std::string>> {
+public:
+  void SetUp() override {
+    sieve = std::get<SieveFunction>(GetParam());
+    sieve_name = std::get<std::string>(GetParam());
+  }
+
+protected:
+  SieveFunction sieve;
+  std::string sieve_name;
+};
+
+TEST_P(SieveFunctionFixture, validate) {
+  const auto begin = std::chrono::system_clock::now();
+  const auto r = sieve(sz);
+  const auto end = std::chrono::system_clock::now();
+
+  const auto ut =
+      std::chrono::duration_cast<std::chrono::microseconds>(end - begin);
+  const auto nt =
+      std::chrono::duration_cast<std::chrono::nanoseconds>(end - begin);
+
+  std::cout << sieve_name << '\n'
+            << "run time: " << ut << '\n'
+            << "prime count: " << r.size() << '\n'
+            << "time per prime: " << nt / r.size() << std::endl;
+
+  EXPECT_EQ(n.size(), r.size());
+  for (size_t i = 0; i < std::min(n.size(), r.size()); ++i) {
+    EXPECT_EQ(n[i], r[i]) << "difference at index " << i;
+  }
+}
+
+INSTANTIATE_TEST_SUITE_P(
+    Eratosthenes, SieveFunctionFixture,
+    Values(std::make_tuple(&sieve1, "sieve1"),
+           std::make_tuple(&sieve2, "sieve2"),
+           std::make_tuple(&sieve3, "sieve3"),
+           std::make_tuple(&segmentedSieve, "segmented sieve"),
+           std::make_tuple(&pthreadSieve, "threaded sieve")));
+
+TEST(benchmark, segmentedSieve) {
+  const auto begin = std::chrono::system_clock::now();
+  const auto r = segmentedSieve(1024 * sz);
+  const auto end = std::chrono::system_clock::now();
+
+  const auto t =
+      std::chrono::duration_cast<std::chrono::microseconds>(end - begin);
+  std::cout << r.size() << " values in " << t << std::endl;
+
+  for (size_t i = 0; i < std::min(n.size(), r.size()); ++i) {
+    EXPECT_EQ(n[i], r[i]) << "difference at index " << i;
+  }
+}
+
+TEST(benchmark, threadedSieve) {
+  const auto begin = std::chrono::system_clock::now();
+  const auto r = pthreadSieve(1024 * sz);
+  const auto end = std::chrono::system_clock::now();
+
+  const auto t =
+      std::chrono::duration_cast<std::chrono::microseconds>(end - begin);
+  std::cout << r.size() << " values in " << t << std::endl;
+
+  for (size_t i = 0; i < std::min(n.size(), r.size()); ++i) {
+    EXPECT_EQ(n[i], r[i]) << "difference at index " << i;
+  }
+}
diff --git a/threaded_sieve.cc b/threaded_sieve.cc
new file mode 100644
index 0000000..3e683c1
--- /dev/null
+++ b/threaded_sieve.cc
@@ -0,0 +1,85 @@
+#include "sieve.hh"
+#include <cmath>
+#include <cstring>
+#include <iostream>
+#include <pthread.h>
+
+struct BlockArgs {
+  BlockArgs(const std::vector<int> &p, size_t block_sz)
+      : s(block_sz), primes(p) {
+    n = (unsigned char *)malloc(s);
+  }
+  ~BlockArgs() { free(n); }
+
+  std::size_t l, h;              // in
+  const size_t s;                // in
+  unsigned char *n;              // in
+  const std::vector<int> primes; // in
+  std::vector<int> r;            // out
+};
+
+static void *worker(void *args) {
+  auto *a = reinterpret_cast<BlockArgs *>(args);
+
+  for (std::size_t l = a->l; l < a->h; l += a->s) {
+    memset(a->n, false, a->s);
+
+    for (const auto prime : a->primes) {
+      // find the first number divisible by prime
+      std::size_t f = (l / prime) * prime;
+      if (f < l)
+        f += prime;
+
+      for (std::size_t i = f; i < l + a->s; i += prime)
+        a->n[i - l] = true;
+    }
+
+    for (size_t i = 0; i < a->s; ++i)
+      if (a->n[i] == false)
+        a->r.emplace_back(i + l);
+  }
+
+  return nullptr;
+}
+
+constexpr auto thread_count = 16;
+
+[[nodiscard]] std::vector<int> pthreadSieve(std::size_t size) {
+  const int s = sqrt(size);
+  const int block_sz = size / thread_count;
+
+  const auto primes = sieve3(s);
+  std::vector<int> r{primes};
+
+  pthread_t threads[thread_count];
+  BlockArgs *args[thread_count]{nullptr};
+  for (int i = 0; i < thread_count; ++i) {
+    args[i] = new BlockArgs(primes, s);
+  }
+
+  {
+    int i = 0; // thread index
+    for (std::size_t l = s; l < size; l += block_sz) {
+
+      // start thread i
+      args[i]->l = l;
+      args[i]->h = std::min(l + block_sz, size);
+      pthread_create(&threads[i], NULL, worker, args[i]);
+      // worker(args[i]);
+
+      ++i;
+    }
+
+    // wait out the threads
+    for (int j = 0; j < i; ++j) {
+      pthread_join(threads[j], NULL);
+      r.reserve(r.size() + args[j]->r.size());
+      r.insert(r.end(), args[j]->r.begin(), args[j]->r.end());
+    }
+  }
+  // cleanup
+  for (int i = 0; i < thread_count; ++i)
+    delete args[i];
+
+  return r;
+}