frobenius_primality_test.sf

#!/usr/bin/ruby

# Simple implementation of the Frobenius primality test.

# See also:
#   https://en.wikipedia.org/wiki/Frobenius_pseudoprime

func lucasUVmod(P, Q, n, m) {

    # Paper "On Lucas Sequences Computation" by Aleksey Koval:
    #   https://file.scirp.org/pdf/IJCNS20101200011_90376712.pdf

    var D  = (P*P - 4*Q)
    var(V1 = 2, V2 = P)
    var(Q1 = 1, Q2 = 1)

    assert(D.is_coprime(m), "P^2 - 4Q must be coprime to #{m}")

    n.as_bin.each {|bit|

        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)
        }
    }

    var Un = mulmod(mulsubmulmod(2, V2, P, V1, m), D.invmod(m), m)

    return (Un, V1)
}

func is_frobenius_pseudoprime(P, Q, n) {
    var D = (P*P - 4*Q)

    gcd(n, 2*Q*D) == 1 || return false

    var (k)    = kronecker(D, n)
    var (U, V) = lucasUVmod(P, Q, n-k, n)

    (U == 0) && ((k==1) ? (V == 2) : (V == (2*Q)%n))
}

func frobenius_primality_test(n) {

    return false if (n <= 1)
    return false if (n.is_power)
    return true  if (n == 2)
    return false if (n.is_even)
    return true  if (n <= 5)

    # Using the parameters (P,Q), where Q = 2 and P is the first
    # odd integer P >= 5 such that kronecker(P^2 - 4*Q, n) = -1.

    var Q = 2
    var P = (5..Inf `by` 2 -> first {|k|
        kronecker(k*k - 4*Q, n) == -1
    })

    is_frobenius_pseudoprime(P, Q, n)
}

var count = 0
var limit = 1e3

for n in (2..limit) {
    if (frobenius_primality_test(n)) {
        ++count
        assert(n.is_prime, n)
    }
    else {
        assert(n.is_composite, n)
    }
}

say "There are #{count} primes <= #{limit}."