Stirling Numbers

Why Stirling Numbers Matters

You have four distinct people and want to split them into two nonempty teams. The teams have no names, and the order of people within a team does not matter. There are seven possible divisions. This is not a binomial coefficient: selecting one team and calling everyone else the second team would describe each division twice, and larger numbers of teams need a more general counting method.

The Stirling numbers of the second kind, written S(n, k), count partitions of n distinct objects into exactly k nonempty, unlabeled blocks. They model grouping tasks, assigning objects to nonempty labeled destinations after a suitable conversion, and changing between polynomial bases.

There is also a different family: the unsigned Stirling numbers of the first kind, written c(n, k), count permutations of n distinct objects with exactly k cycles. For four objects and two cycles, the answer is 11, rather than seven. This lesson distinguishes the two families, explains their closely related recurrences, and implements both with the same one-array dynamic programming pattern. The main worked trace uses S(n, k).

Core Concepts and Terminology

A set partition divides a collection into disjoint nonempty blocks whose union is the whole collection. Block order and the order of elements inside each block do not matter. For example, {1, 2} | {3} and {3} | {2, 1} describe the same partition. The three partitions of {1, 2, 3} into two blocks are {1} | {2, 3}, {2} | {1, 3}, and {3} | {1, 2}, so S(3, 2) = 3.

For the first kind, a permutation is viewed as a collection of disjoint cycles. A cycle (1 2 3) sends 1 to 2, 2 to 3, and 3 back to 1. Rotating a cycle's written starting point does not change it, and the order in which disjoint cycles are written does not matter. Reversing a cycle can change the permutation: (1 2 3) and (1 3 2) are different. Thus cycles contain information that ordinary blocks do not.

For both families, n is a nonnegative integer and k is an integer. Define the count to be zero when k < 0 or k > n. The boundary values are S(0, 0) = c(0, 0) = 1, and S(n, 0) = c(n, 0) = 0 for n > 0. There is one way to divide no objects into no blocks, and one empty permutation with no cycles. Both diagonal values S(n, n) and c(n, n) equal one.

For the second-kind recurrence, inspect the block containing the newest element n. If that block has other elements, remove n: a partition into k blocks remains, and there are k blocks into which n could have been inserted. If n is alone, removing its singleton leaves a partition into k − 1 blocks. The cases are disjoint and exhaustive.

For the unsigned first kind, the newest element either starts its own cycle or is inserted after one of the n − 1 old elements in the existing cycles. Inserting it after a chosen old element preserves the number of cycles and is uniquely reversible. This explains why the coefficient changes from k to n − 1.

How Stirling Numbers Works, Step by Step

A two-dimensional table stores an answer for every smaller pair of indices. Since row i depends only on row i − 1, we can instead use a single array of k + 1 entries. Initialize dp[0] = 1 and all other entries to zero. Columns beyond the current row's diagonal represent impossible counts and remain zero until reached.

For each row i from 1 through n, visit j from min(i, k) down to 1. For the second kind, update dp[j] = j × dp[j] + dp[j − 1]. For the unsigned first kind, update dp[j] = (i − 1) × dp[j] + dp[j − 1]. The coefficient is the only change between the two recurrences.

Descending order ensures that both entries on the right still belong to row i − 1. After the inner loop, set dp[0] = 0: a positive number of objects cannot have zero blocks or cycles. Do not reset dp[0] before the first row's inner loop, because dp[1] needs the old empty-case value one.

At the end, return dp[k]. Invalid k values return zero; n < 0 is rejected. The implementations return one immediately on the diagonal k = n, and return zero immediately when k = 0 and n > 0. Every multiplication and sum is reduced modulo M.

For S(4, 2), the previous stored row is [0, 1, 3]. Updating j = 2 gives 2 × 3 + 1 = 7. Updating j = 1 gives 1 × 1 + 0 = 1. The new row prefix is [0, 1, 7]. For c(4, 2), the previous row prefix is [0, 2, 3], and coefficient i − 1 = 3 gives [0, 6, 11].

flowchart TD
    A["Read n, k, and kind"] --> B{"n negative?"}
    B -->|Yes| X["Reject input"]
    B -->|No| C{"k negative or k greater than n?"}
    C -->|Yes| Z["Return 0"]
    C -->|No| D{"k = n?"}
    D -->|Yes| O["Return 1"]
    D -->|No| E{"k = 0?"}
    E -->|Yes| Z
    E -->|No| F["dp = [1, 0, ..., 0]; i = 1"]
    F --> G{"i at most n?"}
    G -->|No| R["Return dp[k]"]
    G -->|Yes| H["j = min(i, k)"]
    H --> I{"j at least 1?"}
    I -->|Yes| J["weight = i - 1 for first kind; otherwise j"]
    J --> K["dp[j] = (weight times dp[j] + dp[j - 1]) mod M"]
    K --> L["j = j - 1"]
    L --> I
    I -->|No| N["dp[0] = 0; i = i + 1"]
    N --> G

