Merge pull request #337 from chiphogg/chiphogg/prime-wheel

Use wheel factorization for prime numbers

src/core/include/units/bits/prime.h

@@ -0,0 +1,165 @@
+// The MIT License (MIT)
+//
+// Copyright (c) 2018 Mateusz Pusz
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+
+#pragma once
+
+#include <array>
+#include <cassert>
+#include <cstddef>
+#include <numeric>
+#include <optional>
+#include <tuple>
+
+namespace units::detail {
+
+constexpr bool is_prime_by_trial_division(std::size_t n)
+{
+  for (std::size_t f = 2; f * f <= n; f += 1 + (f % 2)) {
+    if (n % f == 0) {
+      return false;
+    }
+  }
+  return true;
+}
+
+// Return the first factor of n, as long as it is either k or n.
+//
+// Precondition: no integer smaller than k evenly divides n.
+constexpr std::optional<std::size_t> first_factor_maybe(std::size_t n, std::size_t k)
+{
+  if (n % k == 0) {
+    return k;
+  }
+  if (k * k > n) {
+    return n;
+  }
+  return std::nullopt;
+}
+
+template<std::size_t N>
+constexpr std::array<std::size_t, N> first_n_primes()
+{
+  std::array<std::size_t, N> primes;
+  primes[0] = 2;
+  for (std::size_t i = 1; i < N; ++i) {
+    primes[i] = primes[i - 1] + 1;
+    while (!is_prime_by_trial_division(primes[i])) {
+      ++primes[i];
+    }
+  }
+  return primes;
+}
+
+template<std::size_t N, typename Callable>
+constexpr void call_for_coprimes_up_to(std::size_t n, const std::array<std::size_t, N>& basis, Callable&& call)
+{
+  for (std::size_t i = 0u; i < n; ++i) {
+    if (std::apply([&i](auto... primes) { return ((i % primes != 0) && ...); }, basis)) {
+      call(i);
+    }
+  }
+}
+
+template<std::size_t N>
+constexpr std::size_t num_coprimes_up_to(std::size_t n, const std::array<std::size_t, N>& basis)
+{
+  std::size_t count = 0u;
+  call_for_coprimes_up_to(n, basis, [&count](auto) { ++count; });
+  return count;
+}
+
+template<std::size_t ResultSize, std::size_t N>
+constexpr auto coprimes_up_to(std::size_t n, const std::array<std::size_t, N>& basis)
+{
+  std::array<std::size_t, ResultSize> coprimes;
+  std::size_t i = 0u;
+
+  call_for_coprimes_up_to(n, basis, [&coprimes, &i](std::size_t cp) { coprimes[i++] = cp; });
+
+  return coprimes;
+}
+
+template<std::size_t N>
+constexpr std::size_t product(const std::array<std::size_t, N>& values)
+{
+  std::size_t product = 1;
+  for (const auto& v : values) {
+    product *= v;
+  }
+  return product;
+}
+
+// A configurable instantiation of the "wheel factorization" algorithm [1] for prime numbers.
+//
+// Instantiate with N to use a "basis" of the first N prime numbers.  Higher values of N use fewer trial divisions, at
+// the cost of additional space.  The amount of space consumed is roughly the total number of numbers that are a) less
+// than the _product_ of the basis elements (first N primes), and b) coprime with every element of the basis.  This
+// means it grows rapidly with N.  Consider this approximate chart:
+//
+//   N | Num coprimes | Trial divisions needed
+//   --+--------------+-----------------------
+//   1 |            1 |                 50.0 %
+//   2 |            2 |                 33.3 %
+//   3 |            8 |                 26.7 %
+//   4 |           48 |                 22.9 %
+//   5 |          480 |                 20.8 %
+//
+// Note the diminishing returns, and the rapidly escalating costs.  Consider this behaviour when choosing the value of N
+// most appropriate for your needs.
+//
+// [1] https://en.wikipedia.org/wiki/Wheel_factorization
+template<std::size_t BasisSize>
+struct wheel_factorizer {
+  static constexpr auto basis = first_n_primes<BasisSize>();
+  static constexpr std::size_t wheel_size = product(basis);
+  static constexpr auto coprimes_in_first_wheel =
+    coprimes_up_to<num_coprimes_up_to(wheel_size, basis)>(wheel_size, basis);
+
+  static constexpr std::size_t find_first_factor(std::size_t n)
+  {
+    for (const auto& p : basis) {
+      if (const auto k = first_factor_maybe(n, p)) {
+        return *k;
+      }
+    }
+
+    for (auto it = std::next(std::begin(coprimes_in_first_wheel)); it != std::end(coprimes_in_first_wheel); ++it) {
+      if (const auto k = first_factor_maybe(n, *it)) {
+        return *k;
+      }
+    }
+
+    for (std::size_t wheel = wheel_size; wheel < n; wheel += wheel_size) {
+      for (const auto& p : coprimes_in_first_wheel) {
+        if (const auto k = first_factor_maybe(n, wheel + p)) {
+          return *k;
+        }
+      }
+    }
+
+    return n;
+  }
+
+  static constexpr bool is_prime(std::size_t n) { return (n > 1) && find_first_factor(n) == n; }
+};
+
+}  // namespace units::detail

