src/HOL/Library/Stirling.thy
author paulson <lp15@cam.ac.uk>
Wed Apr 10 21:29:32 2019 +0100 (4 months ago)
changeset 70113 c8deb8ba6d05
parent 68975 5ce4d117cea7
permissions -rw-r--r--
Fixing the main Homology theory; also moving a lot of sum/prod lemmas into their generic context
     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 [] = [1]"
   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 @ [1]"
   260   by (induct xs arbitrary: y) (simp_all add: zip_with_prev_def)
   261 
   262 lemma stirling_row_code [code]:
   263   "stirling_row 0 = [1]"
   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]] @ [1]"
   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]] @ [1] =
   293       zip_with_prev (\<lambda>a b. a + n * b) 0 (stirling_row n) @ [1]"
   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