Grace/Cumulus/Combinatorics/Stirling.g17

Open in Workbench Download source

Imports

BinomialMore

Declarations

Source

1-- Stirling numbers of the second kind by their triangle recurrence2-- S(n + 1, k + 1) = S(n, k) + (k + 1) · S(n, k + 1): the first and last3-- columns, vanishing above the diagonal, S(n + 1, 2) = 2 ^ n − 1,4-- S(n + 2, n + 1) = (n + 2 choose 2), positivity, and the Bell numbers5-- as row sums6-- Part of Cumulus, the Grace17 mathematics library.7-- Copyright (c) 2026 Grant Passmore.8-- Released under the Apache License 2.0. See LICENSE in the repository root.9-- Authors: Grant Passmore1011module Grace.Cumulus.Combinatorics.Stirling12import BinomialMore1314-- Row 0 of the Stirling triangle: 1 at k = 0, then 0.15def stirlingZero (k : Nat) → Nat16  | zero => 117  | succ(k) => 01819-- The recurrence at column k from the entries above and above-left: the20-- entry at column 0 is 0, and the entry at column succ(j) is the entry at21-- column j of the row above plus succ(j) times the entry at column succ(j).22def stirlingStep (k above left : Nat) → Nat23  | zero, above, left => 024  | succ(j), above, left => add(left, mul(succ(j), above))2526-- The number of partitions of an n-element set into k blocks.27def stirling (n k : Nat) → Nat28  | zero, k => stirlingZero(k)29  | succ(n), k => stirlingStep(k, stirling(n, k), stirling(n, predecessor(k)))30  descends n3132theorem stirling.zero_zero33  shows stirling(0, 0) = 134proof35  compute3637theorem stirling.zero_succ (k : Nat)38  shows stirling(0, k + 1) = 039proof40  have hk : k + 1 = succ(k) := add.right_one with *41  rw hk42  simp [stirling, stirlingZero]4344theorem stirling.succ_zero (n : Nat)45  shows stirling(n + 1, 0) = 046proof47  have hn : n + 1 = succ(n) := add.right_one with *48  rw hn49  simp [stirling, stirlingStep]5051theorem stirling.succ_succ (n k : Nat)52  shows stirling(n + 1, k + 1) = stirling(n, k) + (k + 1) * stirling(n, k + 1)53proof54  have hn : n + 1 = succ(n) := add.right_one with *55  have hk : k + 1 = succ(k) := add.right_one with *56  rw [hn, hk]57  simp [stirling, stirlingStep, predecessor]5859theorem stirling.natural (n : Nat)60  shows ∀ k, Nat(k) → Nat(stirling(n, k))61proof62  induction n63  | zero =>64    intros k kNat65    have split := natural.zero_or_successor at k66    cases split with isZero | isSucc67    . rw [isZero, stirling.zero_zero]68      typecheck69    . obtain ⟨j, jNat, shape⟩ := isSucc70      have hj : succ(j) = j + 171        symm72        exact add.right_one at j73      rw [shape, hj, stirling.zero_succ]74      typecheck75  | m + 1 using ih =>76    intros k kNat77    have split := natural.zero_or_successor at k78    cases split with isZero | isSucc79    . rw [isZero, stirling.succ_zero]80      typecheck81    . obtain ⟨j, jNat, shape⟩ := isSucc82      have hj : succ(j) = j + 183        symm84        exact add.right_one at j85      rw [shape, hj, stirling.succ_succ]86      have first := ih with j, jNat87      have nextNat : Nat(j + 1) by typecheck88      have second := ih with j + 1, nextNat89      typecheck9091theorem stirling.natural_at (n k : Nat)92  shows Nat(stirling(n, k))93proof94  have all := stirling.natural at n95  exact all with k, assumption9697theorem stirling.three_two98  shows stirling(3, 2) = 399proof100  compute101102theorem stirling.four_two103  shows stirling(4, 2) = 7104proof105  compute106107theorem stirling.four_three108  shows stirling(4, 3) = 6109proof110  compute111112theorem stirling.five_two113  shows stirling(5, 2) = 15114proof115  compute116117theorem stirling.five_three118  shows stirling(5, 3) = 25119proof120  compute121122theorem stirling.five_four123  shows stirling(5, 4) = 10124proof125  compute126127-- One block: S(n + 1, 1) = 1.128theorem stirling.one (n : Nat)129  shows stirling(n + 1, 1) = 1130proof131  induction n132  | zero => compute133  | k + 1 using ih =>134    have step := stirling.succ_succ at k + 1, 0135    have one : 0 + 1 = 1 by compute136    rw one at step137    rw [step, stirling.succ_zero, ih]138    compute139140-- Nothing above the diagonal.141theorem stirling.above (n : Nat)142  shows ∀ d, Nat(d) → stirling(n, n + 1 + d) = 0143proof144  induction n145  | zero =>146    intros d dNat147    have shape : 0 + 1 + d = d + 1 by ring148    rw shape149    exact stirling.zero_succ at d150  | k + 1 using ih =>151    intros d dNat152    have shape : k + 1 + 1 + d = (k + 1 + d) + 1 by ring153    rw [shape, stirling.succ_succ]154    have first := ih with d, dNat155    have nextNat : Nat(d + 1) by typecheck156    have second := ih with d + 1, nextNat157    have shape2 : k + 1 + (d + 1) = k + 1 + d + 1 by ring158    rw shape2 at second159    rw [first, second]160    ring161162theorem stirling.above_diagonal (n k : Nat)163  requires below : n < k164  shows stirling(n, k) = 0165proof166  obtain ⟨d, dNat, shape⟩ := below167  have form : k = n + 1 + d168    rw ← shape169    ring170  rw form171  have all := stirling.above at n172  exact all with d, dNat173174-- Singletons: S(n, n) = 1.175theorem stirling.diagonal (n : Nat)176  shows stirling(n, n) = 1177proof178  induction n179  | zero => compute180  | k + 1 using ih =>181    rw [stirling.succ_succ, ih]182    have below : k < k + 1 by arith183    have vanish := stirling.above_diagonal with below184    rw vanish185    ring186187-- Two blocks: S(n + 1, 2) = 2 ^ n − 1.188theorem stirling.two (n : Nat)189  shows stirling(n + 1, 2) + 1 = 2 ^ n190proof191  induction n192  | zero => compute193  | k + 1 using ih =>194    have step := stirling.succ_succ at k + 1, 1195    have two : 1 + 1 = 2 by compute196    rw two at step197    rw [step, stirling.one, pow.add_one]198    have sNat := stirling.natural_at at k + 1, 2199    have powNat := pow.closed at 2, k200    ring [ih]201202-- n + 1 blocks of n + 2 points: S(n + 2, n + 1) = (n + 2 choose 2).203theorem stirling.sub_one (n : Nat)204  shows 2 * stirling(n + 2, n + 1) = (n + 2) * (n + 1)205proof206  induction n207  | zero => compute208  | k + 1 using ih =>209    have shape : k + 1 + 2 = (k + 2) + 1 by ring210    have shape2 : k + 1 + 1 = k + 2 by ring211    rw [shape, stirling.succ_succ, shape2, stirling.diagonal]212    have sNat := stirling.natural_at at k + 2, k + 1213    ring [ih]214215theorem stirling.sub_one_choose (n : Nat)216  shows stirling(n + 2, n + 1) = (n + 2) choose 2217proof218  have first := stirling.sub_one at n219  have second := choose.two at n + 1220  have shape : n + 1 + 1 = n + 2 by ring221  rw shape at second222  rw ← second at first223  have sNat := stirling.natural_at at n + 2, n + 1224  have cNat := choose.natural at n + 2, 2225  have twoPositive : 0 < 2 by arith226  have two := Less.ne_zero_of_positive with twoPositive227  exact mul.left_cancel_positive with two, first228229-- Positivity on and below the diagonal.230theorem stirling.positive (n : Nat)231  shows ∀ k, Nat(k) → k ≤ n → 0 < stirling(n + 1, k + 1)232proof233  induction n234  | zero =>235    intros k kNat le236    have split := natural.zero_or_successor at k237    cases split with isZero | isSucc238    . rw isZero239      have one : 0 + 1 = 1 by compute240      have value : stirling(1, 1) = 1 by compute241      rw [one, value]242      arith243    . obtain ⟨i, iNat, shape⟩ := isSucc244      exfalso245      rw shape at le246      arith247  | j + 1 using ih =>248    intros k kNat le249    have split := natural.zero_or_successor at k250    cases split with isZero | isSucc251    . rw isZero252      have one : 0 + 1 = 1 by compute253      rw [one, stirling.one]254      arith255    . obtain ⟨i, iNat, shape⟩ := isSucc256      have hi : succ(i) = i + 1257        symm258        exact add.right_one at i259      rw [shape, hi, stirling.succ_succ]260      rw [shape, hi] at le261      have earlier : i ≤ j by arith262      have inner := ih with i, iNat, earlier263      have restNat := stirling.natural_at at j + 1, i + 1 + 1264      have productNat : Nat((i + 1 + 1) * stirling(j + 1, i + 1 + 1))265        typecheck266      arith267268theorem stirling.positive_at (n k : Nat)269  requires le : k ≤ n270  shows 0 < stirling(n + 1, k + 1)271proof272  have all := stirling.positive at n273  exact all with k, assumption, le274275-- Partial row sums of the Stirling triangle.276def bellSum (n k : Nat) → Nat277  | n, zero => 0278  | n, succ(k) => bellSum(n, k) + stirling(n, k)279  descends k280281-- The Bell numbers: the number of partitions of an n-element set.282def bell (n : Nat) : Nat := bellSum(n, n + 1)283284theorem bellSum.zero (n : Nat)285  shows bellSum(n, 0) = 0286proof287  simp [bellSum]288289theorem bellSum.succ (n k : Nat)290  shows bellSum(n, k + 1) = bellSum(n, k) + stirling(n, k)291proof292  have hk : k + 1 = succ(k) := add.right_one with *293  rw hk294  simp [bellSum]295296theorem bellSum.natural (n k : Nat)297  shows Nat(bellSum(n, k))298proof299  induction k300  | zero =>301    rw bellSum.zero302    typecheck303  | j + 1 using ih =>304    rw bellSum.succ305    have sNat := stirling.natural_at at n, j306    typecheck307308theorem bell.value (n : Nat)309  shows bell(n) = bellSum(n, n + 1)310proof311  simp [bell]312313theorem bell.natural (n : Nat)314  shows Nat(bell(n))315proof316  rw bell.value317  exact bellSum.natural at n, n + 1318319theorem bell.zero320  shows bell(0) = 1321proof322  compute323324theorem bell.one325  shows bell(1) = 1326proof327  compute328329theorem bell.two330  shows bell(2) = 2331proof332  compute333334theorem bell.three335  shows bell(3) = 5336proof337  compute338339theorem bell.four340  shows bell(4) = 15341proof342  compute343344theorem bell.five345  shows bell(5) = 52346proof347  compute348349theorem bell.ge_one (n : Nat)350  shows 1 ≤ bell(n)351proof352  rw [bell.value, bellSum.succ, stirling.diagonal]353  have sumNat := bellSum.natural at n, n354  arith355356-- The row sum is unchanged beyond the diagonal.357theorem bellSum.beyond (n d : Nat)358  shows bellSum(n, n + 1 + d) = bell(n)359proof360  induction d361  | zero =>362    have shape : n + 1 + 0 = n + 1 by ring363    rw [shape, bell.value]364    rfl365  | j + 1 using ih =>366    have shape : n + 1 + (j + 1) = (n + 1 + j) + 1 by ring367    rw [shape, bellSum.succ, ih]368    have below : n < n + 1 + j by arith369    have vanish := stirling.above_diagonal with below370    rw vanish371    have bNat := bell.natural at n372    ring373