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

A lattice point $(x,y)$ is called inadmissible if $x$, $y$ and $x+y$ are all positive perfect squares. Thus $(9,16)$ is inadmissible, because $9=3^2$, $16=4^2$ and $25=5^2$, while $(0,4)$, $(3,1)$ and $(9,4)$ are not. A path from $(x_1,y_1)$ to $(x_2,y_2)$ uses only unit steps north or east, and it is called admissible if none of its intermediate points is inadmissible. Writing $P(n)$ for the number of admissible paths from $(0,0)$ to $(n,n)$, the statement records $P(5)=252$, $P(16)=596994440$ and $P(1000)\bmod(10^9+7)=341920854$, and asks for $P(10^7)$ modulo $10^9+7$. ^[1]

The whole solution in one sentence: count all $\binom{2n}{n}$ paths, then remove the bad ones by listing the few inadmissible points inside the square and counting, for each of them in turn, the paths that meet it before meeting any other.

Mathematical Approach

1. Counting paths when nothing is forbidden

Let $P=(x,y)$ and $Q=(u,v)$ with $x\le u$ and $y\le v$. A path from $P$ to $Q$ consists of $u-x$ east steps and $v-y$ north steps in some order, and it is determined by which of the $(u-x)+(v-y)$ steps are the east ones. Hence the number of paths is the binomial coefficient ^[2]

$$N(P,Q)=\binom{(u-x)+(v-y)}{u-x}.$$

If $x\gt u$ or $y\gt v$ there is no path at all, since neither step type ever decreases a coordinate. Writing $O=(0,0)$ and $D=(n,n)$, the unrestricted count is $N(O,D)=\binom{2n}{n}$; for $n=3$ that is $\binom{6}{3}=20$.

2. Which points are inadmissible

An inadmissible point is exactly a point $(a^2,b^2)$ with $a,b\ge1$ and $a^2+b^2=c^2$ for some integer $c$, so the inadmissible points correspond one to one with the ordered Pythagorean triples $(a,b,c)$. ^[3] A triple is primitive when its three entries have greatest common divisor $1$; dividing any triple by that divisor leaves a primitive one. Euclid's parametrisation $p=e^2-f^2$, $q=2ef$, $r=e^2+f^2$, with $e\gt f\ge1$ coprime and of opposite parity, yields every primitive triple with its odd leg first, ^[3] and exchanging the legs and taking positive multiples then yields every ordered triple.

A path from $O$ to $D$ stays inside $[0,n]^2$, so only points with $1\le a,b\le\lfloor\sqrt n\rfloor$ can be met. The smallest triple $(3,4,5)$ gives $(9,16)$ and $(16,9)$, and no inadmissible point has both coordinates below $16$; a brute-force count confirms $P(n)=\binom{2n}{n}$ for $n\le15$, so $P(5)=\binom{10}{5}=252$, the first checkpoint. Inside the square of side $400$ there are fourteen:

$$\begin{aligned}&(9,16),\ (16,9),\ (25,144),\\&(36,64),\ (64,36),\ (64,225),\\&(81,144),\ (144,25),\ (144,81),\\&(144,256),\ (225,64),\ (225,400),\\&(256,144),\ (400,225).\end{aligned}$$

3. How many inadmissible points the square contains

The count matters, because the algorithm below is quadratic in it. Fix a leg $a$: since $a^2+b^2=c^2$ is the same as $a^2=(c-b)(c+b)$, every companion $b$ comes from a factorisation $a^2=st$ with $s=c-b$, $t=c+b$, $s\lt t$ and $s\equiv t\pmod 2$, and then $b=(t-s)/2$. For $a=3$ the only candidate is $9=1\cdot9$, giving $b=4$; for $a=12$ the factorisations

$$144=2\cdot72=4\cdot36=6\cdot24=8\cdot18$$

give $b=35,16,9,5$.

