Lucas Test is broken

[JASory]

Lucas test (and by extension the Baillie-PSW test used in "strong_check") flags some primes as composites. A cursory search found no composites that passed (given the rarity of even 2-SPRPs beyond 2^64 this is fully expected) ,however the aforementioned behavior is sufficient to demonstrate that the test does not perform as predicted and may pass more composites than a correct implementation. Below is one example.

N
18446744073710004191 Claim true

Witness 8300848561679324618

Factorization of n-1 2, 5, 1844674407371000419

[mikelodder7]

I'm failing to see where the implementation deviates from the known algorithm. This same use case fails in golangs math/big implementation.

[JASory]

I'm failing to see where the implementation deviates from the known algorithm. This same use case fails in golangs math/big implementation.

Then both are erroneous, it's much more useful to resort to proven results (like the Pratt certificate given above) or statistical measures than referencing some other software. As we regularly see with GMP (arguably the most prestigious mathematical library) they can contain numerous hidden errors. (The Prime and Prejudice paper you cite as the basis being a fine example of how the most popular number-theoretic libraries do not perform as intended by the developers, which makes this issue somewhat ironic).

I'll look at it some more since I am also going to implement a BPSW for integers in 2^64;2^128 in my library. If many people are making this error then it probably needs to be resolved.

[mikelodder7]

I agree with you. I’ll look more into this. In reality I’d love to ECPP instead which is faster but more complex to code

[JASory]

I agree with you. I’ll look more into this. In reality I’d love to ECPP instead which is faster but more complex to code

ECPP is not faster, it's much slower than BPSW. It does however provide a proof (much faster than Lucas' N-1 if that's what you meant).

FYI it also fails primes at a much lower bound as well, I just selected n > 2^64 since that's the machine word size and producing a Pratt certificate for much larger is very difficult. If you want you can use number-theory to rapidly check many smaller primes. It has a rigourously proven deterministic test under 2^64 + 2^45 as well as being able to construct prime-proofs (albeit the slower N-1 variant). (It's how I detected the errors in here, num-prime, and primal).

[mikelodder7]

Yes that’s what meant for comparison. I know ECPP fails at a lower bound which is why I already try the first 2048 primes by trial division

[JASory]

Yes that’s what meant for comparison. I know ECPP fails at a lower bound which is why I already try the first 2048 primes by trial division

?? I'm not aware of Elliptic curve primality failing at certain bounds, and even then trial division would not prevent it. Trial division is simply a method of quickly eliminating composites, even a weak fermat test is many times stronger than trial division by 2048 primes.

You also don't appear to use ECPP, and nothing said appears to apply to the Lucas sequence test either. (The counterexamples I found {if I remember correctly} were far greater than 17881 {2049-th prime}).

[mikelodder7]

Correct this code does trial division followed by Fermat, then miller rabin then Lucas. No ECPP yet. I mention it as another test I’d like to include depending on the size of the prime

[JASory]

Correct this code does trial division followed by Fermat, then miller rabin then Lucas. No ECPP yet. I mention it as another test I’d like to include depending on the size of the prime

I'm curious if you benchmarked to see if the weak fermat is faster than a standard SPRP test? I've experimented it and found it essentially no different, as far as I know it is used in GMP because they have a specialized exponentiation function for it. If the bigint library doesn't also have the function then it probably doesn't actually run any faster. (Purely speculation on my part).

[fjarri]

@JASory, do you know if num-prime still has issues with their Lucas implementation, or have they been fixed? I don't see any recent relevant commits.

I've experimented it and found it essentially no different, as far as I know it is used in GMP

Doesn't GMP use MR basis 2 + Lucas + MR random basis?

[fjarri]

I think I see the reason. In lucas(), after a few steps of looking for P, Q, we are supposed to check if the number is actually a square. The code reads:

if p == 40 && (&n_int * &n_int).sqrt() == n_int {

Which is trivially true. I suppose it is a typo, and should read

let n_sqrt = n_int.sqrt();
if p == 40 && &n_sqrt * &n_sqrt == n_int {