123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140 |
- #!/usr/bin/ruby
- # Daniel "Trizen" Șuteu
- # Date: 04 March 2020
- # https://github.com/trizen
- # High-level implementation of the Baillie-PSW primality.
- # The strong and extra-strong versions.
- # See also:
- # https://oeis.org/A217255 -- Strong Lucas pseudoprimes
- # https://oeis.org/A217719 -- Extra strong Lucas pseudoprimes
- # https://arxiv.org/pdf/2006.14425.pdf -- STRENGTHENING THE BAILLIE-PSW PRIMALITY TEST
- # Find Q for P = 1, such that kronecker(1 - 4Q, n) = -1.
- 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
- }
- return ((1 - D)/4) if (K == -1)
- }
- }
- # Find P >= 3 for Q = 1, such that kronecker(P^2 - 4, n) = -1.
- func findP(N) {
- for P in (3 .. Inf) {
- var D = (P*P - 4)
- var K = kronecker(D, N)
- if (K.is_zero && gcd(D, N).is_between(2, N-1)) {
- return nil
- }
- return P if (K == -1)
- }
- }
- # Strong Lucas pseudoprime test.
- func is_strong_lucas_pseudoprime(n) {
- return false if (n <= 1)
- return true if (n == 2)
- return false if n.is_even
- return false if n.is_power
- var P = 1
- var Q = findQ(n) \\ return false
- var k = valuation(n+1, 2)
- var d = (n+1)>>k
- if (lucasUmod(P, Q, d, n) == 0) {
- return true
- }
- for r in (^k) {
- if (lucasVmod(P, Q, d * 2**r, n) == 0) {
- return true
- }
- }
- return false
- }
- # Extra-strong Lucas pseudoprime test.
- func is_extra_strong_lucas_pseudoprime(n) {
- return false if (n <= 1)
- return true if (n == 2)
- return false if n.is_even
- return false if n.is_power
- var Q = 1
- var P = findP(n) \\ return false
- var k = valuation(n+1, 2)
- var d = (n+1)>>k
- if (lucasUmod(P, Q, d, n) == 0) {
- var V = lucasVmod(P, Q, d, n)
- if (V.is_congruent(2, n) || V.is_congruent(-2, n)) {
- return true
- }
- }
- for r in (^k) {
- if (lucasVmod(P, Q, d * 2**r, n) == 0) {
- return true
- }
- }
- return false
- }
- func strong_Baillie_PSW (n) {
- return false if (n <= 1)
- return true if (n == 2)
- return false if (n.is_even)
- n.is_strong_pseudoprime(2) && is_strong_lucas_pseudoprime(n)
- }
- func extra_strong_Baillie_PSW (n) {
- return false if (n <= 1)
- return true if (n == 2)
- return false if (n.is_even)
- n.is_strong_pseudoprime(2) && is_extra_strong_lucas_pseudoprime(n)
- }
- var strong_lucas_psp = [5459, 5777, 10877, 16109, 18971, 22499, 24569, 25199, 40309, 58519, 75077, 97439, 100127, 113573, 115639, 130139, 155819, 158399, 161027, 162133, 176399, 176471, 189419, 192509, 197801, 224369, 230691, 231703, 243629, 253259, 268349, 288919, 313499, 324899]
- var extra_strong_lucas_psp = [989, 3239, 5777, 10877, 27971, 29681, 30739, 31631, 39059, 72389, 73919, 75077, 100127, 113573, 125249, 137549, 137801, 153931, 155819, 161027, 162133, 189419, 218321, 231703, 249331, 370229, 429479, 430127, 459191, 473891, 480689, 600059, 621781, 632249, 635627]
- assert(strong_lucas_psp.all(is_strong_lucas_pseudoprime))
- assert(extra_strong_lucas_psp.all(is_extra_strong_lucas_pseudoprime))
- assert(!strong_lucas_psp.all(is_extra_strong_lucas_pseudoprime))
- assert(!extra_strong_lucas_psp.all(is_strong_lucas_pseudoprime))
- say 25.by(strong_Baillie_PSW)
- say 25.by(extra_strong_Baillie_PSW)
- say 25.by(is_strong_lucas_pseudoprime)
- say 25.by(is_extra_strong_lucas_pseudoprime)
- assert_eq(strong_Baillie_PSW.grep(1..900), is_strong_lucas_pseudoprime.grep(1..900))
- assert_eq(extra_strong_Baillie_PSW.grep(1..900), is_extra_strong_lucas_pseudoprime.grep(1..900))
|