Worked Example with a Full Trace

Compute S(4, 2). Allocate dp = [1, 0, 0]. The array stores columns j = 0, 1, 2; each completed row represents a new number of objects.

1. Row i = 1: only j = 1 is visited. Set dp[1] = 1 × 0 + 1 = 1, temporarily giving [1, 1, 0]. Now reset dp[0] = 0, giving [0, 1, 0].
2. Row i = 2: visit j = 2 first. Set dp[2] = 2 × 0 + 1 = 1, giving [0, 1, 1]. Next j = 1: dp[1] = 1 × 1 + 0 = 1. Reset dp[0] = 0; the completed prefix is [0, 1, 1].
3. Row i = 3: j = 2 gives dp[2] = 2 × 1 + 1 = 3, so the array becomes [0, 1, 3]. Then j = 1 gives dp[1] = 1 × 1 + 0 = 1. Reset dp[0] = 0; the completed prefix is [0, 1, 3].
4. Row i = 4: j = 2 gives dp[2] = 2 × 3 + 1 = 7, so the array becomes [0, 1, 7]. Then j = 1 gives dp[1] = 1 × 1 + 0 = 1. Reset dp[0] = 0; the completed prefix is [0, 1, 7].

Return dp[2] = 7, agreeing with the seven explicit partitions in the first diagram. An independent check labels the two blocks: there are 2⁴ − 2 = 14 onto assignments to two labeled teams, and each unlabeled partition has 2! labelings, giving 14/2 = 7.

For comparison, apply the first-kind coefficient i − 1 to the same array pattern. Starting from [1, 0, 0], the completed prefixes are [0, 1, 0], [0, 1, 1], [0, 2, 3], and [0, 6, 11]. Thus c(4, 2) = 11. Independently, a permutation on four labels with two cycles has cycle sizes 3 + 1 or 2 + 2. There are 4 × 2 = 8 permutations of the first type and three pairings of the second type, giving 8 + 3 = 11.

Implementation in Python, C++, and Java

Every panel is a complete standalone program. The function computes one query modulo M = 1,000,000,007. Set first_kind or firstKind to false for S(n, k), or true for the unsigned c(n, k). The sample program calls both modes and prints the same queries in the same order. Signed first-kind numbers are discussed below and are not the direct output of the true mode.

All versions reject n < 0 and return zero for k < 0 or k > n. The diagonal and zero-column cases are resolved before allocating. The remaining calculation keeps k + 1 residues, updates right to left, and clears dp[0] after each row. There is no symmetry rule allowing k to be replaced by min(k, n − k), unlike binomial coefficients.

Python uses arbitrary-precision integers. C++ uses long long and reduces the coefficient modulo M before multiplication. Each factor is therefore below M, so a product is at most (M − 1)² = 1,000,000,012,000,000,036, which fits signed 64-bit arithmetic; adding one further residue still fits. Java uses BigInteger for the product, addition, and remainder. C++ and Java take signed 32-bit n and k; feasible sizes are further limited by O(nk) time and O(k) storage.

MOD = 1_000_000_007


def stirling_mod(n: int, k: int, first_kind: bool = False) -> int:
    if n < 0:
        raise ValueError("n must be nonnegative")
    if k < 0 or k > n:
        return 0
    if k == n:
        return 1
    if k == 0:
        return 0
    dp = [0] * (k + 1)
    dp[0] = 1
    for row in range(1, n + 1):
        for j in range(min(row, k), 0, -1):
            weight = ((row - 1) if first_kind else j) % MOD
            dp[j] = (weight * dp[j] + dp[j - 1]) % MOD
        dp[0] = 0
    return dp[k]


def main() -> None:
    samples = [(0, 0), (3, 0), (3, 4), (4, 2), (5, 3),
               (5, 1), (5, 5), (10, 3), (20, 5)]
    print(f"modulus = {MOD}")
    for n, k in samples:
        second = stirling_mod(n, k, False)
        first = stirling_mod(n, k, True)
        print(f"S({n}, {k}) = {second}; c({n}, {k}) = {first}")


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cstddef>
#include <iostream>
#include <stdexcept>
#include <utility>
#include <vector>

constexpr long long MOD = 1'000'000'007LL;

long long stirling_mod(int n, int k, bool first_kind = false) {
    if (n < 0) {
        throw std::invalid_argument("n must be nonnegative");
    }
    if (k < 0 || k > n) {
        return 0;
    }
    if (k == n) {
        return 1;
    }
    if (k == 0) {
        return 0;
    }
    std::vector<long long> dp(static_cast<std::size_t>(k) + 1, 0);
    dp[0] = 1;
    for (long long row = 1; row <= n; ++row) {
        const int upper = static_cast<int>(std::min(row, static_cast<long long>(k)));
        for (int j = upper; j >= 1; --j) {
            const long long weight = (first_kind ? row - 1 : j) % MOD;
            dp[j] = (weight * dp[j] + dp[j - 1]) % MOD;
        }
        dp[0] = 0;
    }
    return dp[k];
}

