Project Euler

Beans in Bowls

Count the number of solutions to x_1 + x_2 +... + x_k = n with 0 <= x_i <= m for all i. Denote this B(n, k, m). Compute the answer modulo a prime p.

Source sync Apr 19, 2026

Problem #0839

Level Level 25

Solved By 410

Languages C++, Python

Answer 150893234438294408

Length 320 words

modular_arithmeticdynamic_programminganalytic_math

The sequence $S_n$ is defined by $S_0 = 290797$ and $S_n = S_{n - 1}^2 \bmod 50515093$ for $n < 0$.

There are $N$ bowls indexed $0,1,\dots ,N-1$. Initially there are $S_n$ beans in bowl $n$.

At each step, the smallest index $n$ is found such that bowl $n$ has strictly more beans than bowl $n+1$. Then one bean is moved from bowl $n$ to bowl $n+1$.

Let $B(N)$ be the number of steps needed to sort the bowls into non-descending order.

For example, $B(5) = 0$, $B(6) = 14263289$ and $B(100)=3284417556$.

Find $B(10^7)$.

Problem 839: Beans in Bowls

Mathematical Foundation

Theorem 1 (Stars and Bars). The number of solutions to $x_1 + \cdots + x_k = n$ with $x_i \ge 0$ (no upper bound) is $\binom{n+k-1}{k-1}$ .

Proof. Represent the $n$ units as stars and the $k-1$ dividers between groups as bars. Each arrangement of $n$ stars and $k-1$ bars in a row of $n+k-1$ positions corresponds bijectively to a solution $(x_1, \ldots, x_k)$ , where $x_i$ is the number of stars between the $(i-1)$ -th and $i$ -th bar. The number of such arrangements is $\binom{n+k-1}{k-1}$ . $\square$

Theorem 2 (Bounded Compositions via Inclusion-Exclusion). The number of solutions to $x_1 + \cdots + x_k = n$ with $0 \le x_i \le m$ for all $i$ is

B(n, k, m) = \sum_{j=0}^{\lfloor n/(m+1) \rfloor} (-1)^j \binom{k}{j} \binom{n - j(m+1) + k - 1}{k - 1}.

Proof. Let $A_i = \{(x_1, \ldots, x_k) : x_1 + \cdots + x_k = n,\; x_j \ge 0 \;\forall j,\; x_i \ge m+1\}$ for $i = 1, \ldots, k$ . We seek $|\overline{A_1 \cup \cdots \cup A_k}|$ , the number of solutions with no variable exceeding $m$ . By inclusion-exclusion:

B(n,k,m) = \sum_{j=0}^{k} (-1)^j \sum_{|T|=j} |A_{i_1} \cap \cdots \cap A_{i_j}|.

For a set $T$ of $j$ indices, $\bigcap_{i \in T} A_i$ requires $x_i \ge m+1$ for all $i \in T$ . Substituting $y_i = x_i - (m+1)$ for $i \in T$ (so $y_i \ge 0$ ) reduces the sum to $n - j(m+1)$ . By Theorem 1, the count is $\binom{n - j(m+1) + k - 1}{k-1}$ when $n - j(m+1) \ge 0$ , and 0 otherwise. Since there are $\binom{k}{j}$ choices of $T$ , the formula follows. The sum terminates at $j = \lfloor n/(m+1) \rfloor$ because the binomial coefficient vanishes when the top argument is negative. $\square$

Lemma 1 (Generating Function). The generating function for a single bounded variable is $g(x) = \frac{1 - x^{m+1}}{1 - x}$ . The answer is the coefficient $[x^n]\, g(x)^k$ .

Proof. We have $g(x) = 1 + x + x^2 + \cdots + x^m$ (geometric sum). The coefficient of $x^n$ in $g(x)^k$ counts the number of ways to choose $x_1, \ldots, x_k \in \{0, 1, \ldots, m\}$ summing to $n$ . Expanding via the binomial theorem:

g(x)^k = (1 - x^{m+1})^k (1-x)^{-k} = \left(\sum_{j=0}^{k} (-1)^j \binom{k}{j} x^{j(m+1)}\right) \left(\sum_{r=0}^{\infty} \binom{r+k-1}{k-1} x^r\right).

