lucas_pseudoprimes_generation_erdos_method.sf

#!/usr/bin/ruby

# Erdos construction method for Lucas D-pseudoprimes, for discriminant D = P^2-4Q:
#   1. Choose an even integer L with many divisors.
#   2. Let P be the set of primes p such that p-kronecker(D,p) divides L and p does not divide L.
#   3. Find a subset S of P such that n = prod(S) satisfies U_n(P,Q) == 0 (mod n) and kronecker(D,n) == -1.

# Alternatively:
#   3. Find a subset S of P such that n = prod(P) / prod(S) satisfies U_n(P,Q) == 0 (mod n) and kronecker(D,n) == -1.

func lucas_pseudoprimes(L, P=1, Q=-1) {

    var D = (P*P - 4*Q)

    var primes = do {
        var divisors = L.divisors

        var A = divisors.map { .dec }.grep {|p|
            p.is_odd && !(p `divides` L) && p.is_prime && (kronecker(D, p) == -1)
        }

        var B = divisors.map { .inc }.grep {|p|
            p.is_odd && !(p `divides` L) && p.is_prime && (kronecker(D, p) == 1)
        }

        A + B -> sort.uniq
    }

    for k in (2 .. primes.len) {
        primes.combinations(k, {|*S|
            var n = S.prod
            var k = kronecker(D, n)     # optionally, check if k == -1
            if (lucasUmod(P, Q, n - k, n) == 0) {
                say n
            }
        })
    }
}

lucas_pseudoprimes(720, 1, -1)

__END__
323
1891
6601
13981
342271
1590841
852841
3348961
9937081
16778881
72881641
10756801
154364221
205534681
609865201
807099601
1438048801
7692170761
921921121
32252538601
222182990161
2051541911881
2217716806743361