Theory Sieve_Primes

theory Sieve_Primes
  imports
    "HOL-Computational_Algebra.Primes"
    "../Num_Class"
    "../HOLCF_Prelude"
(*FIXME: import order matters, if HOLCF_Prelude is not imported last,
constants like "filter" etc. refer to List.thy; I'm afraid however
that the current order only works by accident*)
begin

section ‹The Sieve of Eratosthenes›

declare [[coercion int]]
declare [[coercion_enabled]]

text ‹
  This example proves that the well-known Haskell two-liner that lazily
  calculates the list of all primes does indeed do so. This proof is using
  coinduction.
›

text ‹
  We need to hide some constants again since we imported something from HOL
  not via @{theory "HOLCF-Prelude.HOLCF_Main"}.
›

no_notation
  Rings.divide (infixl "div" 70) and
  Rings.modulo (infixl "mod" 70)

no_notation
  Set.member  ("(:)") and
  Set.member  ("(_/ : _)" [51, 51] 50)

text ‹This is the implementation. We also need a modulus operator.›
fixrec sieve :: "[Integer] → [Integer]" where
  "sieve⋅(p : xs) = p : (sieve⋅(filter⋅(Λ x. neg⋅(eq⋅(mod⋅x⋅p)⋅0))⋅xs))"

fixrec primes :: "[Integer]" where
  "primes = sieve⋅[2..]"

text ‹Simplification rules for modI:›

definition MkI' :: "int ⇒ Integer" where
  "MkI' x = MkI⋅x"

lemma MkI'_simps [simp]:
  shows "MkI' 0 = 0" and "MkI' 1 = 1" and "MkI' (numeral k) = numeral k"
  unfolding MkI'_def zero_Integer_def one_Integer_def numeral_Integer_eq
  by rule+

lemma modI_numeral_numeral [simp]:
  "mod⋅(numeral i)⋅(numeral j) = MkI' (Rings.modulo (numeral i) (numeral j))"
  unfolding numeral_Integer_eq mod_Integer_def MkI'_def by simp

text ‹Some lemmas demonstrating evaluation of our list:›

lemma "primes !! 0 = 2"
  unfolding primes.simps
  apply (simp only: enumFrom_intsFrom_conv)
  apply (subst intsFrom.simps)
  apply simp
  done

lemma "primes !! 1 = 3"
  unfolding primes.simps
  apply (simp only: enumFrom_intsFrom_conv)
  apply (subst intsFrom.simps)
  apply (subst intsFrom.simps)
  apply simp
  done

lemma "primes !! 2 = 5"
  unfolding primes.simps
  apply (simp only: enumFrom_intsFrom_conv)
  apply (subst intsFrom.simps)
  apply (subst intsFrom.simps)
  apply (subst intsFrom.simps)
  apply (subst intsFrom.simps)
  apply simp
  done

lemma "primes !! 3 = 7"
  unfolding primes.simps
  apply (simp only: enumFrom_intsFrom_conv)
  apply (subst intsFrom.simps)
  apply (subst intsFrom.simps)
  apply (subst intsFrom.simps)
  apply (subst intsFrom.simps)
  apply (subst intsFrom.simps)
  apply (subst intsFrom.simps)
  apply (simp del: filter_FF filter_TT)
  (* FIXME: remove these two from the default simpset *)
  done

text ‹Auxiliary lemmas about prime numbers›

lemma find_next_prime_nat:
  fixes n :: nat
  assumes "prime n"
  shows "∃ n'. n' > n ∧ prime n' ∧ (∀k. n < k ⟶ k < n' ⟶ ¬ prime k)"
  using ex_least_nat_le[of "λ k . k > n ∧ prime k"]
  by (metis bigger_prime not_prime_0)


text ‹Simplification for andalso:›

lemma andAlso_Def[simp]: "((Def x) andalso (Def y)) = Def (x ∧ y)"
  by (metis Def_bool2 Def_bool4 andalso_thms(1) andalso_thms(2))

text ‹This defines the bisimulation and proves it to be a list bisimulation.›

