Problem 292: Pythagorean Polygons

Project Euler Problem 292 Solution

EulerSolve provides an optimized solution for Project Euler Problem 292, Pythagorean Polygons, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The goal is to count lattice polygons whose edge lengths are integers and whose perimeter is at most $N$. For the Project Euler target we need $P(120)$. The code does not brute-force polygons directly; instead it builds them from admissible edge directions and uses dynamic programming on displacement and remaining perimeter. Mathematical Approach 1. Which edge directions are allowed? An edge vector $(x,y)$ has integer Euclidean length exactly when $$x^2+y^2=\ell^2$$ for some integer $\ell$. To avoid representing the same direction many times, the code keeps only primitive vectors $$\gcd(x,y)=1,$$ and later allows any positive multiple $g(x,y)$, whose length is $g\ell$. Because no edge of a genuine polygon can exceed half the total perimeter, it is enough to generate primitive vectors with $$\ell \le \left\lfloor \frac{N-1}{2}\right\rfloor.$$ 2. Why directions are processed in polar-angle order The primitive directions are sorted by angle $\theta=\operatorname{atan2}(y,x)$. This gives a canonical way to describe a polygon boundary: after choosing a distinguished lowest vertex as the origin, the boundary edges can be read in nondecreasing angular order. Without this step, the same polygon could be generated many times from different cyclic starting points or different presentations of parallel edges. 3....

Detailed mathematical approach

Problem Summary

A pythagorean polygon is a convex polygon with at least three vertices, no three of them on one line, integer coordinates at every vertex and an integer length for every edge. Translates count once, but a rotated or reflected copy is a new polygon unless it is also a translate. $P(n)$ is the number of pythagorean polygons with perimeter at most $n$. We are given $P(4)=1$, $P(30)=3655$ and $P(60)=891045$, and we must find $P(120)$. ^[1]

The whole solution in one sentence: walked anticlockwise from its lowest vertex, every such polygon is a list of edges whose directions turn steadily left, so we sort the possible directions by angle and use dynamic programming to decide for each direction whether to use it and how long to make that edge, remembering only the current point and the remaining perimeter.

Mathematical Approach

1. The smallest polygons

An edge is a lattice vector $(x,y)$ whose length $\sqrt{x^2+y^2}$ is an integer. None of $1,4,9,16$ is a sum $x^2+y^2$ with $x,y\ne0$, so edges of length $1$ to $4$ are horizontal or vertical; the first slanted edges are the eight vectors $(\pm3,\pm4)$ and $(\pm4,\pm3)$ of length $5$. Every edge is shorter than half the perimeter (Section 4), so up to perimeter $10$ all edges are horizontal or vertical, at most one in each direction (Section 5), and the polygon is an $a\times b$ rectangle with perimeter $2(a+b)$. The $1\times2$ and $2\times1$ rectangles are different, as neither is a translate of the other. Hence $P(4)=1$, the unit square, $P(6)=3$, $P(8)=6$ and $P(10)=10$. At perimeter $12$ a slanted edge has length $5$, so a polygon has at most two of them:

none: the $5$ rectangles with $a+b=6$;
one: the other edges must cancel its components $3$ and $4$ within the remaining length $7$, so they are one edge of length $3$ and one of length $4$, and we get $8$ right triangles, one for each slanted vector;
two opposite ones: the remaining length $2$ must be two opposite unit edges, which gives $4\cdot2=8$ parallelograms;
two others: the remaining axis edges, of total length at most $2$, can cancel their sum only if it is $(\pm1,\pm1)$, which happens for $4$ pairs such as $(3,4)$ and $(-4,-3)$; two perpendicular unit edges complete them to $4$ kites, such as the one with the vertices $(0,0)$, $(1,0)$, $(4,4)$, $(0,1)$.

No perimeter is odd (Section 3), so $P(12)=10+5+8+8+4=35$.

2. Edge directions

