Project Euler Problem 638

Solution

Answer: 18423394

Let the path consist of $a$ right-steps and $b$ up-steps.

If we encode a path by the heights at which right-steps occur, then the area under the path is exactly the sum of those heights.

The key observation is that this is a $q$-binomial coefficient (Gaussian binomial).

For a path $P_{a,b}$,

$$C(a,b,k) = \sum_{P_{a,b}} k^{A(P)}$$

is the generating function for lattice paths weighted by area. This is the classical identity:

$$C(a,b,k) = \binom{a+b}{a}k = \prod{i=1}^{a} \frac{k^{b+i}-1}{k^i-1}.$$

This follows because the area statistic on monotone lattice paths corresponds to inversion counting / Ferrers diagram area, whose generating function is the Gaussian binomial.

For $k=1$, the expression becomes the ordinary binomial coefficient:

$$C(a,b,1)=\binom{a+b}{a},$$

which matches the given example:

$$C(10,10,1)=\binom{20}{10}=184756.$$

We can verify another sample:

$$C(15,10,3) = \prod_{i=1}^{15} \frac{3^{10+i}-1}{3^i-1} \equiv 880419838 \pmod{10^9+7},$$

matching the statement.

Efficient algorithm

We need

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

Since

$$C(n,n,k) = \prod_{i=1}^{n} \frac{k^{n+i}-1}{k^i-1},$$

we evaluate the product iteratively modulo

$$M=1,000,000,007.$$

Because $M$ is prime, divisions are done using modular inverses:

$$x^{-1}\equiv x^{M-2}\pmod M.$$

Python implementation:

MOD = 1_000_000_007

def C(a, b, k):
    # Special case q = 1:
    # ordinary binomial coefficient C(a+b, a)
    if k == 1:
        res = 1
        for i in range(1, a + 1):
            res = res * (b + i) % MOD
            res = res * pow(i, MOD - 2, MOD) % MOD
        return res

    q = k % MOD
    q_to_b = pow(k, b, MOD)

    res = 1
    q_pow_i = 1

    # Product:
    # ∏_{i=1}^{a} (k^(b+i)-1)/(k^i-1)
    for i in range(1, a + 1):
        q_pow_i = q_pow_i * q % MOD

        numerator = (q_to_b * q_pow_i - 1) % MOD
        denominator = (q_pow_i - 1) % MOD

        res = (
            res
            * numerator
            % MOD
            * pow(denominator, MOD - 2, MOD)
            % MOD
        )

    return res

answer = 0

for k in range(1, 8):
    n = 10**k + k
    answer = (answer + C(n, n, k)) % MOD

print(answer)

Walkthrough of the code

C(a, b, k) computes the weighted lattice-path sum.
If k == 1, we compute the ordinary binomial coefficient

$\binom{a+b}{a}$.

Otherwise we use the Gaussian binomial product formula.
q_pow_i stores $k^i$ incrementally.
Each factor

$$\frac{k^{b+i}-1}{k^i-1}$$

is multiplied modulo $10^9+7$, using Fermat inverses.

Sanity checks:

C(2,2,1) = 6
C(2,2,2) = 35
C(10,10,1) = 184756
C(15,10,3) ≡ 880419838 (mod 10^9+7)
C(10000,10000,4) ≡ 395913804 (mod 10^9+7)

All match the problem statement.

The computed values for $k=1,\dots,7$ sum to:

$$18,423,394 \pmod{1,000,000,007}.$$

Answer: 18423394