Add utilities to compute Jacobi symbols #322
Merged
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This is a core need for the Baillie-PSW implementation: we will use
Jacobi symbols to set the parameters for the Strong Lucas test.
The [Jacobi symbol] `(a/n)` is a multiplicative function of two numbers:
`a` being any integer, and `n` any positive odd integer. It only takes
values in `{-1, 0, 1}`. It's defined in terms of another function,
Legendre symbols, which makes it very intimidating to compute... but
actually, computation is straightforward, because the Jacobi symbol has
symmetries and reduction rules that let us skip computing the Legendre
symbols.
I wrote this implementation by reading the linked wikipedia page, and
following the rules. To test it, I wrote a few manual tests for cases
where the right answer was obvious. I also tested against an
implementation which I found on Wikipedia. I decided not to check the
latter test into version control, because I was unsure about the
licensing implications. However, the Jacobi symbol code will still get
robust testing _indirectly_ in the future, because our Baillie-PSW
implementation will depend on it.
Helps #217.
[Jacobi symbol]: https://en.wikipedia.org/wiki/Jacobi_symbol
chiphogg
added
the
release notes: ♻️ lib (refactoring)
geoffviola
approved these changes
Nov 12, 2024
| return a; | ||
| } | ||
|
|
||
| // The conversions `true` -> `1` and `false` -> `0` are guaranteed by the standard. |
There was a problem hiding this comment.
Nit: It might good to comment that this converts true to 1 and false to -1.
chiphogg
added a commit
that referenced
this pull request
Nov 13, 2024
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.
chiphogg
added a commit
that referenced
this pull request
Nov 13, 2024
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=100 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.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This is a core need for the Baillie-PSW implementation: we will use
Jacobi symbols to set the parameters for the Strong Lucas test.
The Jacobi symbol
(a/n)is a multiplicative function of two numbers:abeing any integer, andnany positive odd integer. It only takesvalues in
{-1, 0, 1}. It's defined in terms of another function,Legendre symbols, which makes it very intimidating to compute... but
actually, computation is straightforward, because the Jacobi symbol has
symmetries and reduction rules that let us skip computing the Legendre
symbols.
I wrote this implementation by reading the linked wikipedia page, and
following the rules. To test it, I wrote a few manual tests for cases
where the right answer was obvious. I also tested against an
implementation which I found on Wikipedia. I decided not to check the
latter test into version control, because I was unsure about the
licensing implications. However, the Jacobi symbol code will still get
robust testing indirectly in the future, because our Baillie-PSW
implementation will depend on it.
Helps #217.