src/HOL/Library/Stirling.thy
author haftmann
Fri Mar 22 19:18:08 2019 +0000 (3 months ago)
changeset 69946 494934c30f38
parent 68975 5ce4d117cea7
child 70113 c8deb8ba6d05
permissions -rw-r--r--
improved code equations taken over from AFP
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@68406
   198
      (simp add: algebra_simps sum.distrib sum_distrib_left pochhammer_Suc flip: Suc)
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
hoelzl@63071
   218
subsubsection \<open>Efficient code\<close>
hoelzl@63071
   219
hoelzl@63071
   220
text \<open>
wenzelm@63487
   221
  Naively using the defining equations of the Stirling numbers of the first
wenzelm@63487
   222
  kind to compute them leads to exponential run time due to repeated
wenzelm@63487
   223
  computations. We can use memoisation to compute them row by row without
wenzelm@63487
   224
  repeating computations, at the cost of computing a few unneeded values.
hoelzl@63071
   225
hoelzl@63071
   226
  As a bonus, this is very efficient for applications where an entire row of
wenzelm@63487
   227
  Stirling numbers is needed.
hoelzl@63071
   228
\<close>
hoelzl@63071
   229
wenzelm@63487
   230
definition zip_with_prev :: "('a \<Rightarrow> 'a \<Rightarrow> 'b) \<Rightarrow> 'a \<Rightarrow> 'a list \<Rightarrow> 'b list"
nipkow@66656
   231
  where "zip_with_prev f x xs = map2 f (x # xs) xs"
hoelzl@63071
   232
hoelzl@63071
   233
lemma zip_with_prev_altdef:
hoelzl@63071
   234
  "zip_with_prev f x xs =
wenzelm@63487
   235
    (if xs = [] then [] else f x (hd xs) # [f (xs!i) (xs!(i+1)). i \<leftarrow> [0..<length xs - 1]])"
hoelzl@63071
   236
proof (cases xs)
wenzelm@63487
   237
  case Nil
wenzelm@63487
   238
  then show ?thesis
wenzelm@63487
   239
    by (simp add: zip_with_prev_def)
wenzelm@63487
   240
next
hoelzl@63071
   241
  case (Cons y ys)
wenzelm@63487
   242
  then have "zip_with_prev f x xs = f x (hd xs) # zip_with_prev f y ys"
hoelzl@63071
   243
    by (simp add: zip_with_prev_def)
hoelzl@63071
   244
  also have "zip_with_prev f y ys = map (\<lambda>i. f (xs ! i) (xs ! (i + 1))) [0..<length xs - 1]"
hoelzl@63071
   245
    unfolding Cons
wenzelm@63487
   246
    by (induct ys arbitrary: y)
nipkow@68406
   247
      (simp_all add: zip_with_prev_def upt_conv_Cons flip: map_Suc_upt del: upt_Suc)
wenzelm@63487
   248
  finally show ?thesis
wenzelm@63487
   249
    using Cons by simp
wenzelm@63487
   250
qed
hoelzl@63071
   251
hoelzl@63071
   252
wenzelm@63487
   253
primrec stirling_row_aux
wenzelm@63487
   254
  where
wenzelm@63487
   255
    "stirling_row_aux n y [] = [1]"
wenzelm@63487
   256
  | "stirling_row_aux n y (x#xs) = (y + n * x) # stirling_row_aux n x xs"
hoelzl@63071
   257
hoelzl@63071
   258
lemma stirling_row_aux_correct:
hoelzl@63071
   259
  "stirling_row_aux n y xs = zip_with_prev (\<lambda>a b. a + n * b) y xs @ [1]"
wenzelm@63487
   260
  by (induct xs arbitrary: y) (simp_all add: zip_with_prev_def)
hoelzl@63071
   261
hoelzl@63071
   262
lemma stirling_row_code [code]:
hoelzl@63071
   263
  "stirling_row 0 = [1]"
hoelzl@63071
   264
  "stirling_row (Suc n) = stirling_row_aux n 0 (stirling_row n)"
wenzelm@63487
   265
proof goal_cases
wenzelm@63487
   266
  case 1
wenzelm@63487
   267
  show ?case by (simp add: stirling_row_def)
wenzelm@63487
   268
next
wenzelm@63487
   269
  case 2
hoelzl@63071
   270
  have "stirling_row (Suc n) =
wenzelm@63487
   271
    0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ [1]"
lp15@68975
   272
  proof (rule nth_equalityI, goal_cases length nth)
hoelzl@63071
   273
    case (nth i)
wenzelm@63487
   274
    from nth have "i \<le> Suc n"
wenzelm@63487
   275
      by simp
wenzelm@63487
   276
    then consider "i = 0 \<or> i = Suc n" | "i > 0" "i \<le> n"
wenzelm@63487
   277
      by linarith
wenzelm@63487
   278
    then show ?case
hoelzl@63071
   279
    proof cases
wenzelm@63487
   280
      case 1
wenzelm@63487
   281
      then show ?thesis
wenzelm@63487
   282
        by (auto simp: nth_stirling_row nth_append)
wenzelm@63487
   283
    next
wenzelm@63487
   284
      case 2
wenzelm@63487
   285
      then show ?thesis
wenzelm@63487
   286
        by (cases i) (simp_all add: nth_append nth_stirling_row)
wenzelm@63487
   287
    qed
wenzelm@63487
   288
  next
wenzelm@63487
   289
    case length
wenzelm@63487
   290
    then show ?case by simp
wenzelm@63487
   291
  qed
hoelzl@63071
   292
  also have "0 # [stirling_row n ! i + stirling_row n ! (i+1) * n. i \<leftarrow> [0..<n]] @ [1] =
wenzelm@63487
   293
      zip_with_prev (\<lambda>a b. a + n * b) 0 (stirling_row n) @ [1]"
hoelzl@63071
   294
    by (cases n) (auto simp add: zip_with_prev_altdef stirling_row_def hd_map simp del: upt_Suc)
hoelzl@63071
   295
  also have "\<dots> = stirling_row_aux n 0 (stirling_row n)"
hoelzl@63071
   296
    by (simp add: stirling_row_aux_correct)
wenzelm@63487
   297
  finally show ?case .
wenzelm@63487
   298
qed
hoelzl@63071
   299
hoelzl@63071
   300
lemma stirling_code [code]:
wenzelm@63487
   301
  "stirling n k =
wenzelm@63487
   302
    (if k = 0 then (if n = 0 then 1 else 0)
wenzelm@63487
   303
     else if k > n then 0
wenzelm@63487
   304
     else if k = n then 1
wenzelm@63487
   305
     else stirling_row n ! k)"
hoelzl@63071
   306
  by (simp add: nth_stirling_row)
hoelzl@63071
   307
hoelzl@63071
   308
end