(* find-primes.ec

   Hoare logic example, proving the correctness of a program for
   computing all the primes that are less-than-or-equal to a given
   integer *)

prover quorum=2 ["Alt-Ergo" "Z3"].

require import AllCore IntDiv StdOrder List.
import IntOrder.

(* checking whether the elements of a list of integers appear in a
   strictly ascending order: *)

op nosmt strasc (xs : int list) : bool =  (* smt can't use directly *)
  forall (i j : int),
  0 <= i < j < size xs =>
  nth 0 xs i < nth 0 xs j.

lemma strasc_nil :
  strasc [].
proof. rewrite /strasc /#. qed.

lemma strasc_add (xs : int list, y : int) :
  strasc xs =>
  (xs <> [] => nth 0 xs (size xs - 1) < y) =>
  strasc (rcons xs y).
proof.
rewrite /strasc.
smt(size_rcons nth_rcons).
qed.

lemma strasc_sing (y : int) :
  strasc [y].
proof.
smt(strasc_nil strasc_add).
qed.

(* checking whether an integer is prime: *)

op is_prime (n : int) : bool =
  2 <= n /\
  ! exists (x : int), x %| n /\ 1 < x < n.

lemma has_prime_fact (n : int) :
  2 <= n =>
  exists (m : int), m <= n /\ m %| n /\ is_prime m.
proof.
have H :  (* put in form for using strong induction *)
  forall n,
  0 <= n => 2 <= n =>
  exists m, m <= n /\ m %| n /\ is_prime m.
  elim /sintind => i ge0_i IH ge2_i.
  case (is_prime i) => [prime_i | not_prime_i].
  exists i; smt().
  have [x [#] x_div_i gt1_x lt_x_i] :
    exists (x : int), x %| i /\ 1 < x < i.
    smt().
  have [m [#] m_le_x m_div_x prime_m] :
    exists m, m <= x /\ m %| x /\ is_prime m.
    apply IH => /#.
  exists m; smt(dvdz_trans).
move => ge2_n.
rewrite H /#.
qed.

(* checking whether the elements of xs are all and only the primes
   that are <= n, listed in strictly ascending order: *)

op nosmt primesupto (n : int, xs : int list) : bool =  (* smt can't use *)
  (forall (i : int),
   0 <= i < size xs =>
   nth 0 xs i <= n /\ is_prime (nth 0 xs i)) /\
  (forall (m : int),
   m <= n => is_prime m =>
   (exists (i : int), 0 <= i < size xs /\ nth 0 xs i = m)) /\
  strasc xs.

lemma primesupto_lt2 (n : int) :
  n < 2 => primesupto n [].
proof.
rewrite /primesupto /is_prime.
smt(strasc_nil).
qed.

lemma primesupto_succ_no (n : int, xs : int list) :
  primesupto n xs =>
  (exists (i : int),
   0 <= i < size xs /\ nth 0 xs i %| (n + 1)) =>
  primesupto (n + 1) xs.
proof.
rewrite /primesupto /#.
qed.

lemma primesupto_succ_yes (n : int, xs : int list) :
  primesupto n xs =>
  1 <= n =>
  ! (exists (i : int),
     0 <= i < size xs /\ nth 0 xs i %| (n + 1)) =>
  primesupto (n + 1) (rcons xs (n + 1)).
proof.
rewrite /primesupto.
move => [#] all_in_xs_prime all_prime_in_xs strasc_xs ge1_n good_nplus1.
split => [i [#] ge0_i |].
rewrite size_rcons nth_rcons ltzS => i_le_size_xs.
case (i < size xs) => [i_lt_size_xs /# | not_i_lt_size_xs].
have -> /= : i = size xs.
  smt().
rewrite /is_prime.
split => [/# |].
case (exists x, x %| (n + 1) /\ 1 < x < n + 1) =>
  [[x [#] x_div_n_plus1 gt1_x x_le_n] | //].
rewrite ltzS in x_le_n.
have [m [#] m_le_x m_div_x prime_m] :
  exists m, m <= x /\ m %| x /\ is_prime m.
  rewrite has_prime_fact /#.
smt(dvdz_trans).
split => [m m_le_n_plus1 is_prime_m |].
have [m_le_n | m_eq_n_plus1] : m <= n \/ m = n + 1.
  smt().
have [i [zero_le_i_lt_size_xs <-]] :
  exists i, 0 <= i < size xs /\ nth 0 xs i = m.
  by apply all_prime_in_xs.
exists i; smt(nth_rcons size_rcons).
exists (size xs); smt(nth_rcons size_rcons size_ge0).
rewrite strasc_add //.
case (xs = []) => [// | xs_ne_nil /=].
have /# : 0 < size xs.
  rewrite -(head_behead xs 0) //=.
  smt(size_ge0).
qed.

lemma primesupto2 :
  primesupto 2 [2].
proof.
have -> : [2] = rcons [] 2.
  trivial.
have {1}-> : 2 = 1 + 1.
  trivial.
rewrite primesupto_succ_yes //.
smt(primesupto_lt2).
qed.

module Primes = {
  proc f(n : int) : int list = {
    var m, j : int;
    var prs : int list;
    var found_div : bool;
    if (n < 2) {
      prs <- [];
    }
    else {
      prs <- [2]; m <- 2;
      while (m < n) {
        m <- m + 1;
        found_div <- false; j <- 0;
        while (j < size prs /\ ! found_div) {
          if (nth 0 prs j %| m) {
            found_div <- true;
          }
          j <- j + 1;
        }
        if (! found_div) {
          prs <- rcons prs m;
        }
      }
    }
    return prs;
  }
}.

lemma Primes_f (n_ : int) :
  hoare [Primes.f : n = n_ ==> primesupto n_ res].
proof.

qed.