src/HOL/Library/Stirling.thy
 author wenzelm Wed May 25 11:49:40 2016 +0200 (2016-05-25) changeset 63145 703edebd1d92 parent 63071 3ca3bc795908 child 63487 6e29fb72e659 permissions -rw-r--r--
isabelle update_cartouches -c -t;
1 (* Authors: Amine Chaieb & Florian Haftmann, TU Muenchen
2             with contributions by Lukas Bulwahn and Manuel Eberl*)
4 section \<open>Stirling numbers of first and second kind\<close>
6 theory Stirling
7 imports Binomial
8 begin
10 subsection \<open>Stirling numbers of the second kind\<close>
12 fun Stirling :: "nat \<Rightarrow> nat \<Rightarrow> nat"
13 where
14   "Stirling 0 0 = 1"
15 | "Stirling 0 (Suc k) = 0"
16 | "Stirling (Suc n) 0 = 0"
17 | "Stirling (Suc n) (Suc k) = Suc k * Stirling n (Suc k) + Stirling n k"
19 lemma Stirling_1 [simp]:
20   "Stirling (Suc n) (Suc 0) = 1"
21   by (induct n) simp_all
23 lemma Stirling_less [simp]:
24   "n < k \<Longrightarrow> Stirling n k = 0"
25   by (induct n k rule: Stirling.induct) simp_all
27 lemma Stirling_same [simp]:
28   "Stirling n n = 1"
29   by (induct n) simp_all
31 lemma Stirling_2_2:
32   "Stirling (Suc (Suc n)) (Suc (Suc 0)) = 2 ^ Suc n - 1"
33 proof (induct n)
34   case 0 then show ?case by simp
35 next
36   case (Suc n)
37   have "Stirling (Suc (Suc (Suc n))) (Suc (Suc 0)) =
38     2 * Stirling (Suc (Suc n)) (Suc (Suc 0)) + Stirling (Suc (Suc n)) (Suc 0)" by simp
39   also have "\<dots> = 2 * (2 ^ Suc n - 1) + 1"
40     by (simp only: Suc Stirling_1)
41   also have "\<dots> = 2 ^ Suc (Suc n) - 1"
42   proof -
43     have "(2::nat) ^ Suc n - 1 > 0" by (induct n) simp_all
44     then have "2 * ((2::nat) ^ Suc n - 1) > 0" by simp
45     then have "2 \<le> 2 * ((2::nat) ^ Suc n)" by simp
46     with add_diff_assoc2 [of 2 "2 * 2 ^ Suc n" 1]
47       have "2 * 2 ^ Suc n - 2 + (1::nat) = 2 * 2 ^ Suc n + 1 - 2" .
48     then show ?thesis by (simp add: nat_distrib)
49   qed
50   finally show ?case by simp
51 qed
53 lemma Stirling_2:
54   "Stirling (Suc n) (Suc (Suc 0)) = 2 ^ n - 1"
55   using Stirling_2_2 by (cases n) simp_all
57 subsection \<open>Stirling numbers of the first kind\<close>
59 fun stirling :: "nat \<Rightarrow> nat \<Rightarrow> nat"
60 where
61   "stirling 0 0 = 1"
62 | "stirling 0 (Suc k) = 0"
63 | "stirling (Suc n) 0 = 0"
64 | "stirling (Suc n) (Suc k) = n * stirling n (Suc k) + stirling n k"
66 lemma stirling_0 [simp]: "n > 0 \<Longrightarrow> stirling n 0 = 0"
67   by (cases n) simp_all
69 lemma stirling_less [simp]:
70   "n < k \<Longrightarrow> stirling n k = 0"
71   by (induct n k rule: stirling.induct) simp_all
73 lemma stirling_same [simp]:
74   "stirling n n = 1"
75   by (induct n) simp_all
77 lemma stirling_Suc_n_1:
78   "stirling (Suc n) (Suc 0) = fact n"
79   by (induct n) auto
81 lemma stirling_Suc_n_n:
82   shows "stirling (Suc n) n = Suc n choose 2"
83 by (induct n) (auto simp add: numerals(2))
85 lemma stirling_Suc_n_2:
86   assumes "n \<ge> Suc 0"
87   shows "stirling (Suc n) 2 = (\<Sum>k=1..n. fact n div k)"
88 using assms
89 proof (induct n)
90   case 0 from this show ?case by simp
91 next
92   case (Suc n)
93   show ?case
94   proof (cases n)
95     case 0 from this show ?thesis by (simp add: numerals(2))
96   next
97     case Suc
98     from this have geq1: "Suc 0 \<le> n" by simp
99     have "stirling (Suc (Suc n)) 2 = Suc n * stirling (Suc n) 2 + stirling (Suc n) (Suc 0)"
100       by (simp only: stirling.simps(4)[of "Suc n"] numerals(2))
101     also have "... = Suc n * (\<Sum>k=1..n. fact n div k) + fact n"
102       using Suc.hyps[OF geq1]
103       by (simp only: stirling_Suc_n_1 of_nat_fact of_nat_add of_nat_mult)
104     also have "... = Suc n * (\<Sum>k=1..n. fact n div k) + Suc n * fact n div Suc n"
105       by (metis nat.distinct(1) nonzero_mult_divide_cancel_left)
106     also have "... = (\<Sum>k=1..n. fact (Suc n) div k) + fact (Suc n) div Suc n"
107       by (simp add: setsum_right_distrib div_mult_swap dvd_fact)
108     also have "... = (\<Sum>k=1..Suc n. fact (Suc n) div k)" by simp
109     finally show ?thesis .
110   qed
111 qed
113 lemma of_nat_stirling_Suc_n_2:
114   assumes "n \<ge> Suc 0"
115   shows "(of_nat (stirling (Suc n) 2)::'a::field_char_0) = fact n * (\<Sum>k=1..n. (1 / of_nat k))"
116 using assms
117 proof (induct n)
118   case 0 from this show ?case by simp
119 next
120   case (Suc n)
121   show ?case
122   proof (cases n)
123     case 0 from this show ?thesis by (auto simp add: numerals(2))
124   next
125     case Suc
126     from this have geq1: "Suc 0 \<le> n" by simp
127     have "(of_nat (stirling (Suc (Suc n)) 2)::'a) =
128       of_nat (Suc n * stirling (Suc n) 2 + stirling (Suc n) (Suc 0))"
129       by (simp only: stirling.simps(4)[of "Suc n"] numerals(2))
130     also have "... = of_nat (Suc n) * (fact n * (\<Sum>k = 1..n. 1 / of_nat k)) + fact n"
131       using Suc.hyps[OF geq1]
132       by (simp only: stirling_Suc_n_1 of_nat_fact of_nat_add of_nat_mult)
133     also have "... = fact (Suc n) * (\<Sum>k = 1..n. 1 / of_nat k) + fact (Suc n) * (1 / of_nat (Suc n))"
134       using of_nat_neq_0 by auto
135     also have "... = fact (Suc n) * (\<Sum>k = 1..Suc n. 1 / of_nat k)"
136       by (simp add: distrib_left)
137     finally show ?thesis .
138   qed
139 qed
141 lemma setsum_stirling:
142   "(\<Sum>k\<le>n. stirling n k) = fact n"
143 proof (induct n)
144   case 0
145   from this show ?case by simp
146 next
147   case (Suc n)
148   have "(\<Sum>k\<le>Suc n. stirling (Suc n) k) = stirling (Suc n) 0 + (\<Sum>k\<le>n. stirling (Suc n) (Suc k))"
149     by (simp only: setsum_atMost_Suc_shift)
150   also have "\<dots> = (\<Sum>k\<le>n. stirling (Suc n) (Suc k))" by simp
151   also have "\<dots> = (\<Sum>k\<le>n. n * stirling n (Suc k) + stirling n k)" by simp
152   also have "\<dots> = n * (\<Sum>k\<le>n. stirling n (Suc k)) + (\<Sum>k\<le>n. stirling n k)"
153     by (simp add: setsum.distrib setsum_right_distrib)
154   also have "\<dots> = n * fact n + fact n"
155   proof -
156     have "n * (\<Sum>k\<le>n. stirling n (Suc k)) = n * ((\<Sum>k\<le>Suc n. stirling n k) - stirling n 0)"
157       by (metis add_diff_cancel_left' setsum_atMost_Suc_shift)
158     also have "\<dots> = n * (\<Sum>k\<le>n. stirling n k)" by (cases n) simp+
159     also have "\<dots> = n * fact n" using Suc.hyps by simp
160     finally have "n * (\<Sum>k\<le>n. stirling n (Suc k)) = n * fact n" .
161     moreover have "(\<Sum>k\<le>n. stirling n k) = fact n" using Suc.hyps .
162     ultimately show ?thesis by simp
163   qed
164   also have "\<dots> = fact (Suc n)" by simp
165   finally show ?case .
166 qed
168 lemma stirling_pochhammer:
169   "(\<Sum>k\<le>n. of_nat (stirling n k) * x ^ k) = (pochhammer x n :: 'a :: comm_semiring_1)"
170 proof (induction n)
171   case (Suc n)
172   have "of_nat (n * stirling n 0) = (0 :: 'a)" by (cases n) simp_all
173   hence "(\<Sum>k\<le>Suc n. of_nat (stirling (Suc n) k) * x ^ k) =
174             (of_nat (n * stirling n 0) * x ^ 0 +
175             (\<Sum>i\<le>n. of_nat (n * stirling n (Suc i)) * (x ^ Suc i))) +
176             (\<Sum>i\<le>n. of_nat (stirling n i) * (x ^ Suc i))"
177     by (subst setsum_atMost_Suc_shift) (simp add: setsum.distrib ring_distribs)
178   also have "\<dots> = pochhammer x (Suc n)"
179     by (subst setsum_atMost_Suc_shift [symmetric])
180        (simp add: algebra_simps setsum.distrib setsum_right_distrib pochhammer_Suc Suc [symmetric])
181   finally show ?case .
182 qed simp_all
185 text \<open>A row of the Stirling number triangle\<close>
187 definition stirling_row :: "nat \<Rightarrow> nat list" where
188   "stirling_row n = [stirling n k. k \<leftarrow> [0..<Suc n]]"
190 lemma nth_stirling_row: "k \<le> n \<Longrightarrow> stirling_row n ! k = stirling n k"
191   by (simp add: stirling_row_def del: upt_Suc)
193 lemma length_stirling_row [simp]: "length (stirling_row n) = Suc n"
194   by (simp add: stirling_row_def)
196 lemma stirling_row_nonempty [simp]: "stirling_row n \<noteq> []"
197   using length_stirling_row[of n] by (auto simp del: length_stirling_row)
199 (* TODO Move *)
200 lemma list_ext:
201   assumes "length xs = length ys"
202   assumes "\<And>i. i < length xs \<Longrightarrow> xs ! i = ys ! i"
203   shows   "xs = ys"
204 using assms
205 proof (induction rule: list_induct2)
206   case (Cons x xs y ys)
207   from Cons.prems[of 0] have "x = y" by simp
208   moreover from Cons.prems[of "Suc i" for i] have "xs = ys" by (intro Cons.IH) simp
209   ultimately show ?case by simp
210 qed simp_all
213 subsubsection \<open>Efficient code\<close>
215 text \<open>
216   Naively using the defining equations of the Stirling numbers of the first kind to
217   compute them leads to exponential run time due to repeated computations.
218   We can use memoisation to compute them row by row without repeating computations, at
219   the cost of computing a few unneeded values.
221   As a bonus, this is very efficient for applications where an entire row of
222   Stirling numbers is needed..
223 \<close>
225 definition zip_with_prev :: "('a \<Rightarrow> 'a \<Rightarrow> 'b) \<Rightarrow> 'a \<Rightarrow> 'a list \<Rightarrow> 'b list" where
226   "zip_with_prev f x xs = map (\<lambda>(x,y). f x y) (zip (x # xs) xs)"
228 lemma zip_with_prev_altdef:
229   "zip_with_prev f x xs =
230      (if xs = [] then [] else f x (hd xs) # [f (xs!i) (xs!(i+1)). i \<leftarrow> [0..<length xs - 1]])"
231 proof (cases xs)
232   case (Cons y ys)
233   hence "zip_with_prev f x xs = f x (hd xs) # zip_with_prev f y ys"
234     by (simp add: zip_with_prev_def)
235   also have "zip_with_prev f y ys = map (\<lambda>i. f (xs ! i) (xs ! (i + 1))) [0..<length xs - 1]"
236     unfolding Cons
237     by (induction ys arbitrary: y)
238        (simp_all add: zip_with_prev_def upt_conv_Cons map_Suc_upt [symmetric] del: upt_Suc)
239   finally show ?thesis using Cons by simp
240 qed (simp add: zip_with_prev_def)
243 primrec stirling_row_aux where
244   "stirling_row_aux n y [] = "
245 | "stirling_row_aux n y (x#xs) = (y + n * x) # stirling_row_aux n x xs"
247 lemma stirling_row_aux_correct:
248   "stirling_row_aux n y xs = zip_with_prev (\<lambda>a b. a + n * b) y xs @ "
249   by (induction xs arbitrary: y) (simp_all add: zip_with_prev_def)
251 lemma stirling_row_code [code]:
252   "stirling_row 0 = "
253   "stirling_row (Suc n) = stirling_row_aux n 0 (stirling_row n)"
254 proof -
255   have "stirling_row (Suc n) =
256           0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ "
257   proof (rule list_ext, goal_cases length nth)
258     case (nth i)
259     from nth have "i \<le> Suc n" by simp
260     then consider "i = 0" | j where "i > 0" "i \<le> n" | "i = Suc n" by linarith
261     thus ?case
262     proof cases
263       assume i: "i > 0" "i \<le> n"
264       from this show ?thesis by (cases i) (simp_all add: nth_append nth_stirling_row)
265     qed (simp_all add: nth_stirling_row nth_append)
266   qed simp
267   also have "0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @  =
268                zip_with_prev (\<lambda>a b. a + n * b) 0 (stirling_row n) @ "
269     by (cases n) (auto simp add: zip_with_prev_altdef stirling_row_def hd_map simp del: upt_Suc)
270   also have "\<dots> = stirling_row_aux n 0 (stirling_row n)"
271     by (simp add: stirling_row_aux_correct)
272   finally show "stirling_row (Suc n) = stirling_row_aux n 0 (stirling_row n)" .
273 qed (simp add: stirling_row_def)
275 lemma stirling_code [code]:
276   "stirling n k = (if k = 0 then if n = 0 then 1 else 0
277                    else if k > n then 0 else if k = n then 1
278                    else stirling_row n ! k)"
279   by (simp add: nth_stirling_row)
281 end