src/HOL/Library/Stirling.thy
author wenzelm
Wed May 25 11:49:40 2016 +0200 (2016-05-25)
changeset 63145 703edebd1d92
parent 63071 3ca3bc795908
child 63487 6e29fb72e659
permissions -rw-r--r--
isabelle update_cartouches -c -t;
hoelzl@63071
     1
(* Authors: Amine Chaieb & Florian Haftmann, TU Muenchen
hoelzl@63071
     2
            with contributions by Lukas Bulwahn and Manuel Eberl*)
hoelzl@63071
     3
wenzelm@63145
     4
section \<open>Stirling numbers of first and second kind\<close>
hoelzl@63071
     5
hoelzl@63071
     6
theory Stirling
hoelzl@63071
     7
imports Binomial
hoelzl@63071
     8
begin
hoelzl@63071
     9
wenzelm@63145
    10
subsection \<open>Stirling numbers of the second kind\<close>
hoelzl@63071
    11
hoelzl@63071
    12
fun Stirling :: "nat \<Rightarrow> nat \<Rightarrow> nat"
hoelzl@63071
    13
where
hoelzl@63071
    14
  "Stirling 0 0 = 1"
hoelzl@63071
    15
| "Stirling 0 (Suc k) = 0"
hoelzl@63071
    16
| "Stirling (Suc n) 0 = 0"
hoelzl@63071
    17
| "Stirling (Suc n) (Suc k) = Suc k * Stirling n (Suc k) + Stirling n k"
hoelzl@63071
    18
hoelzl@63071
    19
lemma Stirling_1 [simp]:
hoelzl@63071
    20
  "Stirling (Suc n) (Suc 0) = 1"
hoelzl@63071
    21
  by (induct n) simp_all
hoelzl@63071
    22
hoelzl@63071
    23
lemma Stirling_less [simp]:
hoelzl@63071
    24
  "n < k \<Longrightarrow> Stirling n k = 0"
hoelzl@63071
    25
  by (induct n k rule: Stirling.induct) simp_all
hoelzl@63071
    26
hoelzl@63071
    27
lemma Stirling_same [simp]:
hoelzl@63071
    28
  "Stirling n n = 1"
hoelzl@63071
    29
  by (induct n) simp_all
hoelzl@63071
    30
hoelzl@63071
    31
lemma Stirling_2_2:
hoelzl@63071
    32
  "Stirling (Suc (Suc n)) (Suc (Suc 0)) = 2 ^ Suc n - 1"
hoelzl@63071
    33
proof (induct n)
hoelzl@63071
    34
  case 0 then show ?case by simp
hoelzl@63071
    35
next
hoelzl@63071
    36
  case (Suc n)
hoelzl@63071
    37
  have "Stirling (Suc (Suc (Suc n))) (Suc (Suc 0)) =
hoelzl@63071
    38
    2 * Stirling (Suc (Suc n)) (Suc (Suc 0)) + Stirling (Suc (Suc n)) (Suc 0)" by simp
hoelzl@63071
    39
  also have "\<dots> = 2 * (2 ^ Suc n - 1) + 1"
hoelzl@63071
    40
    by (simp only: Suc Stirling_1)
hoelzl@63071
    41
  also have "\<dots> = 2 ^ Suc (Suc n) - 1"
hoelzl@63071
    42
  proof -
hoelzl@63071
    43
    have "(2::nat) ^ Suc n - 1 > 0" by (induct n) simp_all
hoelzl@63071
    44
    then have "2 * ((2::nat) ^ Suc n - 1) > 0" by simp
hoelzl@63071
    45
    then have "2 \<le> 2 * ((2::nat) ^ Suc n)" by simp
hoelzl@63071
    46
    with add_diff_assoc2 [of 2 "2 * 2 ^ Suc n" 1]
hoelzl@63071
    47
      have "2 * 2 ^ Suc n - 2 + (1::nat) = 2 * 2 ^ Suc n + 1 - 2" .
hoelzl@63071
    48
    then show ?thesis by (simp add: nat_distrib)
hoelzl@63071
    49
  qed
hoelzl@63071
    50
  finally show ?case by simp
hoelzl@63071
    51
qed
hoelzl@63071
    52
hoelzl@63071
    53
lemma Stirling_2:
hoelzl@63071
    54
  "Stirling (Suc n) (Suc (Suc 0)) = 2 ^ n - 1"
hoelzl@63071
    55
  using Stirling_2_2 by (cases n) simp_all
hoelzl@63071
    56
wenzelm@63145
    57
