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