Problem 408: Admissible Paths Through a Grid

Project Euler Problem 408 Solution

EulerSolve provides an optimized solution for Project Euler Problem 408, Admissible Paths Through a Grid, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We count monotone lattice paths from $(0,0)$ to $(n,n)$, where each step is either one unit right or one unit up. A path is admissible if it never visits any forbidden grid point. In the implementations, the forbidden points are exactly the interior square-coordinate points $$F_n=\{(a^2,b^2): 1 \le a,b \le \lfloor \sqrt{n}\rfloor,\ \exists c \in \mathbb{Z}_{>0}: a^2+b^2=c^2\}.$$ The final answer is the number of admissible paths modulo $10^9+7$. Mathematical Approach Step 1: Count paths between two comparable points If $(x_1,y_1)$ and $(x_2,y_2)$ satisfy $x_1 \le x_2$ and $y_1 \le y_2$, then any monotone path from the first point to the second point must use exactly $x_2-x_1$ right moves and $y_2-y_1$ up moves. Hence the number of such paths is $$\binom{(x_2-x_1)+(y_2-y_1)}{x_2-x_1}.$$ In particular, without any forbidden points, the total number of paths from $(0,0)$ to $(n,n)$ is $$\binom{2n}{n}.$$ Step 2: Describe the forbidden set Only points with both coordinates at most $n$ can matter, so it is enough to enumerate $a,b \le \lfloor \sqrt{n}\rfloor$. The condition $a^2+b^2=c^2$ means that $(a,b,c)$ forms a Pythagorean triple, and each such pair produces a forbidden point $(a^2,b^2)$....

Detailed mathematical approach

Problem Summary

We count monotone lattice paths from $(0,0)$ to $(n,n)$, where each step is either one unit right or one unit up. A path is admissible if it never visits any forbidden grid point.

In the implementations, the forbidden points are exactly the interior square-coordinate points

$$F_n=\{(a^2,b^2): 1 \le a,b \le \lfloor \sqrt{n}\rfloor,\ \exists c \in \mathbb{Z}_{>0}: a^2+b^2=c^2\}.$$

The final answer is the number of admissible paths modulo $10^9+7$.

Mathematical Approach

Step 1: Count paths between two comparable points

If $(x_1,y_1)$ and $(x_2,y_2)$ satisfy $x_1 \le x_2$ and $y_1 \le y_2$, then any monotone path from the first point to the second point must use exactly $x_2-x_1$ right moves and $y_2-y_1$ up moves. Hence the number of such paths is

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

In particular, without any forbidden points, the total number of paths from $(0,0)$ to $(n,n)$ is

$$\binom{2n}{n}.$$

Step 2: Describe the forbidden set

Only points with both coordinates at most $n$ can matter, so it is enough to enumerate $a,b \le \lfloor \sqrt{n}\rfloor$. The condition $a^2+b^2=c^2$ means that $(a,b,c)$ forms a Pythagorean triple, and each such pair produces a forbidden point $(a^2,b^2)$.

After sorting the forbidden points lexicographically, write them as

$$P_1=(x_1,y_1),\ P_2=(x_2,y_2),\ \dots,\ P_m=(x_m,y_m).$$

This ordering is compatible with monotone movement: if a path can go from $P_j$ to $P_i$, then necessarily $x_j \le x_i$ and $y_j \le y_i$, which implies $j \lt i$ after lexicographic sorting.

Step 3: Inclusion-exclusion by the first forbidden point

For each forbidden point $P_i$, let $B_i$ be the number of monotone paths from $(0,0)$ to $P_i$ whose first forbidden point is exactly $P_i$. Start from the unrestricted count

$$\binom{x_i+y_i}{x_i}.$$

If a path reaches an earlier forbidden point $P_j$ first and later continues to $P_i$, then the number of such paths is

$$B_j \binom{(x_i-x_j)+(y_i-y_j)}{x_i-x_j},$$

provided $x_j \le x_i$ and $y_j \le y_i$. Therefore

$$B_i=\binom{x_i+y_i}{x_i}-\sum_{\substack{j \lt i \\ x_j \le x_i \\ y_j \le y_i}} B_j \binom{(x_i-x_j)+(y_i-y_j)}{x_i-x_j} \pmod{10^9+7}.$$

This dynamic program is a topological inclusion-exclusion on the partial order induced by coordinate-wise comparison.

Every inadmissible path from $(0,0)$ to $(n,n)$ has a unique first forbidden point, so after all $B_i$ are known we subtract the bad paths by that first hit:

