examples/sieve.mlw

   1
   2 (** Sieve of Eratosthenes
   3
   4     Author: Martin Clochard
   5
   6     See also knuth_prime_numbers.mlw in the gallery.
   7 *)
   8
   9
  10 module Sieve
  11
  12   use int.Int
  13   use array.Array
  14   use ref.Ref
  15   use number.Prime
  16
  17   predicate no_factor_lt (bnd num:int) =
  18     num > 1 /\ forall k l. 1 < l < bnd /\ k > 1 -> num <> k * l
  19
  20   let incr (r:ref int) : unit
  21     ensures { !r = old !r + 1 }
  22   = r := !r + 1
  23
  24   let sieve (n:int) : (m: array bool)
  25     requires { n > 1 }
  26     ensures  { length m = n /\ forall i. 0 <= i < n -> m[i] <-> prime i }
  27   = let t = Array.make n true in
  28     t[0] <- false;
  29     t[1] <- false;
  30     let i = ref 2 in
  31     while !i < n do
  32       invariant { 1 < !i <= n }
  33       invariant { forall j. 0 <= j < n -> t[j] <-> no_factor_lt !i j }
  34       variant { n - !i }
  35       if t[!i] then begin
  36         let s = ref (!i * !i) in
  37         let ghost r = ref !i in
  38         while !s < n do
  39           invariant { 1 < !r <= n /\ !s = !r * !i }
  40           invariant { forall j. 0 <= j < n ->
  41             t[j] <-> (no_factor_lt !i j /\
  42                       forall k. 1 < k < !r -> j <> k * !i) }
  43           variant { n - !r }
  44           t[!s] <- false;
  45           s := !s + !i;
  46           incr r
  47         done;
  48         assert { forall j. 0 <= j < n /\ t[j] ->
  49           (forall k l. 1 < l < !i + 1 -> j = k * l /\ k > 1 ->
  50             (if l = !i then k < !r && false else false) && false) &&
  51           no_factor_lt (!i+1) j }
  52       end else assert { forall j. 0 <= j < n /\ no_factor_lt !i j ->
  53         (forall k l. 1 < l < !i + 1 -> j = k * l /\ k > 1 ->
  54           (if l = !i then (forall k0 l. 1 < l < !i /\ k0 > 1 /\ !i = k0 * l ->
  55               j = (k*k0) * l && false) && false
  56                      else false) && false) && no_factor_lt (!i+1) j };
  57       incr i
  58     done;
  59     assert { forall j. 0 <= j < n /\ no_factor_lt n j -> prime j };
  60     assert { forall j. 0 <= j < n /\ prime j ->
  61       forall k l. 1 < l < n /\ k > 1 -> j = k * l -> l >= j && false };
  62     t
  63
  64 end
  65