Combinations with Repetition and the Stars and Bars Method

In combinatorics and algorithm design, a classic family of counting problems asks: How many ways can we choose $k$ items from $n$ distinct types when repetition is allowed? Equivalently, how many integer solutions exist for a linear equation under various boundary constraints?

This post breaks down the Stars and Bars (also known as the dividers or balls and urns) method across three canonical scenarios:

1. Unbounded repetition ($x_i \ge 0$)
2. At least one of each type ($x_i \ge 1$)
3. Lower and upper bounds ($l_i \le x_i \le u_i$) via Inclusion-Exclusion and Dynamic Programming

1. Combinations with Repetition ($x_i \ge 0$)

Problem Statement

Suppose we have $n$ types of items and want to select $k$ items in total, where:

• Repetition is allowed (we can pick multiple items of the same type, up to $k$), and
• Order does not matter (only the count of each type matters).

The total number of ways to make this selection is:

Equivalence to Non-Negative Integer Solutions

Let $x_i$ denote the number of items chosen from type $i$ (for $i = 1, 2, \dots, n$). Because we choose $k$ items in total and may choose zero or more of each type, this problem is strictly equivalent to counting the number of non-negative integer solutions to the equation:

Why Does $\binom{n+k-1}{k}$ Work?

Imagine representing the $k$ chosen items as $k$ identical stars ($\star$) and using $n - 1$ bars ($|$) as dividers to separate the stars into $n$ bins (one bin per item type):

• Everything to the left of the first bar belongs to $x_1$.
• Everything between the $(i-1)$-th and $i$-th bar belongs to $x_i$.
• Everything to the right of the $(n-1)$-th bar belongs to $x_n$.
• Two adjacent bars ($||$) simply mean the bin between them has $0$ items.

For example, with $n = 4$ types and $k = 5$ items, the configuration:

corresponds to the solution $(x_1, x_2, x_3, x_4) = (2, 0, 3, 0)$, which sums to $5$.

In total, there are $k$ stars and $n - 1$ bars, giving $n + k - 1$ symbols in a row. Choosing the positions of the $k$ stars (or equivalently, the $n - 1$ bars) completely determines the solution:

2. Choosing at Least One of Each Type ($x_i \ge 1$)

Problem Statement

Now suppose we must select at least one item of each of the $n$ types (assuming $k \ge n$). This is equivalent to counting the positive integer solutions to:

The number of solutions is:

Two Ways to Derive This Result

Method A: The Gap / Divider Intuition

Place the $k$ identical items ($\star$) in a line:

There are $k - 1$ internal gaps between adjacent stars. Because no bin can be empty ($x_i \ge 1$), no two bars can occupy the same gap, nor can a bar be placed outside the outer stars. Thus, we simply choose $n - 1$ distinct gaps out of the $k - 1$ available gaps to insert our dividers:

Method B: Variable Substitution

Pre-allocate $1$ item to each of the $n$ variables by defining $y_i = x_i - 1$. Since $x_i \ge 1$, we have $y_i \ge 0$, and the equation transforms into:

Applying the unbounded formula from Section 1 with $k' = k - n$ gives:

3. Bounded Variables ($0 \le x_i \le u_i$)

In many practical problems, each variable is constrained by both a lower bound and an upper bound:

First, any lower bound $x_i \ge l_i$ can be normalized to zero via the substitution $x_i' = x_i - l_i$, which reduces the target sum to $k' = k - \sum_{i=1}^n l_i$ and the upper bounds to $u_i' = u_i - l_i$. Thus, we only need to focus on:

There are two standard techniques to handle upper bounds: Inclusion-Exclusion and Dynamic Programming.

Approach 1: Principle of Inclusion-Exclusion (PIE)

To count solutions satisfying $0 \le x_i \le u_i$, we start with all non-negative solutions ($x_i \ge 0$) and subtract solutions that violate at least one upper bound.

Let $P_i$ be the property that variable $x_i$ violates its upper bound, i.e.,

For any subset of indices $S \subseteq \{1, 2, \dots, n\}$, how many non-negative integer solutions violate all upper bounds in $S$ simultaneously (meaning $x_i \ge u_i + 1$ for all $i \in S$)?

By pre-subtracting $u_i + 1$ from each $x_i$ ($i \in S$), the remaining sum is:

If $k_S \ge 0$, the number of such solutions is $\binom{n + k_S - 1}{n - 1}$ (and $0$ if $k_S < 0$).

By the Principle of Inclusion-Exclusion, the exact number of valid bounded solutions is:

where binomial coefficients with a negative upper index (when $k_S < 0$) are treated as $0$.

Special Case: Uniform Upper Bound ($u_i = u$ for all $i$)

When every variable shares the same upper bound $0 \le x_i \le u$, all subsets $S$ of size $j = |S|$ subtract the same amount $j(u + 1)$. The PIE sum simplifies to an $O(n)$ loop:

Approach 2: Dynamic Programming (DP)

When the upper bounds $u_i$ are distinct and $n$ is too large for an $O(2^n)$ subset enumeration (e.g., $n = 100, k = 1000$), Dynamic Programming is the method of choice.

State Definition

Let $dp[i][s]$ be the number of ways to choose values for the first $i$ variables ($x_1, \dots, x_i$) such that their sum equals $s$.

Transition

For the $i$-th variable, we can pick any value $v \in [0, \min(u_i, s)]$:

Notice that the right-hand side is a contiguous range sum over $dp[i-1][\dots]$ from index $s - \min(u_i, s)$ to $s$. Using prefix sums, each state transition takes $O(1)$ time, reducing the overall time complexity from $O(n k^2)$ to $O(n k)$ and space complexity to $O(k)$.

Python Implementation ($O(nk)$ Time, $O(k)$ Space)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

def count_bounded_solutions(u: list[int], k: int, mod: int = 10**9 + 7) -> int:
    """
    Counts integer solutions to x_1 + x_2 + ... + x_n = k
    subject to 0 <= x_i <= u[i - 1], using prefix-sum optimized DP.
    """
    dp = [0] * (k + 1)
    dp[0] = 1

    for bound in u:
        # Compute prefix sums of the previous DP row
        prefix = [0] * (k + 2)
        for s in range(k + 1):
            prefix[s + 1] = (prefix[s] + dp[s]) % mod

        # Transition in O(1) per sum s:
        # new_dp[s] = sum(dp[s - bound : s + 1])
        new_dp = [0] * (k + 1)
        for s in range(k + 1):
            left = max(0, s - bound)
            new_dp[s] = (prefix[s + 1] - prefix[left]) % mod
        dp = new_dp

    return dp[k]

Summary Cheat Sheet