From 57027fd076504b219f3e886e78535191f6abbb66 Mon Sep 17 00:00:00 2001 From: aqua Date: Sat, 18 Mar 2023 07:17:27 +0200 Subject: Initial commit Sieve of Eratosthenes - simple implementation and optimized versions - Segmented sieve with pthreaded version --- Makefile | 20 +++++++++++ segmented_sieve.cc | 41 +++++++++++++++++++++++ sieve.cc | 58 ++++++++++++++++++++++++++++++++ sieve.hh | 11 +++++++ sieve_utils.cc | 10 ++++++ test.cc | 97 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ threaded_sieve.cc | 85 +++++++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 322 insertions(+) create mode 100644 Makefile create mode 100644 segmented_sieve.cc create mode 100644 sieve.cc create mode 100644 sieve.hh create mode 100644 sieve_utils.cc create mode 100644 test.cc create mode 100644 threaded_sieve.cc 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 +#include + +[[nodiscard]] std::vector segmentedSieve(std::size_t size) { + const size_t s = sqrt(size); + const auto primes = sieve3(s); + std::vector r{primes}; + std::vector 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 +#include + +[[nodiscard]] std::vector sieve1(std::size_t size) { + std::vector n(size, false); + std::vector 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 sieve2(std::size_t size) { + std::vector n(size, false); + std::vector 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 sieve3(std::size_t size) { + std::vector n(size, false); + std::vector 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 + +[[nodiscard]] bool isPrime(std::size_t x); + +[[nodiscard]] std::vector sieve1(std::size_t size); +[[nodiscard]] std::vector sieve2(std::size_t size); +[[nodiscard]] std::vector sieve3(std::size_t size); +[[nodiscard]] std::vector segmentedSieve(std::size_t size); +[[nodiscard]] std::vector 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 + +[[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 +#include +#include +#include + +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::size_t)> SieveFunction; + +class SieveFunctionFixture + : public ::testing::TestWithParam> { +public: + void SetUp() override { + sieve = std::get(GetParam()); + sieve_name = std::get(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(end - begin); + const auto nt = + std::chrono::duration_cast(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(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(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 +#include +#include +#include + +struct BlockArgs { + BlockArgs(const std::vector &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 primes; // in + std::vector r; // out +}; + +static void *worker(void *args) { + auto *a = reinterpret_cast(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 pthreadSieve(std::size_t size) { + const int s = sqrt(size); + const int block_sz = size / thread_count; + + const auto primes = sieve3(s); + std::vector 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; +} -- cgit v1.2.1