summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile20
-rw-r--r--segmented_sieve.cc41
-rw-r--r--sieve.cc58
-rw-r--r--sieve.hh11
-rw-r--r--sieve_utils.cc10
-rw-r--r--test.cc97
-rw-r--r--threaded_sieve.cc85
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;
+}
generated by cgit v1.2.1 (git 2.18.0) at 2024-11-21 12:49:49 +0000