## Knuth's prime numbers

Knuth's algorithm for prime numbers. The Art of Computer Programming, vol 1, page 147.

**Authors:** Jean-Christophe Filliâtre

**Topics:** Array Data Structure

**Tools:** Why3

**References:** The Art of Computer Programming

see also the index (by topic, by tool, by reference, by year)

# Knuth's algorithm for prime numbers

The Art of Computer Programming, vol 1, page 147.

The following code computes a table of the first `m`

prime numbers.
Though there are more efficient ways of computing prime numbers,
the nice thing about this code is that it requires not less than
Bertrand's postulate (for `n > 1`

, there is always a prime `p`

such that
`n < p < 2n`

) to be proved correct.

This program had been proved correct using Coq and an early version of Why back in 2003, by Laurent Théry (INRIA Sophia-Antipolis): Laurent Théry, Proving Pearl: Knuth's Algorithm for Prime Numbers, TPHOLs 2003. Truly a tour de force, this proof includes the full proof of Bertrand's postulate in Coq. Here, we simply focus on the program verification part, posing Bertrand's postulate as a lemma that we do not prove.

Note: Knuth's code is entering the loop where a new prime number is
added, thus resulting into the immediate addition of 3 as `prime[1]`

.
It allows subsequent division tests to start at `prime[1]=3`

, thus
saving the division by `prime[0]=2`

. We did something similar in the
code below, though the use of a while loop looks like we did a
special case for 3 as well.

module PrimeNumbers use int.Int use int.ComputerDivision use number.Parity use number.Divisibility use number.Prime use ref.Refint use map.Map predicate no_prime_in (l u: int) = forall x: int. l < x < u -> not (prime x) predicate first_primes (p: int -> int) (u: int) = p[0] = 2 /\ (* sorted *) (forall i j: int. 0 <= i < j < u -> p[i] < p[j]) /\ (* only primes *) (forall i: int. 0 <= i < u -> prime p[i]) /\ (* all primes *) (forall i: int. 0 <= i < u-1 -> no_prime_in p[i] p[i+1])

`p[0]..p[u-1]`

are the first u prime numbers

lemma exists_prime: forall p: int -> int, u: int. 1 <= u -> first_primes p u -> forall d: int. 2 <= d <= p[u-1] -> prime d -> exists i: int. 0 <= i < u /\ d = p[i] axiom Bertrand_postulate [@W:non_conservative_extension:N] : forall p: int. prime p -> not (no_prime_in p (2*p))

Bertrand's postulate, admitted as an axiom (the label is there to suppress the warning issued by Why3)

use array.Array let prime_numbers (m: int) : array int requires { m >= 2 } ensures { result.length = m } ensures { first_primes result.elts m } = let p = make m 0 in p[0] <- 2; p[1] <- 3; let n = ref 5 in (* candidate for next prime *) for j = 2 to m - 1 do invariant { first_primes p.elts j } invariant { p[j-1] < !n < 2*p[j-1] } invariant { odd !n } invariant { no_prime_in p[j-1] !n } let rec test (k: int) variant { 2*p[j-1] - !n, j - k } requires { 1 <= k < j } requires { first_primes p.elts j } requires { p[j-1] < !n < 2*p[j-1] } requires { odd !n } requires { no_prime_in p[j-1] !n } requires { forall i: int. 0 <= i < k -> not (divides p[i] !n) } ensures { p[j-1] < !n /\ prime !n /\ no_prime_in p[j-1] !n } = if mod !n p[k] = 0 then begin assert { not (prime !n) }; n += 2; test 1 end else if div !n p[k] > p[k] then test (k + 1) else assert { prime !n } in test 1; p[j] <- !n; n += 2 done; p

`prime_numbers m`

returns an array containing the first `m`

prime numbers

exception BenchFailure let bench () raises { BenchFailure -> true } = let t = prime_numbers 100 in if t[0] <> 2 then raise BenchFailure; if t[1] <> 3 then raise BenchFailure; if t[2] <> 5 then raise BenchFailure; if t[3] <> 7 then raise BenchFailure; if t[4] <> 11 then raise BenchFailure; if t[5] <> 13 then raise BenchFailure; if t[6] <> 17 then raise BenchFailure; if t[7] <> 19 then raise BenchFailure; if t[8] <> 23 then raise BenchFailure; if t[9] <> 29 then raise BenchFailure; t end

download ZIP archive