Add the Strong Lucas Probable Prime test

This is more involved than the Miller-Rabin test, but we can tame the
complexity by breaking it down into helper functions, performing tasks
such as:

- Increment the index of the (U_k, V_k) sequence elements by one.
- Double the index of the (U_k, V_k) sequence elements.
- Find an appropriate D parameter.

etc.

With these helpers, the algorithm becomes straightforward (see, for
instance,
https://en.wikipedia.org/wiki/Lucas_pseudoprime#Strong_Lucas_pseudoprimes).
We start by ruling out perfect squares (because if we don't, then the
search for `D` will never terminate).  Then we find our `D`, and
decompose `n + 1` into `s` and `d` parameters (exactly as we did for
Miller-Rabin, except there we used `n - 1`).  At this point, the strong
test is easy: check whether `U_d` is 0, then check `V_d`, as well as `V`
for all successive doublings of the index less than `s`.

A similar testing strategy as for the Miller Rabin gives us sufficient
confidence.

1. Test that we get small primes right.
2. Test that we get known pseudoprimes "correctly wrong".
3. Test some really big primes.

(Remember, a probable prime test must mark every actual prime as
"probably prime".)

Helps #509.

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

@@ -260,6 +260,124 @@ struct NumberDecomposition {
   }
 }
 
