Problem 743: Window into a Matrix

Project Euler Problem 743 Solution

EulerSolve provides an optimized solution for Project Euler Problem 743, Window into a Matrix, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary This problem asks for the quantity $A(k,n)$ arising from the matrix setup, under the condition $k\mid n$, with the final value taken modulo $M=10^9+7$. The concrete target is extremely large, so the essential task is to turn the definition into a summation that can be streamed term by term rather than expanded combinatorially from scratch. The C++, Python, and Java implementations all reveal the same structure: once we write $q=n/k$ and reduce $2^q$ modulo $M$, the whole computation becomes a weighted central-trinomial-type sum whose consecutive terms differ by a simple rational factor. Mathematical Approach Throughout the derivation let $$M=10^9+7,\qquad q=\frac{n}{k},\qquad b=2^q \pmod{M}.$$ The implementations do not build the answer from factorial tables. Instead they exploit a closed form for the summand and a recurrence between neighboring terms. Step 1: Rewrite the target as a weighted constant term Modulo $M$, a convenient reformulation is $$A(k,n)\equiv [z^0]\,(b+z+z^{-1})^k \pmod{M},$$ where $[z^0]$ means the coefficient of $z^0$. To contribute to the constant term, an expansion term must contain the same number of $z$ and $z^{-1}$ factors. If we choose exactly $x$ copies of $z$ and $x$ copies of $z^{-1}$, then the remaining $k-2x$ factors contribute $b$....

Detailed mathematical approach

Problem Summary

This problem asks for the quantity $A(k,n)$ arising from the matrix setup, under the condition $k\mid n$, with the final value taken modulo $M=10^9+7$. The concrete target is extremely large, so the essential task is to turn the definition into a summation that can be streamed term by term rather than expanded combinatorially from scratch.

The C++, Python, and Java implementations all reveal the same structure: once we write $q=n/k$ and reduce $2^q$ modulo $M$, the whole computation becomes a weighted central-trinomial-type sum whose consecutive terms differ by a simple rational factor.

Mathematical Approach

Throughout the derivation let

$$M=10^9+7,\qquad q=\frac{n}{k},\qquad b=2^q \pmod{M}.$$

The implementations do not build the answer from factorial tables. Instead they exploit a closed form for the summand and a recurrence between neighboring terms.

Step 1: Rewrite the target as a weighted constant term

Modulo $M$, a convenient reformulation is

$$A(k,n)\equiv [z^0]\,(b+z+z^{-1})^k \pmod{M},$$

where $[z^0]$ means the coefficient of $z^0$. To contribute to the constant term, an expansion term must contain the same number of $z$ and $z^{-1}$ factors. If we choose exactly $x$ copies of $z$ and $x$ copies of $z^{-1}$, then the remaining $k-2x$ factors contribute $b$.

This gives the explicit sum

$$A(k,n)=\sum_{x=0}^{\lfloor k/2\rfloor}\frac{k!}{x!^2(k-2x)!}\,b^{k-2x}\pmod{M}.$$

An equivalent binomial form is

$$A(k,n)=\sum_{x=0}^{\lfloor k/2\rfloor}\binom{k}{2x}\binom{2x}{x}b^{k-2x}\pmod{M}.$$

This is why the sequence behaves like a weighted version of the central trinomial coefficients.

Step 2: Derive the ratio between consecutive terms

Define

$$T_x=\frac{k!}{x!^2(k-2x)!}\,b^{k-2x}\pmod{M}.$$

Then $A(k,n)=\sum_x T_x$. The key simplification is that $T_x$ can be obtained from $T_{x-1}$ without recomputing any factorials:

$$\frac{T_x}{T_{x-1}}=\frac{(k-2x+2)(k-2x+1)}{x^2}\,b^{-2}.$$

The factorial part collapses because

$$\frac{(x-1)!^2(k-2x+2)!}{x!^2(k-2x)!}=\frac{(k-2x+2)(k-2x+1)}{x^2},$$

and the power of $b$ drops by two at every step. Therefore

$$T_x=T_{x-1}\cdot \frac{(k-2x+2)(k-2x+1)}{x^2}\cdot b^{-2}\pmod{M},$$

starting from

$$T_0=b^k.$$

Step 3: Replace divisions with modular inverses

Because $M$ is prime and $b=2^q$ is nonzero modulo $M$, every required division can be turned into multiplication by an inverse. Fermat's little theorem gives

