Project Euler

Stirling Number Asymptotics

The Bell number B_n = sum_(k=0)^n S(n,k) counts the number of partitions of {1,..., n}. Find B_1000 mod (10^9+7).

Source sync Apr 19, 2026

Problem #0984

Level Level 39

Solved By 89

Languages C++, Python

Answer 465231251

Length 63 words

modular_arithmeticgeometrylinear_algebra

In Western chess, a knight is a piece that moves either two squares horizontally and one square vertically, or one square horizontally and two squares vertically, and is capable of jumping over any intervening pieces.

Chinese chess has a similar piece called the horse, whose moves have an identical displacement as a knight's move; however, a horse, unlike a knight, is unable to jump over intervening pieces.

More specifically, a horse's move consists of two steps: An orthogonal move of one square, followed by a diagonal move by one square in the same direction as the orthogonal move. If the orthogonal square is occupied by another piece, the horse is unable to move in that direction. Specifically the horse in the centre of the above board can move to the squares providing the squares are unoccupied. For example, if was occupied then it could not move to or .

A set of squares on a chessboard is called knight-connected if a knight can travel between any two squares in the set using only legal moves without using any squares not in the set. A set of squares on a chessboard is called horse-disjoint if, when a horse is placed on every square in the set (and no other square), no horse can attack any other. Let be the number of knight-connected, horse-disjoint non-empty subsets of an chessboard. For example, , consisting of the nine singleton sets. You are also given that , and .

Find . Give your answer modulo .

Problem 984: Stirling Number Asymptotics

Mathematical Analysis

Bell Numbers

Theorem (Bell Triangle). $B_n$ can be computed via the Bell triangle:

$a_{0,0} = 1$ .
$a_{n,0} = a_{n-1,n-1}$ .
$a_{n,k} = a_{n,k-1} + a_{n-1,k-1}$ .
$B_n = a_{n,0}$ .

Theorem (Dobinski’s Formula). $B_n = \frac{1}{e} \sum_{k=0}^{\infty} \frac{k^n}{k!}$ .

Growth Rate

Theorem. $\ln B_n \sim n \ln n - n \ln \ln n - n$ as $n \to \infty$ .

For $n = 1000$ : $B_{1000}$ has approximately 2000 digits.

Derivation

Editorial

Use the Bell triangle with two-row rolling array, reducing modulo $10^9 + 7$ at each step.

Complexity Analysis

$O(N^2)$ time, $O(N)$ space.

Answer

\boxed{465231251}

C++ project_euler/problem_984/solution.cpp

#include <bits/stdc++.h>
using namespace std;
const long long MOD = 1e9+7;
int main() {
    int n = 1000;
    vector<long long> row(n+1, 0);
    row[0] = 1;
    for (int i = 1; i <= n; i++) {
        vector<long long> nr(n+1, 0);
        nr[0] = row[i-1];
        for (int j = 1; j <= i; j++)
            nr[j] = (nr[j-1] + row[j-1]) % MOD;
        row = nr;
    }
    cout << row[0] << endl;
    return 0;
}

Python project_euler/problem_984/solution.py

"""
Problem 984: Bell Number B_1000 mod 10^9+7

Compute the 1000th Bell number modulo 10^9 + 7 using the Bell triangle.

Bell numbers B_n count the number of partitions of a set of n elements.
    B_0 = 1, B_1 = 1, B_2 = 2, B_3 = 5, B_4 = 15, B_5 = 52, ...

The Bell triangle recurrence:
    a(0,0) = 1
    a(n,0) = a(n-1, n-1)        (left edge = previous row's last element)
    a(n,k) = a(n,k-1) + a(n-1,k-1)
    B_n = a(n,0)

Key properties:
    - B_n grows super-exponentially: B_n ~ (n / W(n))^n / e^n
    - The Bell triangle modular computation uses O(n) space with rolling rows
    - Dobinski's formula: B_n = (1/e) * sum_{k=0}^{inf} k^n / k!

Answer: computed via Bell triangle mod 10^9+7

Methods:
    - bell_triangle_mod(n, mod): Bell number via triangle with modular arithmetic
    - bell_exact_small(n): Exact Bell numbers for small n (verification)
    - stirling2(n, k): Stirling numbers of second kind (B_n = sum S(n,k))
    - dobinski_approx(n): Approximate B_n via Dobinski's formula
"""

from math import factorial, comb

MOD = 10**9 + 7


def bell_triangle_mod(n, mod):
    """Compute B_n mod m using the Bell triangle with rolling array."""
    if n == 0:
        return 1
    row = [0] * (n + 1)
    row[0] = 1
    for i in range(1, n + 1):
        new_row = [0] * (n + 1)
        new_row[0] = row[i - 1]
        for j in range(1, i + 1):
            new_row[j] = (new_row[j - 1] + row[j - 1]) % mod
        row = new_row
    return row[0]


def bell_exact_small(n_max):
    """Compute exact Bell numbers B_0..B_{n_max} for small n."""
    B = [0] * (n_max + 1)
    B[0] = 1
    # Using Stirling: B_n = sum_{k=0}^{n} S(n,k)
    S = [[0] * (n_max + 1) for _ in range(n_max + 1)]
    S[0][0] = 1
    for nn in range(1, n_max + 1):
        for kk in range(1, nn + 1):
            S[nn][kk] = kk * S[nn - 1][kk] + S[nn - 1][kk - 1]
        B[nn] = sum(S[nn])
    return B


def stirling2_row(n):
    """Compute S(n, k) for k=0..n."""
    if n == 0:
        return [1]
    S = [0] * (n + 1)
    S[0] = 0
    prev = stirling2_row(n - 1) + [0]
    for k in range(1, n + 1):
        S[k] = k * prev[k] + prev[k - 1]
    return S


def dobinski_approx(n, terms=200):
    """Approximate B_n using Dobinski's formula (for small n only)."""
    from math import exp
    total = 0.0
    for k in range(terms):
        total += k ** n / factorial(k) if k > 0 or n == 0 else 1.0 / factorial(0)
    return total / exp(1)

# Verification

# Known Bell numbers
known_bells = [1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975]
exact_bells = bell_exact_small(10)

for i, b in enumerate(known_bells):
    assert exact_bells[i] == b, f"B_{i}: expected {b}, got {exact_bells[i]}"
    assert bell_triangle_mod(i, MOD) == b % MOD, f"Mod mismatch at B_{i}"

# Verify Dobinski for small values
for i in range(8):
    approx = dobinski_approx(i)
    assert abs(approx - known_bells[i]) < 0.5, f"Dobinski approx failed at B_{i}"

answer = bell_triangle_mod(1000, MOD)
print(answer)