## Knuth's prime numbers

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

Auteurs: Jean-Christophe Filliâtre

Catégories: Array Data Structure

Outils: Why3

Références: 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

```