src/HOL/Library/Stirling.thy
 author nipkow Tue Sep 12 20:40:46 2017 +0200 (2017-09-12) changeset 66655 e9be3d6995f9 parent 65552 f533820e7248 child 66656 8f4d252ce2fe permissions -rw-r--r--
introduced zip_with
1 (*  Title:      HOL/Library/Stirling.thy
2     Author:     Amine Chaieb
3     Author:     Florian Haftmann
4     Author:     Lukas Bulwahn
5     Author:     Manuel Eberl
6 *)
8 section \<open>Stirling numbers of first and second kind\<close>
10 theory Stirling
11 imports Main
12 begin
14 subsection \<open>Stirling numbers of the second kind\<close>
16 fun Stirling :: "nat \<Rightarrow> nat \<Rightarrow> nat"
17   where
18     "Stirling 0 0 = 1"
19   | "Stirling 0 (Suc k) = 0"
20   | "Stirling (Suc n) 0 = 0"
21   | "Stirling (Suc n) (Suc k) = Suc k * Stirling n (Suc k) + Stirling n k"
23 lemma Stirling_1 [simp]: "Stirling (Suc n) (Suc 0) = 1"
24   by (induct n) simp_all
26 lemma Stirling_less [simp]: "n < k \<Longrightarrow> Stirling n k = 0"
27   by (induct n k rule: Stirling.induct) simp_all
29 lemma Stirling_same [simp]: "Stirling n n = 1"
30   by (induct n) simp_all
32 lemma Stirling_2_2: "Stirling (Suc (Suc n)) (Suc (Suc 0)) = 2 ^ Suc n - 1"
33 proof (induct n)
34   case 0
35   then show ?case by simp
36 next
37   case (Suc n)
38   have "Stirling (Suc (Suc (Suc n))) (Suc (Suc 0)) =
39       2 * Stirling (Suc (Suc n)) (Suc (Suc 0)) + Stirling (Suc (Suc n)) (Suc 0)"
40     by simp
41   also have "\<dots> = 2 * (2 ^ Suc n - 1) + 1"
42     by (simp only: Suc Stirling_1)
43   also have "\<dots> = 2 ^ Suc (Suc n) - 1"
44   proof -
45     have "(2::nat) ^ Suc n - 1 > 0"
46       by (induct n) simp_all
47     then have "2 * ((2::nat) ^ Suc n - 1) > 0"
48       by simp
49     then have "2 \<le> 2 * ((2::nat) ^ Suc n)"
50       by simp
51     with add_diff_assoc2 [of 2 "2 * 2 ^ Suc n" 1]
52     have "2 * 2 ^ Suc n - 2 + (1::nat) = 2 * 2 ^ Suc n + 1 - 2" .
53     then show ?thesis
54       by (simp add: nat_distrib)
55   qed
56   finally show ?case by simp
57 qed
59 lemma Stirling_2: "Stirling (Suc n) (Suc (Suc 0)) = 2 ^ n - 1"
60   using Stirling_2_2 by (cases n) simp_all
63 subsection \<open>Stirling numbers of the first kind\<close>
65 fun stirling :: "nat \<Rightarrow> nat \<Rightarrow> nat"
66   where
67     "stirling 0 0 = 1"
68   | "stirling 0 (Suc k) = 0"
69   | "stirling (Suc n) 0 = 0"
70   | "stirling (Suc n) (Suc k) = n * stirling n (Suc k) + stirling n k"
72 lemma stirling_0 [simp]: "n > 0 \<Longrightarrow> stirling n 0 = 0"
73   by (cases n) simp_all
75 lemma stirling_less [simp]: "n < k \<Longrightarrow> stirling n k = 0"
76   by (induct n k rule: stirling.induct) simp_all
78 lemma stirling_same [simp]: "stirling n n = 1"
79   by (induct n) simp_all
81 lemma stirling_Suc_n_1: "stirling (Suc n) (Suc 0) = fact n"
82   by (induct n) auto
84 lemma stirling_Suc_n_n: "stirling (Suc n) n = Suc n choose 2"
85   by (induct n) (auto simp add: numerals(2))
87 lemma stirling_Suc_n_2:
88   assumes "n \<ge> Suc 0"
89   shows "stirling (Suc n) 2 = (\<Sum>k=1..n. fact n div k)"
90   using assms
91 proof (induct n)
92   case 0
93   then show ?case by simp
94 next
95   case (Suc n)
96   show ?case
97   proof (cases n)
98     case 0
99     then show ?thesis
100       by (simp add: numerals(2))
101   next
102     case Suc
103     then have geq1: "Suc 0 \<le> n"
104       by simp
105     have "stirling (Suc (Suc n)) 2 = Suc n * stirling (Suc n) 2 + stirling (Suc n) (Suc 0)"
106       by (simp only: stirling.simps(4)[of "Suc n"] numerals(2))
107     also have "\<dots> = Suc n * (\<Sum>k=1..n. fact n div k) + fact n"
108       using Suc.hyps[OF geq1]
109       by (simp only: stirling_Suc_n_1 of_nat_fact of_nat_add of_nat_mult)
110     also have "\<dots> = Suc n * (\<Sum>k=1..n. fact n div k) + Suc n * fact n div Suc n"
111       by (metis nat.distinct(1) nonzero_mult_div_cancel_left)
112     also have "\<dots> = (\<Sum>k=1..n. fact (Suc n) div k) + fact (Suc n) div Suc n"
113       by (simp add: sum_distrib_left div_mult_swap dvd_fact)
114     also have "\<dots> = (\<Sum>k=1..Suc n. fact (Suc n) div k)"
115       by simp
116     finally show ?thesis .
117   qed
118 qed
120 lemma of_nat_stirling_Suc_n_2:
121   assumes "n \<ge> Suc 0"
122   shows "(of_nat (stirling (Suc n) 2)::'a::field_char_0) = fact n * (\<Sum>k=1..n. (1 / of_nat k))"
123   using assms
124 proof (induct n)
125   case 0
126   then show ?case by simp
127 next
128   case (Suc n)
129   show ?case
130   proof (cases n)
131     case 0
132     then show ?thesis
133       by (auto simp add: numerals(2))
134   next
135     case Suc
136     then have geq1: "Suc 0 \<le> n"
137       by simp
138     have "(of_nat (stirling (Suc (Suc n)) 2)::'a) =
139         of_nat (Suc n * stirling (Suc n) 2 + stirling (Suc n) (Suc 0))"
140       by (simp only: stirling.simps(4)[of "Suc n"] numerals(2))
141     also have "\<dots> = of_nat (Suc n) * (fact n * (\<Sum>k = 1..n. 1 / of_nat k)) + fact n"
142       using Suc.hyps[OF geq1]
143       by (simp only: stirling_Suc_n_1 of_nat_fact of_nat_add of_nat_mult)
144     also have "\<dots> = fact (Suc n) * (\<Sum>k = 1..n. 1 / of_nat k) + fact (Suc n) * (1 / of_nat (Suc n))"
145       using of_nat_neq_0 by auto
146     also have "\<dots> = fact (Suc n) * (\<Sum>k = 1..Suc n. 1 / of_nat k)"
147       by (simp add: distrib_left)
148     finally show ?thesis .
149   qed
150 qed
152 lemma sum_stirling: "(\<Sum>k\<le>n. stirling n k) = fact n"
153 proof (induct n)
154   case 0
155   then show ?case by simp
156 next
157   case (Suc n)
158   have "(\<Sum>k\<le>Suc n. stirling (Suc n) k) = stirling (Suc n) 0 + (\<Sum>k\<le>n. stirling (Suc n) (Suc k))"
159     by (simp only: sum_atMost_Suc_shift)
160   also have "\<dots> = (\<Sum>k\<le>n. stirling (Suc n) (Suc k))"
161     by simp
162   also have "\<dots> = (\<Sum>k\<le>n. n * stirling n (Suc k) + stirling n k)"
163     by simp
164   also have "\<dots> = n * (\<Sum>k\<le>n. stirling n (Suc k)) + (\<Sum>k\<le>n. stirling n k)"
165     by (simp add: sum.distrib sum_distrib_left)
166   also have "\<dots> = n * fact n + fact n"
167   proof -
168     have "n * (\<Sum>k\<le>n. stirling n (Suc k)) = n * ((\<Sum>k\<le>Suc n. stirling n k) - stirling n 0)"
169       by (metis add_diff_cancel_left' sum_atMost_Suc_shift)
170     also have "\<dots> = n * (\<Sum>k\<le>n. stirling n k)"
171       by (cases n) simp_all
172     also have "\<dots> = n * fact n"
173       using Suc.hyps by simp
174     finally have "n * (\<Sum>k\<le>n. stirling n (Suc k)) = n * fact n" .
175     moreover have "(\<Sum>k\<le>n. stirling n k) = fact n"
176       using Suc.hyps .
177     ultimately show ?thesis by simp
178   qed
179   also have "\<dots> = fact (Suc n)" by simp
180   finally show ?case .
181 qed
183 lemma stirling_pochhammer:
184   "(\<Sum>k\<le>n. of_nat (stirling n k) * x ^ k) = (pochhammer x n :: 'a::comm_semiring_1)"
185 proof (induct n)
186   case 0
187   then show ?case by simp
188 next
189   case (Suc n)
190   have "of_nat (n * stirling n 0) = (0 :: 'a)" by (cases n) simp_all
191   then have "(\<Sum>k\<le>Suc n. of_nat (stirling (Suc n) k) * x ^ k) =
192       (of_nat (n * stirling n 0) * x ^ 0 +
193       (\<Sum>i\<le>n. of_nat (n * stirling n (Suc i)) * (x ^ Suc i))) +
194       (\<Sum>i\<le>n. of_nat (stirling n i) * (x ^ Suc i))"
195     by (subst sum_atMost_Suc_shift) (simp add: sum.distrib ring_distribs)
196   also have "\<dots> = pochhammer x (Suc n)"
197     by (subst sum_atMost_Suc_shift [symmetric])
198       (simp add: algebra_simps sum.distrib sum_distrib_left pochhammer_Suc Suc [symmetric])
199   finally show ?case .
200 qed
203 text \<open>A row of the Stirling number triangle\<close>
205 definition stirling_row :: "nat \<Rightarrow> nat list"
206   where "stirling_row n = [stirling n k. k \<leftarrow> [0..<Suc n]]"
208 lemma nth_stirling_row: "k \<le> n \<Longrightarrow> stirling_row n ! k = stirling n k"
209   by (simp add: stirling_row_def del: upt_Suc)
211 lemma length_stirling_row [simp]: "length (stirling_row n) = Suc n"
212   by (simp add: stirling_row_def)
214 lemma stirling_row_nonempty [simp]: "stirling_row n \<noteq> []"
215   using length_stirling_row[of n] by (auto simp del: length_stirling_row)
217 (* TODO Move *)
218 lemma list_ext:
219   assumes "length xs = length ys"
220   assumes "\<And>i. i < length xs \<Longrightarrow> xs ! i = ys ! i"
221   shows "xs = ys"
222   using assms
223 proof (induction rule: list_induct2)
224   case Nil
225   then show ?case by simp
226 next
227   case (Cons x xs y ys)
228   from Cons.prems[of 0] have "x = y"
229     by simp
230   moreover from Cons.prems[of "Suc i" for i] have "xs = ys"
231     by (intro Cons.IH) simp
232   ultimately show ?case by simp
233 qed
236 subsubsection \<open>Efficient code\<close>
238 text \<open>
239   Naively using the defining equations of the Stirling numbers of the first
240   kind to compute them leads to exponential run time due to repeated
241   computations. We can use memoisation to compute them row by row without
242   repeating computations, at the cost of computing a few unneeded values.
244   As a bonus, this is very efficient for applications where an entire row of
245   Stirling numbers is needed.
246 \<close>
248 definition zip_with_prev :: "('a \<Rightarrow> 'a \<Rightarrow> 'b) \<Rightarrow> 'a \<Rightarrow> 'a list \<Rightarrow> 'b list"
249   where "zip_with_prev f x xs = zip_with f (x # xs) xs"
251 lemma zip_with_prev_altdef:
252   "zip_with_prev f x xs =
253     (if xs = [] then [] else f x (hd xs) # [f (xs!i) (xs!(i+1)). i \<leftarrow> [0..<length xs - 1]])"
254 proof (cases xs)
255   case Nil
256   then show ?thesis
257     by (simp add: zip_with_prev_def)
258 next
259   case (Cons y ys)
260   then have "zip_with_prev f x xs = f x (hd xs) # zip_with_prev f y ys"
261     by (simp add: zip_with_prev_def)
262   also have "zip_with_prev f y ys = map (\<lambda>i. f (xs ! i) (xs ! (i + 1))) [0..<length xs - 1]"
263     unfolding Cons
264     by (induct ys arbitrary: y)
265       (simp_all add: zip_with_prev_def upt_conv_Cons map_Suc_upt [symmetric] del: upt_Suc)
266   finally show ?thesis
267     using Cons by simp
268 qed
271 primrec stirling_row_aux
272   where
273     "stirling_row_aux n y [] = "
274   | "stirling_row_aux n y (x#xs) = (y + n * x) # stirling_row_aux n x xs"
276 lemma stirling_row_aux_correct:
277   "stirling_row_aux n y xs = zip_with_prev (\<lambda>a b. a + n * b) y xs @ "
278   by (induct xs arbitrary: y) (simp_all add: zip_with_prev_def)
280 lemma stirling_row_code [code]:
281   "stirling_row 0 = "
282   "stirling_row (Suc n) = stirling_row_aux n 0 (stirling_row n)"
283 proof goal_cases
284   case 1
285   show ?case by (simp add: stirling_row_def)
286 next
287   case 2
288   have "stirling_row (Suc n) =
289     0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ "
290   proof (rule list_ext, goal_cases length nth)
291     case (nth i)
292     from nth have "i \<le> Suc n"
293       by simp
294     then consider "i = 0 \<or> i = Suc n" | "i > 0" "i \<le> n"
295       by linarith
296     then show ?case
297     proof cases
298       case 1
299       then show ?thesis
300         by (auto simp: nth_stirling_row nth_append)
301     next
302       case 2
303       then show ?thesis
304         by (cases i) (simp_all add: nth_append nth_stirling_row)
305     qed
306   next
307     case length
308     then show ?case by simp
309   qed
310   also have "0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @  =
311       zip_with_prev (\<lambda>a b. a + n * b) 0 (stirling_row n) @ "
312     by (cases n) (auto simp add: zip_with_prev_altdef stirling_row_def hd_map simp del: upt_Suc)
313   also have "\<dots> = stirling_row_aux n 0 (stirling_row n)"
314     by (simp add: stirling_row_aux_correct)
315   finally show ?case .
316 qed
318 lemma stirling_code [code]:
319   "stirling n k =
320     (if k = 0 then (if n = 0 then 1 else 0)
321      else if k > n then 0
322      else if k = n then 1
323      else stirling_row n ! k)"
324   by (simp add: nth_stirling_row)
326 end