Project Euler

Constrained Sums

Let S(n, k, b) denote the number of solutions to x_1 + x_2 +... + x_k <= n where 0 <= x_m <= b^m for all 1 <= m <= k. Given: S(14, 3, 2) = 135 S(200, 5, 3) = 12949440 S(1000, 10, 5) mod (10^9 + 7)...

Source sync Apr 19, 2026

Problem #0528

Level Level 27

Solved By 363

Languages C++, Python

Answer 779027989

Length 301 words

modular_arithmeticcombinatoricsbrute_force

Let $S(n, k, b)$ represent the number of valid solutions to $x_1 + x_2 + \cdots + x_k \le n$, where $0 \le x_m \le b^m$ for all $1 \le m \le k$.

For example, $S(14,3,2) = 135$, $S(200,5,3) = 12949440$, and $S(1000,10,5) \bmod 1\,000\,000\,007 = 624839075$.

Find $(\sum _{10 \le k \le 15} S(10^k, k, k)) \bmod 1\,000\,000\,007$.

Problem 528: Constrained Sums

Mathematical Foundation

Theorem (Inclusion-Exclusion for Bounded Compositions). Let $n, k \ge 1$ and $u_1, \ldots, u_k \ge 0$ be upper bounds. The number of solutions to $x_1 + \cdots + x_k \le n$ with $0 \le x_m \le u_m$ is

\sum_{T \subseteq \{1, \ldots, k\}} (-1)^{|T|} \binom{n - \sum_{m \in T}(u_m + 1) + k}{k}

where $\binom{a}{k} = 0$ when $a < 0$ .

Proof. Introduce a slack variable $x_0 \ge 0$ to convert the inequality to an equality:

x_0 + x_1 + \cdots + x_k = n, \quad x_0 \ge 0, \; 0 \le x_m \le u_m.

Without upper bounds, the number of non-negative integer solutions to $x_0 + x_1 + \cdots + x_k = n$ is $\binom{n+k}{k}$ by the stars-and-bars theorem.

For each $m \in \{1, \ldots, k\}$ , let $A_m$ denote the set of solutions where $x_m > u_m$ . Substituting $x_m' = x_m - (u_m + 1) \ge 0$ , the number of solutions in $A_m$ is $\binom{n - (u_m + 1) + k}{k}$ .

By inclusion-exclusion:

|A_1^c \cap \cdots \cap A_k^c| = \sum_{T \subseteq \{1, \ldots, k\}} (-1)^{|T|} \binom{n - \sum_{m \in T}(u_m + 1) + k}{k}.

$\square$

Lemma (Exponential Pruning). For $b = k$ and $u_m = k^m$ , the upper bounds grow exponentially. A subset $T \subseteq \{1, \ldots, k\}$ contributes a nonzero term only if

\sum_{m \in T} (k^m + 1) \le n = 10^k.

