Project Euler Problem 408

Let's call a lattice point (x, y) inadmissible if x, y and x+y are all positive perfect squares.

Project Euler Problem 408

Solution

Answer: 299742733

Let

$$T(n)=\binom{2n}{n}$$

be the total number of monotone lattice paths from $(0,0)$ to $(n,n)$.

We must exclude paths that pass through any inadmissible point.


1. Mathematical analysis

A lattice point $(x,y)$ is inadmissible iff

  • $x$ is a positive square,
  • $y$ is a positive square,
  • $x+y$ is also a positive square.

So there exist integers $a,b,c>0$ such that

$$x=a^2,\qquad y=b^2,\qquad a^2+b^2=c^2.$$

Thus inadmissible points correspond exactly to Pythagorean triples.

For $n=10^7$,

$$a,b \le \sqrt{10^7}=3162.$$

Hence the number of inadmissible points is quite small (only 7850).


Counting paths through lattice points

The number of monotone north/east paths from $(x_1,y_1)$ to $(x_2,y_2)$ is

$$\binom{(x_2-x_1)+(y_2-y_1)}{x_2-x_1}.$$

If we sort all inadmissible points by $(x+y,x)$, then every admissible predecessor of a point appears earlier.

Define:

$$F(P)=\text{number of admissible paths from }(0,0)\text{ to }P.$$

For a point $P=(x,y)$,

$$F(P) = \binom{x+y}{x} - \sum_{Q\le P} F(Q)\binom{(x-x_Q)+(y-y_Q)}{x-x_Q},$$

where the sum runs over earlier inadmissible points $Q=(x_Q,y_Q)$ satisfying

$$x_Q\le x,\qquad y_Q\le y.$$

This is standard inclusion–exclusion / DP on a DAG.

Finally,

$$P(n) = \binom{2n}{n} - \sum_Q F(Q) \binom{(n-x_Q)+(n-y_Q)}{n-x_Q}.$$

All computations are modulo

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

Since $M$ is prime and $2n< M$, ordinary factorial/binomial arithmetic modulo $M$ works directly.


2. Python implementation

from math import isqrt
from array import array

MOD = 1_000_000_007
N = 10_000_000

# ------------------------------------------------------------
# Precompute factorials and inverse factorials modulo MOD
# ------------------------------------------------------------

MAXF = 2 * N

fact = array('I', [0]) * (MAXF + 1)
invfact = array('I', [0]) * (MAXF + 1)

fact[0] = 1
for i in range(1, MAXF + 1):
    fact[i] = (fact[i - 1] * i) % MOD

invfact[MAXF] = pow(fact[MAXF], MOD - 2, MOD)

for i in range(MAXF, 0, -1):
    invfact[i - 1] = (invfact[i] * i) % MOD

def C(n, r):
    """Binomial coefficient modulo MOD."""
    if r < 0 or r > n:
        return 0
    return (fact[n] * invfact[r] * invfact[n - r]) % MOD

# ------------------------------------------------------------
# Generate inadmissible points
# ------------------------------------------------------------

limit = isqrt(N)

points = set()

for a in range(1, limit + 1):
    aa = a * a

    for b in range(a + 1, limit + 1):
        bb = b * b

        c2 = aa + bb
        c = isqrt(c2)

        if c * c == c2:
            points.add((aa, bb))
            points.add((bb, aa))

points = sorted(points, key=lambda p: (p[0] + p[1], p[0]))

# ------------------------------------------------------------
# DP over inadmissible points
# ------------------------------------------------------------

k = len(points)

xs = [p[0] for p in points]
ys = [p[1] for p in points]

F = [0] * k

for i in range(k):

    x, y = xs[i], ys[i]

    # Total paths to (x,y)
    val = C(x + y, x)

    # Remove paths already passing through earlier inadmissible points
    for j in range(i):

        xj, yj = xs[j], ys[j]

        if xj <= x and yj <= y:

            ways = C((x - xj) + (y - yj), x - xj)

            val -= F[j] * ways

    F[i] = val % MOD

# ------------------------------------------------------------
# Count admissible paths to (N,N)
# ------------------------------------------------------------

answer = C(2 * N, N)

for j in range(k):

    xj, yj = xs[j], ys[j]

    ways = C((N - xj) + (N - yj), N - xj)

    answer -= F[j] * ways

answer %= MOD

print(answer)

3. Code walkthrough

Factorials

We precompute

$$n! \pmod M$$

and inverse factorials using Fermat:

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

Then

$$\binom{n}{r} = \frac{n!}{r!(n-r)!} \pmod M.$$


Generating inadmissible points

We enumerate all

$$a,b\le \sqrt{N}$$

and test whether

$$a^2+b^2$$

is a square.

Whenever it is, we insert

$$(a^2,b^2)$$

and the symmetric point.

For $N=10^7$, there are only 7850 such points.


DP

For each inadmissible point $P$:

  1. Start with all paths to $P$,
  2. subtract those already forced through earlier inadmissible points.

Thus $F(P)$ counts paths whose first inadmissible point is exactly $P$.

Finally we subtract all paths from $(0,0)$ to $(N,N)$ that pass through any inadmissible point.


4. Verification on the given examples

The program reproduces:

$$P(5)=252$$

$$P(16)=596994440$$

$$P(1000)\bmod 1,000,000,007 = 341920854$$

matching the statement exactly.

The final computation for $N=10^7$ yields:

$$299742733.$$


Answer: 299742733