+// The "D" parameter in the Strong Lucas Probable Prime test.
+//
+// Default construction produces the first value to try, according to Selfridge's parameter selection.
+// Calling `successor()` on this repeatedly will produce the sequence of values to try.
+struct LucasDParameter {
+  uint64_t mag = 5u;
+  bool pos = true;
+
+  friend constexpr int as_int(LucasDParameter p)
+  {
+    int D = static_cast<int>(p.mag);
+    return p.pos ? D : -D;
+  }
+  friend constexpr LucasDParameter successor(LucasDParameter p) { return {.mag = p.mag + 2u, .pos = !p.pos}; }
+};
+
+// The first `D` in the infinite sequence {5, -7, 9, -11, ...} whose Jacobi symbol is -1.  This is the D we
+// want to use for the Strong Lucas Probable Prime test.
+//
+// Precondition: `n` is not a perfect square.
+[[nodiscard]] consteval LucasDParameter find_first_D_with_jacobi_symbol_neg_one(uint64_t n)
+{
+  LucasDParameter D{};
+  while (jacobi_symbol(as_int(D), n) != -1) {
+    D = successor(D);
+  }
+  return D;
+}
+
+// An element of the Lucas sequence, with parameters tuned for the Strong Lucas test.
+//
+// The default values give the first element (i.e., k=1) of the sequence.
+struct LucasSequenceElement {
+  uint64_t u = 1u;
+  uint64_t v = 1u;
+};
+
+// Produce the Lucas element whose index is twice the input element's index.
+[[nodiscard]] consteval LucasSequenceElement double_strong_lucas_index(const LucasSequenceElement& element, uint64_t n,
+                                                                       LucasDParameter D)
+{
+  const auto& [u, v] = element;
+
+  uint64_t v_squared = mul_mod(v, v, n);
+  uint64_t D_u_squared = mul_mod(D.mag, mul_mod(u, u, n), n);
+  uint64_t new_v = D.pos ? add_mod(v_squared, D_u_squared, n) : sub_mod(v_squared, D_u_squared, n);
+  new_v = half_mod_odd(new_v, n);
+
+  return {.u = mul_mod(u, v, n), .v = new_v};
+}
+
+[[nodiscard]] consteval LucasSequenceElement increment_strong_lucas_index(const LucasSequenceElement& element,
+                                                                          uint64_t n, LucasDParameter D)
+{
+  const auto& [u, v] = element;
+
+  const auto new_u = half_mod_odd(add_mod(u, v, n), n);
+
+  const auto D_u = mul_mod(D.mag, u, n);
+  auto new_v = D.pos ? add_mod(v, D_u, n) : sub_mod(v, D_u, n);
+  new_v = half_mod_odd(new_v, n);
+
+  return {.u = new_u, .v = new_v};
+}
+
+[[nodiscard]] consteval LucasSequenceElement find_strong_lucas_element(uint64_t i, uint64_t n, LucasDParameter D)
+{
+  LucasSequenceElement element{};
+
+  bool bits[64] = {};
+  std::size_t n_bits = 0u;
+  while (i > 1u) {
+    bits[n_bits++] = (i & 1u);
+    i >>= 1u;
+  }
+
+  for (auto j = n_bits; j > 0u; --j) {
+    element = double_strong_lucas_index(element, n, D);
+    if (bits[j - 1u]) {
+      element = increment_strong_lucas_index(element, n, D);
+    }
+  }
+
+  return element;
+}
+
+// Perform a Strong Lucas Probable Prime test on `n`.
+//
+// Precondition: (n >= 2).
+// Precondition: (n is odd).
+[[nodiscard]] consteval bool strong_lucas_probable_prime(uint64_t n)
+{
+  MP_UNITS_EXPECTS_DEBUG(n >= 2u);
+  MP_UNITS_EXPECTS_DEBUG(n % 2u == 1u);
+
+  if (is_perfect_square(n)) {
+    return false;
+  }
+
+  const auto D = find_first_D_with_jacobi_symbol_neg_one(n);
+
+  const auto [s, d] = decompose(n + 1u);
+
+  auto element = find_strong_lucas_element(d, n, D);
+  if (element.u == 0u) {
+    return true;
+  }
+
+  for (auto i = 0u; i < s; ++i) {
+    if (element.v == 0u) {
+      return true;
+    }
+    element = double_strong_lucas_index(element, n, D);
+  }
+
+  return false;
+}
+
 [[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)) {

test/static/prime_test.cpp

@@ -165,4 +165,36 @@ static_assert(!is_perfect_square(BIG_SQUARE - 1u));
 static_assert(is_perfect_square(BIG_SQUARE));
 static_assert(!is_perfect_square(BIG_SQUARE + 1u));
 
+// Tests for the Strong Lucas Probable Prime test.
+static_assert(as_int(LucasDParameter{.mag = 5, .pos = true}) == 5);
+static_assert(as_int(LucasDParameter{.mag = 7, .pos = false}) == -7);
+
+static_assert(as_int(LucasDParameter{}) == 5, "First D parameter in the sequence is 5");
+static_assert(as_int(successor(LucasDParameter{})) == -7, "Incrementing adds 2 to the mag, and flips the sign");
+static_assert(as_int(successor(successor(LucasDParameter{}))) == 9);
+static_assert(as_int(successor(successor(successor(LucasDParameter{})))) == -11);
+
+static_assert(strong_lucas_probable_prime(3u), "Known small prime");
+static_assert(strong_lucas_probable_prime(5u), "Known small prime");
+static_assert(strong_lucas_probable_prime(7u), "Known small prime");
+static_assert(!strong_lucas_probable_prime(9u), "Known small composite");
+
+// Test some Miller-Rabin pseudoprimes (https://oeis.org/A001262), which should NOT be marked prime.
+static_assert(!strong_lucas_probable_prime(2047u), "Miller-Rabin pseudoprime");
+static_assert(!strong_lucas_probable_prime(3277u), "Miller-Rabin pseudoprime");
+static_assert(!strong_lucas_probable_prime(486737u), "Miller-Rabin pseudoprime");
+
+// Test some Strong Lucas pseudoprimes (https://oeis.org/A217255).
+static_assert(strong_lucas_probable_prime(5459u), "Strong Lucas pseudoprime");
+static_assert(strong_lucas_probable_prime(5777u), "Strong Lucas pseudoprime");
+static_assert(strong_lucas_probable_prime(10877u), "Strong Lucas pseudoprime");
+static_assert(strong_lucas_probable_prime(324899u), "Strong Lucas pseudoprime");
+
+// Test some actual primes
+static_assert(strong_lucas_probable_prime(225'653'407'801u), "Large known prime");
+static_assert(strong_lucas_probable_prime(334'524'384'739u), "Large known prime");
+static_assert(strong_lucas_probable_prime(9'007'199'254'740'881u), "Large known prime");
+
+static_assert(strong_lucas_probable_prime(18'446'744'073'709'551'557u), "Largest 64-bit prime");
+
 }  // namespace