src/core/include/units/magnitude.h

@@ -23,12 +23,19 @@
 #pragma once
 
 #include <units/bits/external/hacks.h>
+#include <units/bits/prime.h>
 #include <units/ratio.h>
+#include <concepts>
 #include <cstdint>
 #include <numbers>
+#include <optional>
 #include <stdexcept>
 
 namespace units {
+namespace detail {
+// Higher numbers use fewer trial divisions, at the price of more storage space.
+using factorizer = wheel_factorizer<4>;
+}  // namespace detail
 
 /**
  * @brief  Any type which can be used as a basis vector in a BasePower.
@@ -244,17 +251,6 @@ constexpr auto pow(BasePower auto bp, ratio p)
 // A variety of implementation detail helpers.
 namespace detail {
 
-// Find the smallest prime factor of `n`.
-constexpr std::intmax_t find_first_factor(std::intmax_t n)
-{
-  for (std::intmax_t f = 2; f * f <= n; f += 1 + (f % 2)) {
-    if (n % f == 0) {
-      return f;
-    }
-  }
-  return n;
-}
-
 // The exponent of `factor` in the prime factorization of `n`.
 constexpr std::intmax_t multiplicity(std::intmax_t factor, std::intmax_t n)
 {
@@ -278,7 +274,7 @@ constexpr std::intmax_t remove_power(std::intmax_t base, std::intmax_t pow, std:
 }
 
 // A way to check whether a number is prime at compile time.
-constexpr bool is_prime(std::intmax_t n) { return (n > 1) && (find_first_factor(n) == n); }
+constexpr bool is_prime(std::intmax_t n) { return (n >= 0) && factorizer::is_prime(static_cast<std::size_t>(n)); }
 
 constexpr bool is_valid_base_power(const BasePower auto& bp)
 {
@@ -287,6 +283,20 @@ constexpr bool is_valid_base_power(const BasePower auto& bp)
   }
 
   if constexpr (std::is_same_v<decltype(bp.get_base()), std::intmax_t>) {
+    // Some prime numbers are so big, that we can't check their primality without exhausting limits on constexpr steps
+    // and/or iterations.  We can still _perform_ the factorization for these by using the `known_first_factor`
+    // workaround.  However, we can't _check_ that they are prime, because this workaround depends on the input being
+    // usable in a constexpr expression.  This is true for `prime_factorization` (below), where the input `N` is a
+    // template parameter, but is not true for our case, where the input `bp.get_base()` is a function parameter.  (See
+    // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2019/p1045r1.html for some background reading on this
+    // distinction.)
+    //
+    // In our case: we simply give up on excluding every possible ill-formed base power, and settle for catching the
+    // most likely and common mistakes.
+    if (const bool too_big_to_check = (bp.get_base() > 1'000'000'000)) {
+      return true;
+    }
+
     return is_prime(bp.get_base());
   } else {
     return bp.get_base() > 0;
@@ -473,12 +483,30 @@ constexpr auto operator/(Magnitude auto l, Magnitude auto r) { return l * pow<-1
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 // `as_magnitude()` implementation.
 
+// Sometimes we need to give the compiler a "shortcut" when factorizing large numbers (specifically, numbers whose
+// _first factor_ is very large).  If we don't, we can run into limits on the number of constexpr steps or iterations.
+//
+// To provide the first factor for a given number, specialize this variable template.
+//
+// WARNING:  The program behaviour will be undefined if you provide a wrong answer, so check your math!
+template<std::intmax_t N>
+inline constexpr std::optional<std::intmax_t> known_first_factor = std::nullopt;
+
 namespace detail {
 // Helper to perform prime factorization at compile time.
 template<std::intmax_t N>
   requires(N > 0)
 struct prime_factorization {
-  static constexpr std::intmax_t first_base = find_first_factor(N);
+  static constexpr std::intmax_t get_or_compute_first_factor()
+  {
+    if constexpr (known_first_factor<N>.has_value()) {
+      return known_first_factor<N>.value();
+    } else {
+      return static_cast<std::intmax_t>(factorizer::find_first_factor(N));
+    }
+  }
+
+  static constexpr std::intmax_t first_base = get_or_compute_first_factor();
   static constexpr std::intmax_t first_power = multiplicity(first_base, N);
   static constexpr std::intmax_t remainder = remove_power(first_base, first_power, N);
 

test/unit_test/runtime/magnitude_test.cpp

@@ -25,7 +25,13 @@
 #include <units/ratio.h>
 #include <type_traits>
 
-namespace units {
+using namespace units;
+using namespace units::detail;
+
+template<>
+inline constexpr std::optional<std::intmax_t> units::known_first_factor<9223372036854775783> = 9223372036854775783;
+
+namespace {
 
 // A set of non-standard bases for testing purposes.
 struct noninteger_base {
@@ -149,6 +155,17 @@ TEST_CASE("make_ratio performs prime factorization correctly")
     // The failure was due to a prime factor which is larger than 2^31.
     as_magnitude<ratio(16'605'390'666'050, 10'000'000'000'000)>();
   }
+
+  SECTION("Can bypass computing primes by providing known_first_factor<N>")
+  {
+    // Sometimes, even wheel factorization isn't enough to handle the compilers' limits on constexpr steps and/or
+    // iterations.  To work around these cases, we can explicitly provide the correct answer directly to the compiler.
+    //
+    // In this case, we test that we can represent the largest prime that fits in a signed 64-bit int.  The reason this
+    // test can pass is that we have provided the answer, by specializing the `known_first_factor` variable template
+    // above in this file.
+    as_magnitude<9'223'372'036'854'775'783>();
+  }
 }
 
 TEST_CASE("magnitude converts to numerical value")
@@ -334,7 +351,8 @@ TEST_CASE("can distinguish integral, rational, and irrational magnitudes")
   }
 }
 
-namespace detail {
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Detail function tests below.
 
 TEST_CASE("int_power computes integer powers")
 {
@@ -358,16 +376,6 @@ TEST_CASE("int_power computes integer powers")
 
 TEST_CASE("Prime helper functions")
 {
-  SECTION("find_first_factor()")
-  {
-    CHECK(find_first_factor(1) == 1);
-    CHECK(find_first_factor(2) == 2);
-    CHECK(find_first_factor(4) == 2);
-    CHECK(find_first_factor(6) == 2);
-    CHECK(find_first_factor(15) == 3);
-    CHECK(find_first_factor(17) == 17);
-  }
-
   SECTION("multiplicity")
   {
     CHECK(multiplicity(2, 8) == 3);
@@ -514,6 +522,4 @@ TEST_CASE("strictly_increasing")
   }
 }
 
-}  // namespace detail
+}  // namespace
-
-}  // namespace units

test/unit_test/static/CMakeLists.txt

@@ -52,6 +52,7 @@ add_library(unit_tests_static
     kind_test.cpp
     math_test.cpp
     point_origin_test.cpp
+    prime_test.cpp
     ratio_test.cpp
     references_test.cpp
     si_test.cpp

test/unit_test/static/prime_test.cpp

@@ -0,0 +1,76 @@
+// The MIT License (MIT)
+//
+// Copyright (c) 2018 Mateusz Pusz
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+
+#include <units/bits/prime.h>
+#include <type_traits>
+#include <utility>
+
+using namespace units::detail;
+
+namespace {
+
+template<std::size_t BasisSize, std::size_t... Is>
+constexpr bool check_primes(std::index_sequence<Is...>)
+{
+  return ((Is < 2 || wheel_factorizer<BasisSize>::is_prime(Is) == is_prime_by_trial_division(Is)) && ...);
+}
+
+static_assert(check_primes<2>(std::make_index_sequence<122>{}));
+
+// This is the smallest number that can catch the bug where we use only _prime_ numbers in the first wheel, rather than
+// numbers which are _coprime to the basis_.
+//
+// The basis for N = 4 is {2, 3, 5, 7}, so the wheel size is 210.  11 * 11 = 121 is within the first wheel.  It is
+// coprime with every element of the basis, but it is _not_ prime.  If we keep only prime numbers, then we will neglect
+// using numbers of the form (210 * n + 121) as trial divisors, which is a problem if any are prime.  For n = 1, we have
+// a divisor of (210 + 121 = 331), which happens to be prime but will not be used.  Thus, (331 * 331 = 109561) is a
+// composite number which could wrongly appear prime if we skip over 331.
+static_assert(wheel_factorizer<4>::is_prime(109561) == is_prime_by_trial_division(109561));
+
+static_assert(wheel_factorizer<1>::coprimes_in_first_wheel.size() == 1);
+static_assert(wheel_factorizer<2>::coprimes_in_first_wheel.size() == 2);
+static_assert(wheel_factorizer<3>::coprimes_in_first_wheel.size() == 8);
+static_assert(wheel_factorizer<4>::coprimes_in_first_wheel.size() == 48);
+static_assert(wheel_factorizer<5>::coprimes_in_first_wheel.size() == 480);
+
+static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[0] == 1);
+static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[1] == 7);
+static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[2] == 11);
+static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[3] == 13);
+static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[4] == 17);
+static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[5] == 19);
+static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[6] == 23);
+static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[7] == 29);
+
+static_assert(!wheel_factorizer<1>::is_prime(0));
+static_assert(!wheel_factorizer<1>::is_prime(1));
+static_assert(wheel_factorizer<1>::is_prime(2));
+
+static_assert(!wheel_factorizer<2>::is_prime(0));
+static_assert(!wheel_factorizer<2>::is_prime(1));
+static_assert(wheel_factorizer<2>::is_prime(2));
+
+static_assert(!wheel_factorizer<3>::is_prime(0));
+static_assert(!wheel_factorizer<3>::is_prime(1));
+static_assert(wheel_factorizer<3>::is_prime(2));
+
+}  // namespace