src/HOL/Library/Stirling.thy
 author haftmann Fri Mar 22 19:18:08 2019 +0000 (4 months ago) changeset 69946 494934c30f38 parent 68975 5ce4d117cea7 child 70113 c8deb8ba6d05 permissions -rw-r--r--
improved code equations taken over from AFP
```     1 (*  Title:      HOL/Library/Stirling.thy
```
```     2     Author:     Amine Chaieb
```
```     3     Author:     Florian Haftmann
```
```     4     Author:     Lukas Bulwahn
```
```     5     Author:     Manuel Eberl
```
```     6 *)
```
```     7
```
```     8 section \<open>Stirling numbers of first and second kind\<close>
```
```     9
```
```    10 theory Stirling
```
```    11 imports Main
```
```    12 begin
```
```    13
```
```    14 subsection \<open>Stirling numbers of the second kind\<close>
```
```    15
```
```    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"
```
```    22
```
```    23 lemma Stirling_1 [simp]: "Stirling (Suc n) (Suc 0) = 1"
```
```    24   by (induct n) simp_all
```
```    25
```
```    26 lemma Stirling_less [simp]: "n < k \<Longrightarrow> Stirling n k = 0"
```
```    27   by (induct n k rule: Stirling.induct) simp_all
```
```    28
```
```    29 lemma Stirling_same [simp]: "Stirling n n = 1"
```
```    30   by (induct n) simp_all
```
```    31
```
```    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
```
```    58
```
```    59 lemma Stirling_2: "Stirling (Suc n) (Suc (Suc 0)) = 2 ^ n - 1"
```
```    60   using Stirling_2_2 by (cases n) simp_all
```
```    61
```
```    62
```
```    63 subsection \<open>Stirling numbers of the first kind\<close>
```
```    64
```
```    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"
```
```    71
```
```    72 lemma stirling_0 [simp]: "n > 0 \<Longrightarrow> stirling n 0 = 0"
```
```    73   by (cases n) simp_all
```
```    74
```
```    75 lemma stirling_less [simp]: "n < k \<Longrightarrow> stirling n k = 0"
```
```    76   by (induct n k rule: stirling.induct) simp_all
```
```    77
```
```    78 lemma stirling_same [simp]: "stirling n n = 1"
```
```    79   by (induct n) simp_all
```
```    80
```
```    81 lemma stirling_Suc_n_1: "stirling (Suc n) (Suc 0) = fact n"
```
```    82   by (induct n) auto
```
```    83
```
```    84 lemma stirling_Suc_n_n: "stirling (Suc n) n = Suc n choose 2"
```
```    85   by (induct n) (auto simp add: numerals(2))
```
```    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)"
```
```    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
```
```   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))"
```
```   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
```
```   151
```
```   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
```
```   182
```
```   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 flip: Suc)
```
```   199   finally show ?case .
```
```   200 qed
```
```   201
```
```   202
```
```   203 text \<open>A row of the Stirling number triangle\<close>
```
```   204
```
```   205 definition stirling_row :: "nat \<Rightarrow> nat list"
```
```   206   where "stirling_row n = [stirling n k. k \<leftarrow> [0..<Suc n]]"
```
```   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
```
```   218 subsubsection \<open>Efficient code\<close>
```
```   219
```
```   220 text \<open>
```
```   221   Naively using the defining equations of the Stirling numbers of the first
```
```   222   kind to compute them leads to exponential run time due to repeated
```
```   223   computations. We can use memoisation to compute them row by row without
```
```   224   repeating computations, at the cost of computing a few unneeded values.
```
```   225
```
```   226   As a bonus, this is very efficient for applications where an entire row of
```
```   227   Stirling numbers is needed.
```
```   228 \<close>
```
```   229
```
```   230 definition zip_with_prev :: "('a \<Rightarrow> 'a \<Rightarrow> 'b) \<Rightarrow> 'a \<Rightarrow> 'a list \<Rightarrow> 'b list"
```
```   231   where "zip_with_prev f x xs = map2 f (x # xs) xs"
```
```   232
```
```   233 lemma zip_with_prev_altdef:
```
```   234   "zip_with_prev f x xs =
```
```   235     (if xs = [] then [] else f x (hd xs) # [f (xs!i) (xs!(i+1)). i \<leftarrow> [0..<length xs - 1]])"
```
```   236 proof (cases xs)
```
```   237   case Nil
```
```   238   then show ?thesis
```
```   239     by (simp add: zip_with_prev_def)
```
```   240 next
```
```   241   case (Cons y ys)
```
```   242   then have "zip_with_prev f x xs = f x (hd xs) # zip_with_prev f y ys"
```
```   243     by (simp add: zip_with_prev_def)
```
```   244   also have "zip_with_prev f y ys = map (\<lambda>i. f (xs ! i) (xs ! (i + 1))) [0..<length xs - 1]"
```
```   245     unfolding Cons
```
```   246     by (induct ys arbitrary: y)
```
```   247       (simp_all add: zip_with_prev_def upt_conv_Cons flip: map_Suc_upt del: upt_Suc)
```
```   248   finally show ?thesis
```
```   249     using Cons by simp
```
```   250 qed
```
```   251
```
```   252
```
```   253 primrec stirling_row_aux
```
```   254   where
```
```   255     "stirling_row_aux n y [] = "
```
```   256   | "stirling_row_aux n y (x#xs) = (y + n * x) # stirling_row_aux n x xs"
```
```   257
```
```   258 lemma stirling_row_aux_correct:
```
```   259   "stirling_row_aux n y xs = zip_with_prev (\<lambda>a b. a + n * b) y xs @ "
```
```   260   by (induct xs arbitrary: y) (simp_all add: zip_with_prev_def)
```
```   261
```
```   262 lemma stirling_row_code [code]:
```
```   263   "stirling_row 0 = "
```
```   264   "stirling_row (Suc n) = stirling_row_aux n 0 (stirling_row n)"
```
```   265 proof goal_cases
```
```   266   case 1
```
```   267   show ?case by (simp add: stirling_row_def)
```
```   268 next
```
```   269   case 2
```
```   270   have "stirling_row (Suc n) =
```
```   271     0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ "
```
```   272   proof (rule nth_equalityI, goal_cases length nth)
```
```   273     case (nth i)
```
```   274     from nth have "i \<le> Suc n"
```
```   275       by simp
```
```   276     then consider "i = 0 \<or> i = Suc n" | "i > 0" "i \<le> n"
```
```   277       by linarith
```
```   278     then show ?case
```
```   279     proof cases
```
```   280       case 1
```
```   281       then show ?thesis
```
```   282         by (auto simp: nth_stirling_row nth_append)
```
```   283     next
```
```   284       case 2
```
```   285       then show ?thesis
```
```   286         by (cases i) (simp_all add: nth_append nth_stirling_row)
```
```   287     qed
```
```   288   next
```
```   289     case length
```
```   290     then show ?case by simp
```
```   291   qed
```
```   292   also have "0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @  =
```
```   293       zip_with_prev (\<lambda>a b. a + n * b) 0 (stirling_row n) @ "
```
```   294     by (cases n) (auto simp add: zip_with_prev_altdef stirling_row_def hd_map simp del: upt_Suc)
```
```   295   also have "\<dots> = stirling_row_aux n 0 (stirling_row n)"
```
```   296     by (simp add: stirling_row_aux_correct)
```
```   297   finally show ?case .
```
```   298 qed
```
```   299
```
```   300 lemma stirling_code [code]:
```
```   301   "stirling n k =
```
```   302     (if k = 0 then (if n = 0 then 1 else 0)
```
```   303      else if k > n then 0
```
```   304      else if k = n then 1
```
```   305      else stirling_row n ! k)"
```
```   306   by (simp add: nth_stirling_row)
```
```   307
```
```   308 end
```