Write a nonzero lattice vector as $(x,y)=g\,(a,b)$ with $g=\gcd(x,y)$; the primitive vector $(a,b)$ is its direction. Since $x^2+y^2=g^2(a^2+b^2)$, the vector has an integer length $L$ exactly when its direction does: $\sqrt{a^2+b^2}=L/g$ is a fraction whose square is an integer, and such a fraction is itself an integer. The length is then $g\ell$ with $\ell=\sqrt{a^2+b^2}$: $(6,8)=2\,(3,4)$ has length $10$, while $(2,1)$, of length $\sqrt5$, never occurs. The directions of integer length are the four unit vectors and, for every primitive Pythagorean triple $a^2+b^2=c^2$, the eight vectors $(\pm a,\pm b)$ and $(\pm b,\pm a)$ of length $c$. Euclid's formula produces these triples as $(u^2-v^2,\,2uv,\,u^2+v^2)$ with coprime $u\gt v\gt0$ of different parity. ^[2] For $n=120$, Section 4 allows lengths up to $59$; the eight primitive triples with $c\le59$ have $c=5,13,17,25,29,37,41,53$, so there are $4+8\cdot8=68$ usable directions.

3. The perimeter is even

Every integer has the parity of its square, so an edge $(x,y)$ of length $L$ has $L\equiv x^2+y^2\equiv x+y\pmod 2$. Around a polygon the $x$-components add up to $0$, and so do the $y$-components, so the perimeter is even and $P(2k+1)=P(2k)$.

4. How long an edge can be

The other edges of a polygon form a path between the two ends of any edge. That path is not straight, since it passes through a vertex off the line of the edge, so it is longer than the edge: the polygon form of the triangle inequality. ^[5] An edge of length $e$ in a polygon with perimeter $p\le n$ therefore has $2e\lt p\le n$, that is,

$$e\le m=\left\lfloor\frac{n-1}{2}\right\rfloor,$$

which is $59$ for $n=120$. The direction of the edge has a length $\ell\le e$.

5. Convex means sorted by direction

Give every nonzero vector its angle $\theta\in[0,2\pi)$, measured anticlockwise from the positive $x$-axis. Walking anticlockwise around a pythagorean polygon, we turn left at every vertex by an exterior angle strictly between $0$ and $\pi$, since the interior angles of a convex polygon without three aligned vertices are below $\pi$; ^[3] the exterior angles add up to one full turn. ^[4] Start with the edge of smallest angle $\theta_1$ and let $\theta_j$ be $\theta_1$ plus the turns made before the $j$-th edge. Then $\theta_1\lt\theta_2\lt\dots\lt\theta_k\lt\theta_1+2\pi$, and no $\theta_j$ reaches $2\pi$, for $\theta_j-2\pi$ would be the angle of the $j$-th edge, smaller than $\theta_1$. Hence:

(a) Read from the start of its edge of smallest angle, the edges of a pythagorean polygon are lattice vectors with integer lengths, strictly increasing angles in $[0,2\pi)$ and the sum $(0,0)$.

Edges with angles in $[0,\pi]$ never go down and the others do, so the height rises and then falls back to $0$: the start is the lowest vertex, the left one of two. Placed at the origin, the polygon has the partial sums of its list as vertices, so a polygon up to translation is the same thing as its list.

(b) Conversely, take such a list with $k\ge3$ vectors. The gaps between cyclically consecutive angles, $\theta_1+2\pi-\theta_k$ included, are positive and add up to $2\pi$. A gap of at least $\pi$ would put all the vectors into a closed half-plane bounded by a line through the origin; their components across that line would be at least $0$ with sum $0$, hence $0$, so all the vectors would lie on that line, which allows only two distinct angles. So all gaps are below $\pi$. Now follow the walk from the end of an edge $e_j$ back to its start. The later edges, with increasing angles in $(\theta_j,\theta_j+2\pi)$, first move to the left of the line of $e_j$ (angles below $\theta_j+\pi$) and then back, and as the gaps next to $\theta_j$ are below $\pi$, the first one really moves away and the last one really moves back. So every vertex except the ends of $e_j$ lies strictly to the left of that line. As this holds for every edge, every edge is an edge of the convex hull of the vertices, and the walk goes once around this hull: the list is a convex polygon. If a vertex $B$ lay between vertices $A$ and $C$ on a line, the line of the edge leaving $B$ would contain $A$ or $C$ or separate them, contrary to this; so no three vertices are aligned.

For $k=2$ the only lists are $(v,-v)$, a segment walked there and back, which we call a digon. Therefore

$$P(n)=C(n)-1-D(n),$$

where $C(n)$ counts the sorted lists with sum $(0,0)$ and total length at most $n$, the $1$ is the empty list and $D(n)$ counts the digons.

6. Two polygons as sorted lists

