63487
|
1 |
(* Title: HOL/Library/Stirling.thy
|
|
2 |
Author: Amine Chaieb
|
|
3 |
Author: Florian Haftmann
|
|
4 |
Author: Lukas Bulwahn
|
|
5 |
Author: Manuel Eberl
|
|
6 |
*)
|
63071
|
7 |
|
63145
|
8 |
section \<open>Stirling numbers of first and second kind\<close>
|
63071
|
9 |
|
|
10 |
theory Stirling
|
|
11 |
imports Binomial
|
|
12 |
begin
|
|
13 |
|
63145
|
14 |
subsection \<open>Stirling numbers of the second kind\<close>
|
63071
|
15 |
|
|
16 |
fun Stirling :: "nat \<Rightarrow> nat \<Rightarrow> nat"
|
63487
|
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"
|
63071
|
22 |
|
63487
|
23 |
lemma Stirling_1 [simp]: "Stirling (Suc n) (Suc 0) = 1"
|
63071
|
24 |
by (induct n) simp_all
|
|
25 |
|
63487
|
26 |
lemma Stirling_less [simp]: "n < k \<Longrightarrow> Stirling n k = 0"
|
63071
|
27 |
by (induct n k rule: Stirling.induct) simp_all
|
|
28 |
|
63487
|
29 |
lemma Stirling_same [simp]: "Stirling n n = 1"
|
63071
|
30 |
by (induct n) simp_all
|
|
31 |
|
63487
|
32 |
lemma Stirling_2_2: "Stirling (Suc (Suc n)) (Suc (Suc 0)) = 2 ^ Suc n - 1"
|
63071
|
33 |
proof (induct n)
|
63487
|
34 |
case 0
|
|
35 |
then show ?case by simp
|
63071
|
36 |
next
|
|
37 |
case (Suc n)
|
|
38 |
have "Stirling (Suc (Suc (Suc n))) (Suc (Suc 0)) =
|
63487
|
39 |
2 * Stirling (Suc (Suc n)) (Suc (Suc 0)) + Stirling (Suc (Suc n)) (Suc 0)"
|
|
40 |
by simp
|
63071
|
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 -
|
63487
|
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
|
63071
|
51 |
with add_diff_assoc2 [of 2 "2 * 2 ^ Suc n" 1]
|
63487
|
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)
|
63071
|
55 |
qed
|
|
56 |
finally show ?case by simp
|
|
57 |
qed
|
|
58 |
|
63487
|
59 |
lemma Stirling_2: "Stirling (Suc n) (Suc (Suc 0)) = 2 ^ n - 1"
|
63071
|
60 |
using Stirling_2_2 by (cases n) simp_all
|
|
61 |
|
63487
|
62 |
|
63145
|
63 |
subsection \<open>Stirling numbers of the first kind\<close>
|
63071
|
64 |
|
|
65 |
fun stirling :: "nat \<Rightarrow> nat \<Rightarrow> nat"
|
63487
|
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"
|
63071
|
71 |
|
|
72 |
lemma stirling_0 [simp]: "n > 0 \<Longrightarrow> stirling n 0 = 0"
|
|
73 |
by (cases n) simp_all
|
|
74 |
|
63487
|
75 |
lemma stirling_less [simp]: "n < k \<Longrightarrow> stirling n k = 0"
|
63071
|
76 |
by (induct n k rule: stirling.induct) simp_all
|
|
77 |
|
63487
|
78 |
lemma stirling_same [simp]: "stirling n n = 1"
|
63071
|
79 |
by (induct n) simp_all
|
|
80 |
|
63487
|
81 |
lemma stirling_Suc_n_1: "stirling (Suc n) (Suc 0) = fact n"
|
63071
|
82 |
by (induct n) auto
|
|
83 |
|
63487
|
84 |
lemma stirling_Suc_n_n: "stirling (Suc n) n = Suc n choose 2"
|
|
85 |
by (induct n) (auto simp add: numerals(2))
|
63071
|
86 |
|
|
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)"
|
63487
|
90 |
using assms
|
63071
|
91 |
proof (induct n)
|
63487
|
92 |
case 0
|
|
93 |
then show ?case by simp
|
63071
|
94 |
next
|
|
95 |
case (Suc n)
|
|
96 |
show ?case
|
|
97 |
proof (cases n)
|
63487
|
98 |
case 0
|
|
99 |
then show ?thesis
|
|
100 |
by (simp add: numerals(2))
|
63071
|
101 |
next
|
|
102 |
case Suc
|
63487
|
103 |
then have geq1: "Suc 0 \<le> n"
|
|
104 |
by simp
|
63071
|
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))
|
63487
|
107 |
also have "\<dots> = Suc n * (\<Sum>k=1..n. fact n div k) + fact n"
|
63071
|
108 |
using Suc.hyps[OF geq1]
|
|
109 |
by (simp only: stirling_Suc_n_1 of_nat_fact of_nat_add of_nat_mult)
|
63487
|
110 |
also have "\<dots> = Suc n * (\<Sum>k=1..n. fact n div k) + Suc n * fact n div Suc n"
|
64240
|
111 |
by (metis nat.distinct(1) nonzero_mult_div_cancel_left)
|
63487
|
112 |
also have "\<dots> = (\<Sum>k=1..n. fact (Suc n) div k) + fact (Suc n) div Suc n"
|
64267
|
113 |
by (simp add: sum_distrib_left div_mult_swap dvd_fact)
|
63487
|
114 |
also have "\<dots> = (\<Sum>k=1..Suc n. fact (Suc n) div k)"
|
|
115 |
by simp
|
63071
|
116 |
finally show ?thesis .
|
|
117 |
qed
|
|
118 |
qed
|
|
119 |
|
|
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))"
|
63487
|
123 |
using assms
|
63071
|
124 |
proof (induct n)
|
63487
|
125 |
case 0
|
|
126 |
then show ?case by simp
|
63071
|
127 |
next
|
|
128 |
case (Suc n)
|
|
129 |
show ?case
|
|
130 |
proof (cases n)
|
63487
|
131 |
case 0
|
|
132 |
then show ?thesis
|
|
133 |
by (auto simp add: numerals(2))
|
63071
|
134 |
next
|
|
135 |
case Suc
|
63487
|
136 |
then have geq1: "Suc 0 \<le> n"
|
|
137 |
by simp
|
63071
|
138 |
have "(of_nat (stirling (Suc (Suc n)) 2)::'a) =
|
63487
|
139 |
of_nat (Suc n * stirling (Suc n) 2 + stirling (Suc n) (Suc 0))"
|
63071
|
140 |
by (simp only: stirling.simps(4)[of "Suc n"] numerals(2))
|
63487
|
141 |
also have "\<dots> = of_nat (Suc n) * (fact n * (\<Sum>k = 1..n. 1 / of_nat k)) + fact n"
|
63071
|
142 |
using Suc.hyps[OF geq1]
|
|
143 |
by (simp only: stirling_Suc_n_1 of_nat_fact of_nat_add of_nat_mult)
|
63487
|
144 |
also have "\<dots> = fact (Suc n) * (\<Sum>k = 1..n. 1 / of_nat k) + fact (Suc n) * (1 / of_nat (Suc n))"
|
63071
|
145 |
using of_nat_neq_0 by auto
|
63487
|
146 |
also have "\<dots> = fact (Suc n) * (\<Sum>k = 1..Suc n. 1 / of_nat k)"
|
63071
|
147 |
by (simp add: distrib_left)
|
|
148 |
finally show ?thesis .
|
|
149 |
qed
|
|
150 |
qed
|
|
151 |
|
64267
|
152 |
lemma sum_stirling: "(\<Sum>k\<le>n. stirling n k) = fact n"
|
63071
|
153 |
proof (induct n)
|
|
154 |
case 0
|
63487
|
155 |
then show ?case by simp
|
63071
|
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))"
|
64267
|
159 |
by (simp only: sum_atMost_Suc_shift)
|
63487
|
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
|
63071
|
164 |
also have "\<dots> = n * (\<Sum>k\<le>n. stirling n (Suc k)) + (\<Sum>k\<le>n. stirling n k)"
|
64267
|
165 |
by (simp add: sum.distrib sum_distrib_left)
|
63071
|
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)"
|
64267
|
169 |
by (metis add_diff_cancel_left' sum_atMost_Suc_shift)
|
63487
|
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
|
63071
|
174 |
finally have "n * (\<Sum>k\<le>n. stirling n (Suc k)) = n * fact n" .
|
63487
|
175 |
moreover have "(\<Sum>k\<le>n. stirling n k) = fact n"
|
|
176 |
using Suc.hyps .
|
63071
|
177 |
ultimately show ?thesis by simp
|
|
178 |
qed
|
|
179 |
also have "\<dots> = fact (Suc n)" by simp
|
|
180 |
finally show ?case .
|
|
181 |
qed
|
|
182 |
|
|
183 |
lemma stirling_pochhammer:
|
63487
|
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
|
63071
|
189 |
case (Suc n)
|
|
190 |
have "of_nat (n * stirling n 0) = (0 :: 'a)" by (cases n) simp_all
|
63487
|
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))"
|
64267
|
195 |
by (subst sum_atMost_Suc_shift) (simp add: sum.distrib ring_distribs)
|
63071
|
196 |
also have "\<dots> = pochhammer x (Suc n)"
|
64267
|
197 |
by (subst sum_atMost_Suc_shift [symmetric])
|
|
198 |
(simp add: algebra_simps sum.distrib sum_distrib_left pochhammer_Suc Suc [symmetric])
|
63071
|
199 |
finally show ?case .
|
63487
|
200 |
qed
|
63071
|
201 |
|
|
202 |
|
|
203 |
text \<open>A row of the Stirling number triangle\<close>
|
|
204 |
|
63487
|
205 |
definition stirling_row :: "nat \<Rightarrow> nat list"
|
|
206 |
where "stirling_row n = [stirling n k. k \<leftarrow> [0..<Suc n]]"
|
63071
|
207 |
|
|
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)
|
|
210 |
|
|
211 |
lemma length_stirling_row [simp]: "length (stirling_row n) = Suc n"
|
|
212 |
by (simp add: stirling_row_def)
|
|
213 |
|
|
214 |
lemma stirling_row_nonempty [simp]: "stirling_row n \<noteq> []"
|
|
215 |
using length_stirling_row[of n] by (auto simp del: length_stirling_row)
|
|
216 |
|
|
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"
|
63487
|
221 |
shows "xs = ys"
|
|
222 |
using assms
|
63071
|
223 |
proof (induction rule: list_induct2)
|
63487
|
224 |
case Nil
|
|
225 |
then show ?case by simp
|
|
226 |
next
|
63071
|
227 |
case (Cons x xs y ys)
|
63487
|
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
|
63071
|
232 |
ultimately show ?case by simp
|
63487
|
233 |
qed
|
63071
|
234 |
|
|
235 |
|
|
236 |
subsubsection \<open>Efficient code\<close>
|
|
237 |
|
|
238 |
text \<open>
|
63487
|
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.
|
63071
|
243 |
|
|
244 |
As a bonus, this is very efficient for applications where an entire row of
|
63487
|
245 |
Stirling numbers is needed.
|
63071
|
246 |
\<close>
|
|
247 |
|
63487
|
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 = map (\<lambda>(x,y). f x y) (zip (x # xs) xs)"
|
63071
|
250 |
|
|
251 |
lemma zip_with_prev_altdef:
|
|
252 |
"zip_with_prev f x xs =
|
63487
|
253 |
(if xs = [] then [] else f x (hd xs) # [f (xs!i) (xs!(i+1)). i \<leftarrow> [0..<length xs - 1]])"
|
63071
|
254 |
proof (cases xs)
|
63487
|
255 |
case Nil
|
|
256 |
then show ?thesis
|
|
257 |
by (simp add: zip_with_prev_def)
|
|
258 |
next
|
63071
|
259 |
case (Cons y ys)
|
63487
|
260 |
then have "zip_with_prev f x xs = f x (hd xs) # zip_with_prev f y ys"
|
63071
|
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
|
63487
|
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
|
63071
|
269 |
|
|
270 |
|
63487
|
271 |
primrec stirling_row_aux
|
|
272 |
where
|
|
273 |
"stirling_row_aux n y [] = [1]"
|
|
274 |
| "stirling_row_aux n y (x#xs) = (y + n * x) # stirling_row_aux n x xs"
|
63071
|
275 |
|
|
276 |
lemma stirling_row_aux_correct:
|
|
277 |
"stirling_row_aux n y xs = zip_with_prev (\<lambda>a b. a + n * b) y xs @ [1]"
|
63487
|
278 |
by (induct xs arbitrary: y) (simp_all add: zip_with_prev_def)
|
63071
|
279 |
|
|
280 |
lemma stirling_row_code [code]:
|
|
281 |
"stirling_row 0 = [1]"
|
|
282 |
"stirling_row (Suc n) = stirling_row_aux n 0 (stirling_row n)"
|
63487
|
283 |
proof goal_cases
|
|
284 |
case 1
|
|
285 |
show ?case by (simp add: stirling_row_def)
|
|
286 |
next
|
|
287 |
case 2
|
63071
|
288 |
have "stirling_row (Suc n) =
|
63487
|
289 |
0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ [1]"
|
63071
|
290 |
proof (rule list_ext, goal_cases length nth)
|
|
291 |
case (nth i)
|
63487
|
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
|
63071
|
297 |
proof cases
|
63487
|
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
|
63071
|
310 |
also have "0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ [1] =
|
63487
|
311 |
zip_with_prev (\<lambda>a b. a + n * b) 0 (stirling_row n) @ [1]"
|
63071
|
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)
|
63487
|
315 |
finally show ?case .
|
|
316 |
qed
|
63071
|
317 |
|
|
318 |
lemma stirling_code [code]:
|
63487
|
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)"
|
63071
|
324 |
by (simp add: nth_stirling_row)
|
|
325 |
|
|
326 |
end
|