examples/knuth_prime_numbers.mlw

   1
   2 (** {1 Knuth's algorithm for prime numbers}
   3
   4    The Art of Computer Programming, vol 1, page 147.
   5
   6    The following code computes a table of the first `m` prime numbers.
   7    Though there are more efficient ways of computing prime numbers,
   8    the nice thing about this code is that it requires not less than
   9    Bertrand's postulate (for `n > 1`, there is always a prime `p` such that
  10    `n < p < 2n`) to be proved correct.
  11
  12    This program had been proved correct using Coq and an early version of
  13    Why back in 2003, by Laurent Théry (INRIA Sophia-Antipolis): Laurent Théry,
  14    Proving Pearl: Knuth's Algorithm for Prime Numbers, TPHOLs 2003.
  15    Truly a tour de force, this proof includes the full proof of Bertrand's
  16    postulate in Coq. Here, we simply focus on the program verification part,
  17    posing Bertrand's postulate as a lemma that we do not prove.
  18
  19    Note: Knuth's code is entering the loop where a new prime number is
  20    added, thus resulting into the immediate addition of 3 as `prime[1]`.
  21    It allows subsequent division tests to start at `prime[1]=3`, thus
  22    saving the division by `prime[0]=2`. We did something similar in the
  23    code below, though the use of a while loop looks like we did a
  24    special case for 3 as well.  *)
  25
  26 module PrimeNumbers
  27
  28   use int.Int
  29   use int.ComputerDivision
  30   use number.Parity
  31   use number.Divisibility
  32   use number.Prime
  33   use ref.Refint
  34   use map.Map
  35
  36   predicate no_prime_in (l u: int) =
  37     forall x: int. l < x < u -> not (prime x)
  38
  39   (** `p[0]..p[u-1]` are the first u prime numbers *)
  40   predicate first_primes (p: int -> int) (u: int) =
  41     p[0] = 2 /\
  42     (* sorted *)
  43     (forall i j: int. 0 <= i < j < u -> p[i] < p[j]) /\
  44     (* only primes *)
  45     (forall i: int. 0 <= i < u -> prime p[i]) /\
  46     (* all primes *)
  47     (forall i: int. 0 <= i < u-1 -> no_prime_in p[i] p[i+1])
  48
  49   lemma exists_prime:
  50     forall p: int -> int, u: int. 1 <= u -> first_primes p u ->
  51     forall d: int. 2 <= d <= p[u-1] -> prime d ->
  52     exists i: int. 0 <= i < u /\ d = p[i]
  53
  54   (** Bertrand's postulate, admitted as an axiom
  55    (the label is there to suppress the warning issued by Why3) *)
  56   axiom Bertrand_postulate [@W:non_conservative_extension:N] :
  57     forall p: int. prime p -> not (no_prime_in p (2*p))
  58
  59   use array.Array
  60
  61   (** `prime_numbers m` returns an array containing the first `m`
  62   prime numbers *)
  63   let prime_numbers (m: int) : array int
  64     requires { m >= 2 }
  65     ensures  { result.length = m }
  66     ensures  { first_primes result.elts m }
  67   = let p = make m 0 in
  68     p[0] <- 2;
  69     p[1] <- 3;
  70     let n = ref 5 in (* candidate for next prime *)
  71     for j = 2 to m - 1 do
  72       invariant { first_primes p.elts j }
  73       invariant { p[j-1] < !n < 2*p[j-1] }
  74       invariant { odd !n }
  75       invariant { no_prime_in p[j-1] !n }
  76       let rec test (k: int) variant { 2*p[j-1] - !n, j - k }
  77         requires { 1 <= k < j }
  78         requires { first_primes p.elts j }
  79         requires { p[j-1] < !n < 2*p[j-1] }
  80         requires { odd !n }
  81         requires { no_prime_in p[j-1] !n }
  82         requires { forall i: int. 0 <= i < k -> not (divides p[i] !n) }
  83         ensures  { p[j-1] < !n /\ prime !n /\ no_prime_in p[j-1] !n }
  84       = if mod !n p[k] = 0 then begin
  85           assert { not (prime !n) }; n += 2; test 1
  86         end else if div !n p[k] > p[k] then
  87           test (k + 1)
  88         else
  89           assert { prime !n }
  90       in
  91       test 1;
  92       p[j] <- !n;
  93       n += 2
  94     done;
  95     p
  96
  97   exception BenchFailure
  98
  99   let bench () raises { BenchFailure -> true } =
 100      let t = prime_numbers 100 in
 101      if t[0] <> 2 then raise BenchFailure;
 102      if t[1] <> 3 then raise BenchFailure;
 103      if t[2] <> 5 then raise BenchFailure;
 104      if t[3] <> 7 then raise BenchFailure;
 105      if t[4] <> 11 then raise BenchFailure;
 106      if t[5] <> 13 then raise BenchFailure;
 107      if t[6] <> 17 then raise BenchFailure;
 108      if t[7] <> 19 then raise BenchFailure;
 109      if t[8] <> 23 then raise BenchFailure;
 110      if t[9] <> 29 then raise BenchFailure;
 111      t
 112
 113 end
 114