Extracting $[x^n]$ by convolving yields the inclusion-exclusion formula of Theorem 2. $\square$

Lemma 2 (DP Recurrence). Define $f(i, s) =$ number of ways to assign values to bowls $1, \ldots, i$ totaling $s$ . Then:

f(i, s) = \sum_{x=0}^{\min(m, s)} f(i-1, s-x), \quad f(0, 0) = 1, \quad f(0, s) = 0 \text{ for } s > 0.

With prefix-sum optimization, each row $f(i, \cdot)$ is computed from $f(i-1, \cdot)$ in $O(n)$ time.

Proof. The recurrence follows by conditioning on the value of $x_i$ . For the prefix-sum optimization, note that $f(i, s) = \sum_{x=0}^{\min(m,s)} f(i-1, s-x) = F(s) - F(s - m - 1)$ where $F(s) = \sum_{t=0}^{s} f(i-1, t)$ is the prefix sum. Computing $F$ takes $O(n)$ , and each $f(i,s)$ is then $O(1)$ . $\square$

Editorial

Count distributions of n beans into k bowls, each with capacity m. B(n,k,m) = sum_{j=0}^{floor(n/(m+1))} (-1)^j * C(k,j) * C(n-j(m+1)+k-1, k-1). We method: Inclusion-exclusion with modular arithmetic. Finally, precompute factorials and inverse factorials mod p.

Pseudocode

Method: Inclusion-exclusion with modular arithmetic
Precompute factorials and inverse factorials mod p

Complexity Analysis

Time: $O(n + k)$ for factorial precomputation; $O(\min(k, n/(m+1)))$ for the inclusion-exclusion sum. Total: $O(n + k)$ .
Space: $O(n + k)$ for factorial tables.

Answer

\boxed{150893234438294408}

C++ project_euler/problem_839/solution.cpp

#include <bits/stdc++.h>
using namespace std;

/*
 * Problem 839: Beans in Bowls
 *
 * Count solutions to x_1 + ... + x_k = n with 0 <= x_i <= m.
 *
 * B(n,k,m) = sum_{j=0}^{floor(n/(m+1))} (-1)^j C(k,j) C(n-j(m+1)+k-1, k-1)
 *
 * Two methods: inclusion-exclusion and DP.
 */

const long long MOD = 1e9 + 7;
const int MAXFACT = 2000001;

long long fact[MAXFACT], inv_fact[MAXFACT];