The total has an exact formula. Put $L=\lfloor\sqrt n\rfloor$ and let $\mathcal T_L$ hold the primitive triples $(p,q,r)$ with $p\lt q\le L$, each pair of legs taken once; the legs of a primitive triple are never equal, since $2p^2$ is not a square, as Section 4 shows. Every leg pair is a unique multiple $k(p,q)$ of a primitive one, and $kq\le L$ is the binding constraint, so each triple together with each such $k$ gives the two distinct points $((kp)^2,(kq)^2)$ and $((kq)^2,(kp)^2)$, without repetition:

$$m=2\sum_{(p,q,r)\in\mathcal T_L}\left\lfloor\frac Lq\right\rfloor.$$

Enumerating this gives $m=126$, $2068$ and $7850$ at $n=10^4$, $10^6$ and $10^7$, from $18$, $179$ and $566$ primitive triples. Since the larger leg is at least $r/\sqrt2$, Lehmer's density for hypotenuses bounds $|\mathcal T_t|$ by a constant times $t$, ^[4] and partial summation turns the weights $L/q$ into $m=O(\sqrt n\log n)$: a few thousand points against $n=10^7$.

4. Why the two endpoints need no special care

Admissibility bans only intermediate points, so $O$ and $D$ ought to be exempted. Neither needs it: $0$ is not a positive square, so $O$ is admissible for every $n$, and if $D=(n,n)$ were inadmissible we would need $n=a^2$ and $2n=c^2$, hence $c^2=2a^2$ and $c/a=\sqrt2$, which is irrational. ^[5] Every inadmissible point of the square may therefore be forbidden outright — a fact the programs use silently rather than test.

5. Sorting the bad paths by the first point where they go wrong

A path that is not admissible meets at least one forbidden point, and since it visits its points in order it has a well-defined first one. Classifying bad paths by that point splits them into disjoint groups, each factoring into a beginning and a free end. On a tiny board, paths from $(0,0)$ to $(3,3)$ with $(1,1)$ and $(2,2)$ forbidden:

All paths: $\binom{6}{3}=20$.
First forbidden point $(1,1)$: every path reaching $(1,1)$ qualifies, nothing earlier being forbidden, so there are $\binom{2}{1}=2$.
First forbidden point $(2,2)$: of the $\binom{4}{2}=6$ paths to $(2,2)$, those through $(1,1)$ number $2\cdot\binom{2}{1}=4$, leaving $6-4=2$.
Bad paths to $(3,3)$: $2\cdot\binom{4}{2}+2\cdot\binom{2}{1}=16$, so $20-16=4$ survive.

The four survivors, with E for east and N for north, are EEENNN, EENENN, NNENEE and NNNEEE, as a direct grid count confirms. The groups are disjoint: a path through both forbidden points is counted only in the first.

6. The order in which the forbidden points must be processed

Write $P\preceq Q$ when $x_P\le x_Q$ and $y_P\le y_Q$; this is the product, or componentwise, order, ^[6] and by Section 1 a path runs from $P$ to $Q$ precisely when $P\preceq Q$. The recursion of Section 5 needs each forbidden point to come after everything it can be reached from, and lexicographic order does this, being a linear extension of the product order: ^[6] if $P\preceq Q$ and $P\ne Q$, then $x_P\lt x_Q$, or else $x_P=x_Q$ and $y_P\lt y_Q$, and either way $P$ comes first. Sorting by $x$, then $y$, suffices.

7. The recurrence and its proof

Let $P_1,\dots,P_m$ be the inadmissible points of the square in lexicographic order, $P_i=(x_i,y_i)$, and let $B_i$ count the paths from $O$ to $P_i$ whose first inadmissible point is $P_i$ itself. Every path from $O$ to $P_i$ has a first inadmissible point, which lies on the path and so is some $P_j\preceq P_i$; conversely the paths with that first point are a first-hit path to $P_j$, in $B_j$ ways, followed by any of the $N(P_j,P_i)$ paths onward. The groups are disjoint and exhaustive, so

