From 15aa3e84d05c2f6727de1867b1109d5ae2719a6a Mon Sep 17 00:00:00 2001 From: Vesa Karvonen Date: Sat, 14 Sep 2019 01:26:39 +0300 Subject: [PATCH] Added 32-bit and 64-bit `prime_less_than_next_pow_2_or_1` functions --- internals/primes/CMakeLists.txt | 1 + internals/primes/program/main.cpp | 76 ++++++++++ internals/{ => tests}/CMakeLists.txt | 0 .../{ => tests}/testing/defaulted_test.cpp | 0 .../testing/finally_function_test.cpp | 0 .../testing/finally_lambda_test.cpp | 0 .../testing/insertion_sort_test.cpp | 0 internals/tests/testing/primes_test.cpp | 39 ++++++ internals/{ => tests}/testing/ranqd1_test.cpp | 0 provides/CMakeLists.txt | 2 +- provides/include/dumpster_v1/primes.hpp | 3 + provides/include/dumpster_v1/synopsis.hpp | 10 ++ provides/library/primes.cpp | 131 ++++++++++++++++++ 13 files changed, 261 insertions(+), 1 deletion(-) create mode 100644 internals/primes/CMakeLists.txt create mode 100644 internals/primes/program/main.cpp rename internals/{ => tests}/CMakeLists.txt (100%) rename internals/{ => tests}/testing/defaulted_test.cpp (100%) rename internals/{ => tests}/testing/finally_function_test.cpp (100%) rename internals/{ => tests}/testing/finally_lambda_test.cpp (100%) rename internals/{ => tests}/testing/insertion_sort_test.cpp (100%) create mode 100644 internals/tests/testing/primes_test.cpp rename internals/{ => tests}/testing/ranqd1_test.cpp (100%) create mode 100644 provides/include/dumpster_v1/primes.hpp create mode 100644 provides/library/primes.cpp diff --git a/internals/primes/CMakeLists.txt b/internals/primes/CMakeLists.txt new file mode 100644 index 0000000..f1f8d03 --- /dev/null +++ b/internals/primes/CMakeLists.txt @@ -0,0 +1 @@ +add_conventional_executable(primes) diff --git a/internals/primes/program/main.cpp b/internals/primes/program/main.cpp new file mode 100644 index 0000000..65d4aad --- /dev/null +++ b/internals/primes/program/main.cpp @@ -0,0 +1,76 @@ +#include +#include +#include +#include +#include + +// + +static uint64_t bit(int i) { return static_cast(1) << i; } + +static uint64_t bits(int n) { return bit(n) * 2 - 1; } + +// + +static bool is_prime(uint64_t n) { + if (n <= 3) + return n > 1; + else if (n % 2 == 0 || n % 3 == 0) + return false; + else + for (uint64_t i = 5, m = static_cast(sqrt(n) + 1); i <= m; i += 6) + if (n % i == 0 || n % (i + 2) == 0) + return false; + return true; +} + +static std::vector primes_less_than_pow_2(int m) { + std::vector results; + results.reserve(64); + results.push_back(1); + for (int i = 1; i < m; ++i) { + uint64_t n = bits(i); + while (!is_prime(n)) + n -= 2; + results.push_back(n); + } + return results; +} + +// + +static void print_table(int n_bits, + int keep_bits, + std::vector table, + uint64_t de_bruijn_sequence) { + struct entry { + int i; + uint64_t n; + uint64_t v; + }; + + std::vector reindexed; + reindexed.resize(n_bits); + + for (int i = 0; i < n_bits; ++i) { + uint64_t n = bits(i); + int j = static_cast((n * de_bruijn_sequence >> (n_bits - keep_bits)) & + bits(keep_bits - 1)); + uint64_t v = table[i]; + reindexed[j] = entry{i, n, v}; + } + + for (const auto &e : reindexed) + fprintf( + stdout, "0x%016" PRIx64 " - 0x%016" PRIx64 ", // %2d\n", e.n, e.v, e.i); +} + +int main() { + auto primes = primes_less_than_pow_2(64); + + print_table(32, 5, primes, 0x07c4acdd); + fprintf(stdout, "\n"); + print_table(64, 6, primes, 0x03f08a4c6acb9dbd); + + return 0; +} diff --git a/internals/CMakeLists.txt b/internals/tests/CMakeLists.txt similarity index 100% rename from internals/CMakeLists.txt rename to internals/tests/CMakeLists.txt diff --git a/internals/testing/defaulted_test.cpp b/internals/tests/testing/defaulted_test.cpp similarity index 100% rename from internals/testing/defaulted_test.cpp rename to internals/tests/testing/defaulted_test.cpp diff --git a/internals/testing/finally_function_test.cpp b/internals/tests/testing/finally_function_test.cpp similarity index 100% rename from internals/testing/finally_function_test.cpp rename to internals/tests/testing/finally_function_test.cpp diff --git a/internals/testing/finally_lambda_test.cpp b/internals/tests/testing/finally_lambda_test.cpp similarity index 100% rename from internals/testing/finally_lambda_test.cpp rename to internals/tests/testing/finally_lambda_test.cpp diff --git a/internals/testing/insertion_sort_test.cpp b/internals/tests/testing/insertion_sort_test.cpp similarity index 100% rename from internals/testing/insertion_sort_test.cpp rename to internals/tests/testing/insertion_sort_test.cpp diff --git a/internals/tests/testing/primes_test.cpp b/internals/tests/testing/primes_test.cpp new file mode 100644 index 0000000..3505921 --- /dev/null +++ b/internals/tests/testing/primes_test.cpp @@ -0,0 +1,39 @@ +#include "dumpster_v1/primes.hpp" + +#include "testing_v1/test.hpp" + +using namespace testing_v1; + +#define VERIFY_IS_PRIME_OR_1 0 + +#if VERIFY_IS_PRIME_OR_1 +#include + +static bool is_prime(uint64_t n) { + if (n <= 3) + return n > 1; + else if (n % 2 == 0 || n % 3 == 0) + return false; + else + for (uint64_t i = 5, m = static_cast(sqrt(n) + 1); i <= m; i += 6) + if (n % i == 0 || n % (i + 2) == 0) + return false; + return true; +} +#endif + +template static void test() { + for (int i = 0; i < static_cast(sizeof(T) * 8); ++i) { + T n = static_cast(1) << i; + T p = dumpster_v1::prime_less_than_next_pow_2_or_1(n); +#if VERIFY_IS_PRIME_OR_1 + verify(1 == p || is_prime(p)); +#endif + verify(n <= p && p <= n * 2 - 1); + } +} + +auto primes_test = test([]() { + test(); + test(); +}); diff --git a/internals/testing/ranqd1_test.cpp b/internals/tests/testing/ranqd1_test.cpp similarity index 100% rename from internals/testing/ranqd1_test.cpp rename to internals/tests/testing/ranqd1_test.cpp diff --git a/provides/CMakeLists.txt b/provides/CMakeLists.txt index 0af6943..cde9369 100644 --- a/provides/CMakeLists.txt +++ b/provides/CMakeLists.txt @@ -1,2 +1,2 @@ add_conventional_library(dumpster_v1) -target_link_libraries(dumpster_v1 INTERFACE polyfill_v1) +target_link_libraries(dumpster_v1 PUBLIC polyfill_v1) diff --git a/provides/include/dumpster_v1/primes.hpp b/provides/include/dumpster_v1/primes.hpp new file mode 100644 index 0000000..c2e7e64 --- /dev/null +++ b/provides/include/dumpster_v1/primes.hpp @@ -0,0 +1,3 @@ +#pragma once + +#include "dumpster_v1/synopsis.hpp" diff --git a/provides/include/dumpster_v1/synopsis.hpp b/provides/include/dumpster_v1/synopsis.hpp index 742d639..e80cdfb 100644 --- a/provides/include/dumpster_v1/synopsis.hpp +++ b/provides/include/dumpster_v1/synopsis.hpp @@ -64,4 +64,14 @@ void insertion_sort(RandomIt values, Less less, IsSentinel is_sentinel); /// The `ranqd1` generator from Numerical Recipes in C, 2nd Edition. uint32_t ranqd1(uint32_t seed); +// primes.hpp ================================================================== + +/// Returns the largest prime that is less than the given value rounded to the +/// next power of 2 or 1. +uint32_t prime_less_than_next_pow_2_or_1(uint32_t x); + +/// Returns the largest prime that is less than the given value rounded +/// to the next power of 2 or 1. +uint64_t prime_less_than_next_pow_2_or_1(uint64_t x); + } // namespace dumpster_v1 diff --git a/provides/library/primes.cpp b/provides/library/primes.cpp new file mode 100644 index 0000000..a21acb4 --- /dev/null +++ b/provides/library/primes.cpp @@ -0,0 +1,131 @@ +#include "dumpster_v1/primes.hpp" + +// Using de Bruijn Sequences to Index a 1 in a Computer Word +// by Leiserson, Prokop, and Randall + +uint32_t dumpster_v1::prime_less_than_next_pow_2_or_1(uint32_t x) { + alignas(32) static const uint8_t delta[] = { + 0x0000000000000001 - 0x0000000000000001, // 0 + 0x00000000000003ff - 0x00000000000003fd, // 9 + 0x0000000000000003 - 0x0000000000000003, // 1 + 0x00000000000007ff - 0x00000000000007f7, // 10 + 0x0000000000003fff - 0x0000000000003ffd, // 13 + 0x00000000003fffff - 0x00000000003ffffd, // 21 + 0x0000000000000007 - 0x0000000000000007, // 2 + 0x000000003fffffff - 0x000000003fffffdd, // 29 + 0x0000000000000fff - 0x0000000000000ffd, // 11 + 0x0000000000007fff - 0x0000000000007fed, // 14 + 0x000000000001ffff - 0x000000000001ffff, // 16 + 0x000000000007ffff - 0x000000000007ffff, // 18 + 0x00000000007fffff - 0x00000000007ffff1, // 22 + 0x0000000003ffffff - 0x0000000003fffffb, // 25 + 0x000000000000000f - 0x000000000000000d, // 3 + 0x000000007fffffff - 0x000000007fffffff, // 30 + 0x00000000000001ff - 0x00000000000001fd, // 8 + 0x0000000000001fff - 0x0000000000001fff, // 12 + 0x00000000001fffff - 0x00000000001ffff7, // 20 + 0x000000001fffffff - 0x000000001ffffffd, // 28 + 0x000000000000ffff - 0x000000000000fff1, // 15 + 0x000000000003ffff - 0x000000000003fffb, // 17 + 0x0000000001ffffff - 0x0000000001ffffd9, // 24 + 0x00000000000000ff - 0x00000000000000fb, // 7 + 0x00000000000fffff - 0x00000000000ffffd, // 19 + 0x000000000fffffff - 0x000000000fffffc7, // 27 + 0x0000000000ffffff - 0x0000000000fffffd, // 23 + 0x000000000000007f - 0x000000000000007f, // 6 + 0x0000000007ffffff - 0x0000000007ffffd9, // 26 + 0x000000000000003f - 0x000000000000003d, // 5 + 0x000000000000001f - 0x000000000000001f, // 4 + 0x00000000ffffffff - 0x00000000fffffffb, // 31 + }; + + static_assert(sizeof(delta) == 32); + + x |= x >> 1; + x |= x >> 2; + x |= x >> 4; + x |= x >> 8; + x |= x >> 16; + + return x - delta[x * 0x07c4acdd >> (32 - 5)]; +} + +uint64_t dumpster_v1::prime_less_than_next_pow_2_or_1(uint64_t x) { + alignas(64) static const uint8_t delta[] = { + 0x0000000000000001 - 0x0000000000000001, // 0 + 0x0000000000000fff - 0x0000000000000ffd, // 11 + 0x0000000000000003 - 0x0000000000000003, // 1 + 0x0000000000001fff - 0x0000000000001fff, // 12 + 0x000000000001ffff - 0x000000000001ffff, // 16 + 0x000000003fffffff - 0x000000003fffffdd, // 29 + 0x0000000000000007 - 0x0000000000000007, // 2 + 0x0000000000003fff - 0x0000000000003ffd, // 13 + 0x00000000007fffff - 0x00000000007ffff1, // 22 + 0x000000000003ffff - 0x000000000003fffb, // 17 + 0x000003ffffffffff - 0x000003fffffffff5, // 41 + 0x0000000003ffffff - 0x0000000003fffffb, // 25 + 0x000000007fffffff - 0x000000007fffffff, // 30 + 0x0001ffffffffffff - 0x0001ffffffffffaf, // 48 + 0x000000000000000f - 0x000000000000000d, // 3 + 0x3fffffffffffffff - 0x3fffffffffffffc7, // 61 + 0x0000000000007fff - 0x0000000000007fed, // 14 + 0x00000000001fffff - 0x00000000001ffff7, // 20 + 0x0000000000ffffff - 0x0000000000fffffd, // 23 + 0x000000000007ffff - 0x000000000007ffff, // 18 + 0x00000007ffffffff - 0x00000007ffffffe1, // 34 + 0x0000001fffffffff - 0x0000001fffffffe7, // 36 + 0x000007ffffffffff - 0x000007ffffffffc7, // 42 + 0x0000000007ffffff - 0x0000000007ffffd9, // 26 + 0x0000007fffffffff - 0x0000007ffffffff9, // 38 + 0x00000000ffffffff - 0x00000000fffffffb, // 31 + 0x003fffffffffffff - 0x003fffffffffffdf, // 53 + 0x00001fffffffffff - 0x00001fffffffffc9, // 44 + 0x0003ffffffffffff - 0x0003ffffffffffe5, // 49 + 0x01ffffffffffffff - 0x01fffffffffffff3, // 56 + 0x000000000000001f - 0x000000000000001f, // 4 + 0x7fffffffffffffff - 0x7fffffffffffffe7, // 62 + 0x00000000000007ff - 0x00000000000007f7, // 10 + 0x000000000000ffff - 0x000000000000fff1, // 15 + 0x000000001fffffff - 0x000000001ffffffd, // 28 + 0x00000000003fffff - 0x00000000003ffffd, // 21 + 0x000001ffffffffff - 0x000001ffffffffeb, // 40 + 0x0000000001ffffff - 0x0000000001ffffd9, // 24 + 0x0000ffffffffffff - 0x0000ffffffffffc5, // 47 + 0x1fffffffffffffff - 0x1fffffffffffffff, // 60 + 0x00000000000fffff - 0x00000000000ffffd, // 19 + 0x00000003ffffffff - 0x00000003ffffffd7, // 33 + 0x0000000fffffffff - 0x0000000ffffffffb, // 35 + 0x0000003fffffffff - 0x0000003fffffffd3, // 37 + 0x001fffffffffffff - 0x001fffffffffff91, // 52 + 0x00000fffffffffff - 0x00000fffffffffef, // 43 + 0x00ffffffffffffff - 0x00fffffffffffffb, // 55 + 0x00000000000003ff - 0x00000000000003fd, // 9 + 0x000000000fffffff - 0x000000000fffffc7, // 27 + 0x000000ffffffffff - 0x000000ffffffffa9, // 39 + 0x00007fffffffffff - 0x00007fffffffff8d, // 46 + 0x0fffffffffffffff - 0x0fffffffffffffa3, // 59 + 0x00000001ffffffff - 0x00000001fffffff7, // 32 + 0x000fffffffffffff - 0x000fffffffffffd1, // 51 + 0x007fffffffffffff - 0x007fffffffffffc9, // 54 + 0x00000000000001ff - 0x00000000000001fd, // 8 + 0x00003fffffffffff - 0x00003fffffffffeb, // 45 + 0x07ffffffffffffff - 0x07ffffffffffffc9, // 58 + 0x0007ffffffffffff - 0x0007ffffffffff7f, // 50 + 0x00000000000000ff - 0x00000000000000fb, // 7 + 0x03ffffffffffffff - 0x03ffffffffffffe5, // 57 + 0x000000000000007f - 0x000000000000007f, // 6 + 0x000000000000003f - 0x000000000000003d, // 5 + 0xffffffffffffffff - 0xffffffffffffffc5, // 63 + }; + + static_assert(sizeof(delta) == 64); + + x |= x >> 1; + x |= x >> 2; + x |= x >> 4; + x |= x >> 8; + x |= x >> 16; + x |= x >> 32; + + return x - delta[x * 0x03f08a4c6acb9dbd >> (64 - 6)]; +}