move Stirling numbers from AFP/Discrete_Summation
authorhoelzl
Wed May 04 10:19:01 2016 +0200 (2016-05-04)
changeset 630713ca3bc795908
parent 63070 952714a20087
child 63072 eb5d493a9e03
child 63073 413184c7a2a2
move Stirling numbers from AFP/Discrete_Summation
src/HOL/Library/Library.thy
src/HOL/Library/Stirling.thy
     1.1 --- a/src/HOL/Library/Library.thy	Sun May 01 17:26:27 2016 +0200
     1.2 +++ b/src/HOL/Library/Library.thy	Wed May 04 10:19:01 2016 +0200
     1.3 @@ -73,6 +73,7 @@
     1.4    Saturated
     1.5    Set_Algebras
     1.6    State_Monad
     1.7 +  Stirling
     1.8    Stream
     1.9    Sublist
    1.10    Sum_of_Squares
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/src/HOL/Library/Stirling.thy	Wed May 04 10:19:01 2016 +0200
     2.3 @@ -0,0 +1,281 @@
     2.4 +(* Authors: Amine Chaieb & Florian Haftmann, TU Muenchen
     2.5 +            with contributions by Lukas Bulwahn and Manuel Eberl*)
     2.6 +
     2.7 +section {* Stirling numbers of first and second kind *}
     2.8 +
     2.9 +theory Stirling
    2.10 +imports Binomial
    2.11 +begin
    2.12 +
    2.13 +subsection {* Stirling numbers of the second kind *}
    2.14 +
    2.15 +fun Stirling :: "nat \<Rightarrow> nat \<Rightarrow> nat"
    2.16 +where
    2.17 +  "Stirling 0 0 = 1"
    2.18 +| "Stirling 0 (Suc k) = 0"
    2.19 +| "Stirling (Suc n) 0 = 0"
    2.20 +| "Stirling (Suc n) (Suc k) = Suc k * Stirling n (Suc k) + Stirling n k"
    2.21 +
    2.22 +lemma Stirling_1 [simp]:
    2.23 +  "Stirling (Suc n) (Suc 0) = 1"
    2.24 +  by (induct n) simp_all
    2.25 +
    2.26 +lemma Stirling_less [simp]:
    2.27 +  "n < k \<Longrightarrow> Stirling n k = 0"
    2.28 +  by (induct n k rule: Stirling.induct) simp_all
    2.29 +
    2.30 +lemma Stirling_same [simp]:
    2.31 +  "Stirling n n = 1"
    2.32 +  by (induct n) simp_all
    2.33 +
    2.34 +lemma Stirling_2_2:
    2.35 +  "Stirling (Suc (Suc n)) (Suc (Suc 0)) = 2 ^ Suc n - 1"
    2.36 +proof (induct n)
    2.37 +  case 0 then show ?case by simp
    2.38 +next
    2.39 +  case (Suc n)
    2.40 +  have "Stirling (Suc (Suc (Suc n))) (Suc (Suc 0)) =
    2.41 +    2 * Stirling (Suc (Suc n)) (Suc (Suc 0)) + Stirling (Suc (Suc n)) (Suc 0)" by simp
    2.42 +  also have "\<dots> = 2 * (2 ^ Suc n - 1) + 1"
    2.43 +    by (simp only: Suc Stirling_1)
    2.44 +  also have "\<dots> = 2 ^ Suc (Suc n) - 1"
    2.45 +  proof -
    2.46 +    have "(2::nat) ^ Suc n - 1 > 0" by (induct n) simp_all
    2.47 +    then have "2 * ((2::nat) ^ Suc n - 1) > 0" by simp
    2.48 +    then have "2 \<le> 2 * ((2::nat) ^ Suc n)" by simp
    2.49 +    with add_diff_assoc2 [of 2 "2 * 2 ^ Suc n" 1]
    2.50 +      have "2 * 2 ^ Suc n - 2 + (1::nat) = 2 * 2 ^ Suc n + 1 - 2" .
    2.51 +    then show ?thesis by (simp add: nat_distrib)
    2.52 +  qed
    2.53 +  finally show ?case by simp
    2.54 +qed
    2.55 +
    2.56 +lemma Stirling_2:
    2.57 +  "Stirling (Suc n) (Suc (Suc 0)) = 2 ^ n - 1"
    2.58 +  using Stirling_2_2 by (cases n) simp_all
    2.59 +
    2.60 +subsection {* Stirling numbers of the first kind *}
    2.61 +
    2.62 +fun stirling :: "nat \<Rightarrow> nat \<Rightarrow> nat"
    2.63 +where
    2.64 +  "stirling 0 0 = 1"
    2.65 +| "stirling 0 (Suc k) = 0"
    2.66 +| "stirling (Suc n) 0 = 0"
    2.67 +| "stirling (Suc n) (Suc k) = n * stirling n (Suc k) + stirling n k"
    2.68 +
    2.69 +lemma stirling_0 [simp]: "n > 0 \<Longrightarrow> stirling n 0 = 0"
    2.70 +  by (cases n) simp_all
    2.71 +
    2.72 +lemma stirling_less [simp]:
    2.73 +  "n < k \<Longrightarrow> stirling n k = 0"
    2.74 +  by (induct n k rule: stirling.induct) simp_all
    2.75 +
    2.76 +lemma stirling_same [simp]:
    2.77 +  "stirling n n = 1"
    2.78 +  by (induct n) simp_all
    2.79 +
    2.80 +lemma stirling_Suc_n_1:
    2.81 +  "stirling (Suc n) (Suc 0) = fact n"
    2.82 +  by (induct n) auto
    2.83 +
    2.84 +lemma stirling_Suc_n_n:
    2.85 +  shows "stirling (Suc n) n = Suc n choose 2"
    2.86 +by (induct n) (auto simp add: numerals(2))
    2.87 +
    2.88 +lemma stirling_Suc_n_2:
    2.89 +  assumes "n \<ge> Suc 0"
    2.90 +  shows "stirling (Suc n) 2 = (\<Sum>k=1..n. fact n div k)"
    2.91 +using assms
    2.92 +proof (induct n)
    2.93 +  case 0 from this show ?case by simp
    2.94 +next
    2.95 +  case (Suc n)
    2.96 +  show ?case
    2.97 +  proof (cases n)
    2.98 +    case 0 from this show ?thesis by (simp add: numerals(2))
    2.99 +  next
   2.100 +    case Suc
   2.101 +    from this have geq1: "Suc 0 \<le> n" by simp
   2.102 +    have "stirling (Suc (Suc n)) 2 = Suc n * stirling (Suc n) 2 + stirling (Suc n) (Suc 0)"
   2.103 +      by (simp only: stirling.simps(4)[of "Suc n"] numerals(2))
   2.104 +    also have "... = Suc n * (\<Sum>k=1..n. fact n div k) + fact n"
   2.105 +      using Suc.hyps[OF geq1]
   2.106 +      by (simp only: stirling_Suc_n_1 of_nat_fact of_nat_add of_nat_mult)
   2.107 +    also have "... = Suc n * (\<Sum>k=1..n. fact n div k) + Suc n * fact n div Suc n"
   2.108 +      by (metis nat.distinct(1) nonzero_mult_divide_cancel_left)
   2.109 +    also have "... = (\<Sum>k=1..n. fact (Suc n) div k) + fact (Suc n) div Suc n"
   2.110 +      by (simp add: setsum_right_distrib div_mult_swap dvd_fact)
   2.111 +    also have "... = (\<Sum>k=1..Suc n. fact (Suc n) div k)" by simp
   2.112 +    finally show ?thesis .
   2.113 +  qed
   2.114 +qed
   2.115 +
   2.116 +lemma of_nat_stirling_Suc_n_2:
   2.117 +  assumes "n \<ge> Suc 0"
   2.118 +  shows "(of_nat (stirling (Suc n) 2)::'a::field_char_0) = fact n * (\<Sum>k=1..n. (1 / of_nat k))"
   2.119 +using assms
   2.120 +proof (induct n)
   2.121 +  case 0 from this show ?case by simp
   2.122 +next
   2.123 +  case (Suc n)
   2.124 +  show ?case
   2.125 +  proof (cases n)
   2.126 +    case 0 from this show ?thesis by (auto simp add: numerals(2))
   2.127 +  next
   2.128 +    case Suc
   2.129 +    from this have geq1: "Suc 0 \<le> n" by simp
   2.130 +    have "(of_nat (stirling (Suc (Suc n)) 2)::'a) =
   2.131 +      of_nat (Suc n * stirling (Suc n) 2 + stirling (Suc n) (Suc 0))"
   2.132 +      by (simp only: stirling.simps(4)[of "Suc n"] numerals(2))
   2.133 +    also have "... = of_nat (Suc n) * (fact n * (\<Sum>k = 1..n. 1 / of_nat k)) + fact n"
   2.134 +      using Suc.hyps[OF geq1]
   2.135 +      by (simp only: stirling_Suc_n_1 of_nat_fact of_nat_add of_nat_mult)
   2.136 +    also have "... = fact (Suc n) * (\<Sum>k = 1..n. 1 / of_nat k) + fact (Suc n) * (1 / of_nat (Suc n))"
   2.137 +      using of_nat_neq_0 by auto
   2.138 +    also have "... = fact (Suc n) * (\<Sum>k = 1..Suc n. 1 / of_nat k)"
   2.139 +      by (simp add: distrib_left)
   2.140 +    finally show ?thesis .
   2.141 +  qed
   2.142 +qed
   2.143 +
   2.144 +lemma setsum_stirling:
   2.145 +  "(\<Sum>k\<le>n. stirling n k) = fact n"
   2.146 +proof (induct n)
   2.147 +  case 0
   2.148 +  from this show ?case by simp
   2.149 +next
   2.150 +  case (Suc n)
   2.151 +  have "(\<Sum>k\<le>Suc n. stirling (Suc n) k) = stirling (Suc n) 0 + (\<Sum>k\<le>n. stirling (Suc n) (Suc k))"
   2.152 +    by (simp only: setsum_atMost_Suc_shift)
   2.153 +  also have "\<dots> = (\<Sum>k\<le>n. stirling (Suc n) (Suc k))" by simp
   2.154 +  also have "\<dots> = (\<Sum>k\<le>n. n * stirling n (Suc k) + stirling n k)" by simp
   2.155 +  also have "\<dots> = n * (\<Sum>k\<le>n. stirling n (Suc k)) + (\<Sum>k\<le>n. stirling n k)"
   2.156 +    by (simp add: setsum.distrib setsum_right_distrib)
   2.157 +  also have "\<dots> = n * fact n + fact n"
   2.158 +  proof -
   2.159 +    have "n * (\<Sum>k\<le>n. stirling n (Suc k)) = n * ((\<Sum>k\<le>Suc n. stirling n k) - stirling n 0)"
   2.160 +      by (metis add_diff_cancel_left' setsum_atMost_Suc_shift)
   2.161 +    also have "\<dots> = n * (\<Sum>k\<le>n. stirling n k)" by (cases n) simp+
   2.162 +    also have "\<dots> = n * fact n" using Suc.hyps by simp
   2.163 +    finally have "n * (\<Sum>k\<le>n. stirling n (Suc k)) = n * fact n" .
   2.164 +    moreover have "(\<Sum>k\<le>n. stirling n k) = fact n" using Suc.hyps .
   2.165 +    ultimately show ?thesis by simp
   2.166 +  qed
   2.167 +  also have "\<dots> = fact (Suc n)" by simp
   2.168 +  finally show ?case .
   2.169 +qed
   2.170 +
   2.171 +lemma stirling_pochhammer:
   2.172 +  "(\<Sum>k\<le>n. of_nat (stirling n k) * x ^ k) = (pochhammer x n :: 'a :: comm_semiring_1)"
   2.173 +proof (induction n)
   2.174 +  case (Suc n)
   2.175 +  have "of_nat (n * stirling n 0) = (0 :: 'a)" by (cases n) simp_all
   2.176 +  hence "(\<Sum>k\<le>Suc n. of_nat (stirling (Suc n) k) * x ^ k) =
   2.177 +            (of_nat (n * stirling n 0) * x ^ 0 +
   2.178 +            (\<Sum>i\<le>n. of_nat (n * stirling n (Suc i)) * (x ^ Suc i))) +
   2.179 +            (\<Sum>i\<le>n. of_nat (stirling n i) * (x ^ Suc i))"
   2.180 +    by (subst setsum_atMost_Suc_shift) (simp add: setsum.distrib ring_distribs)
   2.181 +  also have "\<dots> = pochhammer x (Suc n)"
   2.182 +    by (subst setsum_atMost_Suc_shift [symmetric])
   2.183 +       (simp add: algebra_simps setsum.distrib setsum_right_distrib pochhammer_Suc Suc [symmetric])
   2.184 +  finally show ?case .
   2.185 +qed simp_all
   2.186 +
   2.187 +
   2.188 +text \<open>A row of the Stirling number triangle\<close>
   2.189 +
   2.190 +definition stirling_row :: "nat \<Rightarrow> nat list" where
   2.191 +  "stirling_row n = [stirling n k. k \<leftarrow> [0..<Suc n]]"
   2.192 +
   2.193 +lemma nth_stirling_row: "k \<le> n \<Longrightarrow> stirling_row n ! k = stirling n k"
   2.194 +  by (simp add: stirling_row_def del: upt_Suc)
   2.195 +
   2.196 +lemma length_stirling_row [simp]: "length (stirling_row n) = Suc n"
   2.197 +  by (simp add: stirling_row_def)
   2.198 +
   2.199 +lemma stirling_row_nonempty [simp]: "stirling_row n \<noteq> []"
   2.200 +  using length_stirling_row[of n] by (auto simp del: length_stirling_row)
   2.201 +
   2.202 +(* TODO Move *)
   2.203 +lemma list_ext:
   2.204 +  assumes "length xs = length ys"
   2.205 +  assumes "\<And>i. i < length xs \<Longrightarrow> xs ! i = ys ! i"
   2.206 +  shows   "xs = ys"
   2.207 +using assms
   2.208 +proof (induction rule: list_induct2)
   2.209 +  case (Cons x xs y ys)
   2.210 +  from Cons.prems[of 0] have "x = y" by simp
   2.211 +  moreover from Cons.prems[of "Suc i" for i] have "xs = ys" by (intro Cons.IH) simp
   2.212 +  ultimately show ?case by simp
   2.213 +qed simp_all
   2.214 +
   2.215 +
   2.216 +subsubsection \<open>Efficient code\<close>
   2.217 +
   2.218 +text \<open>
   2.219 +  Naively using the defining equations of the Stirling numbers of the first kind to
   2.220 +  compute them leads to exponential run time due to repeated computations.
   2.221 +  We can use memoisation to compute them row by row without repeating computations, at
   2.222 +  the cost of computing a few unneeded values.
   2.223 +
   2.224 +  As a bonus, this is very efficient for applications where an entire row of
   2.225 +  Stirling numbers is needed..
   2.226 +\<close>
   2.227 +
   2.228 +definition zip_with_prev :: "('a \<Rightarrow> 'a \<Rightarrow> 'b) \<Rightarrow> 'a \<Rightarrow> 'a list \<Rightarrow> 'b list" where
   2.229 +  "zip_with_prev f x xs = map (\<lambda>(x,y). f x y) (zip (x # xs) xs)"
   2.230 +
   2.231 +lemma zip_with_prev_altdef:
   2.232 +  "zip_with_prev f x xs =
   2.233 +     (if xs = [] then [] else f x (hd xs) # [f (xs!i) (xs!(i+1)). i \<leftarrow> [0..<length xs - 1]])"
   2.234 +proof (cases xs)
   2.235 +  case (Cons y ys)
   2.236 +  hence "zip_with_prev f x xs = f x (hd xs) # zip_with_prev f y ys"
   2.237 +    by (simp add: zip_with_prev_def)
   2.238 +  also have "zip_with_prev f y ys = map (\<lambda>i. f (xs ! i) (xs ! (i + 1))) [0..<length xs - 1]"
   2.239 +    unfolding Cons
   2.240 +    by (induction ys arbitrary: y)
   2.241 +       (simp_all add: zip_with_prev_def upt_conv_Cons map_Suc_upt [symmetric] del: upt_Suc)
   2.242 +  finally show ?thesis using Cons by simp
   2.243 +qed (simp add: zip_with_prev_def)
   2.244 +
   2.245 +
   2.246 +primrec stirling_row_aux where
   2.247 +  "stirling_row_aux n y [] = [1]"
   2.248 +| "stirling_row_aux n y (x#xs) = (y + n * x) # stirling_row_aux n x xs"
   2.249 +
   2.250 +lemma stirling_row_aux_correct:
   2.251 +  "stirling_row_aux n y xs = zip_with_prev (\<lambda>a b. a + n * b) y xs @ [1]"
   2.252 +  by (induction xs arbitrary: y) (simp_all add: zip_with_prev_def)
   2.253 +
   2.254 +lemma stirling_row_code [code]:
   2.255 +  "stirling_row 0 = [1]"
   2.256 +  "stirling_row (Suc n) = stirling_row_aux n 0 (stirling_row n)"
   2.257 +proof -
   2.258 +  have "stirling_row (Suc n) =
   2.259 +          0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ [1]"
   2.260 +  proof (rule list_ext, goal_cases length nth)
   2.261 +    case (nth i)
   2.262 +    from nth have "i \<le> Suc n" by simp
   2.263 +    then consider "i = 0" | j where "i > 0" "i \<le> n" | "i = Suc n" by linarith
   2.264 +    thus ?case
   2.265 +    proof cases
   2.266 +      assume i: "i > 0" "i \<le> n"
   2.267 +      from this show ?thesis by (cases i) (simp_all add: nth_append nth_stirling_row)
   2.268 +    qed (simp_all add: nth_stirling_row nth_append)
   2.269 +  qed simp
   2.270 +  also have "0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ [1] =
   2.271 +               zip_with_prev (\<lambda>a b. a + n * b) 0 (stirling_row n) @ [1]"
   2.272 +    by (cases n) (auto simp add: zip_with_prev_altdef stirling_row_def hd_map simp del: upt_Suc)
   2.273 +  also have "\<dots> = stirling_row_aux n 0 (stirling_row n)"
   2.274 +    by (simp add: stirling_row_aux_correct)
   2.275 +  finally show "stirling_row (Suc n) = stirling_row_aux n 0 (stirling_row n)" .
   2.276 +qed (simp add: stirling_row_def)
   2.277 +
   2.278 +lemma stirling_code [code]:
   2.279 +  "stirling n k = (if k = 0 then if n = 0 then 1 else 0
   2.280 +                   else if k > n then 0 else if k = n then 1
   2.281 +                   else stirling_row n ! k)"
   2.282 +  by (simp add: nth_stirling_row)
   2.283 +
   2.284 +end