Sieve of Eratosthenes

Author: Martin Clochard

See also knuth_prime_numbers.mlw in the gallery.

module Sieve

  use int.Int
  use array.Array
  use ref.Ref
  use number.Prime

  predicate no_factor_lt (bnd num:int) =
    num > 1 /\ forall k l. 1 < l < bnd /\ k > 1 -> num <> k * l

  let incr (r:ref int) : unit
    ensures { !r = old !r + 1 }
  = r := !r + 1

  let sieve (n:int) : (m: array bool)
    requires { n > 1 }
    ensures  { length m = n /\ forall i. 0 <= i < n -> m[i] <-> prime i }
  = let t = Array.make n true in
    t[0] <- false;
    t[1] <- false;
    let i = ref 2 in
    while !i < n do
      invariant { 1 < !i <= n }
      invariant { forall j. 0 <= j < n -> t[j] <-> no_factor_lt !i j }
      variant { n - !i }
      if t[!i] then begin
        let s = ref (!i * !i) in
        let ghost r = ref !i in
        while !s < n do
          invariant { 1 < !r <= n /\ !s = !r * !i }
          invariant { forall j. 0 <= j < n ->
            t[j] <-> (no_factor_lt !i j /\
                      forall k. 1 < k < !r -> j <> k * !i) }
          variant { n - !r }
          t[!s] <- false;
          s := !s + !i;
          incr r
        done;
        assert { forall j. 0 <= j < n /\ t[j] ->
          (forall k l. 1 < l < !i + 1 -> j = k * l /\ k > 1 ->
            (if l = !i then k < !r && false else false) && false) &&
          no_factor_lt (!i+1) j }
      end else assert { forall j. 0 <= j < n /\ no_factor_lt !i j ->
        (forall k l. 1 < l < !i + 1 -> j = k * l /\ k > 1 ->
          (if l = !i then (forall k0 l. 1 < l < !i /\ k0 > 1 /\ !i = k0 * l ->
              j = (k*k0) * l && false) && false
                     else false) && false) && no_factor_lt (!i+1) j };
      incr i
    done;
    assert { forall j. 0 <= j < n /\ no_factor_lt n j -> prime j };
    assert { forall j. 0 <= j < n /\ prime j ->
      forall k l. 1 < l < n /\ k > 1 -> j = k * l -> l >= j && false };
    t

end