The triangle with the vertices $(0,0)$, $(3,0)$, $(0,4)$ has the list $(3,0)$, $(-3,4)$, $(0,-4)$, with the angles $0^\circ$, $\approx126.87^\circ$, $270^\circ$ and the length $3+5+4=12$. The kite of Section 1 has the list $(1,0)$, $(3,4)$, $(-4,-3)$, $(0,-1)$, with the angles $0^\circ$, $\approx53.13^\circ$, $\approx216.87^\circ$, $270^\circ$ and the length $1+5+5+1=12$. Exchanging the first two edges of the kite gives a walk that turns right at $(3,4)$, so it is not convex.

7. Counting sorted lists direction by direction

Sort the directions $d_1,\dots,d_V$ by angle and write $\ell_j$ for the length of $d_j$. A sorted list is a choice of multiplicities $g_j\ge0$: the list consists of the edges $g_jd_j$ with $g_j\ge1$, in the order of $j$. A straight side is the single edge $g\,d_j$, never two edges in one direction, which is how the rule against three aligned vertices enters. This is a counting knapsack problem: the directions are the items, their lengths the weights and the perimeter the capacity. ^[6] Let $f_j(\mathbf p,r)$ count the choices of $g_1,\dots,g_j$ that lead from the origin to the point $\mathbf p$ with the length $n-r$, so that the budget $r$ is left. Then $f_0$ is $1$ at $(\mathbf 0,n)$ and $0$ elsewhere,

$$f_j(\mathbf p,r)=\sum_{g\ge0}f_{j-1}\bigl(\mathbf p-g\,d_j,\;r+g\,\ell_j\bigr),$$

where terms with a budget above $n$ are $0$, and $C(n)$ is the sum of $f_V(\mathbf 0,r)$ over $0\le r\le n$. The program keeps one map from points to counts for each budget. For $d_j$ it visits the budgets $i$ in increasing order and adds the count of every state $(\mathbf p,i)$ to the state $(\mathbf p+g\,d_j,\;i-g\ell_j)$ for $g=1,2,\dots$. The new states have smaller budgets, which were visited already, so a state created for $d_j$ is never extended by $d_j$ again: the map of budget $i$ still holds $f_{j-1}(\cdot,i)$ when it is read and holds $f_j(\cdot,i)$ after the pass. States are never removed, which covers $g_j=0$. For $n=4$ the first direction $(1,0)$ turns $((0,0),4)$ into $((1,0),3)$ and $((2,0),2)$; if budget $3$ were visited after budget $4$, the new state $((1,0),3)$ would be extended to $((2,0),2)$ a second time, counting the side of length $2$ twice.

8. Two tests that never lose a polygon

Before the program stores a state with the point $(x,y)$ and the budget $r$, it checks

$$y\ge0,\qquad x^2+y^2\le r^2.$$

Every beginning of a sorted list that closes within the budget passes both: its height first rises and then falls back to $0$ (Section 5), and its remaining edges, of total length at most $r$, lead back to the origin, while no path is shorter than the straight segment. ^[5] So $C(n)$ is still counted exactly, with much less work: for $n=30$ the tests cut the final number of states from $14800$ to $1924$, and both versions find $3702$ closed lists. The tests also exclude long edges. Before an edge of length $g\ell$ is added at budget $i$, the walk is at most $n-i$ away from the origin, and afterwards at least $g\ell-(n-i)$, so the second test requires $g\ell-(n-i)\le i-g\ell$, that is, $2g\ell\le n$. Finally, every stored point has $|x|\le n$ and $0\le y\le n$.

9. Removing the empty list and the digons

A digon with the edges $g\,v$ and $-g\,v$, where $v$ has the length $\ell$ and an angle below $\pi$, passes both tests exactly when $2g\ell\le n$. Each pair $\{v,-v\}$ thus gives $\lfloor n/(2\ell)\rfloor$ digons, and summing over all directions counts every pair twice:

$$D(n)=\frac12\sum_{v}\left\lfloor\frac{n}{2\ell(v)}\right\rfloor.$$