$$a^{-1}\equiv a^{M-2}\pmod{M}$$

for any nonzero $a$. In particular, the implementations compute $b^{-1}$ as $b^{M-2}$ and then square it once to get the fixed multiplier $b^{-2}$.

They also need $x^{-1}$ for all $1\le x\le \lfloor k/2\rfloor$. Rather than calling fast exponentiation for each $x$, they precompute inverses in linear time using

$$\operatorname{inv}(1)=1,\qquad \operatorname{inv}(i)=M-\left\lfloor\frac{M}{i}\right\rfloor \operatorname{inv}(M\bmod i)\pmod{M}.$$

Then the factor $x^{-2}$ becomes two multiplications by $\operatorname{inv}(x)$.

Step 4: Stream the whole sum in one pass

A direct factorial-based method would require large factorial and inverse-factorial tables up to $k$, which is unnecessary here. The ratio formula reduces the update at step $x$ to five ingredients: the previous term, two descending numerator factors, two copies of $x^{-1}$, and the constant factor $b^{-2}$.

So the computation can march from $x=0$ to $x=\lfloor k/2\rfloor$, updating the current term and adding it to the running total immediately. Nothing from earlier terms needs to be revisited.

Step 5: Worked Example

For the checkpoint $A(3,9)$, we have

$$q=\frac{9}{3}=3,\qquad b=2^3=8.$$

Since $\lfloor 3/2\rfloor=1$, the sum has only two terms:

$$T_0=b^3=8^3=512,$$

$$T_1=\frac{3!}{1!^2\,1!}b^{1}=6\cdot 8=48.$$

Hence

$$A(3,9)=512+48=560,$$

matching the checkpoint used by the implementation. A second checkpoint is

$$A(4,20)=32^4+12\cdot 32^2+6=1{,}060{,}870.$$

How the Code Works

The C++, Python, and Java implementations all follow the same algorithm. First they compute $q=n/k$, then evaluate $b=2^q\pmod{M}$ by fast modular exponentiation, and finally initialize the first summand as $b^k\pmod{M}$. The running answer starts from that same initial term.

Next they compute the inverse of $b$, square it to obtain $b^{-2}$, and build an inverse table for the integers $1$ through $\lfloor k/2\rfloor$. During the main loop, the implementation updates the current summand by multiplying with the next two descending numerator factors, the inverse of the current index twice, and the fixed factor $b^{-2}$. The new summand is then added to the running total modulo $M$.

This approach avoids storing factorial arrays, inverse-factorial arrays, or any large polynomial object. The code keeps only the current summand, the running total, the inverse table, and a few loop variables.

Complexity Analysis

Let $h=\lfloor k/2\rfloor$. The modular exponentiations cost $O(\log q+\log M)$ time, which is tiny compared with the main loop. Filling the inverse table costs $O(h)$ time and $O(h)$ memory, and the streaming summation also costs $O(h)$ time. Therefore the full method runs in $O(k)$ time and uses $O(k)$ memory.

Footnotes and References

Problem page: https://projecteuler.net/problem=743
Central trinomial coefficient: Wikipedia - Central trinomial coefficient
Binomial coefficient: Wikipedia - Binomial coefficient
Modular multiplicative inverse: Wikipedia - Modular multiplicative inverse
Fermat's little theorem: Wikipedia - Fermat's little theorem

Problem 743 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u32 = std::uint32_t;

constexpr u64 kMod = 1'000'000'007ULL;

inline u64 mod_mul(const u64 a, const u64 b) {
    return (a * b) % kMod;
}

u64 mod_pow(u64 base, u64 exp) {
    u64 result = 1ULL;
    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) {
            result = mod_mul(result, base);
        }
        base = mod_mul(base, base);
        exp >>= 1ULL;
    }
    return result;
}