$$N(O,P_i)=\sum_{P_j\preceq P_i}B_j\,N(P_j,P_i).$$

The term $j=i$ is $B_i$ itself, because $N(P_i,P_i)=1$, and by Section 6 every other index with $P_j\preceq P_i$ satisfies $j\lt i$. Hence

$$B_i=N(O,P_i)-\sum_{\substack{j\lt i\\ P_j\preceq P_i}}B_j\,N(P_j,P_i),$$

which determines $B_1,\dots,B_m$ one after another. The same partition at the target point, admissible by Section 4, gives the answer:

$$P(n)=N(O,D)-\sum_{i=1}^{m}B_i\,N(P_i,D).$$

This is a first-passage decomposition, not an alternating inclusion–exclusion over subsets of obstacles: every bad path lies in exactly one first-hit group, and subtracting those groups leaves the admissible paths.

8. Worked example: the checkpoint $P(16)$

For $n=16$ we have $\lfloor\sqrt{16}\rfloor=4$, and among $1\le a,b\le4$ only $(3,4)$ and $(4,3)$ give a square sum, so $m=2$ with $P_1=(9,16)$ and $P_2=(16,9)$. Neither dominates the other, so no path meets both and the sum in Section 7 is empty for both:

$$\begin{aligned}B_1=B_2&=N(O,P_1)=\binom{25}{9}\\&=2042975.\end{aligned}$$

From $P_1$ the only way on to $D=(16,16)$ is seven east steps, so $N(P_1,D)=\binom{7}{7}=1$, and likewise $N(P_2,D)=1$. With $\binom{32}{16}=601080390$ this gives

$$\begin{aligned}P(16)&=601080390-2\cdot2042975\\&=596994440,\end{aligned}$$

the value in the statement. For $n=1000$ the same procedure runs over $m=26$ points and reproduces the third checkpoint, $341920854$.

9. Worked example with a chain: $n=64$

At $n=16$ the first-hit counts need no subtraction between forbidden points. The smallest square that does is $n=64$, where the legs run up to $8$ and the triples $(3,4,5)$ and $(6,8,10)$ give, in lexicographic order,

$$\begin{aligned}P_1&=(9,16),&P_2&=(16,9),\\P_3&=(36,64),&P_4&=(64,36).\end{aligned}$$

Each of $P_1,P_2$ is dominated by each of $P_3,P_4$, so all four possible subtractions occur. Again $B_1=B_2=\binom{25}{9}$, and for $P_3$ the two coordinate differences are $(27,48)$ and $(20,55)$, so

$$B_3=\binom{100}{36}-B_1\binom{75}{27}-B_2\binom{75}{20},$$

which evaluates to $1600605235153583891574863675$; reflection in the diagonal gives $B_4=B_3$, and counting first-hit paths directly on a $64\times64$ grid returns the same four numbers. With $N(P_3,D)=\binom{28}{28}=1$ and $N(P_1,D)=\binom{103}{55}$, the final subtraction leaves a 38-digit $P(64)$; modulo $10^9+7$ it is $566774737$, which the C++ program prints for --n=64.

10. Binomial coefficients modulo $10^9+7$

Only residues modulo the prime $p=10^9+7$ are wanted, which is fortunate, since $\binom{2\cdot10^7}{10^7}$ has about $6020597$ decimal digits. Addition, subtraction and multiplication respect congruences, so the recurrence runs on residues throughout. Division cannot: dividing by a residue $a$ means multiplying by its inverse, the residue $b$ with $ab\equiv1$. Fermat's little theorem gives $a^{p-1}\equiv1\pmod p$ for $a\not\equiv0$, ^[7] hence $a^{-1}\equiv a^{p-2}$. Every factorial index here is at most $2n=2\cdot10^7$, below $p$, so no factorial has a factor divisible by $p$ and all these inverses exist.

