src/HOL/Library/Stirling.thy
author wenzelm
Sat Nov 04 15:24:40 2017 +0100 (21 months ago)
changeset 67003 49850a679c2c
parent 66656 8f4d252ce2fe
child 68406 6beb45f6cf67
permissions -rw-r--r--
more robust sorted_entries;
     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 Suc [symmetric])
   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 (* TODO Move *)
   218 lemma list_ext:
   219   assumes "length xs = length ys"
   220   assumes "\<And>i. i < length xs \<Longrightarrow> xs ! i = ys ! i"
   221   shows "xs = ys"
   222   using assms
   223 proof (induction rule: list_induct2)
   224   case Nil
   225   then show ?case by simp
   226 next
   227   case (Cons x xs y ys)
   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
   232   ultimately show ?case by simp
   233 qed
   234 
   235 
   236 subsubsection \<open>Efficient code\<close>
   237 
   238 text \<open>
   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.
   243 
   244   As a bonus, this is very efficient for applications where an entire row of
   245   Stirling numbers is needed.
   246 \<close>
   247 
   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 = map2 f (x # xs) xs"
   250 
   251 lemma zip_with_prev_altdef:
   252   "zip_with_prev f x xs =
   253     (if xs = [] then [] else f x (hd xs) # [f (xs!i) (xs!(i+1)). i \<leftarrow> [0..<length xs - 1]])"
   254 proof (cases xs)
   255   case Nil
   256   then show ?thesis
   257     by (simp add: zip_with_prev_def)
   258 next
   259   case (Cons y ys)
   260   then have "zip_with_prev f x xs = f x (hd xs) # zip_with_prev f y ys"
   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
   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
   269 
   270 
   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"
   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]"
   278   by (induct xs arbitrary: y) (simp_all add: zip_with_prev_def)
   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)"
   283 proof goal_cases
   284   case 1
   285   show ?case by (simp add: stirling_row_def)
   286 next
   287   case 2
   288   have "stirling_row (Suc n) =
   289     0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ [1]"
   290   proof (rule list_ext, goal_cases length nth)
   291     case (nth i)
   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
   297     proof cases
   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
   310   also have "0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ [1] =
   311       zip_with_prev (\<lambda>a b. a + n * b) 0 (stirling_row n) @ [1]"
   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)
   315   finally show ?case .
   316 qed
   317 
   318 lemma stirling_code [code]:
   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)"
   324   by (simp add: nth_stirling_row)
   325 
   326 end