123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211 |
- #!/usr/bin/ruby
- # Author: Daniel "Trizen" Șuteu
- # Date: 31 October 2023
- # https://github.com/trizen
- # A new primality test, using only the Lucas V sequence.
- # This test is a simplification of the strengthen BPSW test:
- # https://arxiv.org/abs/2006.14425
- define USE_METHOD_A_STAR = false # true to use the A* method in finding (P,Q)
- func partial_lucasVmod_pow2(P, Q, two_val, m, V1=2, V2=P, Q1=1, Q2=1) {
- Q1 = mulmod(Q1, Q2, m)
- Q2 = mulmod(Q1, Q, m)
- V1 = mulsubmulmod(V1, V2, P, Q1, m)
- Q2 = mulmod(Q2, Q1, m)
- two_val.times {
- V1 = mulsubmulmod(V1, V1, 2, Q2, m)
- Q2 = mulmod(Q2, Q2, m)
- }
- return (V1, Q2)
- }
- func partial_lucasVmod(P, Q, bits, m, V1=2, V2=P, Q1=1, Q2=1) {
- for bit in bits {
- Q1 = mulmod(Q1, Q2, m)
- if (bit) {
- Q2 = mulmod(Q1, Q, m)
- V1 = mulsubmulmod(V2, V1, P, Q1, m)
- V2 = mulsubmulmod(V2, V2, 2, Q2, m)
- }
- else {
- Q2 = Q1
- V2 = mulsubmulmod(V2, V1, P, Q1, m)
- V1 = mulsubmulmod(V1, V1, 2, Q2, m)
- }
- }
- return (V1,V2,Q1,Q2)
- }
- func check_lucasV(P,Q,n,m) {
- var b1 = n.bits
- var b2 = n.inc.bits
- var k = 0
- if (b1.len == b2.len) {
- k = b1.range.first{|i| b1[i] != b2[i] }
- }
- var(V1,V2,Q1,Q2) = partial_lucasVmod(P, Q, b1.first(k), m)
- var(V1_a,Q2_a) = partial_lucasVmod_pow2(P, Q, b2.end - k, m, V1, V2, Q1, Q2)
- V1_a.is_congruent(2*Q, m) || return false
- Q2_a.is_congruent(Q*Q, m) || return false
- var(V1_b,_,_,Q2_b) = partial_lucasVmod(P, Q, b1.slice(k), m, V1, V2, Q1, Q2)
- V1_b.is_congruent(P,m) || return false
- Q2_b.is_congruent(Q*kronecker(Q, m), m) || return false
- return true
- }
- func findQ(N) {
- for k in (2 .. Inf) {
- var D = ((-1)**k * (2*k + 1))
- var K = kronecker(D, N)
- if (K.is_zero && gcd(D, N).is_between(2, N-1)) {
- return nil
- }
- elsif ((k == 20) && N.is_square) {
- return nil
- }
- return ((1-D)/4) if (K == -1)
- }
- }
- func findP (N, Q) {
- for P in (2..Inf) {
- var D = (P*P - 4*Q)
- var K = kronecker(D, N)
- if (K == -1) {
- return P
- }
- elsif (K.is_zero && gcd(D, N).is_between(2, N-1)) {
- return nil
- }
- elsif ((P == 20) && N.is_square) {
- return nil
- }
- }
- }
- func is_extra_bfsw_psp(n) {
- n <= 1 && return false
- n == 2 && return true
- n.is_even && return false
- var(P,Q)
- if (USE_METHOD_A_STAR) {
- P = 1
- Q = findQ(n) \\ return false
- if (Q == -1) {
- P = 5
- Q = 5
- }
- }
- else { # this is faster
- Q = -2
- P = findP(n, Q) \\ return false
- }
- check_lucasV(P,Q,n,n)
- }
- say 25.by(is_extra_bfsw_psp)
- assert([913, 150267335403, 430558874533, 14760229232131, 936916995253453].none(is_extra_bfsw_psp))
- for n in (1..1e3) {
- if (is_extra_bfsw_psp(n)) {
- if (!n.is_prime) {
- say "Counter-example: #{n}"
- }
- }
- elsif (n.is_prime) {
- say "Missed-prime: #{n}"
- }
- }
- __END__
- Inspired by the paper "Strengthening the Baillie-PSW primality test", I propose a simplified test based on Lucas V-pseudoprimes, that requires computing only the Lucas V sequence, making it faster than the full BPSW test, while being about as strong.
- The first observation was that none of the 5 vpsp terms < 10^15 satisfy:
- Q^(n+1) == Q^2 (mod n)
- This gives us a simple test:
- V_{n+1}(P,Q) == 2*Q (mod n)
- Q^(n+1) == Q^2 (mod n)
- where (P,Q) are selected using Method A*.
- At very little additional computational cost (on average), we can make the test even stronger, by also checking:
- V_n(P,Q) == P (mod n)
- Notice that also none of the 5 vpsp terms < 10^15 satisfy the above congruence.
- The trick for computing V_n with very little additional computational cost (on average), is to compute the partial value of the Lucas V sequence, using the most significant overlapping bits of n and n+1.
- First we compute:
- V_d(P,Q) mod n
- where d is the "most significant overlapping binary part" of n and n+1.
- For example, if n = 43, we have:
- n = 101011_2
- n+1 = 101100_2
- The most significant overlapping bits of n and n+1 are: "101", therefore d = 101_2 = 5.
- From V_d(P,Q) mod n, we compute V_{n+1}(P,Q) mod n, using the remaining bits of n+1: "100".
- Notice that the remaining bits of n+1 always form a power of two, allowing us to optimize the computation of V_{n+1}(P,Q) mod n.
- At this stage, we check the necessary congruences trying to return early:
- V_{n+1}(P,Q) == 2*Q (mod n)
- Q^(n+1) == Q^2 (mod n)
- If the number passed the above congruences, we compute V_n(P,Q) mod n from V_d(P,Q) mod n, using the remaining bits of n: "011", then we check:
- V_n(P,Q) == P (mod n)
- Q^((n+1)/2) == Q*(Q|n) (mod n)
- Finally, we return true if the number satisfied all the congruences, indicating that it is probably prime.
- There are no known counter-examples to the presented test.
- Remarks:
- - For numbers of the form n = 4*x + 1, only the last last two bits differ from n and n+1, therefore only two extra steps in the "partial_lucasVmod()" function are needed to also compute V_n(P,Q) mod n, which is very cheap.
- - On the other hand, for numbers of the form n = 2^k - 1, all the bits of n and n+1 are different, which makes the computation of V_n(P,Q) quite expensive. But we can use the Lucas-Lehmer test for such numbers.
- - Numbers of the form x*2^k - 1, with x < 2^k, also take longer to check, but we can use the Lucas-Lehmer-Riesel (LLR) test for those.
- Optimization ideas:
- - To ensure that the test is always fast, we can skip the computation of V_n(P,Q) if the length of the remaining bits of n is too large (e.g. larger than the number of bits of d). This bounds the running time of the test to: 1.5 * (the cost of computing V_n(P,Q) mod n), while still having no known counter-examples.
- - In the selection of parameters (P,Q), we can start with Q = -2 and finding the first P >= 2 that satisfies jacobi(P^2 - 4*Q, n) = -1. The reason being that it is faster for computers to multiply by powers of two, and thus it makes the computation of the Lucas V sequence a bit faster, since |Q| is a power of two and, most of the time, P is also a power of 2.
- - In a general-purpose "is_prime(n)" function, for performance reasons, we should also do a little bit of trial-division (or gcd with primorials) and then a strong pseudoprime test to base 2, trying to return early if possible.
|