Since $k^k + 1 \le 10^k + 1$ and $\sum_{m \in T} k^m \ge k^{\max T}$ , the constraint forces $|T|$ to be small. Specifically, if $\max T = k$ , then $k^k \le 10^k$ , so $T$ can contain at most one large element. The total number of contributing subsets is at most $O(k \cdot 2^{k'})$ where $k'$ is the number of “small” indices, yielding a manageable enumeration.

Proof. For any $T$ with $\sum_{m \in T}(k^m + 1) > n$ , the binomial coefficient $\binom{n - \sum(k^m+1) + k}{k} = 0$ since the upper argument is negative. The exponential growth of $k^m$ means that including the index $m$ in $T$ “uses up” at least $k^m$ of the budget $n$ . Since $\sum_{m=1}^{k} k^m = k(k^k - 1)/(k-1)$ , which is comparable to $k^{k+1}/(k-1)$ , only subsets with controlled total can contribute. For practical $k \le 15$ , brute-force enumeration of all $2^k$ subsets with early termination suffices. $\square$

Lemma (Modular Binomial Coefficient for Large $n$ , Small $k$ ). For $n$ potentially as large as $10^{15} + 15$ and $k \le 15$ , the binomial coefficient is

\binom{n}{k} = \frac{n(n-1)\cdots(n-k+1)}{k!} \bmod p

computed as a product of $k$ terms modulo $p$ , multiplied by $(k!)^{-1} \bmod p$ .

Proof. Since $k \le 15$ and $p = 10^9 + 7 > k$ , we have $\gcd(k!, p) = 1$ , so the modular inverse of $k!$ exists. The product $n(n-1)\cdots(n-k+1)$ involves $k$ terms, each reduced modulo $p$ . $\square$

Editorial

S(n,k,b) = number of solutions to x1+x2+…+xk <= n with 0 <= xm <= b^m. Find (sum_{k=10}^{15} S(10^k, k, k)) mod 10^9+7. Uses inclusion-exclusion with modular arithmetic.

Pseudocode

upper bounds: u[m] = k^m for m = 1..k
Enumerate all subsets T of {1, ..., k}
if bit m of T is set
if valid

Complexity Analysis

Time: For each $k$ , we enumerate at most $2^k$ subsets (with pruning reducing the effective count). For each subset, the binomial coefficient computation is $O(k)$ . Total: $O\bigl(\sum_{k=10}^{15} 2^k \cdot k\bigr) = O(2^{15} \cdot 15) = O(491520)$ .
Space: $O(k)$ for the current subset and intermediate values.

Answer

\boxed{779027989}

C++ project_euler/problem_528/solution.cpp

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

/*
 * Project Euler Problem 528: Constrained Sums
 *
 * S(n,k,b) = number of solutions to x1+x2+...+xk <= n with 0 <= xm <= b^m.
 * Find (sum_{k=10}^{15} S(10^k, k, k)) mod 10^9+7.
 *
 * Uses inclusion-exclusion with modular arithmetic.
 */

typedef long long ll;
typedef __int128 lll;

const ll MOD = 1000000007;

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

ll mod_inv(ll a, ll mod) {
    return power(a, mod - 2, mod);
}

// Compute C(N, k) mod p where N can be very large but k is small
// C(N, k) = N*(N-1)*...*(N-k+1) / k!
ll binom_large_n(ll N, int k) {
    if (N < 0 || N < k) return 0;
    if (k == 0) return 1;

    // Compute numerator mod MOD
    ll num = 1;
    for (int i = 0; i < k; i++) {
        ll term = ((N - i) % MOD + MOD) % MOD;
        num = num * term % MOD;
    }

    // Compute k! mod MOD
    ll den = 1;
    for (int i = 1; i <= k; i++) {
        den = den * i % MOD;
    }

    return num % MOD * mod_inv(den, MOD) % MOD;
}

// Compute S(n, k, b) mod MOD using inclusion-exclusion
// The upper bounds are b^1, b^2, ..., b^k
// We need subsets T of {1,...,k} where sum of (b^m + 1) for m in T <= n
ll compute_S(ll n, int k, int b) {
    // Precompute b^m for m = 1..k (as regular integers, might overflow for large b,k)
    // For b=k<=15, b^k can be up to 15^15 ~ 4.37e17, fits in long long
    vector<ll> bpow(k + 1);
    bpow[0] = 1;
    for (int m = 1; m <= k; m++) {
        // Check overflow
        if (bpow[m-1] > 1e18 / b) {
            bpow[m] = (ll)2e18; // sentinel: too large
        } else {
            bpow[m] = bpow[m-1] * b;
        }
    }

    ll result = 0;

    // Enumerate subsets of {1,...,k} using bitmask
    for (int mask = 0; mask < (1 << k); mask++) {
        int bits = __builtin_popcount(mask);
        ll subtract = 0;
        bool overflow = false;

        for (int m = 0; m < k; m++) {
            if (mask & (1 << m)) {
                ll val = bpow[m + 1] + 1;
                if (subtract > 1e18 - val) {
                    overflow = true;
                    break;
                }
                subtract += val;
            }
        }

        if (overflow || subtract > n) continue;

        ll remaining = n - subtract;
        ll coeff = binom_large_n(remaining + k, k);

        if (bits % 2 == 0) {
            result = (result + coeff) % MOD;
        } else {
            result = (result - coeff + MOD) % MOD;
        }
    }

    return result;
}

int main() {
    // Verify small cases
    cout << "S(14,3,2) = " << compute_S(14, 3, 2) << endl;      // 135
    cout << "S(200,5,3) = " << compute_S(200, 5, 3) << endl;    // 12949440

    ll total = 0;
    for (int k = 10; k <= 15; k++) {
        // n = 10^k
        ll n = 1;
        for (int i = 0; i < k; i++) n *= 10;

        ll s = compute_S(n, k, k);
        cout << "S(10^" << k << "," << k << "," << k << ") mod 10^9+7 = " << s << endl;
        total = (total + s) % MOD;
    }

    cout << "Answer: " << total << endl;

    return 0;
}

Python project_euler/problem_528/solution.py

"""
Project Euler Problem 528: Constrained Sums

S(n,k,b) = number of solutions to x1+x2+...+xk <= n with 0 <= xm <= b^m.
Find (sum_{k=10}^{15} S(10^k, k, k)) mod 10^9+7.

Uses inclusion-exclusion with modular arithmetic.
"""

MOD = 10**9 + 7

def mod_inv(a, m=MOD):
    return pow(a, m - 2, m)

def binom_large_n(N, k):
    """Compute C(N, k) mod MOD where N can be very large but k is small."""
    if N < 0 or N < k:
        return 0
    if k == 0:
        return 1

    num = 1
    for i in range(k):
        num = num * ((N - i) % MOD) % MOD

    den = 1
    for i in range(1, k + 1):
        den = den * i % MOD

    return num * mod_inv(den) % MOD

def compute_S(n, k, b):
    """Compute S(n, k, b) mod MOD using inclusion-exclusion."""
    # Precompute b^m for m = 1..k
    bpow = [1] * (k + 1)
    for m in range(1, k + 1):
        bpow[m] = bpow[m-1] * b

    result = 0

    # Enumerate subsets of {1,...,k}
    for mask in range(1 << k):
        bits = bin(mask).count('1')
        subtract = 0
        valid = True

        for m in range(k):
            if mask & (1 << m):
                subtract += bpow[m + 1] + 1
                if subtract > n:
                    valid = False
                    break

        if not valid:
            continue

        remaining = n - subtract
        coeff = binom_large_n(remaining + k, k)

        if bits % 2 == 0:
            result = (result + coeff) % MOD
        else:
            result = (result - coeff + MOD) % MOD

    return result

# Verify small cases
print(f"S(14,3,2) = {compute_S(14, 3, 2)}")      # 135
print(f"S(200,5,3) = {compute_S(200, 5, 3)}")    # 12949440

total = 0
for k in range(10, 16):
    n = 10**k
    s = compute_S(n, k, k)
    print(f"S(10^{k},{k},{k}) mod 10^9+7 = {s}")
    total = (total + s) % MOD

print(f"Answer: {total}")