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