Sieve of Eratosthenes
authorhaftmann
Sun Feb 17 21:29:30 2013 +0100 (2013-02-17)
changeset 511733cbb4e95a565
parent 51172 16eb76ca1e4a
child 51174 071674018df9
Sieve of Eratosthenes
CONTRIBUTORS
src/HOL/Codegenerator_Test/Candidates.thy
src/HOL/Divides.thy
src/HOL/List.thy
src/HOL/Nat.thy
src/HOL/Number_Theory/Eratosthenes.thy
src/HOL/Number_Theory/Number_Theory.thy
src/HOL/Product_Type.thy
src/HOL/Set.thy
     1.1 --- a/CONTRIBUTORS	Sun Feb 17 20:45:49 2013 +0100
     1.2 +++ b/CONTRIBUTORS	Sun Feb 17 21:29:30 2013 +0100
     1.3 @@ -6,10 +6,13 @@
     1.4  Contributions to this Isabelle version
     1.5  --------------------------------------
     1.6  
     1.7 -* 2013: Florian Haftmann, TUM
     1.8 +* Feb. 2013: Florian Haftmann, TUM
     1.9    Reworking and consolidation of code generation for target
    1.10    language numerals.
    1.11  
    1.12 +* Feb. 2013: Florian Haftmann, TUM
    1.13 +  Sieve of Eratosthenes.
    1.14 +
    1.15  
    1.16  Contributions to Isabelle2013
    1.17  -----------------------------
     2.1 --- a/src/HOL/Codegenerator_Test/Candidates.thy	Sun Feb 17 20:45:49 2013 +0100
     2.2 +++ b/src/HOL/Codegenerator_Test/Candidates.thy	Sun Feb 17 21:29:30 2013 +0100
     2.3 @@ -8,7 +8,7 @@
     2.4    Complex_Main
     2.5    "~~/src/HOL/Library/Library"
     2.6    "~~/src/HOL/Library/Sublist_Order"
     2.7 -  "~~/src/HOL/Number_Theory/Primes"
     2.8 +  "~~/src/HOL/Number_Theory/Eratosthenes"
     2.9    "~~/src/HOL/ex/Records"
    2.10  begin
    2.11  
     3.1 --- a/src/HOL/Divides.thy	Sun Feb 17 20:45:49 2013 +0100
     3.2 +++ b/src/HOL/Divides.thy	Sun Feb 17 21:29:30 2013 +0100
     3.3 @@ -740,6 +740,10 @@
     3.4    shows "m mod n < (n::nat)"
     3.5    using assms divmod_nat_rel [of m n] unfolding divmod_nat_rel_def by auto
     3.6  
     3.7 +lemma mod_Suc_le_divisor [simp]:
     3.8 +  "m mod Suc n \<le> n"
     3.9 +  using mod_less_divisor [of "Suc n" m] by arith
    3.10 +
    3.11  lemma mod_less_eq_dividend [simp]:
    3.12    fixes m n :: nat
    3.13    shows "m mod n \<le> m"
     4.1 --- a/src/HOL/List.thy	Sun Feb 17 20:45:49 2013 +0100
     4.2 +++ b/src/HOL/List.thy	Sun Feb 17 21:29:30 2013 +0100
     4.3 @@ -196,6 +196,9 @@
     4.4  abbreviation length :: "'a list \<Rightarrow> nat" where
     4.5  "length \<equiv> size"
     4.6  
     4.7 +definition enumerate :: "nat \<Rightarrow> 'a list \<Rightarrow> (nat \<times> 'a) list" where
     4.8 +enumerate_eq_zip: "enumerate n xs = zip [n..<n + length xs] xs"
     4.9 +
    4.10  primrec rotate1 :: "'a list \<Rightarrow> 'a list" where
    4.11  "rotate1 [] = []" |
    4.12  "rotate1 (x # xs) = xs @ [x]"
    4.13 @@ -245,6 +248,7 @@
    4.14  @{lemma "foldl f x [a,b,c] = f (f (f x a) b) c" by simp}\\
    4.15  @{lemma "zip [a,b,c] [x,y,z] = [(a,x),(b,y),(c,z)]" by simp}\\
    4.16  @{lemma "zip [a,b] [x,y,z] = [(a,x),(b,y)]" by simp}\\
    4.17 +@{lemma "enumerate 3 [a,b,c] = [(3,a),(4,b),(5,c)]" by normalization}\\
    4.18  @{lemma "List.product [a,b] [c,d] = [(a, c), (a, d), (b, c), (b, d)]" by simp}\\
    4.19  @{lemma "splice [a,b,c] [x,y,z] = [a,x,b,y,c,z]" by simp}\\
    4.20  @{lemma "splice [a,b,c,d] [x,y] = [a,x,b,y,c,d]" by simp}\\
    4.21 @@ -2479,6 +2483,20 @@
    4.22    "length xs = length ys \<Longrightarrow> zip xs ys = zs \<longleftrightarrow> map fst zs = xs \<and> map snd zs = ys"
    4.23    by (auto simp add: zip_map_fst_snd)
    4.24  
    4.25 +lemma in_set_zip:
    4.26 +  "p \<in> set (zip xs ys) \<longleftrightarrow> (\<exists>n. xs ! n = fst p \<and> ys ! n = snd p
    4.27 +    \<and> n < length xs \<and> n < length ys)"
    4.28 +  by (cases p) (auto simp add: set_zip)
    4.29 +
    4.30 +lemma pair_list_eqI:
    4.31 +  assumes "map fst xs = map fst ys" and "map snd xs = map snd ys"
    4.32 +  shows "xs = ys"
    4.33 +proof -
    4.34 +  from assms(1) have "length xs = length ys" by (rule map_eq_imp_length_eq)
    4.35 +  from this assms show ?thesis
    4.36 +    by (induct xs ys rule: list_induct2) (simp_all add: prod_eqI)
    4.37 +qed
    4.38 +
    4.39  
    4.40  subsubsection {* @{const list_all2} *}
    4.41  
    4.42 @@ -3880,6 +3898,57 @@
    4.43  qed
    4.44  
    4.45  
    4.46 +subsubsection {* @{const enumerate} *}
    4.47 +
    4.48 +lemma enumerate_simps [simp, code]:
    4.49 +  "enumerate n [] = []"
    4.50 +  "enumerate n (x # xs) = (n, x) # enumerate (Suc n) xs"
    4.51 +  apply (auto simp add: enumerate_eq_zip not_le)
    4.52 +  apply (cases "n < n + length xs")
    4.53 +  apply (auto simp add: upt_conv_Cons)
    4.54 +  done
    4.55 +
    4.56 +lemma length_enumerate [simp]:
    4.57 +  "length (enumerate n xs) = length xs"
    4.58 +  by (simp add: enumerate_eq_zip)
    4.59 +
    4.60 +lemma map_fst_enumerate [simp]:
    4.61 +  "map fst (enumerate n xs) = [n..<n + length xs]"
    4.62 +  by (simp add: enumerate_eq_zip)
    4.63 +
    4.64 +lemma map_snd_enumerate [simp]:
    4.65 +  "map snd (enumerate n xs) = xs"
    4.66 +  by (simp add: enumerate_eq_zip)
    4.67 +  
    4.68 +lemma in_set_enumerate_eq:
    4.69 +  "p \<in> set (enumerate n xs) \<longleftrightarrow> n \<le> fst p \<and> fst p < length xs + n \<and> nth xs (fst p - n) = snd p"
    4.70 +proof -
    4.71 +  { fix m
    4.72 +    assume "n \<le> m"
    4.73 +    moreover assume "m < length xs + n"
    4.74 +    ultimately have "[n..<n + length xs] ! (m - n) = m \<and>
    4.75 +      xs ! (m - n) = xs ! (m - n) \<and> m - n < length xs" by auto
    4.76 +    then have "\<exists>q. [n..<n + length xs] ! q = m \<and>
    4.77 +        xs ! q = xs ! (m - n) \<and> q < length xs" ..
    4.78 +  } then show ?thesis by (cases p) (auto simp add: enumerate_eq_zip in_set_zip)
    4.79 +qed
    4.80 +
    4.81 +lemma nth_enumerate_eq:
    4.82 +  assumes "m < length xs"
    4.83 +  shows "enumerate n xs ! m = (n + m, xs ! m)"
    4.84 +  using assms by (simp add: enumerate_eq_zip)
    4.85 +
    4.86 +lemma enumerate_replicate_eq:
    4.87 +  "enumerate n (replicate m a) = map (\<lambda>q. (q, a)) [n..<n + m]"
    4.88 +  by (rule pair_list_eqI)
    4.89 +    (simp_all add: enumerate_eq_zip comp_def map_replicate_const)
    4.90 +
    4.91 +lemma enumerate_Suc_eq:
    4.92 +  "enumerate (Suc n) xs = map (apfst Suc) (enumerate n xs)"
    4.93 +  by (rule pair_list_eqI)
    4.94 +    (simp_all add: not_le, simp del: map_map [simp del] add: map_Suc_upt map_map [symmetric])
    4.95 +
    4.96 +
    4.97  subsubsection {* @{const rotate1} and @{const rotate} *}
    4.98  
    4.99  lemma rotate0[simp]: "rotate 0 = id"
     5.1 --- a/src/HOL/Nat.thy	Sun Feb 17 20:45:49 2013 +0100
     5.2 +++ b/src/HOL/Nat.thy	Sun Feb 17 21:29:30 2013 +0100
     5.3 @@ -1587,6 +1587,12 @@
     5.4  lemma less_diff_conv: "(i < j-k) = (i+k < (j::nat))"
     5.5  by arith
     5.6  
     5.7 +lemma less_diff_conv2:
     5.8 +  fixes j k i :: nat
     5.9 +  assumes "k \<le> j"
    5.10 +  shows "j - k < i \<longleftrightarrow> j < i + k"
    5.11 +  using assms by arith
    5.12 +
    5.13  lemma le_diff_conv: "(j-k \<le> (i::nat)) = (j \<le> i+k)"
    5.14  by arith
    5.15  
    5.16 @@ -1801,6 +1807,74 @@
    5.17    shows "0 < m \<Longrightarrow> m < n \<Longrightarrow> \<not> n dvd m"
    5.18  by (auto elim!: dvdE) (auto simp add: gr0_conv_Suc)
    5.19  
    5.20 +lemma dvd_plusE:
    5.21 +  fixes m n q :: nat
    5.22 +  assumes "m dvd n + q" "m dvd n"
    5.23 +  obtains "m dvd q"
    5.24 +proof (cases "m = 0")
    5.25 +  case True with assms that show thesis by simp
    5.26 +next
    5.27 +  case False then have "m > 0" by simp
    5.28 +  from assms obtain r s where "n = m * r" and "n + q = m * s" by (blast elim: dvdE)
    5.29 +  then have *: "m * r + q = m * s" by simp
    5.30 +  show thesis proof (cases "r \<le> s")
    5.31 +    case False then have "s < r" by (simp add: not_le)
    5.32 +    with * have "m * r + q - m * s = m * s - m * s" by simp
    5.33 +    then have "m * r + q - m * s = 0" by simp
    5.34 +    with `m > 0` `s < r` have "m * r - m * s + q = 0" by simp
    5.35 +    then have "m * (r - s) + q = 0" by auto
    5.36 +    then have "m * (r - s) = 0" by simp
    5.37 +    then have "m = 0 \<or> r - s = 0" by simp
    5.38 +    with `s < r` have "m = 0" by arith
    5.39 +    with `m > 0` show thesis by auto
    5.40 +  next
    5.41 +    case True with * have "m * r + q - m * r = m * s - m * r" by simp
    5.42 +    with `m > 0` `r \<le> s` have "m * r - m * r + q = m * s - m * r" by simp
    5.43 +    then have "q = m * (s - r)" by (simp add: diff_mult_distrib2)
    5.44 +    with assms that show thesis by (auto intro: dvdI)
    5.45 +  qed
    5.46 +qed
    5.47 +
    5.48 +lemma dvd_plus_eq_right:
    5.49 +  fixes m n q :: nat
    5.50 +  assumes "m dvd n"
    5.51 +  shows "m dvd n + q \<longleftrightarrow> m dvd q"
    5.52 +  using assms by (auto elim: dvd_plusE)
    5.53 +
    5.54 +lemma dvd_plus_eq_left:
    5.55 +  fixes m n q :: nat
    5.56 +  assumes "m dvd q"
    5.57 +  shows "m dvd n + q \<longleftrightarrow> m dvd n"
    5.58 +  using assms by (simp add: dvd_plus_eq_right add_commute [of n])
    5.59 +
    5.60 +lemma less_dvd_minus:
    5.61 +  fixes m n :: nat
    5.62 +  assumes "m < n"
    5.63 +  shows "m dvd n \<longleftrightarrow> m dvd (n - m)"
    5.64 +proof -
    5.65 +  from assms have "n = m + (n - m)" by arith
    5.66 +  then obtain q where "n = m + q" ..
    5.67 +  then show ?thesis by (simp add: dvd_reduce add_commute [of m])
    5.68 +qed
    5.69 +
    5.70 +lemma dvd_minus_self:
    5.71 +  fixes m n :: nat
    5.72 +  shows "m dvd n - m \<longleftrightarrow> n < m \<or> m dvd n"
    5.73 +  by (cases "n < m") (auto elim!: dvdE simp add: not_less le_imp_diff_is_add)
    5.74 +
    5.75 +lemma dvd_minus_add:
    5.76 +  fixes m n q r :: nat
    5.77 +  assumes "q \<le> n" "q \<le> r * m"
    5.78 +  shows "m dvd n - q \<longleftrightarrow> m dvd n + (r * m - q)"
    5.79 +proof -
    5.80 +  have "m dvd n - q \<longleftrightarrow> m dvd r * m + (n - q)"
    5.81 +    by (auto elim: dvd_plusE)
    5.82 +  also with assms have "\<dots> \<longleftrightarrow> m dvd r * m + n - q" by simp
    5.83 +  also with assms have "\<dots> \<longleftrightarrow> m dvd (r * m - q) + n" by simp
    5.84 +  also have "\<dots> \<longleftrightarrow> m dvd n + (r * m - q)" by (simp add: add_commute)
    5.85 +  finally show ?thesis .
    5.86 +qed
    5.87 +
    5.88  
    5.89  subsection {* aliasses *}
    5.90  
     6.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     6.2 +++ b/src/HOL/Number_Theory/Eratosthenes.thy	Sun Feb 17 21:29:30 2013 +0100
     6.3 @@ -0,0 +1,276 @@
     6.4 +(*  Title:      HOL/Number_Theory/Eratosthenes.thy
     6.5 +    Author:     Florian Haftmann, TU Muenchen
     6.6 +*)
     6.7 +
     6.8 +header {* The sieve of Eratosthenes *}
     6.9 +
    6.10 +theory Eratosthenes
    6.11 +imports Primes
    6.12 +begin
    6.13 +
    6.14 +subsection {* Preliminary: strict divisibility *}
    6.15 +
    6.16 +context dvd
    6.17 +begin
    6.18 +
    6.19 +abbreviation dvd_strict :: "'a \<Rightarrow> 'a \<Rightarrow> bool" (infixl "dvd'_strict" 50)
    6.20 +where
    6.21 +  "b dvd_strict a \<equiv> b dvd a \<and> \<not> a dvd b"
    6.22 +
    6.23 +end
    6.24 +
    6.25 +subsection {* Main corpus *}
    6.26 +
    6.27 +text {* The sieve is modelled as a list of booleans, where @{const False} means \emph{marked out}. *}
    6.28 +
    6.29 +type_synonym marks = "bool list"
    6.30 +
    6.31 +definition numbers_of_marks :: "nat \<Rightarrow> marks \<Rightarrow> nat set"
    6.32 +where
    6.33 +  "numbers_of_marks n bs = fst ` {x \<in> set (enumerate n bs). snd x}"
    6.34 +
    6.35 +lemma numbers_of_marks_simps [simp, code]:
    6.36 +  "numbers_of_marks n [] = {}"
    6.37 +  "numbers_of_marks n (True # bs) = insert n (numbers_of_marks (Suc n) bs)"
    6.38 +  "numbers_of_marks n (False # bs) = numbers_of_marks (Suc n) bs"
    6.39 +  by (auto simp add: numbers_of_marks_def intro!: image_eqI)
    6.40 +
    6.41 +lemma numbers_of_marks_Suc:
    6.42 +  "numbers_of_marks (Suc n) bs = Suc ` numbers_of_marks n bs"
    6.43 +  by (auto simp add: numbers_of_marks_def enumerate_Suc_eq image_iff Bex_def)
    6.44 +
    6.45 +lemma numbers_of_marks_replicate_False [simp]:
    6.46 +  "numbers_of_marks n (replicate m False) = {}"
    6.47 +  by (auto simp add: numbers_of_marks_def enumerate_replicate_eq)
    6.48 +
    6.49 +lemma numbers_of_marks_replicate_True [simp]:
    6.50 +  "numbers_of_marks n (replicate m True) = {n..<n+m}"
    6.51 +  by (auto simp add: numbers_of_marks_def enumerate_replicate_eq image_def)
    6.52 +
    6.53 +lemma in_numbers_of_marks_eq:
    6.54 +  "m \<in> numbers_of_marks n bs \<longleftrightarrow> m \<in> {n..<n + length bs} \<and> bs ! (m - n)"
    6.55 +  by (simp add: numbers_of_marks_def in_set_enumerate_eq image_iff add_commute)
    6.56 +
    6.57 +
    6.58 +text {* Marking out multiples in a sieve  *}
    6.59 + 
    6.60 +definition mark_out :: "nat \<Rightarrow> marks \<Rightarrow> marks"
    6.61 +where
    6.62 +  "mark_out n bs = map (\<lambda>(q, b). b \<and> \<not> Suc n dvd Suc (Suc q)) (enumerate n bs)"
    6.63 +
    6.64 +lemma mark_out_Nil [simp]:
    6.65 +  "mark_out n [] = []"
    6.66 +  by (simp add: mark_out_def)
    6.67 +  
    6.68 +lemma length_mark_out [simp]:
    6.69 +  "length (mark_out n bs) = length bs"
    6.70 +  by (simp add: mark_out_def)
    6.71 +
    6.72 +lemma numbers_of_marks_mark_out:
    6.73 +  "numbers_of_marks n (mark_out m bs) = {q \<in> numbers_of_marks n bs. \<not> Suc m dvd Suc q - n}"
    6.74 +  by (auto simp add: numbers_of_marks_def mark_out_def in_set_enumerate_eq image_iff
    6.75 +    nth_enumerate_eq less_dvd_minus)
    6.76 +
    6.77 +
    6.78 +text {* Auxiliary operation for efficient implementation  *}
    6.79 +
    6.80 +definition mark_out_aux :: "nat \<Rightarrow> nat \<Rightarrow> marks \<Rightarrow> marks"
    6.81 +where
    6.82 +  "mark_out_aux n m bs =
    6.83 +    map (\<lambda>(q, b). b \<and> (q < m + n \<or> \<not> Suc n dvd Suc (Suc q) + (n - m mod Suc n))) (enumerate n bs)"
    6.84 +
    6.85 +lemma mark_out_code [code]:
    6.86 +  "mark_out n bs = mark_out_aux n n bs"
    6.87 +proof -
    6.88 +  { fix a
    6.89 +    assume A: "Suc n dvd Suc (Suc a)"
    6.90 +      and B: "a < n + n"
    6.91 +      and C: "n \<le> a"
    6.92 +    have False
    6.93 +    proof (cases "n = 0")
    6.94 +      case True with A B C show False by simp
    6.95 +    next
    6.96 +      def m \<equiv> "Suc n" then have "m > 0" by simp
    6.97 +      case False then have "n > 0" by simp
    6.98 +      from A obtain q where q: "Suc (Suc a) = Suc n * q" by (rule dvdE)
    6.99 +      have "q > 0"
   6.100 +      proof (rule ccontr)
   6.101 +        assume "\<not> q > 0"
   6.102 +        with q show False by simp
   6.103 +      qed
   6.104 +      with `n > 0` have "Suc n * q \<ge> 2" by (auto simp add: gr0_conv_Suc)
   6.105 +      with q have a: "a = Suc n * q - 2" by simp
   6.106 +      with B have "q + n * q < n + n + 2"
   6.107 +        by auto
   6.108 +      then have "m * q < m * 2" by (simp add: m_def)
   6.109 +      with `m > 0` have "q < 2" by simp
   6.110 +      with `q > 0` have "q = 1" by simp
   6.111 +      with a have "a = n - 1" by simp
   6.112 +      with `n > 0` C show False by simp
   6.113 +    qed
   6.114 +  } note aux = this 
   6.115 +  show ?thesis
   6.116 +    by (auto simp add: mark_out_def mark_out_aux_def in_set_enumerate_eq intro: aux)
   6.117 +qed
   6.118 +
   6.119 +lemma mark_out_aux_simps [simp, code]:
   6.120 +  "mark_out_aux n m [] = []" (is ?thesis1)
   6.121 +  "mark_out_aux n 0 (b # bs) = False # mark_out_aux n n bs" (is ?thesis2)
   6.122 +  "mark_out_aux n (Suc m) (b # bs) = b # mark_out_aux n m bs" (is ?thesis3)
   6.123 +proof -
   6.124 +  show ?thesis1
   6.125 +    by (simp add: mark_out_aux_def)
   6.126 +  show ?thesis2
   6.127 +    by (auto simp add: mark_out_code [symmetric] mark_out_aux_def mark_out_def
   6.128 +      enumerate_Suc_eq in_set_enumerate_eq less_dvd_minus)
   6.129 +  { def v \<equiv> "Suc m" and w \<equiv> "Suc n"
   6.130 +    fix q
   6.131 +    assume "m + n \<le> q"
   6.132 +    then obtain r where q: "q = m + n + r" by (auto simp add: le_iff_add)
   6.133 +    { fix u
   6.134 +      from w_def have "u mod w < w" by simp
   6.135 +      then have "u + (w - u mod w) = w + (u - u mod w)"
   6.136 +        by simp
   6.137 +      then have "u + (w - u mod w) = w + u div w * w"
   6.138 +        by (simp add: div_mod_equality' [symmetric])
   6.139 +    }
   6.140 +    then have "w dvd v + w + r + (w - v mod w) \<longleftrightarrow> w dvd m + w + r + (w - m mod w)"
   6.141 +      by (simp add: add_assoc add_left_commute [of m] add_left_commute [of v]
   6.142 +        dvd_plus_eq_left dvd_plus_eq_right)
   6.143 +    moreover from q have "Suc q = m + w + r" by (simp add: w_def)
   6.144 +    moreover from q have "Suc (Suc q) = v + w + r" by (simp add: v_def w_def)
   6.145 +    ultimately have "w dvd Suc (Suc (q + (w - v mod w))) \<longleftrightarrow> w dvd Suc (q + (w - m mod w))"
   6.146 +      by (simp only: add_Suc [symmetric])
   6.147 +    then have "Suc n dvd Suc (Suc (Suc (q + n) - Suc m mod Suc n)) \<longleftrightarrow>
   6.148 +      Suc n dvd Suc (Suc (q + n - m mod Suc n))"
   6.149 +      by (simp add: v_def w_def Suc_diff_le trans_le_add2)
   6.150 +  }
   6.151 +  then show ?thesis3
   6.152 +    by (auto simp add: mark_out_aux_def
   6.153 +      enumerate_Suc_eq in_set_enumerate_eq not_less)
   6.154 +qed
   6.155 +
   6.156 +
   6.157 +text {* Main entry point to sieve *}
   6.158 +
   6.159 +fun sieve :: "nat \<Rightarrow> marks \<Rightarrow> marks"
   6.160 +where
   6.161 +  "sieve n [] = []"
   6.162 +| "sieve n (False # bs) = False # sieve (Suc n) bs"
   6.163 +| "sieve n (True # bs) = True # sieve (Suc n) (mark_out n bs)"
   6.164 +
   6.165 +text {*
   6.166 +  There are the following possible optimisations here:
   6.167 +
   6.168 +  \begin{itemize}
   6.169 +
   6.170 +    \item @{const sieve} can abort as soon as @{term n} is too big to let
   6.171 +      @{const mark_out} have any effect.
   6.172 +
   6.173 +    \item Search for further primes can be given up as soon as the search
   6.174 +      position exceeds the square root of the maximum candidate.
   6.175 +
   6.176 +  \end{itemize}
   6.177 +
   6.178 +  This is left as an constructive exercise to the reader.
   6.179 +*}
   6.180 +
   6.181 +lemma numbers_of_marks_sieve:
   6.182 +  "numbers_of_marks (Suc n) (sieve n bs) =
   6.183 +    {q \<in> numbers_of_marks (Suc n) bs. \<forall>m \<in> numbers_of_marks (Suc n) bs. \<not> m dvd_strict q}"
   6.184 +proof (induct n bs rule: sieve.induct)
   6.185 +  case 1 show ?case by simp
   6.186 +next
   6.187 +  case 2 then show ?case by simp
   6.188 +next
   6.189 +  case (3 n bs)
   6.190 +  have aux: "\<And>M n. n \<in> Suc ` M \<longleftrightarrow> n > 0 \<and> n - 1 \<in> M"
   6.191 +  proof
   6.192 +    fix M and n
   6.193 +    assume "n \<in> Suc ` M" then show "n > 0 \<and> n - 1 \<in> M" by auto
   6.194 +  next
   6.195 +    fix M and n :: nat
   6.196 +    assume "n > 0 \<and> n - 1 \<in> M"
   6.197 +    then have "n > 0" and "n - 1 \<in> M" by auto
   6.198 +    then have "Suc (n - 1) \<in> Suc ` M" by blast
   6.199 +    with `n > 0` show "n \<in> Suc ` M" by simp
   6.200 +  qed
   6.201 +  { fix m :: nat
   6.202 +    assume "Suc (Suc n) \<le> m" and "m dvd Suc n"
   6.203 +    from `m dvd Suc n` obtain q where "Suc n = m * q" ..
   6.204 +    with `Suc (Suc n) \<le> m` have "Suc (m * q) \<le> m" by simp
   6.205 +    then have "m * q < m" by arith
   6.206 +    then have "q = 0" by simp
   6.207 +    with `Suc n = m * q` have False by simp
   6.208 +  } note aux1 = this
   6.209 +  { fix m q :: nat
   6.210 +    assume "\<forall>q>0. 1 < q \<longrightarrow> Suc n < q \<longrightarrow> q \<le> Suc (n + length bs)
   6.211 +      \<longrightarrow> bs ! (q - Suc (Suc n)) \<longrightarrow> \<not> Suc n dvd q \<longrightarrow> q dvd m \<longrightarrow> m dvd q"
   6.212 +    then have *: "\<And>q. Suc n < q \<Longrightarrow> q \<le> Suc (n + length bs)
   6.213 +      \<Longrightarrow> bs ! (q - Suc (Suc n)) \<Longrightarrow> \<not> Suc n dvd q \<Longrightarrow> q dvd m \<Longrightarrow> m dvd q"
   6.214 +      by auto
   6.215 +    assume "\<not> Suc n dvd m" and "q dvd m"
   6.216 +    then have "\<not> Suc n dvd q" by (auto elim: dvdE)
   6.217 +    moreover assume "Suc n < q" and "q \<le> Suc (n + length bs)"
   6.218 +      and "bs ! (q - Suc (Suc n))"
   6.219 +    moreover note `q dvd m`
   6.220 +    ultimately have "m dvd q" by (auto intro: *)
   6.221 +  } note aux2 = this
   6.222 +  from 3 show ?case
   6.223 +    apply (simp_all add: numbers_of_marks_mark_out numbers_of_marks_Suc Compr_image_eq inj_image_eq_iff
   6.224 +      in_numbers_of_marks_eq Ball_def imp_conjL aux)
   6.225 +    apply safe
   6.226 +    apply (simp_all add: less_diff_conv2 le_diff_conv2 dvd_minus_self not_less)
   6.227 +    apply (clarsimp dest!: aux1)
   6.228 +    apply (simp add: Suc_le_eq less_Suc_eq_le)
   6.229 +    apply (rule aux2) apply (clarsimp dest!: aux1)+
   6.230 +    done
   6.231 +qed
   6.232 +
   6.233 +
   6.234 +text {* Relation the sieve algorithm to actual primes *}
   6.235 +
   6.236 +definition primes_upto :: "nat \<Rightarrow> nat set"
   6.237 +where
   6.238 +  "primes_upto n = {m. m \<le> n \<and> prime m}"
   6.239 +
   6.240 +lemma in_primes_upto:
   6.241 +  "m \<in> primes_upto n \<longleftrightarrow> m \<le> n \<and> prime m"
   6.242 +  by (simp add: primes_upto_def)
   6.243 +
   6.244 +lemma primes_upto_sieve [code]:
   6.245 +  "primes_upto n = numbers_of_marks 2 (sieve 1 (replicate (n - 1) True))"
   6.246 +proof (cases "n > 1")
   6.247 +  case False then have "n = 0 \<or> n = 1" by arith
   6.248 +  then show ?thesis
   6.249 +    by (auto simp add: numbers_of_marks_sieve One_nat_def numeral_2_eq_2 primes_upto_def dest: prime_gt_Suc_0_nat)
   6.250 +next
   6.251 +  { fix m q
   6.252 +    assume "Suc (Suc 0) \<le> q"
   6.253 +      and "q < Suc n"
   6.254 +      and "m dvd q"
   6.255 +    then have "m < Suc n" by (auto dest: dvd_imp_le)
   6.256 +    assume *: "\<forall>m\<in>{Suc (Suc 0)..<Suc n}. m dvd q \<longrightarrow> q dvd m"
   6.257 +      and "m dvd q" and "m \<noteq> 1"
   6.258 +    have "m = q" proof (cases "m = 0")
   6.259 +      case True with `m dvd q` show ?thesis by simp
   6.260 +    next
   6.261 +      case False with `m \<noteq> 1` have "Suc (Suc 0) \<le> m" by arith
   6.262 +      with `m < Suc n` * `m dvd q` have "q dvd m" by simp
   6.263 +      with `m dvd q` show ?thesis by (simp add: dvd.eq_iff)
   6.264 +    qed
   6.265 +  }
   6.266 +  then have aux: "\<And>m q. Suc (Suc 0) \<le> q \<Longrightarrow>
   6.267 +    q < Suc n \<Longrightarrow>
   6.268 +    m dvd q \<Longrightarrow>
   6.269 +    \<forall>m\<in>{Suc (Suc 0)..<Suc n}. m dvd q \<longrightarrow> q dvd m \<Longrightarrow>
   6.270 +    m dvd q \<Longrightarrow> m \<noteq> q \<Longrightarrow> m = 1" by auto
   6.271 +  case True then show ?thesis
   6.272 +    apply (auto simp add: numbers_of_marks_sieve One_nat_def numeral_2_eq_2 primes_upto_def dest: prime_gt_Suc_0_nat)
   6.273 +    apply (simp add: prime_nat_def dvd_def)
   6.274 +    apply (auto simp add: prime_nat_def aux)
   6.275 +    done
   6.276 +qed
   6.277 +
   6.278 +end
   6.279 +
     7.1 --- a/src/HOL/Number_Theory/Number_Theory.thy	Sun Feb 17 20:45:49 2013 +0100
     7.2 +++ b/src/HOL/Number_Theory/Number_Theory.thy	Sun Feb 17 21:29:30 2013 +0100
     7.3 @@ -2,7 +2,8 @@
     7.4  header {* Comprehensive number theory *}
     7.5  
     7.6  theory Number_Theory
     7.7 -imports Fib Residues
     7.8 +imports Fib Residues Eratosthenes
     7.9  begin
    7.10  
    7.11  end
    7.12 +
     8.1 --- a/src/HOL/Product_Type.thy	Sun Feb 17 20:45:49 2013 +0100
     8.2 +++ b/src/HOL/Product_Type.thy	Sun Feb 17 21:29:30 2013 +0100
     8.3 @@ -835,18 +835,34 @@
     8.4    "fst (apfst f x) = f (fst x)"
     8.5    by (cases x) simp
     8.6  
     8.7 +lemma fst_comp_apfst [simp]:
     8.8 +  "fst \<circ> apfst f = f \<circ> fst"
     8.9 +  by (simp add: fun_eq_iff)
    8.10 +
    8.11  lemma fst_apsnd [simp]:
    8.12    "fst (apsnd f x) = fst x"
    8.13    by (cases x) simp
    8.14  
    8.15 +lemma fst_comp_apsnd [simp]:
    8.16 +  "fst \<circ> apsnd f = fst"
    8.17 +  by (simp add: fun_eq_iff)
    8.18 +
    8.19  lemma snd_apfst [simp]:
    8.20    "snd (apfst f x) = snd x"
    8.21    by (cases x) simp
    8.22  
    8.23 +lemma snd_comp_apfst [simp]:
    8.24 +  "snd \<circ> apfst f = snd"
    8.25 +  by (simp add: fun_eq_iff)
    8.26 +
    8.27  lemma snd_apsnd [simp]:
    8.28    "snd (apsnd f x) = f (snd x)"
    8.29    by (cases x) simp
    8.30  
    8.31 +lemma snd_comp_apsnd [simp]:
    8.32 +  "snd \<circ> apsnd f = f \<circ> snd"
    8.33 +  by (simp add: fun_eq_iff)
    8.34 +
    8.35  lemma apfst_compose:
    8.36    "apfst f (apfst g x) = apfst (f \<circ> g) x"
    8.37    by (cases x) simp
     9.1 --- a/src/HOL/Set.thy	Sun Feb 17 20:45:49 2013 +0100
     9.2 +++ b/src/HOL/Set.thy	Sun Feb 17 21:29:30 2013 +0100
     9.3 @@ -908,6 +908,10 @@
     9.4    -- {* The eta-expansion gives variable-name preservation. *}
     9.5    by (unfold image_def) blast
     9.6  
     9.7 +lemma Compr_image_eq:
     9.8 +  "{x \<in> f ` A. P x} = f ` {x \<in> A. P (f x)}"
     9.9 +  by auto
    9.10 +
    9.11  lemma image_Un: "f`(A Un B) = f`A Un f`B"
    9.12    by blast
    9.13