Store $F_k\equiv k!$ and $I_k\equiv(k!)^{-1}$. Build $F_k$ forwards, obtain $I_{2n}$ by one modular exponentiation, and run $I_{k-1}\equiv kI_k$ backwards, an identity that follows from $k!=k\,(k-1)!$. Then every query is

$$\binom sr\equiv F_s\,I_r\,I_{s-r}\pmod p,$$

two modular multiplications. ^[8] A small check modulo $7$: $5!\equiv1$, the inverses of $2!$ and $3!$ are $4$ and $6$, and $1\cdot4\cdot6\equiv3$, matching $\binom52=10\equiv3$.

How the Code Works

The C++ program ^[9] reads --n= and --skip-checkpoints in parse_arguments and requires $n\ge1$; its helper uses std::stoi, which accepts a numeric prefix even when text follows, so --n=64junk also runs $n=64$, and nothing enforces the bound $2n\lt p$ of Section 10. Unless checkpoints are skipped, main calls P_mod for $5$, $16$ and $1000$ first, and each call uses the binary-exponentiation routine mod_pow once for its inverse-factorial seed, four times in all.

P_mod builds the tables of Section 10, then calls inadmissible_points for Section 2. That routine sets lim to $\lfloor\sqrt n\rfloor$ and lim2 to $\lfloor\sqrt{2n}\rfloor$, marks a byte array is_square of length $2n+1$ at every $c^2$ with $c\le$ lim2, and keeps $(a^2,b^2)$ for all $a,b\le$ lim with is_square[a*a + b*b] set. The pair loop looks up at most index $2\lfloor\sqrt n\rfloor^2$, the marking loop writes at most $\lfloor\sqrt{2n}\rfloor^2$, and both are at most $2n$, so every access stays inside the array. The lexicographic sort is what Section 6 requires; the following std::unique removes nothing, because $(a,b)\mapsto(a^2,b^2)$ is injective.

The tables are fact up to $2n$ and inv_fact, seeded at inv_fact[2n] with mod_pow(fact[2n], MOD - 2) and swept downwards, so nCr costs two multiplications and returns $0$ outside range. The first loop over pts fills ways[i], the $B_i$: from nCr(xi + yi, xi) it subtracts ways[j] * paths for every earlier point with xj <= xi && yj <= yi, adding MOD when the running value turns negative. After the sort xj <= xi can never fail, so only the second coordinate does real work. The second loop is the final subtraction, and run_checkpoints tests $P(5)$, $P(16)$ and $P(1000)$ beforehand.

Python and Java use the same method with $n$ fixed at $10^7$ and no checkpoints. Python builds its tables before the points, as C++ does, Java the points first; Python keeps the squares in the set is_sq, C++ in a byte array, Java in the boolean array isSquare. Python's bound comes from the exact math.isqrt, C++'s from truncating std::sqrt of a long double, Java's from Math.sqrt of a double; a separate comparison gives identical $\lfloor\sqrt n\rfloor$ and $\lfloor\sqrt{2n}\rfloor$ at $n=5,16,1000,10^6,10^7$. C++ and Java hold each partial result in 64 bits, where a product of two residues below $10^9+7$ stays under $2^{63}$; Python's integers are unbounded but costly, its two factorial lists holding $20000001$ references each, $40000002$ in all, with the integer objects stored separately.

Complexity and Verification

Marking the squares and enumerating the $L^2\le n$ candidate pairs is $O(n)$ time and memory, and so is building the two factorial arrays. The recurrence compares every point with every earlier one, costing $\Theta(m^2)$ time and $\Theta(m)$ extra memory, with $m=O(\sqrt n\log n)$ by Section 3; the sort adds only $O(m\log m)$. In total the method is $O(n+m^2)$ time and $O(n)$ memory, and at $n=10^7$ the two terms are comparable: with $m=7850$ the inner loop runs $m(m-1)/2=30807325$ times, of which $19362376$ pass the dominance test, against $4\cdot10^7$ factorial-table updates, $2\cdot10^7$ in each of the forward and backward passes. That quadratic term comes straight from the nested loops as written: each eligible earlier point carries a binomial weight depending on both coordinate differences, so a single running total cannot supply the sum. The small value of $m$ is what makes the direct pairwise evaluation practical.

