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.
The Strong Lucas test coming in the next PR will already be complicated
enough. It'll be convenient, and less distracting, if we already have
functions for certain operations we'll need.
One thing we'll need to do is detect inputs that are perfect squares.
Fortunately, this is pretty easy to do robustly and quickly, with
Newton's method. We _don't_ want to use `std::sqrt`, because that takes
us into the floating point domain for no good reason, which could give
us wrong answers for larger integers.
The other thing we need is Jacobi symbols. These are a lot more
obscure, but thankfully, still resonably straightforward to compute.
The Wikipedia page (https://en.wikipedia.org/wiki/Jacobi_symbol) has a
good explanation, and in particular, good instructions for computing
values.
With these utilities in place, the Strong Lucas code should be easier to
review.
This can mark a number as either "probably prime", or "definitely
composite". The first parameter is the base, and the second is the
number to test.
Future PRs will build up the Strong Lucas test which complements this,
and then form the Baillie-PSW test by composing the two.
Helps #506.
The prime-testing techniques we will use (Miller-Rabin, Strong Lucas)
all make heavy usage of modular arithmetic. Therefore, we lay those
foundations here, adding utilities to perform the basic arithmetic
operations robustly.
Since these are internal-only helper functions, we don't bother checking
the preconditions, although we state them clearly in the contract
comment for each utility. After C++26, we could add contracts for
these.
Helps #509.