subsection \<open>Stirling numbers of the first kind\<close>
hoelzl@63071
    58
hoelzl@63071
    59
fun stirling :: "nat \<Rightarrow> nat \<Rightarrow> nat"
hoelzl@63071
    60
where
hoelzl@63071
    61
  "stirling 0 0 = 1"
hoelzl@63071
    62
| "stirling 0 (Suc k) = 0"
hoelzl@63071
    63
| "stirling (Suc n) 0 = 0"
hoelzl@63071
    64
| "stirling (Suc n) (Suc k) = n * stirling n (Suc k) + stirling n k"
hoelzl@63071
    65
hoelzl@63071
    66
lemma stirling_0 [simp]: "n > 0 \<Longrightarrow> stirling n 0 = 0"
hoelzl@63071
    67
  by (cases n) simp_all
hoelzl@63071
    68
hoelzl@63071
    69
lemma stirling_less [simp]:
hoelzl@63071
    70
  "n < k \<Longrightarrow> stirling n k = 0"
hoelzl@63071
    71
  by (induct n k rule: stirling.induct) simp_all
hoelzl@63071
    72
hoelzl@63071
    73
lemma stirling_same [simp]:
hoelzl@63071
    74
  "stirling n n = 1"
hoelzl@63071
    75
  by (induct n) simp_all
hoelzl@63071
    76
hoelzl@63071
    77
lemma stirling_Suc_n_1:
hoelzl@63071
    78
  "stirling (Suc n) (Suc 0) = fact n"
hoelzl@63071
    79
  by (induct n) auto
hoelzl@63071
    80
hoelzl@63071
    81
lemma stirling_Suc_n_n:
hoelzl@63071
    82
  shows "stirling (Suc n) n = Suc n choose 2"
hoelzl@63071
    83
by (induct n) (auto simp add: numerals(2))
hoelzl@63071
    84
hoelzl@63071
    85
lemma stirling_Suc_n_2:
hoelzl@63071
    86
  assumes "n \<ge> Suc 0"
hoelzl@63071
    87
  shows "stirling (Suc n) 2 = (\<Sum>k=1..n. fact n div k)"
hoelzl@63071
    88
using assms
hoelzl@63071
    89
proof (induct n)
hoelzl@63071
    90
  case 0 from this show ?case by simp
hoelzl@63071
    91
next
hoelzl@63071
    92
  case (Suc n)
hoelzl@63071
    93
  show ?case
hoelzl@63071
    94
  proof (cases n)
hoelzl@63071
    95
    case 0 from this show ?thesis by (simp add: numerals(2))
hoelzl@63071
    96
  next
hoelzl@63071
    97
    case Suc
hoelzl@63071
    98
    from this have geq1: "Suc 0 \<le> n" by simp
hoelzl@63071
    99
    have "stirling (Suc (Suc n)) 2 = Suc n * stirling (Suc n) 2 + stirling (Suc n) (Suc 0)"
hoelzl@63071
   100
      by (simp only: stirling.simps(4)[of "Suc n"] numerals(2))
hoelzl@63071
   101
    also have "... = Suc n * (\<Sum>k=1..n. fact n div k) + fact n"
hoelzl@63071
   102
      using Suc.hyps[OF geq1]
hoelzl@63071
   103
      by (simp only: stirling_Suc_n_1 of_nat_fact of_nat_add of_nat_mult)
hoelzl@63071
   104
    also have "... = Suc n * (\<Sum>k=1..n. fact n div k) + Suc n * fact n div Suc n"
hoelzl@63071
   105
      by (metis nat.distinct(1) nonzero_mult_divide_cancel_left)
hoelzl@63071
   106
    also have "... = (\<Sum>k=1..n. fact (Suc n) div k) + fact (Suc n) div Suc n"
hoelzl@63071
   107
      by (simp add: setsum_right_distrib div_mult_swap dvd_fact)
hoelzl@63071
   108
    also have "... = (\<Sum>k=1..Suc n. fact (Suc n) div k)" by simp
hoelzl@63071
   109
    finally show ?thesis .
hoelzl@63071
   110
  qed
hoelzl@63071
   111
qed
hoelzl@63071
   112
hoelzl@63071
   113
lemma of_nat_stirling_Suc_n_2:
hoelzl@63071
   114
  assumes "n \<ge> Suc 0"
hoelzl@63071
   115
  shows "(of_nat (stirling (Suc n) 2)::'a::field_char_0) = fact n * (\<Sum>k=1..n. (1 / of_nat k))"
hoelzl@63071
   116
using assms
hoelzl@63071
   117