$$A(n)=\binom{2n}{n}-\sum_{i=1}^{m} B_i \binom{(n-x_i)+(n-y_i)}{n-x_i} \pmod{10^9+7}.$$

This is exactly the recurrence used in the implementations.

Step 4: Worked examples from the checkpoints

For $n=5$, there are no forbidden points because the only possible square coordinates are $1$ and $4$, and none of the sums $1+1$, $1+4$, $4+1$, $4+4$ is a square. Hence

$$A(5)=\binom{10}{5}=252.$$

For $n=16$, the forbidden points are $(9,16)$ and $(16,9)$, coming from $3^2+4^2=5^2$ and $4^2+3^2=5^2$. The unrestricted count is

$$\binom{32}{16}=601080390.$$

Each forbidden point receives

$$\binom{9+16}{9}=\binom{25}{9}=2042975$$

paths from the origin. Neither forbidden point dominates the other, so no path can pass through both. Therefore

$$A(16)=\binom{32}{16}-2\binom{25}{9}=596994440,$$

which matches the checkpoint value used by the implementations.

Step 5: Fast modular binomial coefficients

All binomial coefficients are evaluated modulo the prime

$$M=10^9+7.$$

The implementations precompute factorials and inverse factorials up to $2n$, so every combination query becomes

$$\binom{N}{R}\equiv N! \cdot (R!)^{-1} \cdot ((N-R)!)^{-1} \pmod{M}.$$

The inverse factorials are derived using Fermat's little theorem:

$$a^{-1}\equiv a^{M-2}\pmod{M}, \qquad a \not\equiv 0 \pmod{M}.$$

After this preprocessing, each path-count query is $O(1)$.

How the Code Works

The C++, Python, and Java implementations follow the same structure. They first build a fast square-membership structure up to $2n$, enumerate every forbidden point $(a^2,b^2)$ inside the $n \times n$ grid, sort the resulting list, and remove duplicates. Next they precompute factorials and inverse factorials so that any binomial coefficient needed by the path formulas can be answered in constant time modulo $10^9+7$.

Once those tables are ready, the implementation scans the forbidden points in sorted order and computes the number of paths whose first forbidden point is each obstacle. The final admissible count is obtained by subtracting all bad-path contributions from the unrestricted total $\binom{2n}{n}$.

Complexity Analysis

Let $m$ be the number of forbidden points. Enumerating candidate square pairs takes $O(n)$ time because both $a$ and $b$ range up to $\lfloor \sqrt{n}\rfloor$. Precomputing factorials and inverse factorials up to $2n$ also takes $O(n)$ time and $O(n)$ memory. The dynamic program compares each forbidden point with all earlier ones, so it costs $O(m^2)$ time and $O(m)$ additional memory. Overall, the algorithm runs in $O(n+m^2)$ time and uses $O(n+m)$ memory.

Footnotes and References

Problem page: https://projecteuler.net/problem=408
Lattice path counting: Wikipedia — Lattice path
Binomial coefficients: Wikipedia — Binomial coefficient
Fermat's little theorem: Wikipedia — Fermat's little theorem
Pythagorean triples: Wikipedia — Pythagorean triple

Problem 408 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <utility>
#include <vector>
#include <cmath>
#include <functional>