definition prim_bisim:
  "prim_bisim x1 x2 = (∃ n . prime n  ∧
    x1 = sieve⋅(filter⋅(Λ (MkI⋅i). Def ((∀d. d > 1 ⟶ d < n ⟶ ¬ (d dvd i))))⋅[MkI⋅n..]) ∧
    x2 = filter⋅(Λ (MkI⋅i). Def (prime (nat ¦i¦)))⋅[MkI⋅n..])"

lemma prim_bisim_is_bisim: "list_bisim prim_bisim"
proof -
  {
    fix xs ys
    assume "prim_bisim xs ys"
    then obtain n :: nat where
      "prime n" and
      "n > 1" and
      "xs = sieve⋅(filter⋅(Λ (MkI⋅i). Def ((∀d::int. d > 1 ⟶ d < n ⟶ ¬ (d dvd i))))⋅[MkI⋅n..])" (is "_ = sieve⋅?xs") and
      "ys = filter⋅(Λ (MkI⋅i). Def (prime (nat ¦i¦)))⋅[MkI⋅n..]"
    proof -
      assume a1: "⋀n. ⟦prime n; 1 < n; xs = sieve⋅ (Data_List.filter⋅ (Λ (MkI⋅i). Def (∀d>1. d < int n ⟶ ¬ d dvd i))⋅ [MkI⋅(int n)..]); ys = Data_List.filter⋅ (Λ (MkI⋅i). Def (prime (nat ¦i¦)))⋅ [MkI⋅(int n)..]⟧ ⟹ thesis"
      obtain ii where
        f2: "∀is isa. (¬ prim_bisim is isa ∨ prime (ii is isa) ∧ is = sieve⋅ (Data_List.filter⋅ (Λ (MkI⋅i). Def (∀ia>1. ia < ii is isa ⟶ ¬ ia dvd i))⋅ [MkI⋅(ii is isa)..]) ∧ isa = Data_List.filter⋅ (Λ (MkI⋅i). Def (prime (nat ¦i¦)))⋅ [MkI⋅(ii is isa)..]) ∧ ((∀i. ¬ prime i ∨ is ≠ sieve⋅ (Data_List.filter⋅ (Λ (MkI⋅ia). Def (∀ib>1. ib < i ⟶ ¬ ib dvd ia))⋅ [MkI⋅i..]) ∨ isa ≠ Data_List.filter⋅ (Λ (MkI⋅i). Def (prime (nat ¦i¦)))⋅ [MkI⋅i..]) ∨ prim_bisim is isa)"
        unfolding prim_bisim by (atomize_elim, (subst choice_iff[symmetric])+, blast)
      then have f3: "prime (ii xs ys)"
        using ‹prim_bisim xs ys› by blast
      then obtain nn :: "int ⇒ nat" where
        f4: "int (nn (ii xs ys)) = ii xs ys"
        by (metis (no_types) prime_gt_0_int zero_less_imp_eq_int)
      then have "prime (nn (ii xs ys))"
        using f3 by (metis (no_types) prime_nat_int_transfer)
      then show ?thesis
        using f4 f2 a1 ‹prim_bisim xs ys› prime_gt_1_nat by presburger
    qed

    obtain n' where "n' > n" and "prime n'" and not_prime: "∀k. n < k ⟶ k < n' ⟶ ¬ prime k"
      using find_next_prime_nat[OF ‹prime n›] by auto

    {
      fix k :: int
      assume "n < k" and "k < n'"

      have  "k > 1" using ‹n < k› ‹n > 1› by auto
      then obtain k' :: nat where "k = int k'" by (cases k) auto
      then obtain p where "prime p" and "p dvd k"
        using ‹k > 1› ‹k = int k'›
        by (metis (full_types) less_numeral_extra(4) of_nat_1 of_nat_dvd_iff prime_factor_nat prime_nat_int_transfer)  
      then have "p < n'" using ‹k < n'›  ‹k > 1›
        using zdvd_imp_le [of p k] by simp
      then have "p ≤ n" using ‹prime p› not_prime
        using not_le prime_gt_0_int zero_less_imp_eq_int
        by (metis of_nat_less_iff prime_nat_int_transfer) 
      then have "∃d::int>1. d ≤ n ∧ d dvd k" using ‹p dvd k› ‹prime p› of_nat_le_iff prime_gt_1_nat
          prime_gt_1_int by auto
    }
    then have between_have_divisors: "⋀k::int. n < k ⟹ k < n' ⟹ ∃d::int>1. d ≤ n ∧ d dvd k".

    {
      fix i
      {
        assume small: "∀d::int>1. d ≤ n ⟶ ¬ d dvd i"
        fix d
        assume "1 < d" and "d dvd i" and "d < n'"
        with small have "d > n" by auto

        obtain d'::int where "d' > 1" and "d' ≤ n" and "d' dvd d" using between_have_divisors[OF ‹n < d› ‹d < n'›] by auto
        with ‹d dvd i› small have False by (metis (full_types) dvd_trans)
      }
      then have "(∀d::int. d > 1 ⟶ d ≤ n ⟶ ¬ (d dvd i)) = (∀d::int. d > 1 ⟶ d < n' ⟶ ¬ (d dvd i))"
        using ‹n' > n› by auto
    }
    then have between_not_relvant:  "⋀ i. (∀d::int. d > 1 ⟶ d ≤ n ⟶ ¬ (d dvd i)) = (∀d::int. d > 1 ⟶ d < n' ⟶ ¬ (d dvd i))" .

    from ‹prime n›
    have "∀d::int >1. d < n ⟶ ¬ d dvd n"
      unfolding prime_int_altdef using int_one_le_iff_zero_less le_less
      by (simp add: prime_int_not_dvd)
    then
    obtain xs' where "?xs = (MkI⋅n) : xs'" and "xs' = filter⋅(Λ (MkI⋅i). Def ((∀d::int. d > 1 ⟶ d < n ⟶ ¬ (d dvd i))))⋅[MkI⋅(n+1)..]"
      by (subst (asm) intsFrom.simps[unfolded enumFrom_intsFrom_conv[symmetric]], simp add: one_Integer_def TT_def[symmetric] add.commute)

    {
      have "filter⋅(Λ x. neg⋅(eq⋅(mod⋅x⋅(MkI⋅n))⋅0))⋅(filter⋅(Λ (MkI⋅i). Def ((∀d::int. d > 1 ⟶ d < n ⟶ ¬ (d dvd i))))⋅[MkI⋅(n+1)..])
            = filter⋅(Λ (MkI⋅x). Def (¬ n dvd x))⋅(filter⋅(Λ (MkI⋅i). Def ((∀d::int. d > 1 ⟶ d < n ⟶ ¬ (d dvd i))))⋅[MkI⋅(n+1)..])"
        apply (rule filter_cong[rule_format])
        apply (rename_tac x, case_tac x)
         apply (simp)
        apply (auto simp add: zero_Integer_def)
         apply (rule FF_def)
        apply (simp add: TT_def)
        by (metis dvd_eq_mod_eq_0)
      also have "... = filter⋅(Λ (MkI⋅i). Def ((¬ n dvd i) ∧ (∀d::int. d > 1 ⟶ d < n ⟶ ¬ (d dvd i))))⋅[MkI⋅(n+1)..]"
        by (auto intro!: filter_cong[rule_format])
      also have "... = filter⋅(Λ (MkI⋅i). Def ((∀d::int. d > 1 ⟶ d ≤ n ⟶ ¬ (d dvd i))))⋅[MkI⋅(n+1)..]"
        apply (rule filter_cong[rule_format])
        apply (rename_tac x, case_tac x)
        using ‹n > 1›
         apply auto
        done
      also have "... = filter⋅(Λ (MkI⋅i). Def ((∀d::int. d > 1 ⟶ d ≤ n ⟶ ¬ (d dvd i))))⋅[MkI⋅(int n+1)..]"
        by (metis (no_types, lifting) of_nat_1 of_nat_add)
      also have "... = filter⋅(Λ (MkI⋅i). Def ((∀d::int. d > 1 ⟶ d ≤ n ⟶ ¬ (d dvd i))))⋅[MkI⋅n'..]"
        apply (rule filter_fast_forward[of n n'])  using ‹n' > n›
         apply (auto simp add: between_have_divisors)
        done
      also have "... = filter⋅(Λ (MkI⋅i). Def ((∀d::int. d > 1 ⟶ d < n' ⟶ ¬ (d dvd i))))⋅[MkI⋅n'..]"
        by (auto intro: filter_cong[rule_format] simp add: between_not_relvant)
      also note calculation
    }
    note tmp = this

    {
      have "xs = sieve⋅?xs" by fact
      also have "... = sieve⋅((MkI⋅n) : xs')" using ‹?xs = _› by simp
      also have "... = sieve⋅((MkI⋅n) : filter⋅(Λ (MkI⋅i). Def ((∀d::int. d > 1 ⟶ d < n ⟶ ¬ (d dvd i))))⋅[MkI⋅(n+1)..])" using ‹xs' = _› by simp
      also have "... = (MkI⋅n) : sieve⋅(filter⋅(Λ (MkI⋅i). Def ((∀d::int. d > 1 ⟶ d < n' ⟶ ¬ (d dvd i))))⋅[MkI⋅n'..])"  using tmp by simp
      also note calculation
    }
    moreover
    {
      have "ys = filter⋅(Λ (MkI⋅i). Def (prime (nat ¦i¦)))⋅[MkI⋅n..]" by fact
      also have "... = (MkI⋅n) : filter⋅(Λ (MkI⋅i). Def (prime (nat ¦i¦)))⋅[MkI⋅(int n+1)..]"
        using ‹prime n›
        by (subst intsFrom.simps[unfolded enumFrom_intsFrom_conv[symmetric]])(simp add: one_Integer_def TT_def[symmetric] add.commute)
      also have "... = (MkI⋅n) : filter⋅(Λ (MkI⋅i). Def (prime (nat ¦i¦)))⋅[MkI⋅n'..]"
        apply (subst filter_fast_forward[of n n']) using ‹n' > n› and not_prime
          apply auto
        apply (metis (full_types) ‹⋀k. ⟦int n < k; k < int n'⟧ ⟹ ∃d>1. d ≤ int n ∧ d dvd k› le_less not_le prime_gt_0_int prime_int_not_dvd zdvd_imp_le) 
        done
      also note calculation
    }
    moreover
    note ‹prime n'›
    ultimately
    have "∃ p xs' ys'. xs = p : xs' ∧ ys = p : ys' ∧ prim_bisim xs' ys'" unfolding prim_bisim
      using prime_nat_int_transfer by blast
  }
  then show ?thesis unfolding list.bisim_def by metis
qed

text ‹Now we apply coinduction:›

lemma sieve_produces_primes:
  fixes n :: nat
  assumes "prime n"
  shows "sieve⋅(filter⋅(Λ (MkI⋅i). Def ((∀d::int. d > 1 ⟶ d < n ⟶ ¬ (d dvd i))))⋅[MkI⋅n..])
    = filter⋅(Λ (MkI⋅i). Def (prime (nat ¦i¦)))⋅[MkI⋅n..]"
  using assms
  apply -
  apply (rule list.coinduct[OF prim_bisim_is_bisim], auto simp add: prim_bisim)
  using prime_nat_int_transfer by blast

text ‹And finally show the correctness of primes.›

theorem primes:
  shows "primes = filter⋅(Λ (MkI⋅i). Def (prime (nat ¦i¦)))⋅[MkI⋅2..]"
proof -
  have "primes = sieve⋅[2 ..]" by (rule primes.simps)
  also have "... = sieve⋅(filter⋅(Λ (MkI⋅i). Def ((∀d::int. d > 1 ⟶ d < (int 2) ⟶ ¬ (d dvd i))))⋅[MkI⋅(int 2)..])"
    unfolding numeral_Integer_eq
    by (subst filter_TT, auto)
  also have "... = filter⋅(Λ (MkI⋅i). Def (prime (nat ¦i¦)))⋅[MkI⋅(int 2)..]"
    by (rule sieve_produces_primes[OF two_is_prime_nat])
  also have "... = filter⋅(Λ (MkI⋅i). Def (prime (nat ¦i¦)))⋅[MkI⋅2..]"
    by simp
  finally show ?thesis .
qed

end