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(modxp)0))xs))"

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

text ‹Simplification rules for modI:›

definition MkI' :: "int  Integer" where
  "MkI' x = MkIx"

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(Λ (MkIi). Def ((d. d > 1  d < n  ¬ (d dvd i))))[MkIn..]) 
    x2 = filter(Λ (MkIi). Def (prime (nat ¦i¦)))[MkIn..])"

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(Λ (MkIi). Def ((d::int. d > 1  d < n  ¬ (d dvd i))))[MkIn..])" (is "_ = sieve?xs") and
      "ys = filter(Λ (MkIi). Def (prime (nat ¦i¦)))[MkIn..]"
    proof -
      assume a1: "n. prime n; 1 < n; xs = sieve (Data_List.filter (Λ (MkIi). Def (d>1. d < int n  ¬ d dvd i)) [MkI(int n)..]); ys = Data_List.filter (Λ (MkIi). 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 (Λ (MkIi). Def (ia>1. ia < ii is isa  ¬ ia dvd i)) [MkI(ii is isa)..])  isa = Data_List.filter (Λ (MkIi). Def (prime (nat ¦i¦))) [MkI(ii is isa)..])  ((i. ¬ prime i  is  sieve (Data_List.filter (Λ (MkIia). Def (ib>1. ib < i  ¬ ib dvd ia)) [MkIi..])  isa  Data_List.filter (Λ (MkIi). Def (prime (nat ¦i¦))) [MkIi..])  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 = (MkIn) : xs'" and "xs' = filter(Λ (MkIi). 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(modx(MkIn))0))(filter(Λ (MkIi). Def ((d::int. d > 1  d < n  ¬ (d dvd i))))[MkI(n+1)..])
            = filter(Λ (MkIx). Def (¬ n dvd x))(filter(Λ (MkIi). 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(Λ (MkIi). 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(Λ (MkIi). 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(Λ (MkIi). 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(Λ (MkIi). Def ((d::int. d > 1  d  n  ¬ (d dvd i))))[MkIn'..]"
        apply (rule filter_fast_forward[of n n'])  using n' > n
         apply (auto simp add: between_have_divisors)
        done
      also have "... = filter(Λ (MkIi). Def ((d::int. d > 1  d < n'  ¬ (d dvd i))))[MkIn'..]"
        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((MkIn) : xs')" using ?xs = _ by simp
      also have "... = sieve((MkIn) : filter(Λ (MkIi). Def ((d::int. d > 1  d < n  ¬ (d dvd i))))[MkI(n+1)..])" using xs' = _ by simp
      also have "... = (MkIn) : sieve(filter(Λ (MkIi). Def ((d::int. d > 1  d < n'  ¬ (d dvd i))))[MkIn'..])"  using tmp by simp
      also note calculation
    }
    moreover
    {
      have "ys = filter(Λ (MkIi). Def (prime (nat ¦i¦)))[MkIn..]" by fact
      also have "... = (MkIn) : filter(Λ (MkIi). 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 "... = (MkIn) : filter(Λ (MkIi). Def (prime (nat ¦i¦)))[MkIn'..]"
        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(Λ (MkIi). Def ((d::int. d > 1  d < n  ¬ (d dvd i))))[MkIn..])
    = filter(Λ (MkIi). Def (prime (nat ¦i¦)))[MkIn..]"
  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(Λ (MkIi). Def (prime (nat ¦i¦)))[MkI2..]"
proof -
  have "primes = sieve[2 ..]" by (rule primes.simps)
  also have "... = sieve(filter(Λ (MkIi). 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(Λ (MkIi). Def (prime (nat ¦i¦)))[MkI(int 2)..]"
    by (rule sieve_produces_primes[OF two_is_prime_nat])
  also have "... = filter(Λ (MkIi). Def (prime (nat ¦i¦)))[MkI2..]"
    by simp
  finally show ?thesis .
qed

end