proof (induct n)
hoelzl@63071
   118
  case 0 from this show ?case by simp
hoelzl@63071
   119
next
hoelzl@63071
   120
  case (Suc n)
hoelzl@63071
   121
  show ?case
hoelzl@63071
   122
  proof (cases n)
hoelzl@63071
   123
    case 0 from this show ?thesis by (auto simp add: numerals(2))
hoelzl@63071
   124
  next
hoelzl@63071
   125
    case Suc
hoelzl@63071
   126
    from this have geq1: "Suc 0 \<le> n" by simp
hoelzl@63071
   127
    have "(of_nat (stirling (Suc (Suc n)) 2)::'a) =
hoelzl@63071
   128
      of_nat (Suc n * stirling (Suc n) 2 + stirling (Suc n) (Suc 0))"
hoelzl@63071
   129
      by (simp only: stirling.simps(4)[of "Suc n"] numerals(2))
hoelzl@63071
   130
    also have "... = of_nat (Suc n) * (fact n * (\<Sum>k = 1..n. 1 / of_nat k)) + fact n"
hoelzl@63071
   131
      using Suc.hyps[OF geq1]
hoelzl@63071
   132
      by (simp only: stirling_Suc_n_1 of_nat_fact of_nat_add of_nat_mult)
hoelzl@63071
   133
    also have "... = fact (Suc n) * (\<Sum>k = 1..n. 1 / of_nat k) + fact (Suc n) * (1 / of_nat (Suc n))"
hoelzl@63071
   134
      using of_nat_neq_0 by auto
hoelzl@63071
   135
    also have "... = fact (Suc n) * (\<Sum>k = 1..Suc n. 1 / of_nat k)"
hoelzl@63071
   136
      by (simp add: distrib_left)
hoelzl@63071
   137
    finally show ?thesis .
hoelzl@63071
   138
  qed
hoelzl@63071
   139
qed
hoelzl@63071
   140
hoelzl@63071
   141
lemma setsum_stirling:
hoelzl@63071
   142
  "(\<Sum>k\<le>n. stirling n k) = fact n"
hoelzl@63071
   143
proof (induct n)
hoelzl@63071
   144
  case 0
hoelzl@63071
   145
  from this show ?case by simp
hoelzl@63071
   146
next
hoelzl@63071
   147
  case (Suc n)
hoelzl@63071
   148
  have "(\<Sum>k\<le>Suc n. stirling (Suc n) k) = stirling (Suc n) 0 + (\<Sum>k\<le>n. stirling (Suc n) (Suc k))"
hoelzl@63071
   149
    by (simp only: setsum_atMost_Suc_shift)
hoelzl@63071
   150
  also have "\<dots> = (\<Sum>k\<le>n. stirling (Suc n) (Suc k))" by simp
hoelzl@63071
   151
  also have "\<dots> = (\<Sum>k\<le>n. n * stirling n (Suc k) + stirling n k)" by simp
hoelzl@63071
   152
  also have "\<dots> = n * (\<Sum>k\<le>n. stirling n (Suc k)) + (\<Sum>k\<le>n. stirling n k)"
hoelzl@63071
   153
    by (simp add: setsum.distrib setsum_right_distrib)
hoelzl@63071
   154
  also have "\<dots> = n * fact n + fact n"
hoelzl@63071
   155
  proof -
hoelzl@63071
   156
    have "n * (\<Sum>k\<le>n. stirling n (Suc k)) = n * ((\<Sum>k\<le>Suc n. stirling n k) - stirling n 0)"