int main() {
    const std::vector<std::pair<int, int>> samples = {
        {0, 0}, {3, 0}, {3, 4}, {4, 2}, {5, 3},
        {5, 1}, {5, 5}, {10, 3}, {20, 5}
    };
    std::cout << "modulus = " << MOD << '\n';
    for (const auto& sample : samples) {
        const int n = sample.first;
        const int k = sample.second;
        std::cout << "S(" << n << ", " << k << ") = "
                  << stirling_mod(n, k, false)
                  << "; c(" << n << ", " << k << ") = "
                  << stirling_mod(n, k, true) << '\n';
    }
}

import java.math.BigInteger;
import java.util.Arrays;

public class Main {
    private static final BigInteger MOD = BigInteger.valueOf(1_000_000_007L);

    public static BigInteger stirlingMod(int n, int k, boolean firstKind) {
        if (n < 0) {
            throw new IllegalArgumentException("n must be nonnegative");
        }
        if (k < 0 || k > n) {
            return BigInteger.ZERO;
        }
        if (k == n) {
            return BigInteger.ONE;
        }
        if (k == 0) {
            return BigInteger.ZERO;
        }
        BigInteger[] dp = new BigInteger[Math.addExact(k, 1)];
        Arrays.fill(dp, BigInteger.ZERO);
        dp[0] = BigInteger.ONE;
        for (long row = 1; row <= n; ++row) {
            int upper = (int) Math.min(row, (long) k);
            for (int j = upper; j >= 1; --j) {
                long coefficient = firstKind ? row - 1 : j;
                BigInteger weight = BigInteger.valueOf(coefficient).mod(MOD);
                dp[j] = weight.multiply(dp[j]).add(dp[j - 1]).mod(MOD);
            }
            dp[0] = BigInteger.ZERO;
        }
        return dp[k];
    }

    public static void main(String[] args) {
        int[][] samples = {
            {0, 0}, {3, 0}, {3, 4}, {4, 2}, {5, 3},
            {5, 1}, {5, 5}, {10, 3}, {20, 5}
        };
        System.out.println("modulus = " + MOD);
        for (int[] sample : samples) {
            int n = sample[0];
            int k = sample[1];
            System.out.println("S(" + n + ", " + k + ") = "
                + stirlingMod(n, k, false) + "; c(" + n + ", " + k + ") = "
                + stirlingMod(n, k, true));
        }
    }
}

Run Python with Python 3.9 or later, compile C++ with C++14 or later, and save Java as Main.java before compiling with javac Main.java. Each program prints these exact ten lines:

• modulus = 1000000007
• S(0, 0) = 1; c(0, 0) = 1
• S(3, 0) = 0; c(3, 0) = 0
• S(3, 4) = 0; c(3, 4) = 0
• S(4, 2) = 7; c(4, 2) = 11
• S(5, 3) = 25; c(5, 3) = 35
• S(5, 1) = 1; c(5, 1) = 24
• S(5, 5) = 1; c(5, 5) = 1
• S(10, 3) = 9330; c(10, 3) = 1172700
• S(20, 5) = 206085257; c(20, 5) = 745534512

The last line contains remainders: the full values are S(20, 5) = 749206090500 and c(20, 5) = 371384787345228000. To obtain full counts in Python or Java, remove the coefficient reduction as well as the result reduction and retain arbitrary-precision arithmetic. C++ requires arbitrary-precision arithmetic or a separate proof that every intermediate value fits. Removing only the final remainder operation while leaving a reduced coefficient would not produce general full counts.

Complexity, Correctness, and Invariants

For a nonboundary query with 1 ≤ k < n, the number of updates is min(1, k) + min(2, k) + … + min(n, k) = nk − k(k − 1)/2. Both modes therefore take Θ(nk) time and Θ(k + 1) auxiliary storage under fixed-modulus, unit-cost arithmetic. Computing both kinds separately doubles the work but does not change the asymptotic bounds.

For fixed n and k in that range, best, average, and worst cases have the same update count: the loop schedule does not depend on the values stored. The implemented invalid-k, diagonal, and zero-column shortcuts take Θ(1) time and storage. Across k values for fixed n, the nonboundary running time ranges from Θ(n) when k = 1 to Θ(n²) when k is proportional to n. No input distribution is needed for these deterministic bounds.

For unrestricted full counts, the entries and coefficients require a separate bit-cost analysis. Integer multiplication and addition no longer have constant cost as their operands grow. The bounds above count fixed-modulus arithmetic operations and stored residues, rather than the cost of manipulating arbitrary-size full counts.

