Detailed mathematical approach

Problem Summary

The quantity to compute is \(B(k,n)\) modulo \(M=10^9+7\). The C++, Python, and Java implementations do not enumerate matrix states directly. Instead, they rewrite the target value as a coefficient extraction problem built from factorial-based series. The difficult part is a family of sums involving sixteenth powers of binomial coefficients, and the implementations evaluate all of those sums at once by fast polynomial convolution.

For the target instance, \(k\) is as large as \(100000\), so a direct quadratic summation would be far too slow. The entire strategy is therefore organized around turning the problem into two convolutions of length about \(k\), followed by a final modular reconstruction.

Mathematical Approach

Write

$q=\left\lfloor\frac{n}{k}\right\rfloor,\qquad M=10^9+7,\qquad \mu\equiv 2^q-2 \pmod{M}.$

The implemented formula can be derived cleanly in the following steps.

Step 1: Reduce the exponential parameter

All dependence on \(n\) enters through the single scalar \(\mu\). Because \(M\) is prime and \(2\not\equiv 0\pmod{M}\), Fermat's little theorem gives

$2^q \equiv 2^{\,q\bmod (M-1)} \pmod{M}.$

So even when \(q\) is enormous, the implementation only needs the reduced exponent \(q\bmod(M-1)\). After this reduction, the rest of the computation depends only on \(k\) and on the already-computed value of \(\mu\).

Step 2: Convert the inner binomial-power sum into a convolution

Define the sequence

$c_i=(i!)^{-16}\pmod{M}\qquad (0\le i\le k).$

Now form its self-convolution:

$d_s=\sum_{i=0}^{s} c_i\,c_{s-i} \qquad (0\le s\le k).$

Multiplying by \((s!)^{16}\) transforms this into

$a_s=(s!)^{16}d_s=(s!)^{16}\sum_{i=0}^{s}\frac{1}{(i!)^{16}((s-i)!)^{16}}=\sum_{i=0}^{s}\left(\frac{s!}{i!(s-i)!}\right)^{16}.$

Therefore

$\boxed{a_s=\sum_{i=0}^{s}\binom{s}{i}^{16}}.$

This identity is the key simplification: instead of evaluating each \(a_s\) separately, one self-convolution computes every value \(a_0,a_1,\dots,a_k\) in one shot.

Step 3: Rewrite \(B(k,n)\) as one more coefficient extraction

After the quantities \(a_s\) are known, define

$g_s=\frac{a_s}{s!},\qquad e_j=\frac{\mu^j}{j!}.$

Convolving these two sequences gives

$f_t=\sum_{s=0}^{t} g_s\,e_{t-s}=\sum_{s=0}^{t}\frac{a_s}{s!}\frac{\mu^{\,t-s}}{(t-s)!}.$

The final answer is the \(t=k\) case multiplied by \(k!\):

$B(k,n)=k!\,f_k.$

Expanding the factorials yields the explicit summation formula

$\boxed{B(k,n)=\sum_{s=0}^{k}\binom{k}{s}a_s\,\mu^{\,k-s}\pmod{M}.}$

Equivalently, if we package the coefficients into an exponential generating function, then

$\frac{B(k,n)}{k!}=\left[x^k\right]\left(\sum_{s\ge 0}a_s\frac{x^s}{s!}\right)e^{\mu x}.$

So the whole problem becomes a coefficient lookup after a second convolution.

Step 4: Use NTT and CRT for both convolutions

A naive evaluation of all \(a_s\) from the formula \(\sum_{i=0}^{s}\binom{s}{i}^{16}\) would cost

$\sum_{s=0}^{k} O(s)=O(k^2),$

which is too slow when \(k=100000\). Both required convolutions therefore use the Number Theoretic Transform. The target modulus \(M=10^9+7\) is not suitable for an NTT of the needed lengths, so the implementation performs each convolution under three NTT-friendly primes:

$998244353,\qquad 1004535809,\qquad 469762049.$

Each coefficient is first computed modulo these three primes and then reconstructed by the Chinese Remainder Theorem. The product of the three auxiliary moduli is far larger than the raw coefficient sizes that occur here, so the reconstruction is exact before the final reduction modulo \(M\).

Step 5: Worked example

The small checkpoint \(B(2,4)\) illustrates every layer of the formula. Here

$q=\left\lfloor\frac{4}{2}\right\rfloor=2,\qquad \mu=2^2-2=2.$

Next compute the first three values of \(a_s\):

$a_0=\binom{0}{0}^{16}=1,$

$a_1=\binom{1}{0}^{16}+\binom{1}{1}^{16}=1+1=2,$

$a_2=\binom{2}{0}^{16}+\binom{2}{1}^{16}+\binom{2}{2}^{16}=1+2^{16}+1=65538.$

Substituting into the closed form gives

$\begin{aligned} B(2,4) &=\binom{2}{0}a_0\mu^2+\binom{2}{1}a_1\mu+\binom{2}{2}a_2\\ &=1\cdot 1\cdot 4+2\cdot 2\cdot 2+1\cdot 65538\\ &=4+8+65538\\ &=65550. \end{aligned}$

This matches the small-value check used by the implementation.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. First they precompute factorials and inverse factorials modulo \(M\) up to \(k\). That makes it cheap to evaluate \((i!)^{-16}\), \((s!)^{16}\), and the factorial denominators that appear in the two generating series.

Next they compute \(\mu=2^{\lfloor n/k\rfloor}-2\pmod{M}\) with fast modular exponentiation, reducing the exponent modulo \(M-1\) before the power is taken. Then they build the sequence \(c_i=(i!)^{-16}\), convolve it with itself, and rescale the result to recover every value \(a_s=\sum_{i=0}^{s}\binom{s}{i}^{16}\).

After that, the implementation forms one series from \(a_s/s!\) and a second series from \(\mu^j/j!\), performs a second convolution, and multiplies the coefficient of degree \(k\) by \(k!\). The C++ version additionally parallelizes the three auxiliary-modulus convolutions and the CRT recombination, while the Python and Java versions execute the same mathematical steps sequentially.

Complexity Analysis

Precomputing factorials, inverse factorials, and the powers needed for the exponential series costs \(O(k)\) time and \(O(k)\) memory. Each NTT-based convolution costs \(O(k\log k)\) time. There are two logical convolutions, each evaluated under a fixed set of three auxiliary primes, so the asymptotic running time remains \(O(k\log k)\) and the memory usage remains \(O(k)\). Parallel execution in the C++ implementation improves wall-clock time but does not change the asymptotic order.

Footnotes and References

1. Problem page: Project Euler 767
2. Binomial coefficients: Wikipedia - Binomial coefficient
3. Generating functions: Wikipedia - Generating function
4. Number Theoretic Transform: cp-algorithms - Fast Fourier transform and NTT
5. Chinese Remainder Theorem: Wikipedia - Chinese remainder theorem