u64 A(const u64 k, const u64 n) {
    assert(n % k == 0ULL);

    const u64 q = n / k;
    const u64 b = mod_pow(2ULL, q);

    u64 term = mod_pow(b, k);
    u64 answer = term;

    const u64 inv_b = mod_pow(b, kMod - 2ULL);
    const u64 inv_b2 = mod_mul(inv_b, inv_b);

    const u64 half = k / 2ULL;
    std::vector<u32> inverses(static_cast<std::size_t>(half + 1ULL), 0U);
    if (half >= 1ULL) {
        inverses[1] = 1U;
    }
    for (u64 i = 2ULL; i <= half; ++i) {
        inverses[static_cast<std::size_t>(i)] = static_cast<u32>(
            kMod - mod_mul(kMod / i, inverses[static_cast<std::size_t>(kMod % i)]));
    }

    u64 num1 = k % kMod;
    u64 num2 = (k + kMod - 1ULL) % kMod;
    for (u64 x = 1ULL; x <= half; ++x) {
        const u64 inv_x = inverses[static_cast<std::size_t>(x)];

        term = mod_mul(term, num1);
        term = mod_mul(term, num2);
        term = mod_mul(term, inv_x);
        term = mod_mul(term, inv_x);
        term = mod_mul(term, inv_b2);

        answer += term;
        if (answer >= kMod) {
            answer -= kMod;
        }

        num1 += kMod - 2ULL;
        if (num1 >= kMod) {
            num1 -= kMod;
        }
        num2 += kMod - 2ULL;
        if (num2 >= kMod) {
            num2 -= kMod;
        }
    }

    return answer;
}

}  // namespace

int main() {
    assert(A(3ULL, 9ULL) == 560ULL);
    assert(A(4ULL, 20ULL) == 1'060'870ULL);

    std::cout << A(100'000'000ULL, 10'000'000'000'000'000ULL) << '\n';
    return 0;
}

Python

def solve():
    MOD = 1_000_000_007
    k = 100_000_000
    n = 10_000_000_000_000_000

    def mod_pow(base, exp, mod=MOD):
        r = 1
        base %= mod
        while exp > 0:
            if exp & 1: r = r * base % mod
            base = base * base % mod
            exp >>= 1
        return r

    def mod_mul(a, b): return a * b % MOD

    q = n // k
    b = mod_pow(2, q)
    term = mod_pow(b, k)
    answer = term

    inv_b = mod_pow(b, MOD - 2)
    inv_b2 = mod_mul(inv_b, inv_b)

    half = k // 2
    inverses = [0] * (half + 1)
    if half >= 1: inverses[1] = 1
    for i in range(2, half + 1):
        inverses[i] = MOD - (MOD // i) * inverses[MOD % i] % MOD

    num1 = k % MOD
    num2 = (k - 1) % MOD
    for x in range(1, half + 1):
        inv_x = inverses[x]
        term = mod_mul(term, num1)
        term = mod_mul(term, num2)
        term = mod_mul(term, inv_x)
        term = mod_mul(term, inv_x)
        term = mod_mul(term, inv_b2)
        answer = (answer + term) % MOD
        num1 = (num1 + MOD - 2) % MOD
        num2 = (num2 + MOD - 2) % MOD

    return str(answer)

if __name__ == '__main__':
    print(solve())

Java

public class Euler743 {
    static final long kMod = 1000000007L;

    static long modMul(long a, long b) {
        return (a * b) % kMod;
    }

    static long modPow(long base, long exp) {
        long result = 1L;
        while (exp > 0L) {
            if ((exp & 1L) != 0L) {
                result = modMul(result, base);
            }
            base = modMul(base, base);
            exp >>= 1L;
        }
        return result;
    }

    public static String solve() {
        long k = 100000000L;
        long n = 10000000000000000L;

        long q = n / k;
        long b = modPow(2L, q);

        long term = modPow(b, k);
        long answer = term;

        long invB = modPow(b, kMod - 2L);
        long invB2 = modMul(invB, invB);

        int half = (int) (k / 2L);
        int[] inverses = new int[half + 1];
        if (half >= 1) {
            inverses[1] = 1;
        }
        for (int i = 2; i <= half; ++i) {
            inverses[i] = (int) (kMod - modMul(kMod / i, inverses[(int) (kMod % i)]));
        }

        long num1 = k % kMod;
        long num2 = (k + kMod - 1L) % kMod;

        for (int x = 1; x <= half; ++x) {
            long invX = inverses[x];

            term = modMul(term, num1);
            term = modMul(term, num2);
            term = modMul(term, invX);
            term = modMul(term, invX);
            term = modMul(term, invB2);

            answer += term;
            if (answer >= kMod) {
                answer -= kMod;
            }

            num1 += kMod - 2L;
            if (num1 >= kMod)
                num1 -= kMod;

            num2 += kMod - 2L;
            if (num2 >= kMod)
                num2 -= kMod;
        }

        return Long.toString(answer);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}