Correctness combines the combinatorial recurrence with an array invariant. Initially dp[0] = 1 correctly represents row zero, and higher entries are zero. Assume the array represents row i − 1 before a new row begins. At an update of column j, descending order keeps dp[j] and dp[j − 1] unchanged from that old row. The selected recurrence coefficient therefore produces the correct entry for row i. Processed columns belong to the new row, while the operands for unprocessed columns still belong to the old row. Finally, resetting dp[0] = 0 supplies the new row's boundary value. Induction proves the returned entry, modulo M.

The recurrence itself is correct because removing the newest element gives a unique predecessor and a recoverable insertion choice. For partitions, the two cases are joining a block and creating a singleton. For permutations, they are insertion into a cycle and creation of a new singleton cycle. Different insertion choices yield different counted objects, so the coefficients multiply counts without duplicating objects.

For the formula row, k! × S(n, k) equals the sum over i = 0…k of (−1)^(k − i) × binomial(k, i) × iⁿ. It counts onto functions to k labeled destinations by excluding missing destinations. For n = k = 0, the counting convention is 0⁰ = 1. Division by k! is exact over the integers, but cannot be replaced by ordinary division of a modular remainder.

Variants, Pitfalls, and Optimizations

Name the family before choosing the coefficient. Second-kind S(n, k) uses k, while unsigned first-kind c(n, k) uses n − 1. Signed first-kind numbers are s(n, k) = (−1)^(n − k) × c(n, k). If signed values are required modulo M, normalize a negative residue: when n − k is odd, use (M − c_mod) mod M. These families are also distinct from Stirling's approximation for factorials.

Distinguish labeled destinations from unlabeled blocks. S(n, k) counts unlabeled nonempty groups. Assigning names to those k groups gives k! × S(n, k) onto assignments to k labeled destinations. If destinations may be empty, the assignment count is kⁿ instead. Block size restrictions, capacities, or distinguished groups generally require a modified recurrence or additional DP state.

Preserve old-row values. A left-to-right update can read a newly computed dp[j − 1], corrupting the recurrence. Do not assume a small example will reveal the bug: for the second kind, column one stays equal to one, so errors can remain hidden when only two columns are stored. The example below reaches column three, where the corruption is visible.

Keep the boundary timing correct. The initial dp[0] = 1 must survive until row one has been computed, then become zero. Leaving it equal to one introduces nonexistent positive-size partitions or permutations with zero groups. Clearing it before row one's inner loop destroys the only initial source of nonzero counts.

Choose storage for the query pattern. A full triangle supports many queries at the cost of quadratic storage. If only one row is needed, compute through column n and sum second-kind entries to obtain the Bell number B(n), which counts all set partitions. Summing unsigned first-kind entries gives n!, because every permutation has one definite number of cycles. Unlike binomial coefficients, Stirling rows are not symmetric in k, so truncating to min(k, n − k) is invalid.

Choose arithmetic for the output. The recurrence uses no division and applies modulo prime or composite M when the base values are reduced appropriately. Formulas involving division by k! need additional care: for a prime modulus p and k < p, k! has a modular inverse; if k ≥ p, that simple inverse step fails. Never divide a residue as though it were a full integer count.

Practice Problems and Real-World Use

Start by stating exactly what identifies an object and exactly what identifies a group. Are people distinct? Do groups have names? Does order within a group matter? Is that order cyclic? Those questions decide whether the answer is a second-kind number, an unsigned first-kind number, a factorial multiple, or a different counting model.

For grouping applications, Stirling numbers count the possible configurations rather than evaluate their quality. For example, S(n, k) counts candidate partitions of n distinct records into k nonempty unlabeled clusters, but does not choose useful clusters or model similarity constraints. For permutation applications, the first-kind numbers classify all permutations by cycle count.

For the polynomial example, define the falling factorial x-underlined-k as x(x − 1)…(x − k + 1), with the degree-zero case equal to one. Then xⁿ is the sum of S(n, k) × x-underlined-k over k = 0…n. For n = 4, this reads x⁴ = x + 7x(x − 1) + 6x(x − 1)(x − 2) + x(x − 1)(x − 2)(x − 3). The unsigned first-kind numbers are the coefficients of the rising factorial x(x + 1)…(x + n − 1) in ordinary powers of x.

Verify small second-kind rows by generating partitions in a canonical representation, so reordering blocks does not create duplicates. Verify small first-kind rows by enumerating permutations and counting their cycles. Also check row sums against Bell numbers and factorials, the diagonal values, the zero column, impossible k, and the differing answers S(4, 2) = 7 versus c(4, 2) = 11. Use full counts or suitable exact arithmetic for probabilities; modular residues cannot be substituted for numerators of ordinary probabilities.