Any set of directions that contains all those of length at most $m$, and $-v$ with every $v$, gives the same $P(n)$ when used in both the table and this sum: polygons use only the short directions (Section 4), and each extra pair $\{v,-v\}$ adds the same digons to $C(n)$ and $D(n)$. For $n=4$ the directions are the four unit vectors. The $C(4)=6$ closed lists are the empty list, two digons of length $2$, two of length $4$ and the unit square, and $D(4)=\tfrac12\cdot4\cdot2=4$ gives $P(4)=6-1-4=1$. For $n=12$ the directions of length at most $m=5$ are the $4$ unit vectors and the $8$ slanted vectors of length $5$, so $D(12)=\tfrac12(4\cdot6+8\cdot1)=16$, and $C(12)=52$ gives $52-1-16=35$, the hand count of Section 1. In the same way $C(30)=3702$ and $D(30)=46$ give $3655$, and $C(60)=891150$ and $D(60)=104$ give $891045$.

How the Code Works

The C++ program ^[8] follows Sections 2 to 9. init_vectors scans every $(x,y)$ with $|x|,|y|\le m$ and keeps it when gcd2 returns $1$ and $x^2+y^2$ is a perfect square, tested by rounding the floating-point square root. A Direction stores the length, the coordinates and the angle from atan2, whose values lie in $(-\pi,\pi]$; ^[7] negative angles get $2\pi$ added. They are sorted by angle; the tie-breaks on the coordinates never act. The square is symmetric and holds every direction of length at most $m$, as Section 9 requires, and some longer ones: for $n=120$ it holds $84$ vectors, the $68$ of Section 2 and $16$ of the lengths $65$ and $73$ from $(33,56,65)$ and $(48,55,73)$, which are longer than $n/2$, never produce a state (Section 8) and add $0$ to $D(n)$.

count_polygons returns $0$ for $n\lt3$. Otherwise arr holds $n+1$ hash maps of type unordered_map<int, u64>, and arr[r] maps a packed point to its count with budget $r$. The key of $(x,y)$ is $(x+n)(n+1)+y$, unique because $0\le y\le n$ and $|x|\le n$; arr[n] starts with the origin, key $n(n+1)$, and count $1$. The loops are those of Section 7: for each direction, each budget i from its length up to $n$ with a nonempty map, and each g up to $\lfloor i/\ell\rfloor$, add_transitions decodes the keys of arr[i], moves each point by $g\,d$, applies the two tests with the new budget nrl $=i-g\ell$ and adds the count to arr[nrl]. Finally diagonal is $D(n)$, and the result is the sum of the origin counts over all maps minus diagonal minus $1$.

main reads --n=, by default $120$, and --skip-checkpoints; without the latter, run_checkpoints first recomputes $P(4)$, $P(30)$ and $P(60)$ and throws an exception if one of them is wrong. The Python and Java files contain the same three functions, in Java named initVectors, addTransitions and countPolygons, and no checkpoints; their entry points use $n=120$. The counts are unsigned 64-bit integers in C++, long values in Java's HashMap<Integer, Long> and unbounded integers in Python's dictionaries; all of them stay below $2^{40}$.

Complexity and Verification

A stored state $(x,y,r)$ satisfies $y\ge0$ and $x^2+y^2\le r^2$ by the tests, $x^2+y^2\le(n-r)^2$ because the walk so far has the length $n-r$, and $x+y\equiv n-r\pmod2$ by Section 3. A half-disc of radius $R$ holds about $\pi R^2/2$ lattice points, half of them with the right parity, so there are about

$$\sum_{r=0}^{n}\frac{\pi}{4}\min(r,n-r)^2\approx\frac{\pi n^3}{48},$$

such triples, $113\,097$ for $n=120$, and the memory is $O(n^3)$ entries. A check for every $n$ from $3$ to $120$ shows that the tables reach every one of these triples, $114\,919$ for $n=120$; the largest map has $2855$ entries, and the C++ program needs about 7 MB. Scanning the square costs $O(n^2\log n)$. A direction of length $\ell$ reads each entry at most $\lfloor n/\ell\rfloor$ times, so the table work is at most the number of states times $B$, the sum of $\lfloor n/\ell\rfloor$ over the generated directions; for $n=120$, $B=936$ and the bound is about $1.1\cdot10^8$ reads. In fact $25\,927\,148$ reads happen, $6\,761\,639$ of which pass the tests, since each direction sees only the states created before it. On our machine the C++ program, checkpoints included, needs about $0.1$ s, Java $0.3$ s and Python $3$ s.

For this article we verified:

