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
download ZIP archive