sidef-scripts/Math/Baillie-PSW_high-level.sf

#!/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))