long long power(long long base, long long exp, long long mod) {
    long long result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

void precompute_factorials(int n) {
    fact[0] = 1;
    for (int i = 1; i <= n; i++)
        fact[i] = fact[i - 1] * i % MOD;
    inv_fact[n] = power(fact[n], MOD - 2, MOD);
    for (int i = n - 1; i >= 0; i--)
        inv_fact[i] = inv_fact[i + 1] * (i + 1) % MOD;
}

long long comb(int n, int r) {
    if (r < 0 || r > n) return 0;
    return fact[n] % MOD * inv_fact[r] % MOD * inv_fact[n - r] % MOD;
}

// Method 1: Inclusion-exclusion
long long solve_ie(int n, int k, int m) {
    long long result = 0;
    for (int j = 0; j <= k; j++) {
        long long rem = n - (long long)j * (m + 1);
        if (rem < 0) break;
        long long term = comb(k, j) * comb(rem + k - 1, k - 1) % MOD;
        if (j % 2 == 0)
            result = (result + term) % MOD;
        else
            result = (result - term + MOD) % MOD;
    }
    return result;
}

// Method 2: DP with prefix sums
long long solve_dp(int n, int k, int m) {
    vector<long long> dp(n + 1, 0);
    dp[0] = 1;
    for (int bowl = 0; bowl < k; bowl++) {
        vector<long long> ndp(n + 1, 0);
        vector<long long> prefix(n + 2, 0);
        for (int s = 0; s <= n; s++)
            prefix[s + 1] = (prefix[s] + dp[s]) % MOD;
        for (int s = 0; s <= n; s++) {
            int lo = max(0, s - m);
            ndp[s] = (prefix[s + 1] - prefix[lo] + MOD) % MOD;
        }
        dp = ndp;
    }
    return dp[n];
}

int main() {
    precompute_factorials(MAXFACT - 1);

    // Verify on small cases
    assert(solve_ie(3, 2, 2) == 2);
    assert(solve_ie(4, 3, 2) == 6);

    // Cross-check methods on small inputs
    for (int n = 0; n <= 20; n++)
        for (int k = 1; k <= 5; k++)
            for (int m = 1; m <= 5; m++)
                assert(solve_ie(n, k, m) == solve_dp(n, k, m));

    // Solve main problem
    long long ans = solve_ie(1000, 100, 20);
    cout << ans << endl;
    return 0;
}

Python project_euler/problem_839/solution.py

"""
Problem 839: Beans in Bowls

Count distributions of n beans into k bowls, each with capacity m.
B(n,k,m) = sum_{j=0}^{floor(n/(m+1))} (-1)^j * C(k,j) * C(n-j(m+1)+k-1, k-1)
"""

from math import comb
from functools import lru_cache

MOD = 10**9 + 7

# --- Method 1: Inclusion-Exclusion ---
def beans_inclusion_exclusion(n: int, k: int, m: int):
    """Inclusion-exclusion formula for bounded compositions."""
    result = 0
    for j in range(k + 1):
        rem = n - j * (m + 1)
        if rem < 0:
            break
        term = comb(k, j) * comb(rem + k - 1, k - 1)
        if j % 2 == 0:
            result += term
        else:
            result -= term
    return result

# --- Method 2: Dynamic Programming ---
def beans_dp(n: int, k: int, m: int):
    """DP approach with prefix-sum optimization."""
    # dp[s] = number of ways to fill first i bowls to sum s
    dp = [0] * (n + 1)
    dp[0] = 1
    for _ in range(k):
        new_dp = [0] * (n + 1)
        prefix = [0] * (n + 2)
        for s in range(n + 1):
            prefix[s + 1] = prefix[s] + dp[s]
        for s in range(n + 1):
            lo = max(0, s - m)
            new_dp[s] = prefix[s + 1] - prefix[lo]
        dp = new_dp
    return dp[n]

# --- Method 3: Generating function (polynomial multiplication) ---
def beans_genfunc(n: int, k: int, m: int):
    """Multiply the polynomial (1 + x + ... + x^m) k times."""
    poly = [1] * (m + 1)  # single bowl GF
    result = [1]  # identity polynomial
    for _ in range(k):
        new_result = [0] * (len(result) + len(poly) - 1)
        for i, a in enumerate(result):
            for j, b in enumerate(poly):
                if i + j <= n:
                    new_result[i + j] += a * b
        result = new_result[:n + 1]
    return result[n] if n < len(result) else 0

# --- Verify all methods agree ---
for nn in range(0, 20):
    for kk in range(1, 6):
        for mm in range(1, 8):
            a = beans_inclusion_exclusion(nn, kk, mm)
            b = beans_dp(nn, kk, mm)
            c = beans_genfunc(nn, kk, mm)
            assert a == b == c, f"Mismatch at ({nn},{kk},{mm}): IE={a}, DP={b}, GF={c}"

print("All verification passed!")

# Modular version for large inputs
def modinv(a, m=MOD):
    return pow(a, m - 2, m)

def comb_mod(n, r, mod=MOD):
    if r < 0 or r > n:
        return 0
    if r == 0:
        return 1
    num = 1
    den = 1
    for i in range(r):
        num = num * ((n - i) % mod) % mod
        den = den * ((i + 1) % mod) % mod
    return num * modinv(den, mod) % mod

def beans_mod(n, k, m, mod=MOD):
    result = 0
    for j in range(k + 1):
        rem = n - j * (m + 1)
        if rem < 0:
            break
        term = comb_mod(k, j, mod) * comb_mod(rem + k - 1, k - 1, mod) % mod
        if j % 2 == 0:
            result = (result + term) % mod
        else:
            result = (result - term + mod) % mod
    return result

# Example computation
answer = beans_mod(1000, 100, 20)
print(answer)