sieve.mlw 1.93 KB
Newer Older
1

2 3 4 5 6 7 8 9
(** Sieve of Eratosthenes

    Author: Martin Clochard

    See also knuth_prime_numbers.mlw in the gallery.
*)


10
module Sieve
11

12 13 14 15
  use import int.Int
  use import array.Array
  use import ref.Ref
  use import number.Prime
16

17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
  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) : array bool
    requires { n > 1 }
    returns { m -> 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
36 37 38 39
        let s = ref (!i * !i) in
        let ghost r = ref !i in
        while !s < n do
          invariant { 1 < !r <= n /\ !s = !r * !i }
40 41 42 43
          invariant { forall j. 0 <= j < n ->
            t[j] <-> (no_factor_lt !i j /\
                      forall k. 1 < k < !r -> j <> k * !i) }
          variant { n - !r }
44 45 46
          t[!s] <- false;
          s := !s + !i;
          incr r
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
        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
63

64 65
end