The two factorial arrays hold $20000001$ 32-bit values each, $160000008$ bytes together, and the byte table adds $20000001$ more: about $180$ megabytes before the point vectors and allocator overhead. Compiled with -O2 and run at a size just below the requested one, the program took $0.42$ seconds and $182$ megabytes of resident memory on one machine, matching that calculation; elapsed time of course depends on build and machine.

These checks were run for this article. An exact brute force that counts paths on the grid with the forbidden points blocked, using unbounded integers and no binomial coefficients, agrees with the recurrence of Section 7 for every $n$ from $0$ to $250$; in particular it gives $P(n)=\binom{2n}{n}$ for $n\le15$ and the 38-digit $P(64)$. The counts $B_i$ at $n=64$ were obtained directly as well, by blocking all forbidden points except the target one, and matched Section 9. The three checkpoints hold, in the program's own pass and in an independent re-implementation. The factorisation characterisation of Section 3 was checked against a direct search for all $a\le60$, and every quoted $m$ was obtained twice, once by enumerating all pairs of legs and once from the primitive-triple formula.

Footnotes and References

The references below supply the background facts; the first-passage recurrence, its proof, the exact count of forbidden points and the worked examples are derived in this article.

Project Euler 408 — Admissible Paths Through a Grid. The official statement calls a lattice point inadmissible when x, y and x + y are all positive perfect squares, and a path of unit north and east steps admissible when no intermediate point is inadmissible; P(n) counts those from (0, 0) to (n, n). The checkpoints are P(5) = 252, P(16) = 596994440 and P(1000) mod 1000000007 = 341920854, and P(10000000) modulo 1000000007 is asked for.
Wikipedia — Lattice path. A north-east lattice path has steps (0, 1) and (1, 0), and the number of such paths from (0, 0) to (n, k) is C(n + k, n), with the example C(5, 2) = 10. This is the count used throughout as N(P, Q).
Wikipedia — Pythagorean triple. Euclid's formula a = m² − n², b = 2mn, c = m² + n² with m > n > 0 generates the triples, primitive exactly when m and n are coprime and of opposite parity, and every triple is a positive multiple of a primitive one — the facts Sections 2 and 3 use.
MathWorld — Pythagorean Triple. The page states that Lehmer (1900) proved the number of primitive triples with hypotenuse less than N, divided by N, tends to 1/(2π) ≈ 0.1591549. Section 3 needs only that this count is at most a constant times N.
Wikipedia — Square root of 2. Several proofs are given that the square root of two is irrational, among them infinite descent on a coprime fraction and a parity count of the factors of 2 in a² = 2b², the fact Sections 3 and 4 need.
Wikipedia — Product order. The page defines the product, or componentwise, order on a Cartesian product and states that the lexicographic combination of two total orders is a linear extension of their product order — why the sort of Section 6 is valid.
Wikipedia — Fermat's little theorem. For a prime p and any integer a, the number a^p − a is a multiple of p, and when a is not divisible by p this is equivalent to a^(p−1) ≡ 1 (mod p). The form a^(p−2) used in Section 10 follows; the page does not spell it out.
cp-algorithms — Binomial Coefficients. The page answers binomial queries modulo a prime in constant time as factorial[n] · inverse_factorial[k] · inverse_factorial[n − k], and gives (x!)^(−1) ≡ ((x−1)!)^(−1) · x^(−1), which rearranges into the backward sweep used in Section 10.
C++ — Euler408.cpp. The linked immutable revision defines parse_arguments, mod_pow, inadmissible_points, P_mod and run_checkpoints, and accepts --n= and --skip-checkpoints. It is the implementation record for this article.

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());
    }
}