namespace {

using i64 = long long;
using u64 = std::uint64_t;
constexpr int MOD = 1000000007;

struct Options {
    int n = 10000000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        value = std::stoi(tail);
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1;
}

int mod_pow(int base, i64 exp) {
    i64 result = 1;
    i64 cur = base % MOD;
    i64 e = exp;
    while (e > 0) {
        if (e & 1LL) {
            result = (result * cur) % MOD;
        }
        cur = (cur * cur) % MOD;
        e >>= 1LL;
    }
    return static_cast<int>(result);
}

std::vector<std::pair<int, int>> inadmissible_points(int n) {
    const int lim = static_cast<int>(std::sqrt(static_cast<long double>(n)));
    const int lim2 = static_cast<int>(std::sqrt(static_cast<long double>(2LL * n)));

    std::vector<std::uint8_t> is_square(static_cast<std::size_t>(2 * n + 1), 0U);
    for (int c = 1; c <= lim2; ++c) {
        is_square[static_cast<std::size_t>(c * c)] = 1U;
    }

    std::vector<std::pair<int, int>> pts;
    pts.reserve(8000);
    for (int a = 1; a <= lim; ++a) {
        const int a2 = a * a;
        for (int b = 1; b <= lim; ++b) {
            const int b2 = b * b;
            if (is_square[static_cast<std::size_t>(a2 + b2)] != 0U) {
                pts.emplace_back(a2, b2);
            }
        }
    }

    std::sort(pts.begin(), pts.end());
    pts.erase(std::unique(pts.begin(), pts.end()), pts.end());
    return pts;
}

int P_mod(const int n) {
    const int maxv = 2 * n;
    std::vector<int> fact(static_cast<std::size_t>(maxv + 1), 1);
    for (int i = 1; i <= maxv; ++i) {
        fact[static_cast<std::size_t>(i)] =
            static_cast<int>((1LL * fact[static_cast<std::size_t>(i - 1)] * i) % MOD);
    }

    std::vector<int> inv_fact(static_cast<std::size_t>(maxv + 1), 1);
    inv_fact[static_cast<std::size_t>(maxv)] = mod_pow(fact[static_cast<std::size_t>(maxv)], MOD - 2);
    for (int i = maxv; i >= 1; --i) {
        inv_fact[static_cast<std::size_t>(i - 1)] =
            static_cast<int>((1LL * inv_fact[static_cast<std::size_t>(i)] * i) % MOD);
    }

    auto nCr = [&](int nn, int rr) -> int {
        if (rr < 0 || rr > nn) {
            return 0;
        }
        return static_cast<int>(
            1LL * fact[static_cast<std::size_t>(nn)] * inv_fact[static_cast<std::size_t>(rr)] % MOD *
            inv_fact[static_cast<std::size_t>(nn - rr)] % MOD);
    };

    std::vector<std::pair<int, int>> pts = inadmissible_points(n);
    const int m = static_cast<int>(pts.size());
    std::vector<int> ways(static_cast<std::size_t>(m), 0);

    for (int i = 0; i < m; ++i) {
        const int xi = pts[static_cast<std::size_t>(i)].first;
        const int yi = pts[static_cast<std::size_t>(i)].second;
        i64 w = nCr(xi + yi, xi);

        for (int j = 0; j < i; ++j) {
            const int xj = pts[static_cast<std::size_t>(j)].first;
            const int yj = pts[static_cast<std::size_t>(j)].second;
            if (xj <= xi && yj <= yi) {
                const int paths = nCr((xi - xj) + (yi - yj), xi - xj);
                w -= 1LL * ways[static_cast<std::size_t>(j)] * paths % MOD;
                if (w < 0) {
                    w += MOD;
                }
            }
        }

        ways[static_cast<std::size_t>(i)] = static_cast<int>(w % MOD);
    }

    i64 ans = nCr(2 * n, n);
    for (int i = 0; i < m; ++i) {
        const int x = pts[static_cast<std::size_t>(i)].first;
        const int y = pts[static_cast<std::size_t>(i)].second;
        const int tail_paths = nCr((n - x) + (n - y), n - x);
        ans -= 1LL * ways[static_cast<std::size_t>(i)] * tail_paths % MOD;
        if (ans < 0) {
            ans += MOD;
        }
    }

    return static_cast<int>(ans % MOD);
}

bool run_checkpoints() {
    if (P_mod(5) != 252) {
        std::cerr << "Checkpoint failed: P(5)\n";
        return false;
    }
    if (P_mod(16) != 596994440) {
        std::cerr << "Checkpoint failed: P(16)\n";
        return false;
    }
    if (P_mod(1000) != 341920854) {
        std::cerr << "Checkpoint failed: P(1000)\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << P_mod(options.n) << '\n';
    return 0;
}

Python

import math

def solve():
    N = 10_000_000
    MOD = 1_000_000_007

    def mod_pow(base, exp, mod):
        result = 1
        cur = base % mod
        e = exp
        while e > 0:
            if e & 1: result = result * cur % mod
            cur = cur * cur % mod
            e >>= 1
        return result

    maxv = 2 * N
    fact = [1] * (maxv + 1)
    for i in range(1, maxv + 1):
        fact[i] = fact[i-1] * i % MOD
    inv_fact = [1] * (maxv + 1)
    inv_fact[maxv] = mod_pow(fact[maxv], MOD - 2, MOD)
    for i in range(maxv, 0, -1):
        inv_fact[i-1] = inv_fact[i] * i % MOD

    def nCr(n, r):
        if r < 0 or r > n: return 0
        return fact[n] * inv_fact[r] % MOD * inv_fact[n-r] % MOD

    # Find inadmissible points (a², b²) where a²+b² is a perfect square
    lim = int(math.isqrt(N))
    lim2 = int(math.isqrt(2 * N))
    is_sq = set(c*c for c in range(1, lim2+1))

    pts = set()
    for a in range(1, lim+1):
        a2 = a * a
        for b in range(1, lim+1):
            b2 = b * b
            if a2 + b2 in is_sq:
                pts.add((a2, b2))
    pts = sorted(pts)

    m = len(pts)
    ways = [0] * m
    for i in range(m):
        xi, yi = pts[i]
        w = nCr(xi + yi, xi)
        for j in range(i):
            xj, yj = pts[j]
            if xj <= xi and yj <= yi:
                paths = nCr((xi-xj)+(yi-yj), xi-xj)
                w = (w - ways[j] * paths) % MOD
        ways[i] = w % MOD

    ans = nCr(2*N, N)
    for i in range(m):
        x, y = pts[i]
        tail = nCr((N-x)+(N-y), N-x)
        ans = (ans - ways[i] * tail) % MOD

    return str(ans % MOD)

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

Java

import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

public class Euler408 {
    private static final int MOD = 1000000007;

    private static int modPow(int base, long exp) {
        long result = 1;
        long cur = base % MOD;
        long e = exp;
        while (e > 0) {
            if ((e & 1L) != 0) {
                result = (result * cur) % MOD;
            }
            cur = (cur * cur) % MOD;
            e >>= 1L;
        }
        return (int) result;
    }

    private static class Point implements Comparable<Point> {
        int x, y;

        Point(int x, int y) {
            this.x = x;
            this.y = y;
        }

        @Override
        public int compareTo(Point o) {
            if (this.x != o.x)
                return Integer.compare(this.x, o.x);
            return Integer.compare(this.y, o.y);
        }

        @Override
        public boolean equals(Object obj) {
            if (!(obj instanceof Point))
                return false;
            Point o = (Point) obj;
            return this.x == o.x && this.y == o.y;
        }
    }

    public static String solve() {
        int n = 10000000;
        int lim = (int) Math.sqrt(n);
        int lim2 = (int) Math.sqrt(2L * n);

        boolean[] isSquare = new boolean[2 * n + 1];
        for (int c = 1; c <= lim2; ++c) {
            isSquare[c * c] = true;
        }

        List<Point> pts = new ArrayList<>();
        for (int a = 1; a <= lim; ++a) {
            int a2 = a * a;
            for (int b = 1; b <= lim; ++b) {
                int b2 = b * b;
                if (isSquare[a2 + b2]) {
                    pts.add(new Point(a2, b2));
                }
            }
        }

        Collections.sort(pts);
        List<Point> uniquePts = new ArrayList<>();
        for (Point p : pts) {
            if (uniquePts.isEmpty() || !uniquePts.get(uniquePts.size() - 1).equals(p)) {
                uniquePts.add(p);
            }
        }
        pts = uniquePts;

        int maxv = 2 * n;
        int[] fact = new int[maxv + 1];
        fact[0] = 1;
        for (int i = 1; i <= maxv; ++i) {
            fact[i] = (int) ((1L * fact[i - 1] * i) % MOD);
        }

        int[] invFact = new int[maxv + 1];
        invFact[maxv] = modPow(fact[maxv], MOD - 2);
        for (int i = maxv; i >= 1; --i) {
            invFact[i - 1] = (int) ((1L * invFact[i] * i) % MOD);
        }

        int m = pts.size();
        int[] ways = new int[m];

        for (int i = 0; i < m; ++i) {
            int xi = pts.get(i).x;
            int yi = pts.get(i).y;
            long w = nCr(xi + yi, xi, fact, invFact);

            for (int j = 0; j < i; ++j) {
                int xj = pts.get(j).x;
                int yj = pts.get(j).y;
                if (xj <= xi && yj <= yi) {
                    long paths = nCr((xi - xj) + (yi - yj), xi - xj, fact, invFact);
                    w = (w - 1L * ways[j] * paths) % MOD;
                    if (w < 0)
                        w += MOD;
                }
            }
            ways[i] = (int) w;
        }

        long ans = nCr(2 * n, n, fact, invFact);
        for (int i = 0; i < m; ++i) {
            int x = pts.get(i).x;
            int y = pts.get(i).y;
            long tailPaths = nCr((n - x) + (n - y), n - x, fact, invFact);
            ans = (ans - 1L * ways[i] * tailPaths) % MOD;
            if (ans < 0)
                ans += MOD;
        }

        return String.valueOf(ans);
    }

    private static long nCr(int nn, int rr, int[] fact, int[] invFact) {
        if (rr < 0 || rr > nn)
            return 0;
        return 1L * fact[nn] * invFact[rr] % MOD * invFact[nn - rr] % MOD;
    }

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