Use Baillie-PSW for prime factorization

Formerly, `is_prime` depended on `find_first_factor`.  Now, we reverse
that dependency!  This is good, because while factoring is generally
hard enough that we've built the foundations of global banking on that
difficulty, `is_prime` can be done in polynomial time --- and now is,
because we're using `baillie_psw`.  We have a `static_assert` to make
sure it's restricted to 64 bits, but even this could probably be
removed, because there aren't any known counterexamples of _any_ size.

For efficiency reasons, when factoring a number, it's common to do trial
division before moving on to the "fast" primality checkers.  We hardcode
a list of the "first N primes", taking N=1000 for now.  (We could also
compute them at compile time, but this turned out to have a very large
impact on compilation times.)  It should be easy to adjust the size of
this list if we decide later: there are tests to make sure that it
contains only primes, has them in order, and doesn't skip any.

The new `is_prime` is indeed very fast and accurate.  It now correctly
handles all of the "problem numbers" mentioned in #217, as well as the
(much much larger) largest 64-bit prime.

One more tiny fix: we had switched to use `std::abs` in #322, but this
doesn't actually work: `std::abs` won't be `constexpr` compatible until
C++23!  So as part of this change, we switch back to something that will
work.

Fixes #217.

au/code/au/utility/factoring.hh

@@ -16,18 +16,61 @@
 
 #include <cstdint>
 
+#include "au/utility/probable_primes.hh"
+
 namespace au {
 namespace detail {
 
+// Check whether a number is prime.
+constexpr bool is_prime(std::uintmax_t n) {
+    static_assert(sizeof(std::uintmax_t) <= sizeof(std::uint64_t),
+                  "Baillie-PSW only strictly guaranteed for 64-bit numbers");
+
+    return baillie_psw(n) == PrimeResult::PROBABLY_PRIME;
+}
+
+template <typename T = void>
+struct FirstPrimesImpl {
+    static constexpr uint16_t values[] = {
+        2,   3,   5,   7,   11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,  53,  59,
+        61,  67,  71,  73,  79,  83,  89,  97,  101, 103, 107, 109, 113, 127, 131, 137, 139,
+        149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233,
+        239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337,
+        347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439,
+        443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541};
+    static constexpr std::size_t N = sizeof(values) / sizeof(values[0]);
+};
+template <typename T>
+constexpr uint16_t FirstPrimesImpl<T>::values[];
+template <typename T>
+constexpr std::size_t FirstPrimesImpl<T>::N;
+using FirstPrimes = FirstPrimesImpl<>;
+
 // Find the smallest factor which divides n.
 //
 // Undefined unless (n > 1).
 constexpr std::uintmax_t find_first_factor(std::uintmax_t n) {
-    if (n % 2u == 0u) {
-        return 2u;
+    const auto &first_primes = FirstPrimes::values;
+    const auto &n_primes = FirstPrimes::N;
+
+    // First, do trial division against the first N primes.
+    for (const auto &p : first_primes) {
+        if (n % p == 0u) {
+            return p;
+        }
+
+        if (p * p > n) {
+            return n;
+        }
     }
 
-    std::uintmax_t factor = 3u;
+    // If we got this far, and haven't found a factor nor terminated, do a fast primality check.
+    if (is_prime(n)) {
+        return n;
+    }
+
+    // If we're here, we know `n` is composite, so continue with trial division for all odd numbers.
+    std::uintmax_t factor = first_primes[n_primes - 1u] + 2u;
     while (factor * factor <= n) {
         if (n % factor == 0u) {
             return factor;
@@ -38,9 +81,6 @@ constexpr std::uintmax_t find_first_factor(std::uintmax_t n) {
     return n;
 }
 
-// Check whether a number is prime.
-constexpr bool is_prime(std::uintmax_t n) { return (n > 1) && (find_first_factor(n) == n); }
-
 // Find the largest power of `factor` which divides `n`.
 //
 // Undefined unless n > 0, and factor > 1.

au/code/au/utility/probable_primes.hh

@@ -179,7 +179,7 @@ constexpr int jacobi_symbol(int64_t raw_a, uint64_t n) {
     // Starting conditions: transform `a` to strictly non-negative values, setting `result` to the
     // sign we pick up from this operation (if any).
     int result = bool_sign((raw_a >= 0) || (n % 4u == 1u));
-    auto a = static_cast<uint64_t>(std::abs(raw_a)) % n;
+    auto a = static_cast<uint64_t>(raw_a * bool_sign(raw_a >= 0)) % n;
 
     // Delegate to an implementation which can only handle positive numbers.
     return jacobi_symbol_positive_numerator(a, n, result);

au/code/au/utility/test/factoring_test.cc

@@ -22,6 +22,20 @@ namespace {
 std::uintmax_t cube(std::uintmax_t n) { return n * n * n; }
 }  // namespace
 
+TEST(FirstPrimes, HasOnlyPrimesInOrderAndDoesntSkipAny) {
+    const auto &first_primes = FirstPrimes::values;
+    const auto &n_primes = FirstPrimes::N;
+    auto i_prime = 0u;
+    for (auto i = 2u; i_prime < n_primes; ++i) {
+        if (i == first_primes[i_prime]) {
+            EXPECT_TRUE(is_prime(i)) << i;
+            ++i_prime;
+        } else {
+            EXPECT_FALSE(is_prime(i)) << i;
+        }
+    }
+}
+
 TEST(FindFirstFactor, ReturnsInputForPrimes) {
     EXPECT_EQ(find_first_factor(2u), 2u);
     EXPECT_EQ(find_first_factor(3u), 3u);
@@ -60,6 +74,17 @@ TEST(IsPrime, FindsPrimes) {
     EXPECT_FALSE(is_prime(196962u));
 }
 
+TEST(IsPrime, CanHandleVeryLargePrimes) {
+    for (const auto &p : {
+             225'653'407'801ull,
+             334'524'384'739ull,
+             9'007'199'254'740'881ull,
+             18'446'744'073'709'551'557ull,
+         }) {
+        EXPECT_TRUE(is_prime(p)) << p;
+    }
+}
+
 TEST(Multiplicity, CountsFactors) {
     constexpr std::uintmax_t n = (2u * 2u * 2u) * (3u) * (5u * 5u);
     EXPECT_EQ(multiplicity(2u, n), 3u);