Problem 292: Pythagorean Polygons
View on Project EulerProject 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
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. DP state: current displacement and remaining perimeter
The state stored by the program is
$$ (X,Y,r), $$
where \((X,Y)\) is the current displacement from the chosen origin and \(r\) is the remaining perimeter budget. If we use the current primitive direction \((dx,dy)\) with primitive length \(\ell\) and multiplicity \(g\ge 1\), then the transition is
$$ (X,Y,r)\longrightarrow (X+g\,dx,\;Y+g\,dy,\;r-g\ell). $$
In other words, the DP is not guessing whole polygons at once; it is appending one angularly ordered edge block at a time.
4. The two pruning rules in the code
The transition is accepted only if two geometric filters hold.
Upper-half-plane rule. The code requires
$$Y_{\text{new}}\ge 0.$$
This chooses the lowest vertex of the polygon as the origin and forbids the partial walk from going below it. That gives one canonical representative instead of counting vertical translations or cyclic rotations of the same boundary.
Return-to-origin feasibility. The code also requires
$$X_{\text{new}}^2+Y_{\text{new}}^2\le r_{\text{new}}^2.$$
By the triangle inequality, if the Euclidean distance back to the origin already exceeds the remaining perimeter, closure is impossible. This single test removes a huge number of hopeless states.
5. Why returning to the origin is almost enough
After all directions have been processed, every state with displacement \((0,0)\) corresponds to a closed walk built from admissible integer-length edges. So the DP sums all states with origin key over every remaining budget level.
However, this raw count still includes two non-polygonal artifacts:
The empty walk. The initial state at the origin contributes one trivial closure and must be removed.
Two-edge backtracks. For any primitive direction of length \(\ell\), choosing an edge of length \(g\ell\) and then the opposite edge of the same length yields a closed degenerate segment of perimeter \(2g\ell\), not a genuine polygon. The code counts how many such objects exist via
$$\sum_{\text{primitive }v}\left\lfloor\frac{N}{2\ell(v)}\right\rfloor,$$
then divides by \(2\) because each direction and its opposite both appear in the vector list.
That is exactly why the final answer is computed as
$$\text{answer}=\Bigl(\text{all closed states}\Bigr)-1-\text{diagonal}.$$
6. Small checks
The implementation validates itself on
$$P(4)=1,\qquad P(30)=3655,\qquad P(60)=891045.$$
For \(N=4\), the only polygon is the unit square, so the first checkpoint is easy to visualize. The program then computes
$$P(120)=3600060866.$$
How the Code Works
The function init_vectors generates all primitive integer vectors with integer norm up to \((N-1)/2\), then sorts them by angle. The array arr[r] is a hash map of reachable displacements at remaining budget \(r\). It starts from the origin at budget \(N\). For each direction, the code tries every admissible multiplicity \(g\), applies the two pruning rules, and adds the resulting state into the hash map for the smaller remaining budget. Finally it sums all origin states and subtracts the empty path and the degenerate two-edge corrections.
Complexity Analysis
If \(V\) is the number of primitive directions and \(S_r\) is the number of reachable states at remaining budget \(r\), then the running time is roughly
$$O\!\left(\sum_{v\in V}\sum_r S_r\,M_{v,r}\right),$$
where \(M_{v,r}\) is the number of multiplicities \(g\) that still fit inside budget \(r\). Memory usage is the total size of the hash maps \(\texttt{arr}[0],\dots,\texttt{arr}[N]\). The essential gain comes from the geometric pruning, especially the distance-to-origin bound.
Further Reading
- Problem page: https://projecteuler.net/problem=292
- Pythagorean triples: https://en.wikipedia.org/wiki/Pythagorean_triple
- Lattice polygons: https://en.wikipedia.org/wiki/Lattice_polygon
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());
}
}