The C++ checkpoints $P(4)=1$, $P(30)=3655$ and $P(60)=891045$ pass.
An independent search builds polygons vertex by vertex from their lowest, leftmost vertex, with any integer-length lattice vector as the next edge, checks every closed walk explicitly (strict left turns, one full turn, convexity, no three vertices on a line) and merges translates. Without directions, sorting or the digon correction, it agrees with the program for every $n\le60$, for example $P(20)=379$ and $P(40)=26613$, and finds only even perimeters.
A Python transcription of the tables reproduces all the counts above, and the C++, Python and Java programs print the same value for $n=120$.

Footnotes and References

The references supply the background facts; the correspondence between polygons and sorted lists, the two tests and all counts are derived above.

Project Euler 292 — Pythagorean Polygons. The official statement defines pythagorean polygons and P(n), the number of them with perimeter at most n, counting translates once, gives P(4) = 1, P(30) = 3655 and P(60) = 891045, and asks for P(120).
Wikipedia — Pythagorean triple. Euclid's formula a = m² − n², b = 2mn, c = m² + n², with coprime m > n > 0 of different parity, gives every primitive triple up to exchanging a and b; the article lists the 16 primitive triples with entries up to 100.
Wikipedia — Convex polygon. Every interior angle of a convex polygon is at most 180 degrees; it is strictly convex if no line contains more than two of its vertices, the condition that no three vertices are aligned.
Wikipedia — Internal and external angles. The external angles of a simple polygon, also called turning angles, add up to a full turn of 360 degrees.
Wikipedia — Triangle inequality. No polygonal path between two points is shorter than the straight line between them, and any side of a polygon is shorter than the sum of the other sides.
Wikipedia — Knapsack problem. The variants differ in how many copies of each item may be taken; the dynamic programming solution considers the items one at a time, with a table indexed by the capacity.
Wikipedia — atan2. atan2(y, x) is the angle between the positive x-axis and the ray from the origin to (x, y), with values in (−π, π].
C++ — Euler292.cpp. The linked immutable revision contains init_vectors, add_transitions, count_polygons and run_checkpoints; it is the implementation record for this article.

Problem 292 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Options {
    int perimeter_limit = 120;
    bool run_checkpoints = true;
    unsigned requested_threads = 0U;
};

struct Direction {
    int length = 0;
    int dx = 0;
    int dy = 0;
    double angle = 0.0;
};

