Project Euler Problem 528

Solution

Answer: 779027989

Let

$$S(n,k,b)=#\left{(x_1,\dots,x_k)\in \mathbb Z_{\ge 0}^k: x_1+\cdots+x_k\le n,;;0\le x_m\le b^m \right}.$$

We need

$$\left(\sum_{10\le k\le 15} S(10^k,k,k)\right)\bmod 1{,}000{,}000{,}007.$$

1. Mathematical analysis

We want to count bounded integer solutions to

$$x_1+\cdots+x_k\le n$$

with

$$0\le x_m\le b^m.$$

A standard trick converts the inequality into an equality.

Introduce a slack variable

$$x_{k+1}=n-(x_1+\cdots+x_k)\ge 0.$$

Then counting solutions to

$$x_1+\cdots+x_k\le n$$

is equivalent to counting nonnegative solutions of

$$x_1+\cdots+x_k+x_{k+1}=n$$

with upper bounds only on $x_1,\dots,x_k$.

Step 1: Unrestricted count

Without upper bounds, the number of nonnegative solutions to

$$x_1+\cdots+x_{k+1}=n$$

is the classic stars-and-bars result:

$$\binom{n+k}{k}.$$

Step 2: Inclusion–exclusion

We must enforce

$$x_m\le b^m.$$

Let

$$u_m=b^m.$$

For a subset $A\subseteq {1,\dots,k}$, suppose every variable in $A$ violates its bound:

$$x_i\ge u_i+1 \qquad (i\in A).$$

Set

$$y_i=x_i-(u_i+1)$$

for $i\in A$. Then the equation becomes

$$\sum y_i + \sum_{i\notin A}x_i + x_{k+1} = n-\sum_{i\in A}(u_i+1).$$

Hence the number of such solutions is

$$\binom{ n-\sum_{i\in A}(u_i+1)+k }{k},$$

interpreted as $0$ if the top value is $<k$.

By inclusion–exclusion:

$$S(n,k,b) = \sum_{A\subseteq {1,\dots,k}} (-1)^{|A|} \binom{ n-\sum_{i\in A}(b^i+1)+k }{k}.$$

Since $k\le 15$, only $2^{15}=32768$ subsets are needed, making this extremely fast.

2. Python implementation

from math import comb

MOD = 1_000_000_007

def S(n, k, b):
    # Upper bounds u_i = b^i
    bounds = [b ** i for i in range(1, k + 1)]

    total = 0

    # Inclusion-exclusion over subsets
    for mask in range(1 << k):
        shift = 0
        bits = 0

        for i in range(k):
            if (mask >> i) & 1:
                shift += bounds[i] + 1
                bits += 1

        top = n - shift + k

        # Binomial coefficient contributes only if top >= k
        ways = comb(top, k) if top >= k else 0

        if bits % 2:
            total -= ways
        else:
            total += ways

    return total

# Verify examples
print(S(14, 3, 2))                  # 135
print(S(200, 5, 3))                 # 12949440
print(S(1000, 10, 5) % MOD)         # 624839075

# Compute Euler 528 answer
answer = sum(S(10**k, k, k) for k in range(10, 16)) % MOD

print(answer)

3. Code walkthrough

`bounds = [b ** i for i in range(1, k + 1)]`

Creates:

$$[b^1,b^2,\dots,b^k].$$

These are the upper limits for each variable.

`for mask in range(1 << k)`

Iterates through all subsets of variables using a bitmask.

Each subset represents the variables that violate their upper bounds.

`shift += bounds[i] + 1`

If variable $x_i$ is forced to exceed its limit, we subtract off:

$$b^i+1$$

because we substitute

$$x_i = y_i + (b^i+1).$$

`top = n - shift + k`

After shifting, the remaining equation has

$$k+1$$

nonnegative variables summing to a smaller total.

Stars-and-bars gives:

$$\binom{\text{top}}{k}.$$

Inclusion–exclusion sign

Odd-sized subsets subtract; even-sized subsets add:

if bits % 2:
    total -= ways
else:
    total += ways

4. Verification against examples

The code reproduces all given values:

$$S(14,3,2)=135$$

$$S(200,5,3)=12949440$$

$$S(1000,10,5)\bmod 1{,}000{,}000{,}007 = 624839075$$

So the formula and implementation are consistent with the problem statement.

Computing

$$\sum_{k=10}^{15} S(10^k,k,k) \pmod{1{,}000{,}000{,}007}$$

gives:

Answer: 779027989