Replace old factoring with Baillie-PSW

src/core/include/mp-units/ext/prime.h

@@ -401,30 +401,6 @@ struct LucasSequenceElement {
   return strong_lucas_probable_prime(n);
 }
 
-[[nodiscard]] consteval bool is_prime_by_trial_division(std::uintmax_t n)
-{
-  for (std::uintmax_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.
-[[nodiscard]] constexpr std::optional<std::uintmax_t> first_factor_maybe(std::uintmax_t n, std::uintmax_t k)
-{
-  if (n % k == 0) {
-    return k;
-  }
-  if (k * k > n) {
-    return n;
-  }
-  return std::nullopt;
-}
-
 template<std::size_t N>
 [[nodiscard]] consteval std::array<std::uintmax_t, N> first_n_primes()
 {
@@ -432,42 +408,13 @@ template<std::size_t N>
   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])) {
+    while (!baillie_psw_probable_prime(primes[i])) {
       ++primes[i];
     }
   }
   return primes;
 }
 
-template<std::size_t N, typename Callable>
-consteval void call_for_coprimes_up_to(std::uintmax_t n, const std::array<std::uintmax_t, N>& basis, Callable call)
-{
-  for (std::uintmax_t i = 0u; i < n; ++i) {
-    if (std::apply([&i](auto... primes) { return ((i % primes != 0) && ...); }, basis)) {
-      call(i);
-    }
-  }
-}
-
-template<std::size_t N>
-[[nodiscard]] consteval std::size_t num_coprimes_up_to(std::uintmax_t n, const std::array<std::uintmax_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>
-[[nodiscard]] consteval auto coprimes_up_to(std::size_t n, const std::array<std::uintmax_t, N>& basis)
-{
-  std::array<std::uintmax_t, ResultSize> coprimes{};
-  std::size_t i = 0u;
-
-  call_for_coprimes_up_to(n, basis, [&coprimes, &i](std::uintmax_t cp) { coprimes[i++] = cp; });
-
-  return coprimes;
-}
-
 template<std::size_t N>
 [[nodiscard]] consteval std::uintmax_t product(const std::array<std::uintmax_t, N>& values)
 {
@@ -494,48 +441,34 @@ constexpr auto get_first_of(const Rng& rng, UnaryFunction f)
   return get_first_of(begin(rng), end(rng), f);
 }
 
-// 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 auto wheel_size = product(basis);
-  static constexpr auto coprimes_in_first_wheel =
-    coprimes_up_to<num_coprimes_up_to(wheel_size, basis)>(wheel_size, basis);
+[[nodiscard]] consteval std::uintmax_t find_first_factor(std::uintmax_t n)
+{
+  constexpr auto first_100_primes = first_n_primes<100>();
 
-  [[nodiscard]] static consteval std::uintmax_t find_first_factor(std::uintmax_t n)
-  {
-    if (const auto k = detail::get_first_of(basis, [&](auto p) { return first_factor_maybe(n, p); })) return *k;
+  for (const auto& p : first_100_primes) {
+    if (n % p == 0u) {
+      return p;
+    }
+    if (p * p > n) {
+      return n;
+    }
+  }
 
-    if (const auto k = detail::get_first_of(std::next(begin(coprimes_in_first_wheel)), end(coprimes_in_first_wheel),
-                                            [&](auto p) { return first_factor_maybe(n, p); }))
-      return *k;
-
-    for (std::size_t wheel = wheel_size; wheel < n; wheel += wheel_size)
-      if (const auto k =
-            detail::get_first_of(coprimes_in_first_wheel, [&](auto p) { return first_factor_maybe(n, wheel + p); }))
-        return *k;
+  // If we got this far, and we haven't found a factor, and haven't terminated... maybe it's prime?  Do a fast check.
+  if (baillie_psw_probable_prime(n)) {
     return n;
   }
 
-  [[nodiscard]] static consteval bool is_prime(std::size_t n) { return (n > 1) && find_first_factor(n) == n; }
-};
+  // If we're here, we know `n` is composite, so continue with trial division for all odd numbers.
+  std::uintmax_t factor = first_100_primes.back() + 2u;
+  while (factor * factor <= n) {
+    if (n % factor == 0u) {
+      return factor;
+    }
+    factor += 2u;
+  }
+
+  return n;  // Technically unreachable.
+}
 
 }  // namespace mp_units::detail

src/core/include/mp-units/framework/magnitude.h

@@ -51,13 +51,6 @@ import std;
 
 namespace mp_units {
 
-namespace detail {
-
-// Higher numbers use fewer trial divisions, at the price of more storage space.
-using factorizer = wheel_factorizer<4>;
-
-}  // namespace detail
-
 #if defined MP_UNITS_COMP_CLANG || MP_UNITS_COMP_CLANG < 18
 
 MP_UNITS_EXPORT template<symbol_text Symbol>
@@ -629,14 +622,10 @@ using common_magnitude_type = decltype(_common_magnitude_type_impl(M));
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 // `mag()` 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!
 MP_UNITS_EXPORT template<std::intmax_t N>
-constexpr std::optional<std::intmax_t> known_first_factor = std::nullopt;
+[[deprecated("`known_first_factor` is no longer necessary and can simply be removed")]]
+constexpr std::optional<std::intmax_t>
+  known_first_factor = std::nullopt;
 
 namespace detail {
 
@@ -646,12 +635,7 @@ template<std::intmax_t N>
 struct prime_factorization {
   [[nodiscard]] static consteval std::intmax_t get_or_compute_first_factor()
   {
-    constexpr auto opt = known_first_factor<N>;
-    if constexpr (opt.has_value()) {
-      return opt.value();  // NOLINT(bugprone-unchecked-optional-access)
-    } else {
-      return static_cast<std::intmax_t>(factorizer::find_first_factor(N));
-    }
+    return static_cast<std::intmax_t>(find_first_factor(N));
   }
 
   static constexpr std::intmax_t first_base = get_or_compute_first_factor();

src/systems/include/mp-units/systems/hep.h

@@ -35,9 +35,6 @@
 MP_UNITS_EXPORT
 namespace mp_units {
 
-template<>
-constexpr std::optional<std::intmax_t> known_first_factor<334'524'384'739> = 334'524'384'739;
-
 namespace hep {
 
 // energy

test/static/magnitude_test.cpp

@@ -31,9 +31,6 @@ import mp_units;
 using namespace units;
 using namespace units::detail;
 
-template<>
-constexpr std::optional<std::intmax_t> units::known_first_factor<9223372036854775783> = 9223372036854775783;
-
 namespace {
 
 // A set of non-standard bases for testing purposes.
@@ -182,20 +179,6 @@ static_assert(std::is_same_v<decltype(get_base(power_v<mag_2, 5, 8>{})), mag_2_>
 //     mag_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.
-//     mag<9'223'372'036'854'775'783>();
-//   }
-// }
-
 // TEST_CASE("magnitude converts to numerical value")
 // {
 //   SECTION("Positive integer powers of integer bases give integer values")

test/static/prime_test.cpp

@@ -35,51 +35,6 @@ namespace {
 
 inline constexpr auto MAX_U64 = std::numeric_limits<std::uint64_t>::max();
 
-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(109'561) == is_prime_by_trial_division(109'561));
-
-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));
-
 // Modular arithmetic.
 static_assert(add_mod(1u, 2u, 5u) == 3u);
 static_assert(add_mod(4u, 4u, 5u) == 3u);