int gcd2(int x, int y) {
    x = std::abs(x);
    y = std::abs(y);
    if (y == 0) return x;
    if (x == 0) return y;
    while (y != 0) {
        const int r = x % y;
        x = y;
        y = r;
    }
    return x;
}

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

    long long parsed = 0;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10LL + static_cast<long long>(c - '0');
        if (parsed > static_cast<long long>(std::numeric_limits<int>::max())) {
            return false;
        }
    }

    value = static_cast<int>(parsed);
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const std::string& prefix,
                                 unsigned& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10ULL + static_cast<u64>(c - '0');
        if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
            return false;
        }
    }

    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const 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.perimeter_limit)) {
            continue;
        }
        if (parse_unsigned_after_prefix(arg, "--threads=", options.requested_threads)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

std::vector<Direction> init_vectors(const int n) {
    std::vector<Direction> vectors;
    vectors.reserve(static_cast<std::size_t>(4 * n * n + 16));

    for (int y = -n; y <= n; ++y) {
        for (int x = -n; x <= n; ++x) {
            const int t = x * x + y * y;
            const int t2 = static_cast<int>(std::sqrt(static_cast<double>(t)) + 0.5);
            if (gcd2(x, y) != 1 || t2 * t2 != t) {
                continue;
            }

            double ang = std::atan2(static_cast<double>(y), static_cast<double>(x));
            if (ang < 0.0) {
                ang += 6.2831853071795862;
            }
            vectors.push_back(Direction{t2, x, y, ang});
        }
    }

    std::sort(vectors.begin(), vectors.end(),
              [](const Direction& a, const Direction& b) {
                  if (a.angle != b.angle) {
                      return a.angle < b.angle;
                  }
                  if (a.dx != b.dx) {
                      return a.dx < b.dx;
                  }
                  return a.dy < b.dy;
              });

    return vectors;
}

void add_transitions(const std::unordered_map<int, u64>& src,
                     const int dx,
                     const int dy,
                     const int nrl,
                     const int n,
                     std::unordered_map<int, u64>& dst) {
    const int base = n + 1;
    const int nrl2 = nrl * nrl;

    for (const auto& kv : src) {
        const int packed = kv.first;
        const u64 count = kv.second;

        const int cy = packed % base;
        const int cx = packed / base - n;
        const int nx = cx + dx;
        const int ny = cy + dy;
        const int np = nx * nx + ny * ny;

        if (ny >= 0 && np <= nrl2) {
            const int key = (nx + n) * base + ny;
            dst[key] += count;
        }
    }
}

u64 count_polygons(const int n) {
    if (n < 3) {
        return 0ULL;
    }

    const std::vector<Direction> vectors = init_vectors((n - 1) / 2);

    std::vector<std::unordered_map<int, u64>> arr(static_cast<std::size_t>(n + 1));
    for (auto& h : arr) {
        h.reserve(8);
    }

    const int origin_key = n * (n + 1);
    arr[static_cast<std::size_t>(n)][origin_key] = 1ULL;

    for (const Direction& v : vectors) {
        for (int i = v.length; i <= n; ++i) {
            const auto& src = arr[static_cast<std::size_t>(i)];
            if (src.empty()) {
                continue;
            }
            for (int g = 1; g <= i / v.length; ++g) {
                const int nrl = i - v.length * g;
                add_transitions(src, v.dx * g, v.dy * g, nrl, n, arr[static_cast<std::size_t>(nrl)]);
            }
        }
    }

    long long diagonal = 0;
    for (const Direction& v : vectors) {
        diagonal += n / (2 * v.length);
    }
    diagonal /= 2;

    long long answer = -diagonal - 1;
    for (int i = 0; i <= n; ++i) {
        const auto it = arr[static_cast<std::size_t>(i)].find(origin_key);
        if (it != arr[static_cast<std::size_t>(i)].end()) {
            answer += static_cast<long long>(it->second);
        }
    }

    return static_cast<u64>(answer);
}

void run_checkpoints() {
    struct Checkpoint {
        int n;
        u64 expected;
    };

    const std::vector<Checkpoint> checkpoints = {
        {4, 1ULL},
        {30, 3655ULL},
        {60, 891045ULL},
    };

    for (const Checkpoint& cp : checkpoints) {
        const u64 got = count_polygons(cp.n);
        if (got != cp.expected) {
            throw std::runtime_error("Checkpoint failed for P(" +
                                     std::to_string(cp.n) + "): got " +
                                     std::to_string(got) + ", expected " +
                                     std::to_string(cp.expected));
        }
    }
}

}  // namespace

