Detailed mathematical approach

Problem Summary

For integers \(n\), \(k\), and \(b\), define \(S(n,k,b)\) as the number of \(k\)-tuples \((x_1,\dots,x_k)\) satisfying

$0\le x_m\le b^m \quad (1\le m\le k),\qquad \sum_{m=1}^{k} x_m \le n.$

The overall task is to evaluate

$\sum_{k=10}^{15} S(10^k,k,k)\pmod{10^9+7}.$

The bounds \(b^m\) grow quickly, while \(n\) can be as large as \(10^{15}\), so a direct dynamic program on the sum would be far too large. The successful approach is to turn the inequality into an equality, apply inclusion-exclusion to the upper bounds, and exploit the fact that \(k\le 15\) is tiny.

Mathematical Approach

The key point is that the number of variables is small even though the target sum is huge. That makes subset enumeration feasible, while any method indexed directly by \(n\) would be wasteful.

Step 1: Convert the Inequality into an Equality

Introduce a slack variable \(s\ge 0\) and rewrite the condition as

$x_1+\cdots+x_k+s=n.$

Now we are counting nonnegative integer solutions in \(k+1\) variables, with upper bounds only on \(x_1,\dots,x_k\).

If those upper bounds were absent, stars and bars would give

$\#\{(x_1,\dots,x_k,s)\in \mathbb{Z}_{\ge 0}^{k+1}:x_1+\cdots+x_k+s=n\}=\binom{n+k}{k}.$

So the whole problem reduces to correcting this unrestricted count for the coordinates that exceed their allowed caps.

Step 2: Apply Inclusion-Exclusion to the Upper Bounds

For each index \(m\), let \(E_m\) be the event that the upper bound is violated:

$E_m=\{x_m\ge b^m+1\}.$

Fix a subset \(A\subseteq\{1,\dots,k\}\). If every event in \(A\) occurs, then for each \(m\in A\) we can offset

$x_m=y_m+(b^m+1),\qquad y_m\ge 0.$

After this change of variables, the equation becomes

$\sum_{m\notin A} x_m+\sum_{m\in A} y_m+s=n-\sum_{m\in A}(b^m+1).$

Therefore the number of solutions contributing to the intersection of all events in \(A\) is

$\binom{n-\sum_{m\in A}(b^m+1)+k}{k},$

with the usual convention that the term is \(0\) when the upper argument is negative or smaller than \(k\).

Inclusion-exclusion now yields the exact formula

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

This formula is already efficient enough conceptually because there are only \(2^k\) subsets.

Step 3: Evaluate the Binomial Coefficients Modulo \(10^9+7\)

Let \(p=10^9+7\). The program always uses \(k\le 15\), so in particular \(k<p\).

Write any upper argument as

$N=N_0+N_1p+N_2p^2+\cdots,\qquad 0\le N_0<p.$

Since \(k\) has only one nonzero base-\(p\) digit, Lucas' theorem gives

$\binom{N}{k}\equiv \binom{N_0}{k}\pmod{p}.$

So only the residue \(N\bmod p\) matters. If \(N_0<k\), then \(\binom{N}{k}\equiv 0\pmod p\). Otherwise we can compute

$\binom{N}{k}\equiv \frac{N_0(N_0-1)\cdots(N_0-k+1)}{k!}\pmod p.$

Because \(k\) is tiny, the numerator needs only \(k\) modular multiplications, and \(1/k!\) can be precomputed once using Fermat's little theorem.

Step 4: Precompute the Subset Shifts

For each index \(m\), define its violation cost

$c_m=b^m+1.$

For a subset \(A\), the inclusion-exclusion formula depends only on the total offset

$\sigma(A)=\sum_{m\in A} c_m.$

Rather than recomputing that sum from scratch for every subset, one can build all \(2^k\) values incrementally. If \(A\neq \varnothing\) and we remove one chosen element \(r\in A\), then

$\sigma(A)=\sigma(A\setminus\{r\})+c_r.$

This turns the subset enumeration into a simple table build followed by a signed accumulation. The mathematics stays the same; the precomputation only saves repeated work.

Step 5: Worked Example with \(S(14,3,2)\)

Here the bounds are

$x_1\le 2,\qquad x_2\le 4,\qquad x_3\le 8,\qquad x_1+x_2+x_3\le 14.$

After adding slack \(s\), the unrestricted count is

$\binom{14+3}{3}=\binom{17}{3}=680.$

The single-violation offsets are \(3\), \(5\), and \(9\), so the three single-event terms are

$\binom{14}{3}=364,\qquad \binom{12}{3}=220,\qquad \binom{8}{3}=56.$

The double intersections use offsets \(8\), \(12\), and \(14\), giving

$\binom{9}{3}=84,\qquad \binom{5}{3}=10,\qquad \binom{3}{3}=1.$

The triple intersection uses offset \(17\), which leaves a negative remainder, so that term is \(0\).

Hence

$S(14,3,2)=680-(364+220+56)+(84+10+1)=135.$

This is one of the checkpoint values matched by the implementation.

Step 6: Assemble the Final Euler Sum

Once \(S(n,k,b)\) can be evaluated quickly, the complete answer is just

$\sum_{k=10}^{15} S(10^k,k,k)\pmod{10^9+7}.$

There are only six queries, and each of them uses \(k\le 15\), so the subset-based method remains easily fast enough.

How the Code Works

The C++, Python, and Java implementations follow the same structure. First they precompute factorials and inverse factorials up to \(15\), which is enough for every binomial denominator that appears. Then, for each required value of \(k\), they build the list \(b^1+1,b^2+1,\dots,b^k+1\).

Next, they generate the total offset for every subset of \(\{1,\dots,k\}\) by reusing the previously computed total of a smaller subset. For each subset they check whether that offset already exceeds \(n\); if so, the corresponding inclusion-exclusion term contributes nothing. Otherwise they form the upper argument \(N=(n-\sigma(A))+k\), reduce \(N\) modulo \(10^9+7\), evaluate the binomial coefficient with a short falling product, and add or subtract the result according to the parity of \(|A|\).

After one query \(S(n,k,b)\) is finished, the program moves on to the next \(k\) and finally sums the six query results modulo \(10^9+7\).

Complexity Analysis

For one triple \((n,k,b)\), generating the violation costs takes \(O(k)\) time. Building all subset offsets takes \(O(2^k)\) time and \(O(2^k)\) memory. Evaluating the inclusion-exclusion sum needs \(2^k\) binomial terms, and each term costs \(O(k)\) modular operations because \(k\) is computed by a short falling product.

Therefore one query costs \(O(2^k\cdot k)\) time and \(O(2^k)\) memory. Since the full problem uses only \(k=10,\dots,15\), the total work is just a small constant multiple of that bound.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=528
2. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
3. Stars and bars: Wikipedia — Stars and bars
4. Lucas' theorem: Wikipedia — Lucas' theorem
5. Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse