src/HOL/Library/Stirling.thy
 author hoelzl Wed May 04 10:19:01 2016 +0200 (2016-05-04) changeset 63071 3ca3bc795908 child 63145 703edebd1d92 permissions -rw-r--r--
move Stirling numbers from AFP/Discrete_Summation
```     1 (* Authors: Amine Chaieb & Florian Haftmann, TU Muenchen
```
```     2             with contributions by Lukas Bulwahn and Manuel Eberl*)
```
```     3
```
```     4 section {* Stirling numbers of first and second kind *}
```
```     5
```
```     6 theory Stirling
```
```     7 imports Binomial
```
```     8 begin
```
```     9
```
```    10 subsection {* Stirling numbers of the second kind *}
```
```    11
```
```    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"
```
```    18
```
```    19 lemma Stirling_1 [simp]:
```
```    20   "Stirling (Suc n) (Suc 0) = 1"
```
```    21   by (induct n) simp_all
```
```    22
```
```    23 lemma Stirling_less [simp]:
```
```    24   "n < k \<Longrightarrow> Stirling n k = 0"
```
```    25   by (induct n k rule: Stirling.induct) simp_all
```
```    26
```
```    27 lemma Stirling_same [simp]:
```
```    28   "Stirling n n = 1"
```
```    29   by (induct n) simp_all
```
```    30
```
```    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
```
```    52
```
```    53 lemma Stirling_2:
```
```    54   "Stirling (Suc n) (Suc (Suc 0)) = 2 ^ n - 1"
```
```    55   using Stirling_2_2 by (cases n) simp_all
```
```    56
```
```    57 subsection {* Stirling numbers of the first kind *}
```
```    58
```
```    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"
```
```    65
```
```    66 lemma stirling_0 [simp]: "n > 0 \<Longrightarrow> stirling n 0 = 0"
```
```    67   by (cases n) simp_all
```
```    68
```
```    69 lemma stirling_less [simp]:
```
```    70   "n < k \<Longrightarrow> stirling n k = 0"
```
```    71   by (induct n k rule: stirling.induct) simp_all
```
```    72
```
```    73 lemma stirling_same [simp]:
```
```    74   "stirling n n = 1"
```
```    75   by (induct n) simp_all
```
```    76
```
```    77 lemma stirling_Suc_n_1:
```
```    78   "stirling (Suc n) (Suc 0) = fact n"
```
```    79   by (induct n) auto
```
```    80
```
```    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))
```
```    84
```
```    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
```
```   112
```
```   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
```
```   140
```
```   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
```
```   167
```
```   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
```
```   183
```
```   184
```
```   185 text \<open>A row of the Stirling number triangle\<close>
```
```   186
```
```   187 definition stirling_row :: "nat \<Rightarrow> nat list" where
```
```   188   "stirling_row n = [stirling n k. k \<leftarrow> [0..<Suc n]]"
```
```   189
```
```   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)
```
```   192
```
```   193 lemma length_stirling_row [simp]: "length (stirling_row n) = Suc n"
```
```   194   by (simp add: stirling_row_def)
```
```   195
```
```   196 lemma stirling_row_nonempty [simp]: "stirling_row n \<noteq> []"
```
```   197   using length_stirling_row[of n] by (auto simp del: length_stirling_row)
```
```   198
```
```   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
```
```   211
```
```   212
```
```   213 subsubsection \<open>Efficient code\<close>
```
```   214
```
```   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.
```
```   220
```
```   221   As a bonus, this is very efficient for applications where an entire row of
```
```   222   Stirling numbers is needed..
```
```   223 \<close>
```
```   224
```
```   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)"
```
```   227
```
```   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)
```
```   241
```
```   242
```
```   243 primrec stirling_row_aux where
```
```   244   "stirling_row_aux n y [] = [1]"
```
```   245 | "stirling_row_aux n y (x#xs) = (y + n * x) # stirling_row_aux n x xs"
```
```   246
```
```   247 lemma stirling_row_aux_correct:
```
```   248   "stirling_row_aux n y xs = zip_with_prev (\<lambda>a b. a + n * b) y xs @ [1]"
```
```   249   by (induction xs arbitrary: y) (simp_all add: zip_with_prev_def)
```
```   250
```
```   251 lemma stirling_row_code [code]:
```
```   252   "stirling_row 0 = [1]"
```
```   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]] @ [1]"
```
```   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]] @ [1] =
```
```   268                zip_with_prev (\<lambda>a b. a + n * b) 0 (stirling_row n) @ [1]"
```
```   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)
```
```   274
```
```   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)
```
```   280
```
```   281 end
```