hoelzl@63071
   157
      by (metis add_diff_cancel_left' setsum_atMost_Suc_shift)
hoelzl@63071
   158
    also have "\<dots> = n * (\<Sum>k\<le>n. stirling n k)" by (cases n) simp+
hoelzl@63071
   159
    also have "\<dots> = n * fact n" using Suc.hyps by simp
hoelzl@63071
   160
    finally have "n * (\<Sum>k\<le>n. stirling n (Suc k)) = n * fact n" .
hoelzl@63071
   161
    moreover have "(\<Sum>k\<le>n. stirling n k) = fact n" using Suc.hyps .
hoelzl@63071
   162
    ultimately show ?thesis by simp
hoelzl@63071
   163
  qed
hoelzl@63071
   164
  also have "\<dots> = fact (Suc n)" by simp
hoelzl@63071
   165
  finally show ?case .
hoelzl@63071
   166
qed
hoelzl@63071
   167
hoelzl@63071
   168
lemma stirling_pochhammer:
hoelzl@63071
   169
  "(\<Sum>k\<le>n. of_nat (stirling n k) * x ^ k) = (pochhammer x n :: 'a :: comm_semiring_1)"
hoelzl@63071
   170
proof (induction n)
hoelzl@63071
   171
  case (Suc n)
hoelzl@63071
   172
  have "of_nat (n * stirling n 0) = (0 :: 'a)" by (cases n) simp_all
hoelzl@63071
   173
  hence "(\<Sum>k\<le>Suc n. of_nat (stirling (Suc n) k) * x ^ k) =
hoelzl@63071
   174
            (of_nat (n * stirling n 0) * x ^ 0 +
hoelzl@63071
   175
            (\<Sum>i\<le>n. of_nat (n * stirling n (Suc i)) * (x ^ Suc i))) +
hoelzl@63071
   176
            (\<Sum>i\<le>n. of_nat (stirling n i) * (x ^ Suc i))"
hoelzl@63071
   177
    by (subst setsum_atMost_Suc_shift) (simp add: setsum.distrib ring_distribs)
hoelzl@63071
   178
  also have "\<dots> = pochhammer x (Suc n)"
hoelzl@63071
   179
    by (subst setsum_atMost_Suc_shift [symmetric])
hoelzl@63071
   180
       (simp add: algebra_simps setsum.distrib setsum_right_distrib pochhammer_Suc Suc [symmetric])
hoelzl@63071
   181
  finally show ?case .
hoelzl@63071
   182
qed simp_all
hoelzl@63071
   183
hoelzl@63071
   184
hoelzl@63071
   185
text \<open>A row of the Stirling number triangle\<close>
hoelzl@63071
   186
hoelzl@63071
   187
definition stirling_row :: "nat \<Rightarrow> nat list" where
hoelzl@63071
   188
  "stirling_row n = [stirling n k. k \<leftarrow> [0..<Suc n]]"
hoelzl@63071
   189
hoelzl@63071
   190
lemma nth_stirling_row: "k \<le> n \<Longrightarrow> stirling_row n ! k = stirling n k"
hoelzl@63071
   191
  by (simp add: stirling_row_def del: upt_Suc)
hoelzl@63071
   192
hoelzl@63071
   193
lemma length_stirling_row [simp]: "length (stirling_row n) = Suc n"
hoelzl@63071
   194
  by (simp add: stirling_row_def)
hoelzl@63071
   195
hoelzl@63071
   196
lemma stirling_row_nonempty [simp]: "stirling_row n \<noteq> []"
hoelzl@63071
   197
  using length_stirling_row[of n] by (auto simp del: length_stirling_row)
hoelzl@63071
   198
hoelzl@63071
   199
(* TODO Move *)
hoelzl@63071
   200
lemma list_ext:
hoelzl@63071
   201
  assumes "length xs = length ys"
hoelzl@63071
   202
  assumes "\<And>i. i < length xs \<Longrightarrow> xs ! i = ys ! i"
hoelzl@63071
   203
  shows   "xs = ys"
hoelzl@63071
   204
using assms
hoelzl@63071
   205
proof (induction rule: list_induct2)
hoelzl@63071
   206
  case (Cons x xs y ys)
hoelzl@63071
   207
  from Cons.prems[of 0] have "x = y" by simp
hoelzl@63071
   208
  moreover from Cons.prems[of "Suc i" for i] have "xs = ys" by (intro Cons.IH) simp
hoelzl@63071
   209
  ultimately show ?case by simp
hoelzl@63071
   210
qed simp_all
hoelzl@63071
   211
hoelzl@63071
   212
hoelzl@63071
   213
subsubsection \<open>Efficient code\<close>
hoelzl@63071
   214
hoelzl@63071
   215
text \<open>
hoelzl@63071
   216
  Naively using the defining equations of the Stirling numbers of the first kind to
hoelzl@63071
   217
  compute them leads to exponential run time due to repeated computations.
hoelzl@63071
   218
  We can use memoisation to compute them row by row without repeating computations, at
hoelzl@63071
   219
  the cost of computing a few unneeded values.
hoelzl@63071
   220
hoelzl@63071
   221
  As a bonus, this is very efficient for applications where an entire row of
hoelzl@63071
   222
  Stirling numbers is needed..
hoelzl@63071
   223
\<close>
hoelzl@63071
   224
hoelzl@63071
   225
definition zip_with_prev :: "('a \<Rightarrow> 'a \<Rightarrow> 'b) \<Rightarrow> 'a \<Rightarrow> 'a list \<Rightarrow> 'b list" where
hoelzl@63071
   226
  "zip_with_prev f x xs = map (\<lambda>(x,y). f x y) (zip (x # xs) xs)"
hoelzl@63071
   227
hoelzl@63071
   228
lemma zip_with_prev_altdef:
hoelzl@63071
   229
  "zip_with_prev f x xs =
hoelzl@63071
   230
     (if xs = [] then [] else f x (hd xs) # [f (xs!i) (xs!(i+1)). i \<leftarrow> [0..<length xs - 1]])"
hoelzl@63071
   231
proof (cases xs)
hoelzl@63071
   232
  case (Cons y ys)
hoelzl@63071
   233
  hence "zip_with_prev f x xs = f x (hd xs) # zip_with_prev f y ys"
hoelzl@63071
   234
    by (simp add: zip_with_prev_def)
hoelzl@63071
   235
  also have "zip_with_prev f y ys = map (\<lambda>i. f (xs ! i) (xs ! (i + 1))) [0..<length xs - 1]"
hoelzl@63071
   236
    unfolding Cons
hoelzl@63071
   237
    by (induction ys arbitrary: y)
hoelzl@63071
   238
       (simp_all add: zip_with_prev_def upt_conv_Cons map_Suc_upt [symmetric] del: upt_Suc)
hoelzl@63071
   239
  finally show ?thesis using Cons by simp
hoelzl@63071
   240
qed (simp add: zip_with_prev_def)
hoelzl@63071
   241
hoelzl@63071
   242
hoelzl@63071
   243
primrec stirling_row_aux where
hoelzl@63071
   244
  "stirling_row_aux n y [] = [1]"
hoelzl@63071
   245
| "stirling_row_aux n y (x#xs) = (y + n * x) # stirling_row_aux n x xs"
hoelzl@63071
   246
hoelzl@63071
   247
lemma stirling_row_aux_correct:
hoelzl@63071
   248
  "stirling_row_aux n y xs = zip_with_prev (\<lambda>a b. a + n * b) y xs @ [1]"
hoelzl@63071
   249
  by (induction xs arbitrary: y) (simp_all add: zip_with_prev_def)
hoelzl@63071
   250
hoelzl@63071
   251
lemma stirling_row_code [code]:
hoelzl@63071
   252
  "stirling_row 0 = [1]"
hoelzl@63071
   253
  "stirling_row (Suc n) = stirling_row_aux n 0 (stirling_row n)"
hoelzl@63071
   254
proof -
hoelzl@63071
   255
  have "stirling_row (Suc n) =
hoelzl@63071
   256
          0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ [1]"
hoelzl@63071
   257
  proof (rule list_ext, goal_cases length nth)
hoelzl@63071
   258
    case (nth i)
hoelzl@63071
   259
    from nth have "i \<le> Suc n" by simp
hoelzl@63071
   260
    then consider "i = 0" | j where "i > 0" "i \<le> n" | "i = Suc n" by linarith
hoelzl@63071
   261
    thus ?case
hoelzl@63071
   262
    proof cases
hoelzl@63071
   263
      assume i: "i > 0" "i \<le> n"
hoelzl@63071
   264
      from this show ?thesis by (cases i) (simp_all add: nth_append nth_stirling_row)
hoelzl@63071
   265
    qed (simp_all add: nth_stirling_row nth_append)
hoelzl@63071
   266
  qed simp
hoelzl@63071
   267
  also have "0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ [1] =
hoelzl@63071
   268
               zip_with_prev (\<lambda>a b. a + n * b) 0 (stirling_row n) @ [1]"
hoelzl@63071
   269
    by (cases n) (auto simp add: zip_with_prev_altdef stirling_row_def hd_map simp del: upt_Suc)
hoelzl@63071
   270
  also have "\<dots> = stirling_row_aux n 0 (stirling_row n)"
hoelzl@63071
   271
    by (simp add: stirling_row_aux_correct)
hoelzl@63071
   272
  finally show "stirling_row (Suc n) = stirling_row_aux n 0 (stirling_row n)" .
hoelzl@63071
   273
qed (simp add: stirling_row_def)
hoelzl@63071
   274
hoelzl@63071
   275
lemma stirling_code [code]:
hoelzl@63071
   276
  "stirling n k = (if k = 0 then if n = 0 then 1 else 0
hoelzl@63071
   277
                   else if k > n then 0 else if k = n then 1
hoelzl@63071
   278
                   else stirling_row n ! k)"
hoelzl@63071
   279
  by (simp add: nth_stirling_row)
hoelzl@63071
   280
hoelzl@63071
   281
end