int main(int argc, char** argv) {
    try {
        Options options;
        if (!parse_arguments(argc, argv, options)) {
            return 1;
        }

        if (options.perimeter_limit < 0) {
            throw std::invalid_argument("--n must be non-negative");
        }

        if (options.run_checkpoints) {
            run_checkpoints();
        }

        const u64 answer = count_polygons(options.perimeter_limit);
        std::cout << answer << '\n';
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

import math

class Direction:
    def __init__(self, length, dx, dy, angle):
        self.length = length
        self.dx = dx
        self.dy = dy
        self.angle = angle

def gcd2(x, y):
    x = abs(x)
    y = abs(y)
    if y == 0: return x
    if x == 0: return y
    while y != 0:
        x, y = y, x % y
    return x

def init_vectors(n):
    vectors = []
    for y in range(-n, n + 1):
        for x in range(-n, n + 1):
            t = x * x + y * y
            t2 = int(math.sqrt(t) + 0.5)
            if gcd2(x, y) != 1 or t2 * t2 != t:
                continue
            
            ang = math.atan2(y, x)
            if ang < 0.0:
                ang += 6.2831853071795862
            vectors.append(Direction(t2, x, y, ang))
            
    vectors.sort(key=lambda d: (d.angle, d.dx, d.dy))
    return vectors

def add_transitions(src, dx, dy, nrl, n, dst):
    base = n + 1
    nrl2 = nrl * nrl
    
    for packed, count in src.items():
        cy = packed % base
        cx = packed // base - n
        nx = cx + dx
        ny = cy + dy
        np = nx * nx + ny * ny
        
        if ny >= 0 and np <= nrl2:
            key = (nx + n) * base + ny
            dst[key] = dst.get(key, 0) + count

def count_polygons(n):
    if n < 3:
        return 0
        
    vectors = init_vectors((n - 1) // 2)
    
    arr = [{} for _ in range(n + 1)]
    origin_key = n * (n + 1)
    arr[n][origin_key] = 1
    
    for v in vectors:
        for i in range(v.length, n + 1):
            src = arr[i]
            if not src:
                continue
            for g in range(1, i // v.length + 1):
                nrl = i - v.length * g
                add_transitions(src, v.dx * g, v.dy * g, nrl, n, arr[nrl])
                
    diagonal = 0
    for v in vectors:
        diagonal += n // (2 * v.length)
    diagonal //= 2
    
    answer = -diagonal - 1
    for i in range(n + 1):
        if origin_key in arr[i]:
            answer += arr[i][origin_key]
            
    return answer

def solve(perimeter_limit=120):
    return str(count_polygons(perimeter_limit))

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

Java

import java.util.*;

public class Euler292 {
    static class Direction implements Comparable<Direction> {
        int length;
        int dx;
        int dy;
        double angle;

        Direction(int length, int dx, int dy, double angle) {
            this.length = length;
            this.dx = dx;
            this.dy = dy;
            this.angle = angle;
        }

        @Override
        public int compareTo(Direction b) {
            if (this.angle != b.angle) {
                return Double.compare(this.angle, b.angle);
            }
            if (this.dx != b.dx) {
                return Integer.compare(this.dx, b.dx);
            }
            return Integer.compare(this.dy, b.dy);
        }
    }

    static int gcd2(int x, int y) {
        x = Math.abs(x);
        y = Math.abs(y);
        if (y == 0)
            return x;
        if (x == 0)
            return y;
        while (y != 0) {
            int r = x % y;
            x = y;
            y = r;
        }
        return x;
    }

    static List<Direction> initVectors(int n) {
        List<Direction> vectors = new ArrayList<>();
        for (int y = -n; y <= n; ++y) {
            for (int x = -n; x <= n; ++x) {
                int t = x * x + y * y;
                int t2 = (int) Math.round(Math.sqrt(t));
                if (gcd2(x, y) != 1 || t2 * t2 != t) {
                    continue;
                }
                double ang = Math.atan2(y, x);
                if (ang < 0.0) {
                    ang += 6.2831853071795862;
                }
                vectors.add(new Direction(t2, x, y, ang));
            }
        }
        Collections.sort(vectors);
        return vectors;
    }

    static void addTransitions(Map<Integer, Long> src, int dx, int dy, int nrl, int n, Map<Integer, Long> dst) {
        int base = n + 1;
        int nrl2 = nrl * nrl;

        for (Map.Entry<Integer, Long> entry : src.entrySet()) {
            int packed = entry.getKey();
            long count = entry.getValue();

            int cy = packed % base;
            int cx = packed / base - n;
            int nx = cx + dx;
            int ny = cy + dy;
            int np = nx * nx + ny * ny;

            if (ny >= 0 && np <= nrl2) {
                int key = (nx + n) * base + ny;
                dst.put(key, dst.getOrDefault(key, 0L) + count);
            }
        }
    }

    static long countPolygons(int n) {
        if (n < 3)
            return 0;

        List<Direction> vectors = initVectors((n - 1) / 2);

        List<Map<Integer, Long>> arr = new ArrayList<>();
        for (int i = 0; i <= n; ++i) {
            arr.add(new HashMap<>());
        }

        int originKey = n * (n + 1);
        arr.get(n).put(originKey, 1L);

        for (Direction v : vectors) {
            for (int i = v.length; i <= n; ++i) {
                Map<Integer, Long> src = arr.get(i);
                if (src.isEmpty())
                    continue;

                for (int g = 1; g <= i / v.length; ++g) {
                    int nrl = i - v.length * g;
                    addTransitions(src, v.dx * g, v.dy * g, nrl, n, arr.get(nrl));
                }
            }
        }

        long diagonal = 0;
        for (Direction v : vectors) {
            diagonal += n / (2 * v.length);
        }
        diagonal /= 2;

        long answer = -diagonal - 1;
        for (int i = 0; i <= n; ++i) {
            if (arr.get(i).containsKey(originKey)) {
                answer += arr.get(i).get(originKey);
            }
        }

        return answer;
    }

    public static String solve() {
        return String.valueOf(countPolygons(120));
    }

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