pratt.dfy

// *** Basics *****************************************************************

// A nonzero number's divisors are the numbers that divide it evenly.
//
// Also called "factors". The divisors of 12 are 1, 2, 3, 4, 6, and 12.
function divisors(k: nat): set<nat>
    requires k > 0
{
    set n: int | 0 < n <= k && k % n == 0
}

// A number is prime if it has exactly two divisors: itself and 1.
predicate isPrime(n: nat)
{
    1 < n && divisors(n) == {1, n}
}

method testPrime() {
    assert !isPrime(1);
    assert isPrime(2);
    assert isPrime(3);
    assert 2 in divisors(4);
    assert !isPrime(4);
}

// Integer exponentiation.
function ipow(base: nat, exponent: nat): nat
{
    if exponent == 0 then 1 else base * ipow(base, exponent - 1)
}

function product(s: seq<nat>): nat {
  if s == [] then 1 else s[0] * product(s[1..])
}

predicate isPrimeFactorization(s: seq<nat>, n: nat) {
    0 < n &&
        product(s) == n &&
        forall k :: k in s ==> isPrime(k)
}

predicate isPrimeFactorOf(p: nat, n: nat) {
    isPrime(p) && n % p == 0
}

lemma primeFactorIsPresent(s: seq<nat>, n: nat, p: nat)
    requires isPrimeFactorization(s, n)
    requires isPrimeFactorOf(p, n)
    ensures p in s


// *** Map ********************************************************************

function seqmap<A, B>(f: A -> B, a: seq<A>): (b: seq<B>)
    ensures |a| == |b| && forall i :: 0 <= i < |a| ==> f(a[i]) == b[i]
{
    if a == [] then [] else [f(a[0])] + seqmap(f, a[1..])
}


// *** Theory *****************************************************************

// --- Part 1

predicate coprime(a: nat, b: nat)
    requires a != 0 && b != 0
{
    !exists k | 1 < k < b :: a % k == 0 && b % k == 0
}

lemma lucas_part_1(n: nat, a: nat)
    requires 1 < n && 0 < a < n
    requires ipow(a, n - 1) % n == 1
    ensures coprime(n, a)

// --- Part 2

function units(n: nat): set<nat>
    requires 1 < n
{
    set k: nat | 0 < k < n && coprime(k, n) :: k
}

lemma lucas_part_2(n: nat, a: nat)
    requires 1 < n && 0 < a < n
    requires coprime(n, a)
    requires forall q :: isPrimeFactorOf(q, n - 1) ==> ipow(a, (n - 1) / q) % n != 1
    ensures |units(n)| == n - 1

// --- Part 3: If the order of the group of units of n is n - 1, then n is prime.

function range(a: nat, b: nat): set<nat>
    requires a <= b
{
    set k: nat | a <= k < b
}

lemma range_size(a: nat, b: nat)
    requires a <= b
    ensures |range(a, b)| == b - a
{
    if a == b {
        // Trivial.
    } else {
        range_size(a, b - 1);
        assert range(a, b) == range(a, b - 1) + {b - 1};
        assert |range(a, b)| == |range(a, b - 1)| + 1;
    }
}

lemma full_finite_subset_is_equal(ss: set<nat>, s: set<nat>)
    requires ss <= s
    requires |ss| >= |s|
    ensures ss == s
{
    assert forall k: nat | k in ss :: k in s; // definition of <= on sets
    assert forall k: nat | k in s :: k in ss by {
        forall k: nat | k in s ensures k in ss {
            if k !in ss {
                full_finite_subset_is_equal(ss, s - {k});
                assert |ss| < |s|;
            }
        }
    }
}

lemma lucas_part_3(n: nat)
    requires 1 < n
    requires |units(n)| == n - 1
    ensures isPrime(n)
{
    assert units(n) == range(1, n) by {
        assert forall k: nat | k in units(n) :: 1 <= k < n;
        assert units(n) <= range(1, n);
        range_size(1, n);
        full_finite_subset_is_equal(units(n), range(1, n));
    }
    forall k: nat | 0 < k < n ensures coprime(k, n) {
        assert k in units(n);
    }
    assert divisors(n) == {1, n} by {
        assert n % 1 == 0;
        assert n % n == 0;
        assert forall k: nat | 1 < k < n :: coprime(k, n) ==> k !in divisors(n);
        assert forall k: nat | 1 < k < n :: k !in divisors(n);
    }
}

// --- Put the pieces together

// Lucas primality test: Suppose we have an integer a such that:
// *   a^(n − 1) ≡ 1 (mod n),
// *   for every prime factor q of n − 1, it is not the case that
//     a^((n − 1)/q) ≡ 1 (mod n).
//
// Then n is prime.
//
lemma lucas(n: nat, a: nat)
    requires 1 < n && 0 < a < n
    requires ipow(a, n - 1) % n == 1
    requires forall q :: isPrimeFactorOf(q, n - 1) ==> ipow(a, (n - 1) / q) % n != 1
    ensures isPrime(n)
{
    lucas_part_1(n, a);
    lucas_part_2(n, a);
    lucas_part_3(n);
}


// *** Pratt primality certificates *******************************************

// A Pratt certificate is proof that a number n is prime.
//
// It has three fields:
// *   n, the prime number
// *   a, a number such that a^(n - 1) ≡ 1 (mod n)
// *   f, the prime factorization of n - 1.
//
// To complete the proof, the prime factors in f are themselves given as Certs.
// Eventually this bottoms out in `Cert(n: 2, a: 1, f: [])`, the prime factorization of 1 being empty.
///
// The Cert `c` only certifies that `c.n` is prime if `valid(c)`.
//
datatype Cert = Cert(n: nat, a: nat, f: seq<Cert>)

function factors(cert: Cert): seq<nat> {
    seqmap((subcert: Cert) => subcert.n, cert.f)
}

function getSubcert(cert: Cert, q: nat): (subcert: Cert)
    requires q in factors(cert)
    ensures subcert in cert.f
{
    var i :| 0 <= i < |factors(cert)| && factors(cert)[i] == q;
    assert ((subcert: Cert) => subcert.n)(cert.f[i]) == q;
    cert.f[i]
}

predicate valid(cert: Cert) {
    var Cert(n, a, f) := cert;
    1 < n &&
        0 < a < n &&
        ipow(a, n - 1) % n == 1 &&
        product(factors(cert)) == n - 1 &&
        forall subcert :: subcert in f ==> (
            valid(subcert) &&
            ipow(a, (n - 1) / subcert.n) % n != 1
        )
}

method testValidCert() {
    assert valid(Cert(2, 1, []));
}

lemma certified(cert: Cert)
    requires valid(cert)
    ensures isPrime(cert.n)
{
    var Cert(n, a, f) := cert;
    forall subcert | subcert in f
        ensures isPrime(subcert.n)
    {
        certified(subcert);
    }
    assert isPrimeFactorization(factors(cert), cert.n - 1) by {
        assert 0 < cert.n - 1;
        assert product(factors(cert)) == cert.n - 1;
        forall f | f in factors(cert) ensures isPrime(f) {
            certified(getSubcert(cert, f));
        }
    }
    forall q | isPrimeFactorOf(q, n - 1)
        ensures ipow(a, (n - 1) / q) % n != 1
    {
        primeFactorIsPresent(factors(cert), cert.n - 1, q);
        assert exists subcert :: subcert in f && subcert.n == q;
    }
    lucas(n, a);
}