Problem 594: Rhombus Tilings
View on Project EulerProject Euler Problem 594 Solution
EulerSolve provides an optimized solution for Project Euler Problem 594, Rhombus Tilings, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Let \(O_{a,b}\) be the equal-angled octagon whose side lengths, read cyclically, are \(a,b,a,b,a,b,a,b\). We must compute \(t(O_{4,2})\), where \(t(O_{a,b})\) counts tilings of this octagon by unit-edge rhombi and unit squares. The implementations do not use a closed formula; instead they count tilings recursively by describing every remaining subregion only through its boundary. Mathematical Approach The key observation is that every tile edge lies in one of eight directions separated by \(45^\circ\). After scaling coordinates by \(2\), those directions become $$ (2,0),\ (\sqrt{2},\sqrt{2}),\ (0,2),\ (-\sqrt{2},\sqrt{2}),\ (-2,0),\ (-\sqrt{2},-\sqrt{2}),\ (0,-2),\ (\sqrt{2},-\sqrt{2}). $$ Every coordinate can therefore be written as \(r+s\sqrt{2}\), which matches the algebraic representation used by the implementations. Step 1: Encode the Region by an Eight-Direction Boundary Word The initial octagon is encoded by the cyclic word $$ w_{a,b}=(\underbrace{0,\dots,0}_{a},\underbrace{1,\dots,1}_{b},\underbrace{2,\dots,2}_{a},\underbrace{3,\dots,3}_{b},\underbrace{4,\dots,4}_{a},\underbrace{5,\dots,5}_{b},\underbrace{6,\dots,6}_{a},\underbrace{7,\dots,7}_{b}). $$ More generally, any remaining subregion is represented by the cyclic sequence of its unit boundary steps. If \(w\) is such a word, let \(T(w)\) denote the number of tilings of the region enclosed by \(w\)....
Detailed mathematical approach
Problem Summary
Let \(O_{a,b}\) be the equal-angled octagon whose side lengths, read cyclically, are \(a,b,a,b,a,b,a,b\). We must compute \(t(O_{4,2})\), where \(t(O_{a,b})\) counts tilings of this octagon by unit-edge rhombi and unit squares. The implementations do not use a closed formula; instead they count tilings recursively by describing every remaining subregion only through its boundary.
Mathematical Approach
The key observation is that every tile edge lies in one of eight directions separated by \(45^\circ\). After scaling coordinates by \(2\), those directions become
$$ (2,0),\ (\sqrt{2},\sqrt{2}),\ (0,2),\ (-\sqrt{2},\sqrt{2}),\ (-2,0),\ (-\sqrt{2},-\sqrt{2}),\ (0,-2),\ (\sqrt{2},-\sqrt{2}). $$
Every coordinate can therefore be written as \(r+s\sqrt{2}\), which matches the algebraic representation used by the implementations.
Step 1: Encode the Region by an Eight-Direction Boundary Word
The initial octagon is encoded by the cyclic word
$$ w_{a,b}=(\underbrace{0,\dots,0}_{a},\underbrace{1,\dots,1}_{b},\underbrace{2,\dots,2}_{a},\underbrace{3,\dots,3}_{b},\underbrace{4,\dots,4}_{a},\underbrace{5,\dots,5}_{b},\underbrace{6,\dots,6}_{a},\underbrace{7,\dots,7}_{b}). $$
More generally, any remaining subregion is represented by the cyclic sequence of its unit boundary steps. If \(w\) is such a word, let \(T(w)\) denote the number of tilings of the region enclosed by \(w\).
Step 2: Choose a Canonical Boundary Edge
To avoid double counting, the recursion always selects the lowest boundary vertex, breaking ties by taking the leftmost one. Let that vertex be \(p\), and let the outgoing boundary edge from \(p\) have direction \(d\). In any valid tiling, exactly one tile uses that exposed boundary edge, so the recursion only needs to enumerate the possible tiles attached there.
Step 3: Only Three Local Tiles Are Possible
Let \(u\) be the unit vector in direction \(d\). The second side of the tile must turn into the interior by \(45^\circ\), \(90^\circ\), or \(135^\circ\). Therefore the only candidates use the unit vector in direction \(d+\delta \pmod 8\) with \(\delta\in\{1,2,3\}\), and their corners are
$$ p,\qquad p+u,\qquad p+u+v,\qquad p+v. $$
The case \(\delta=2\) gives a square; the cases \(\delta=1\) and \(\delta=3\) give the two rhombus orientations. No other convex unit-edge square or rhombus can share the chosen boundary edge and still lie inside the region.
Step 4: Enforce Geometric Validity
A candidate tile is accepted only when its two interior corners lie inside or on the current polygon and none of its four edges crosses an existing boundary edge, except for shared endpoints or the shared boundary edge itself. Cross products, segment tests, and vertex coordinates all stay inside the algebraic lattice generated by \(1\) and \(\sqrt{2}\), so the geometric tests follow the same exact combinatorial geometry as the tiles.
Step 5: Update the Boundary by Symmetric Difference
Write \(E(P)\) for the set of unit edges on the current boundary and \(E(Q)\) for the boundary of the chosen tile. When the tile is removed from the remaining region, shared edges disappear and newly exposed tile edges appear, so
$$ E(P\setminus Q)=E(P)\triangle E(Q), $$
where \(\triangle\) denotes symmetric difference. The resulting edge set is rebuilt into one or more counterclockwise boundary cycles by a deterministic walk: start from the lowest-leftmost unused edge and always take the smallest left turn. If the region splits into cycles \(w_1,\dots,w_r\), then the subproblems are independent and
$$ T(w)=\sum_{Q} \prod_{i=1}^{r} T(w_i), $$
with the empty boundary contributing \(1\).
Worked Example: \(O_{1,1}\)
For \(a=b=1\), the initial boundary word is
$$ w_{1,1}=(0,1,2,3,4,5,6,7). $$
The canonical starting vertex is the lowest-leftmost corner, and its outgoing boundary edge points east. The recursion therefore tests exactly three first tiles: one rhombus turning northeast, one square turning north, and one rhombus turning northwest. Summing all valid descendants gives
$$ t(O_{1,1})=8, $$
which matches one of the small checkpoint values used to validate the recurrence.
How the Code Works
The C++, Python, and Java implementations all follow the same recurrence. They build the initial boundary word for \(O_{4,2}\), expand it into vertices, choose the canonical lowest-leftmost start, try up to three candidate tiles, reject candidates that fail point-in-polygon or segment-intersection tests, and update the boundary through symmetric difference of unit edges. Every reconstructed cycle is oriented counterclockwise before memo lookup so that equivalent subregions share the same cache entry. The count is stored in arbitrary-precision integers, and the C++ implementation also checks small benchmark values such as \(t(O_{1,1})=8\), \(t(O_{2,1})=76\), and \(t(O_{3,2})=456572\).
Complexity Analysis
If a state has boundary length \(m\), then building its vertices, scanning for the canonical start, testing the three candidate tiles, checking segment intersections, and reconstructing the next boundary cycles each take linear time in \(m\). Thus one visited state costs \(O(m)\) time and \(O(m)\) temporary memory. The full search is still exponential in the worst case because the number of distinct boundary states can grow exponentially, but memoization, immediate geometric rejection, and decomposition into independent cycles reduce the practical search space sharply.
Footnotes and References
- Problem page: https://projecteuler.net/problem=594
- Tiling problem: Wikipedia — Tiling problem
- Quadratic integer: Wikipedia — Quadratic integer
- Symmetric difference: Wikipedia — Symmetric difference
- Point in polygon: Wikipedia — Point in polygon
Problem 594 source code
C++
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <future>
#include <iostream>
#include <map>
#include <string>
#include <thread>
#include <unordered_map>
#include <utility>
#include <vector>
#include <boost/multiprecision/cpp_int.hpp>
using namespace std;
using boost::multiprecision::cpp_int;
namespace {
constexpr long double kSqrt2 = 1.4142135623730950488016887242097L;
struct Num {
long long a = 0;
long long b = 0; // value = a + b*sqrt(2)
Num() = default;
Num(long long a_, long long b_) : a(a_), b(b_) {}
Num operator+(const Num& other) const { return {a + other.a, b + other.b}; }
Num operator-(const Num& other) const { return {a - other.a, b - other.b}; }
Num operator*(const Num& other) const {
// (a + b*sqrt(2))(c + d*sqrt(2)) = (ac + 2bd) + (ad + bc)*sqrt(2)
return {a * other.a + 2 * b * other.b, a * other.b + b * other.a};
}
bool operator==(const Num& other) const { return a == other.a && b == other.b; }
bool is_zero() const { return a == 0 && b == 0; }
long double value() const { return static_cast<long double>(a) + static_cast<long double>(b) * kSqrt2; }
};
struct Point {
Num x;
Num y;
bool operator==(const Point& other) const { return x == other.x && y == other.y; }
};
struct PointKey {
long long xa = 0;
long long xb = 0;
long long ya = 0;
long long yb = 0;
bool operator==(const PointKey& other) const {
return xa == other.xa && xb == other.xb && ya == other.ya && yb == other.yb;
}
};
PointKey key_of(const Point& p) { return {p.x.a, p.x.b, p.y.a, p.y.b}; }
Point point_from_key(const PointKey& k) { return {{k.xa, k.xb}, {k.ya, k.yb}}; }
struct PointKeyHash {
size_t operator()(const PointKey& k) const {
size_t h = std::hash<long long>{}(k.xa);
h ^= std::hash<long long>{}(k.xb) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
h ^= std::hash<long long>{}(k.ya) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
h ^= std::hash<long long>{}(k.yb) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
return h;
}
};
struct EdgeKey {
PointKey a;
PointKey b;
bool operator==(const EdgeKey& other) const { return a == other.a && b == other.b; }
};
struct EdgeKeyHash {
size_t operator()(const EdgeKey& k) const {
size_t h = PointKeyHash{}(k.a);
h ^= PointKeyHash{}(k.b) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
return h;
}
};
struct EdgeVal {
Point from;
Point to;
};
int sign_num(const Num& v) {
if (v.is_zero()) return 0;
return v.value() < 0 ? -1 : 1;
}
int cmp_num(const Num& a, const Num& b) { return sign_num(a - b); }
bool point_less_geom(const Point& a, const Point& b) {
int cy = cmp_num(a.y, b.y);
if (cy != 0) return cy < 0;
return cmp_num(a.x, b.x) < 0;
}
bool pointkey_less(const PointKey& a, const PointKey& b) {
if (a.xa != b.xa) return a.xa < b.xa;
if (a.xb != b.xb) return a.xb < b.xb;
if (a.ya != b.ya) return a.ya < b.ya;
return a.yb < b.yb;
}
Point operator+(const Point& a, const Point& b) { return {a.x + b.x, a.y + b.y}; }
Point operator-(const Point& a, const Point& b) { return {a.x - b.x, a.y - b.y}; }
Num cross(const Point& a, const Point& b, const Point& c) {
Point ab = b - a;
Point ac = c - a;
Num term1 = ab.x * ac.y;
Num term2 = ab.y * ac.x;
return term1 - term2;
}
Num min_num(const Num& a, const Num& b) { return cmp_num(a, b) <= 0 ? a : b; }
Num max_num(const Num& a, const Num& b) { return cmp_num(a, b) >= 0 ? a : b; }
bool on_segment(const Point& a, const Point& b, const Point& p) {
if (!cross(a, b, p).is_zero()) return false;
if (cmp_num(p.x, min_num(a.x, b.x)) < 0 || cmp_num(p.x, max_num(a.x, b.x)) > 0) return false;
if (cmp_num(p.y, min_num(a.y, b.y)) < 0 || cmp_num(p.y, max_num(a.y, b.y)) > 0) return false;
return true;
}
int orient(const Point& a, const Point& b, const Point& c) { return sign_num(cross(a, b, c)); }
bool seg_intersect(const Point& a, const Point& b, const Point& c, const Point& d) {
int o1 = orient(a, b, c);
int o2 = orient(a, b, d);
int o3 = orient(c, d, a);
int o4 = orient(c, d, b);
if (o1 == 0 && on_segment(a, b, c)) return true;
if (o2 == 0 && on_segment(a, b, d)) return true;
if (o3 == 0 && on_segment(c, d, a)) return true;
if (o4 == 0 && on_segment(c, d, b)) return true;
return (o1 != o2 && o3 != o4);
}
bool point_in_polygon(const Point& p, const vector<Point>& poly) {
long double px = p.x.value();
long double py = p.y.value();
bool inside = false;
size_t n = poly.size();
for (size_t i = 0; i < n; ++i) {
const Point& a = poly[i];
const Point& b = poly[(i + 1) % n];
if (on_segment(a, b, p)) return true;
long double ay = a.y.value();
long double by = b.y.value();
if ((ay > py) != (by > py)) {
long double ax = a.x.value();
long double bx = b.x.value();
long double xin = ax + (bx - ax) * (py - ay) / (by - ay);
if (xin > px) inside = !inside;
}
}
return inside;
}
EdgeKey make_edge_key(const Point& a, const Point& b) {
PointKey ka = key_of(a);
PointKey kb = key_of(b);
if (pointkey_less(ka, kb)) return {ka, kb};
return {kb, ka};
}
vector<Point> build_vertices(const vector<uint8_t>& steps, const array<Point, 8>& dirs) {
vector<Point> pts;
pts.reserve(steps.size() + 1);
Point cur{{0, 0}, {0, 0}};
pts.push_back(cur);
for (uint8_t d : steps) {
cur = cur + dirs[d];
pts.push_back(cur);
}
return pts;
}
long double signed_area(const vector<uint8_t>& steps, const array<Point, 8>& dirs) {
if (steps.empty()) return 0.0L;
vector<Point> pts = build_vertices(steps, dirs);
// drop repeated last point
pts.pop_back();
long double area = 0.0L;
for (size_t i = 0; i < pts.size(); ++i) {
const Point& a = pts[i];
const Point& b = pts[(i + 1) % pts.size()];
area += a.x.value() * b.y.value() - b.x.value() * a.y.value();
}
return 0.5L * area;
}
vector<uint8_t> to_ccw(const vector<uint8_t>& steps) {
vector<uint8_t> out;
out.reserve(steps.size());
for (size_t i = steps.size(); i > 0; --i) {
out.push_back(static_cast<uint8_t>((steps[i - 1] + 4) & 7));
}
return out;
}
unordered_map<EdgeKey, EdgeVal, EdgeKeyHash> edges_from_steps(
const vector<uint8_t>& steps, const array<Point, 8>& dirs) {
unordered_map<EdgeKey, EdgeVal, EdgeKeyHash> edges;
edges.reserve(steps.size() * 2);
Point cur{{0, 0}, {0, 0}};
for (uint8_t d : steps) {
Point nxt = cur + dirs[d];
edges.emplace(make_edge_key(cur, nxt), EdgeVal{cur, nxt});
cur = nxt;
}
return edges;
}
bool extract_cycles_from_edges(const unordered_map<EdgeKey, EdgeVal, EdgeKeyHash>& edge_map,
const array<Point, 8>& dirs,
vector<vector<uint8_t>>& cycles) {
// Walk each directed boundary cycle by keeping the interior on the left
// (smallest CCW turn from the incoming direction), even when components touch.
cycles.clear();
if (edge_map.empty()) return true;
struct EdgeInfo {
PointKey from;
PointKey to;
int dir = -1;
bool used = false;
};
vector<EdgeInfo> edges;
edges.reserve(edge_map.size());
unordered_map<PointKey, vector<int>, PointKeyHash> outgoing;
outgoing.reserve(edge_map.size() * 2);
for (const auto& kv : edge_map) {
const EdgeVal& e = kv.second;
Point diff = e.to - e.from;
int dir_idx = -1;
for (int i = 0; i < 8; ++i) {
if (diff.x == dirs[i].x && diff.y == dirs[i].y) {
dir_idx = i;
break;
}
}
if (dir_idx < 0) return false;
int id = static_cast<int>(edges.size());
edges.push_back({key_of(e.from), key_of(e.to), dir_idx, false});
outgoing[edges.back().from].push_back(id);
}
auto pick_start = [&]() -> int {
bool init = false;
Point min_point;
int best_id = -1;
for (int i = 0; i < static_cast<int>(edges.size()); ++i) {
if (edges[i].used) continue;
Point p = point_from_key(edges[i].from);
if (!init || point_less_geom(p, min_point)) {
min_point = p;
best_id = i;
init = true;
}
}
return best_id;
};
size_t edges_used = 0;
while (edges_used < edges.size()) {
int start_id = pick_start();
if (start_id < 0) return false;
int cur_id = start_id;
vector<uint8_t> steps;
while (true) {
EdgeInfo& cur = edges[cur_id];
if (cur.used) return false;
cur.used = true;
++edges_used;
steps.push_back(static_cast<uint8_t>(cur.dir));
PointKey head = cur.to;
auto it_out = outgoing.find(head);
if (it_out == outgoing.end()) return false;
int best_next = -1;
int best_delta = 9;
for (int cand_id : it_out->second) {
if (edges[cand_id].used && cand_id != start_id) continue;
int delta = (edges[cand_id].dir - cur.dir + 8) & 7;
if (delta < best_delta) {
best_delta = delta;
best_next = cand_id;
}
}
if (best_next < 0) return false;
if (best_next == start_id) {
break;
}
cur_id = best_next;
}
cycles.push_back(std::move(steps));
}
return edges_used == edges.size();
}
string steps_key(const vector<uint8_t>& steps) {
string key;
key.resize(steps.size());
for (size_t i = 0; i < steps.size(); ++i) key[i] = static_cast<char>(steps[i]);
return key;
}
struct Solver {
array<Point, 8> dirs;
map<string, cpp_int> memo;
Solver() {
// Directions are scaled by 2 to avoid halves: (x, y) = (a + b*sqrt(2), c + d*sqrt(2))
dirs[0] = {{2, 0}, {0, 0}}; // E
dirs[1] = {{0, 1}, {0, 1}}; // NE
dirs[2] = {{0, 0}, {2, 0}}; // N
dirs[3] = {{0, -1}, {0, 1}}; // NW
dirs[4] = {{-2, 0}, {0, 0}}; // W
dirs[5] = {{0, -1}, {0, -1}}; // SW
dirs[6] = {{0, 0}, {-2, 0}}; // S
dirs[7] = {{0, 1}, {0, -1}}; // SE
}
cpp_int dfs(const string& key) {
auto it = memo.find(key);
if (it != memo.end()) return it->second;
if (key.empty()) return cpp_int(1);
vector<uint8_t> steps;
steps.reserve(key.size());
for (unsigned char c : key) steps.push_back(static_cast<uint8_t>(c));
vector<Point> pts = build_vertices(steps, dirs);
size_t n = steps.size();
// Find the lowest (y, then x) vertex.
size_t min_idx = 0;
Point min_point = pts[0];
for (size_t i = 0; i < n; ++i) {
if (point_less_geom(pts[i], min_point)) {
min_point = pts[i];
min_idx = i;
}
}
const Point p = min_point;
const uint8_t d = steps[min_idx];
vector<pair<Point, Point>> boundary_edges;
boundary_edges.reserve(n);
for (size_t i = 0; i < n; ++i) {
boundary_edges.emplace_back(pts[i], pts[i + 1]);
}
cpp_int total = 0;
auto base_edges = edges_from_steps(steps, dirs);
for (int delta = 1; delta <= 3; ++delta) {
uint8_t vd = static_cast<uint8_t>((d + delta) & 7);
Point u = dirs[d];
Point v = dirs[vd];
Point p1 = p;
Point p2 = p + u;
Point p3 = p2 + v;
Point p4 = p + v;
if (!point_in_polygon(p4, pts) || !point_in_polygon(p3, pts)) continue;
bool ok = true;
array<pair<Point, Point>, 4> tile_edges = {{
{p1, p2}, {p2, p3}, {p3, p4}, {p4, p1},
}};
for (const auto& te : tile_edges) {
for (const auto& be : boundary_edges) {
const Point& a = te.first;
const Point& b = te.second;
const Point& c = be.first;
const Point& d2 = be.second;
if ((a == c && b == d2) || (a == d2 && b == c)) continue;
if (a == c || a == d2 || b == c || b == d2) continue;
if (seg_intersect(a, b, c, d2)) {
ok = false;
break;
}
}
if (!ok) break;
}
if (!ok) continue;
auto edges = base_edges;
array<pair<Point, Point>, 4> cw_edges = {{
{p1, p4}, {p4, p3}, {p3, p2}, {p2, p1},
}};
for (const auto& e : cw_edges) {
EdgeKey ek = make_edge_key(e.first, e.second);
auto it_edge = edges.find(ek);
if (it_edge != edges.end()) {
edges.erase(it_edge);
} else {
edges.emplace(ek, EdgeVal{e.first, e.second});
}
}
vector<vector<uint8_t>> cycles;
if (!extract_cycles_from_edges(edges, dirs, cycles)) continue;
if (cycles.empty()) {
total += 1;
continue;
}
cpp_int ways = 1;
for (auto& cycle : cycles) {
if (signed_area(cycle, dirs) < 0) cycle = to_ccw(cycle);
ways *= dfs(steps_key(cycle));
}
total += ways;
}
memo.emplace(key, total);
return total;
}
cpp_int solve(int a, int b) {
vector<uint8_t> steps;
steps.reserve(8 * (a + b));
array<int, 8> lens = {a, b, a, b, a, b, a, b};
for (int i = 0; i < 8; ++i) {
for (int j = 0; j < lens[i]; ++j) steps.push_back(static_cast<uint8_t>(i));
}
return dfs(steps_key(steps));
}
};
bool run_validation(unsigned threads) {
struct Case {
int a;
int b;
const char* expected;
};
const Case cases[] = {
{1, 1, "8"},
{2, 1, "76"},
{3, 2, "456572"},
};
bool ok = true;
for (const auto& c : cases) {
Solver solver;
cpp_int got = solver.solve(c.a, c.b);
std::string got_str = got.convert_to<std::string>();
if (got_str != c.expected) {
cerr << "Validation failed for O_" << c.a << "," << c.b
<< ": got " << got_str << ", expected " << c.expected << "\n";
ok = false;
}
}
if (ok) cerr << "Validation checkpoints passed.\n";
(void)threads;
return ok;
}
cpp_int solve_parallel(int a, int b, unsigned threads) {
Solver base;
vector<uint8_t> steps;
steps.reserve(8 * (a + b));
array<int, 8> lens = {a, b, a, b, a, b, a, b};
for (int i = 0; i < 8; ++i) {
for (int j = 0; j < lens[i]; ++j) steps.push_back(static_cast<uint8_t>(i));
}
// Build initial state and split on the first edge placements.
vector<Point> pts = build_vertices(steps, base.dirs);
size_t n = steps.size();
size_t min_idx = 0;
Point min_point = pts[0];
for (size_t i = 0; i < n; ++i) {
if (point_less_geom(pts[i], min_point)) {
min_point = pts[i];
min_idx = i;
}
}
Point p = min_point;
uint8_t d = steps[min_idx];
vector<pair<Point, Point>> boundary_edges;
boundary_edges.reserve(n);
for (size_t i = 0; i < n; ++i) boundary_edges.emplace_back(pts[i], pts[i + 1]);
auto base_edges = edges_from_steps(steps, base.dirs);
vector<string> sub_tasks;
vector<cpp_int> fixed_tasks;
for (int delta = 1; delta <= 3; ++delta) {
uint8_t vd = static_cast<uint8_t>((d + delta) & 7);
Point u = base.dirs[d];
Point v = base.dirs[vd];
Point p1 = p;
Point p2 = p + u;
Point p3 = p2 + v;
Point p4 = p + v;
if (!point_in_polygon(p4, pts) || !point_in_polygon(p3, pts)) continue;
bool ok = true;
array<pair<Point, Point>, 4> tile_edges = {{
{p1, p2}, {p2, p3}, {p3, p4}, {p4, p1},
}};
for (const auto& te : tile_edges) {
for (const auto& be : boundary_edges) {
const Point& a0 = te.first;
const Point& b0 = te.second;
const Point& c0 = be.first;
const Point& d0 = be.second;
if ((a0 == c0 && b0 == d0) || (a0 == d0 && b0 == c0)) continue;
if (a0 == c0 || a0 == d0 || b0 == c0 || b0 == d0) continue;
if (seg_intersect(a0, b0, c0, d0)) {
ok = false;
break;
}
}
if (!ok) break;
}
if (!ok) continue;
auto edges = base_edges;
array<pair<Point, Point>, 4> cw_edges = {{
{p1, p4}, {p4, p3}, {p3, p2}, {p2, p1},
}};
for (const auto& e : cw_edges) {
EdgeKey ek = make_edge_key(e.first, e.second);
auto it_edge = edges.find(ek);
if (it_edge != edges.end()) {
edges.erase(it_edge);
} else {
edges.emplace(ek, EdgeVal{e.first, e.second});
}
}
vector<vector<uint8_t>> cycles;
if (!extract_cycles_from_edges(edges, base.dirs, cycles)) continue;
if (cycles.empty()) {
fixed_tasks.push_back(1);
continue;
}
if (cycles.size() == 1) {
if (signed_area(cycles[0], base.dirs) < 0) cycles[0] = to_ccw(cycles[0]);
sub_tasks.push_back(steps_key(cycles[0]));
continue;
}
// Multiple components: solve independently and multiply here.
cpp_int ways = 1;
for (auto& cycle : cycles) {
if (signed_area(cycle, base.dirs) < 0) cycle = to_ccw(cycle);
Solver solver;
ways *= solver.dfs(steps_key(cycle));
}
fixed_tasks.push_back(ways);
}
if (threads <= 1 || sub_tasks.size() <= 1) {
Solver solver;
cpp_int total = 0;
for (const auto& val : fixed_tasks) total += val;
for (const auto& key : sub_tasks) total += solver.dfs(key);
return total;
}
vector<future<cpp_int>> futures;
futures.reserve(sub_tasks.size());
for (const auto& key : sub_tasks) {
futures.push_back(std::async(std::launch::async, [key]() {
Solver solver;
return solver.dfs(key);
}));
}
cpp_int total = 0;
for (const auto& val : fixed_tasks) total += val;
for (auto& f : futures) total += f.get();
return total;
}
} // namespace
int main(int argc, char** argv) {
ios::sync_with_stdio(false);
cin.tie(nullptr);
int a = 4;
int b = 2;
unsigned threads = thread::hardware_concurrency();
if (threads == 0) threads = 1;
threads = min(threads, 4u);
bool validate = true;
// Optional CLI: ./a.out [a] [b] [threads] [validate(0/1)]
if (argc >= 2) a = stoi(argv[1]);
if (argc >= 3) b = stoi(argv[2]);
if (argc >= 4) threads = static_cast<unsigned>(stoul(argv[3]));
if (argc >= 5) validate = (stoi(argv[4]) != 0);
if (validate && !run_validation(min(threads, 2u))) return 1;
cpp_int answer = solve_parallel(a, b, threads);
cout << answer << "\n";
return 0;
}
Python
import math
def solve():
SQRT2 = math.sqrt(2)
class Num:
__slots__=('a','b')
def __init__(self,a=0,b=0): self.a=a; self.b=b
def __add__(s,o): return Num(s.a+o.a,s.b+o.b)
def __sub__(s,o): return Num(s.a-o.a,s.b-o.b)
def __mul__(s,o): return Num(s.a*o.a+2*s.b*o.b,s.a*o.b+s.b*o.a)
def __eq__(s,o): return s.a==o.a and s.b==o.b
def __hash__(s): return hash((s.a,s.b))
def val(s): return s.a+s.b*SQRT2
def zero(s): return s.a==0 and s.b==0
class Pt:
__slots__=('x','y')
def __init__(s,x=None,y=None): s.x=x or Num(); s.y=y or Num()
def __add__(s,o): return Pt(s.x+o.x,s.y+o.y)
def __sub__(s,o): return Pt(s.x-o.x,s.y-o.y)
def __eq__(s,o): return s.x==o.x and s.y==o.y
def __hash__(s): return hash((s.x.a,s.x.b,s.y.a,s.y.b))
def key(s): return (s.x.a,s.x.b,s.y.a,s.y.b)
def sign_num(v):
if v.zero(): return 0
return -1 if v.val()<0 else 1
def cross(a,b,c):
ab=b-a; ac=c-a
return ab.x*ac.y-ab.y*ac.x
def on_seg(a,b,p):
cr=cross(a,b,p)
if not cr.zero(): return False
def mn(x,y): return x if sign_num(x-y)<=0 else y
def mx(x,y): return x if sign_num(x-y)>=0 else y
if sign_num(p.x-mn(a.x,b.x))<0 or sign_num(p.x-mx(a.x,b.x))>0: return False
if sign_num(p.y-mn(a.y,b.y))<0 or sign_num(p.y-mx(a.y,b.y))>0: return False
return True
def seg_ix(a,b,c,d):
o1=sign_num(cross(a,b,c)); o2=sign_num(cross(a,b,d))
o3=sign_num(cross(c,d,a)); o4=sign_num(cross(c,d,b))
if o1==0 and on_seg(a,b,c): return True
if o2==0 and on_seg(a,b,d): return True
if o3==0 and on_seg(c,d,a): return True
if o4==0 and on_seg(c,d,b): return True
return o1!=o2 and o3!=o4
def pip(p,poly):
px,py=p.x.val(),p.y.val(); inside=False; n=len(poly)
for i in range(n):
a,b=poly[i],poly[(i+1)%n]
if on_seg(a,b,p): return True
ay,by=a.y.val(),b.y.val()
if (ay>py)!=(by>py):
ax,bx=a.x.val(),b.x.val()
xin=ax+(bx-ax)*(py-ay)/(by-ay)
if xin>px: inside=not inside
return inside
dirs=[Pt(Num(2,0),Num(0,0)),Pt(Num(0,1),Num(0,1)),Pt(Num(0,0),Num(2,0)),
Pt(Num(0,-1),Num(0,1)),Pt(Num(-2,0),Num(0,0)),Pt(Num(0,-1),Num(0,-1)),
Pt(Num(0,0),Num(-2,0)),Pt(Num(0,1),Num(0,-1))]
def build_verts(steps):
pts=[Pt()]; cur=Pt()
for d in steps: cur=cur+dirs[d]; pts.append(cur)
return pts
def sa(steps):
if not steps: return 0.0
pts=build_verts(steps)[:-1]
area=0.0
for i in range(len(pts)):
a,b=pts[i],pts[(i+1)%len(pts)]
area+=a.x.val()*b.y.val()-b.x.val()*a.y.val()
return 0.5*area
def to_ccw(steps):
return [((s+4)&7) for s in reversed(steps)]
def make_ek(a,b):
ka,kb=a.key(),b.key()
return (ka,kb) if ka<kb else (kb,ka)
def edges_from(steps):
edges={}; cur=Pt()
for d in steps:
nxt=cur+dirs[d]; ek=make_ek(cur,nxt)
edges[ek]=(cur,nxt); cur=nxt
return edges
def extract_cycles(edge_map):
if not edge_map: return []
edges=[]; outgoing={}
for ek,(fr,to) in edge_map.items():
diff=to-fr; di=-1
for i in range(8):
if diff.x==dirs[i].x and diff.y==dirs[i].y: di=i; break
if di<0: return None
eid=len(edges); edges.append([fr.key(),to.key(),di,False])
outgoing.setdefault(fr.key(),[]).append(eid)
cycles=[]; used=0
def pick():
best=None; bi=-1
for i,(fk,_,_,u) in enumerate(edges):
if u: continue
fp=Pt(Num(fk[0],fk[1]),Num(fk[2],fk[3]))
if best is None or (fp.y.val(),fp.x.val())<(best.y.val(),best.x.val()):
best=fp; bi=i
return bi
while used<len(edges):
si=pick()
if si<0: return None
ci=si; steps=[]
while True:
e=edges[ci]
if e[3]: return None
e[3]=True; used+=1; steps.append(e[2])
head=e[1]; outs=outgoing.get(head,[])
bn=-1; bd=9
for cid in outs:
if edges[cid][3] and cid!=si: continue
delta=(edges[cid][2]-e[2]+8)&7
if delta<bd: bd=delta; bn=cid
if bn<0: return None
if bn==si: break
ci=bn
cycles.append(steps)
return cycles
memo={}
def dfs(key):
if key in memo: return memo[key]
if not key: return 1
steps=list(key); pts=build_verts(steps); n=len(steps)
mi=0; mp=pts[0]
for i in range(n):
if (pts[i].y.val(),pts[i].x.val())<(mp.y.val(),mp.x.val()): mp=pts[i]; mi=i
p=mp; d=steps[mi]
be=[(pts[i],pts[i+1]) for i in range(n)]
base_e=edges_from(steps)
total=0
for delta in range(1,4):
vd=(d+delta)&7; u=dirs[d]; v=dirs[vd]
p1=p; p2=p+u; p3=p2+v; p4=p+v
if not pip(p4,pts) or not pip(p3,pts): continue
te=[(p1,p2),(p2,p3),(p3,p4),(p4,p1)]; ok=True
for ta,tb in te:
for ca,cb in be:
if (ta==ca and tb==cb) or (ta==cb and tb==ca): continue
if ta==ca or ta==cb or tb==ca or tb==cb: continue
if seg_ix(ta,tb,ca,cb): ok=False; break
if not ok: break
if not ok: continue
edges=dict(base_e)
cw=[(p1,p4),(p4,p3),(p3,p2),(p2,p1)]
for ea,eb in cw:
ek=make_ek(ea,eb)
if ek in edges: del edges[ek]
else: edges[ek]=(ea,eb)
cycles=extract_cycles(edges)
if cycles is None: continue
if not cycles: total+=1; continue
ways=1
for cy in cycles:
if sa(cy)<0: cy=to_ccw(cy)
ways*=dfs(tuple(cy))
total+=ways
memo[key]=total
return total
a,b=4,2
steps=[]
lens=[a,b,a,b,a,b,a,b]
for i in range(8):
for _ in range(lens[i]): steps.append(i)
return str(dfs(tuple(steps)))
if __name__=='__main__':
print(solve())
Java
import java.util.*;
import java.math.BigInteger;
public class Euler594 {
static final double kSqrt2 = 1.4142135623730951;
static class Num {
long a, b;
Num(long a, long b) {
this.a = a;
this.b = b;
}
Num add(Num o) {
return new Num(a + o.a, b + o.b);
}
Num sub(Num o) {
return new Num(a - o.a, b - o.b);
}
Num mul(Num o) {
return new Num(a * o.a + 2 * b * o.b, a * o.b + b * o.a);
}
boolean equals(Num o) {
return a == o.a && b == o.b;
}
boolean isZero() {
return a == 0 && b == 0;
}
double value() {
return (double) a + (double) b * kSqrt2;
}
}
static class Point {
Num x, y;
Point(Num x, Num y) {
this.x = x;
this.y = y;
}
Point add(Point o) {
return new Point(x.add(o.x), y.add(o.y));
}
Point sub(Point o) {
return new Point(x.sub(o.x), y.sub(o.y));
}
boolean equals(Point o) {
return x.equals(o.x) && y.equals(o.y);
}
}
static class PointKey {
long xa, xb, ya, yb;
PointKey(long xa, long xb, long ya, long yb) {
this.xa = xa;
this.xb = xb;
this.ya = ya;
this.yb = yb;
}
@Override
public int hashCode() {
int h = Long.hashCode(xa);
h ^= Long.hashCode(xb) + 0x9e3779b9 + (h << 6) + (h >> 2);
h ^= Long.hashCode(ya) + 0x9e3779b9 + (h << 6) + (h >> 2);
h ^= Long.hashCode(yb) + 0x9e3779b9 + (h << 6) + (h >> 2);
return h;
}
@Override
public boolean equals(Object obj) {
PointKey o = (PointKey) obj;
return xa == o.xa && xb == o.xb && ya == o.ya && yb == o.yb;
}
}
static PointKey keyOf(Point p) {
return new PointKey(p.x.a, p.x.b, p.y.a, p.y.b);
}
static Point pointFromKey(PointKey k) {
return new Point(new Num(k.xa, k.xb), new Num(k.ya, k.yb));
}
static class EdgeKey {
PointKey a, b;
EdgeKey(PointKey a, PointKey b) {
this.a = a;
this.b = b;
}
@Override
public int hashCode() {
int h = a.hashCode();
h ^= b.hashCode() + 0x9e3779b9 + (h << 6) + (h >> 2);
return h;
}
@Override
public boolean equals(Object obj) {
EdgeKey o = (EdgeKey) obj;
return a.equals(o.a) && b.equals(o.b);
}
}
static class EdgeVal {
Point from, to;
EdgeVal(Point from, Point to) {
this.from = from;
this.to = to;
}
}
static int signNum(Num v) {
if (v.isZero())
return 0;
return v.value() < 0 ? -1 : 1;
}
static int cmpNum(Num a, Num b) {
return signNum(a.sub(b));
}
static boolean pointLessGeom(Point a, Point b) {
int cy = cmpNum(a.y, b.y);
if (cy != 0)
return cy < 0;
return cmpNum(a.x, b.x) < 0;
}
static boolean pointkeyLess(PointKey a, PointKey b) {
if (a.xa != b.xa)
return a.xa < b.xa;
if (a.xb != b.xb)
return a.xb < b.xb;
if (a.ya != b.ya)
return a.ya < b.ya;
return a.yb < b.yb;
}
static EdgeKey makeEdgeKey(Point a, Point b) {
PointKey ka = keyOf(a);
PointKey kb = keyOf(b);
if (pointkeyLess(ka, kb))
return new EdgeKey(ka, kb);
return new EdgeKey(kb, ka);
}
static Num cross(Point a, Point b, Point c) {
Point ab = b.sub(a);
Point ac = c.sub(a);
Num term1 = ab.x.mul(ac.y);
Num term2 = ab.y.mul(ac.x);
return term1.sub(term2);
}
static Num minNum(Num a, Num b) {
return cmpNum(a, b) <= 0 ? a : b;
}
static Num maxNum(Num a, Num b) {
return cmpNum(a, b) >= 0 ? a : b;
}
static boolean onSegment(Point a, Point b, Point p) {
if (!cross(a, b, p).isZero())
return false;
if (cmpNum(p.x, minNum(a.x, b.x)) < 0 || cmpNum(p.x, maxNum(a.x, b.x)) > 0)
return false;
if (cmpNum(p.y, minNum(a.y, b.y)) < 0 || cmpNum(p.y, maxNum(a.y, b.y)) > 0)
return false;
return true;
}
static int orient(Point a, Point b, Point c) {
return signNum(cross(a, b, c));
}
static boolean segIntersect(Point a, Point b, Point c, Point d) {
int o1 = orient(a, b, c);
int o2 = orient(a, b, d);
int o3 = orient(c, d, a);
int o4 = orient(c, d, b);
if (o1 == 0 && onSegment(a, b, c))
return true;
if (o2 == 0 && onSegment(a, b, d))
return true;
if (o3 == 0 && onSegment(c, d, a))
return true;
if (o4 == 0 && onSegment(c, d, b))
return true;
return (o1 != o2 && o3 != o4);
}
static boolean pointInPolygon(Point p, List<Point> poly) {
double px = p.x.value();
double py = p.y.value();
boolean inside = false;
int n = poly.size();
for (int i = 0; i < n; i++) {
Point a = poly.get(i);
Point b = poly.get((i + 1) % n);
if (onSegment(a, b, p))
return true;
double ay = a.y.value();
double by = b.y.value();
if ((ay > py) != (by > py)) {
double ax = a.x.value();
double bx = b.x.value();
double xin = ax + (bx - ax) * (py - ay) / (by - ay);
if (xin > px)
inside = !inside;
}
}
return inside;
}
static List<Point> buildVertices(List<Integer> steps, Point[] dirs) {
List<Point> pts = new ArrayList<>();
Point cur = new Point(new Num(0, 0), new Num(0, 0));
pts.add(cur);
for (int d : steps) {
cur = cur.add(dirs[d]);
pts.add(cur);
}
return pts;
}
static double signedArea(List<Integer> steps, Point[] dirs) {
if (steps.isEmpty())
return 0.0;
List<Point> pts = buildVertices(steps, dirs);
pts.remove(pts.size() - 1);
double area = 0.0;
int n = pts.size();
for (int i = 0; i < n; i++) {
Point a = pts.get(i);
Point b = pts.get((i + 1) % n);
area += a.x.value() * b.y.value() - b.x.value() * a.y.value();
}
return 0.5 * area;
}
static List<Integer> toCcw(List<Integer> steps) {
List<Integer> out = new ArrayList<>(steps.size());
for (int i = steps.size(); i > 0; i--) {
out.add((steps.get(i - 1) + 4) & 7);
}
return out;
}
static Map<EdgeKey, EdgeVal> edgesFromSteps(List<Integer> steps, Point[] dirs) {
Map<EdgeKey, EdgeVal> edges = new HashMap<>();
Point cur = new Point(new Num(0, 0), new Num(0, 0));
for (int d : steps) {
Point nxt = cur.add(dirs[d]);
edges.put(makeEdgeKey(cur, nxt), new EdgeVal(cur, nxt));
cur = nxt;
}
return edges;
}
static class EdgeInfo {
PointKey from, to;
int dir;
boolean used;
EdgeInfo(PointKey from, PointKey to, int dir) {
this.from = from;
this.to = to;
this.dir = dir;
this.used = false;
}
}
static boolean extractCyclesFromEdges(Map<EdgeKey, EdgeVal> edgeMap, Point[] dirs, List<List<Integer>> cycles) {
cycles.clear();
if (edgeMap.isEmpty())
return true;
List<EdgeInfo> edges = new ArrayList<>(edgeMap.size());
Map<PointKey, List<Integer>> outgoing = new HashMap<>();
for (EdgeVal e : edgeMap.values()) {
Point diff = e.to.sub(e.from);
int dirIdx = -1;
for (int i = 0; i < 8; i++) {
if (diff.equals(dirs[i])) {
dirIdx = i;
break;
}
}
if (dirIdx < 0)
return false;
int id = edges.size();
PointKey kFrom = keyOf(e.from);
PointKey kTo = keyOf(e.to);
edges.add(new EdgeInfo(kFrom, kTo, dirIdx));
outgoing.computeIfAbsent(kFrom, k -> new ArrayList<>()).add(id);
}
int edgesUsed = 0;
while (edgesUsed < edges.size()) {
int startId = -1;
boolean init = false;
Point minPoint = null;
for (int i = 0; i < edges.size(); i++) {
if (edges.get(i).used)
continue;
Point p = pointFromKey(edges.get(i).from);
if (!init || pointLessGeom(p, minPoint)) {
minPoint = p;
startId = i;
init = true;
}
}
if (startId < 0)
return false;
int curId = startId;
List<Integer> steps = new ArrayList<>();
while (true) {
EdgeInfo cur = edges.get(curId);
if (cur.used)
return false;
cur.used = true;
edgesUsed++;
steps.add(cur.dir);
List<Integer> outList = outgoing.get(cur.to);
if (outList == null)
return false;
int bestNext = -1;
int bestDelta = 9;
for (int candId : outList) {
if (edges.get(candId).used && candId != startId)
continue;
int delta = (edges.get(candId).dir - cur.dir + 8) & 7;
if (delta < bestDelta) {
bestDelta = delta;
bestNext = candId;
}
}
if (bestNext < 0)
return false;
if (bestNext == startId)
break;
curId = bestNext;
}
cycles.add(steps);
}
return edgesUsed == edges.size();
}
static String stepsKey(List<Integer> steps) {
StringBuilder sb = new StringBuilder(steps.size());
for (int c : steps)
sb.append((char) c);
return sb.toString();
}
static class Pair<K, V> {
K first;
V second;
Pair(K first, V second) {
this.first = first;
this.second = second;
}
}
static class Solver {
Point[] dirs = new Point[8];
Map<String, BigInteger> memo = new HashMap<>();
Solver() {
dirs[0] = new Point(new Num(2, 0), new Num(0, 0));
dirs[1] = new Point(new Num(0, 1), new Num(0, 1));
dirs[2] = new Point(new Num(0, 0), new Num(2, 0));
dirs[3] = new Point(new Num(0, -1), new Num(0, 1));
dirs[4] = new Point(new Num(-2, 0), new Num(0, 0));
dirs[5] = new Point(new Num(0, -1), new Num(0, -1));
dirs[6] = new Point(new Num(0, 0), new Num(-2, 0));
dirs[7] = new Point(new Num(0, 1), new Num(0, -1));
}
BigInteger dfs(String key) {
if (memo.containsKey(key))
return memo.get(key);
if (key.isEmpty())
return BigInteger.ONE;
List<Integer> steps = new ArrayList<>(key.length());
for (int i = 0; i < key.length(); i++) {
steps.add((int) key.charAt(i));
}
List<Point> pts = buildVertices(steps, dirs);
int n = steps.size();
int minIdx = 0;
Point minPoint = pts.get(0);
for (int i = 0; i < n; i++) {
if (pointLessGeom(pts.get(i), minPoint)) {
minPoint = pts.get(i);
minIdx = i;
}
}
Point p = minPoint;
int d = steps.get(minIdx);
List<Pair<Point, Point>> bounds = new ArrayList<>(n);
for (int i = 0; i < n; i++) {
bounds.add(new Pair<>(pts.get(i), pts.get(i + 1)));
}
BigInteger total = BigInteger.ZERO;
Map<EdgeKey, EdgeVal> baseEdges = edgesFromSteps(steps, dirs);
for (int delta = 1; delta <= 3; delta++) {
int vd = (d + delta) & 7;
Point u = dirs[d];
Point v = dirs[vd];
Point p1 = p;
Point p2 = p.add(u);
Point p3 = p2.add(v);
Point p4 = p.add(v);
if (!pointInPolygon(p4, pts) || !pointInPolygon(p3, pts))
continue;
boolean ok = true;
Pair<Point, Point>[] tileEdges = new Pair[] {
new Pair<>(p1, p2), new Pair<>(p2, p3), new Pair<>(p3, p4), new Pair<>(p4, p1)
};
for (Pair<Point, Point> te : tileEdges) {
for (Pair<Point, Point> be : bounds) {
Point a = te.first;
Point bX = te.second;
Point c = be.first;
Point d2 = be.second;
if ((a.equals(c) && bX.equals(d2)) || (a.equals(d2) && bX.equals(c)))
continue;
if (a.equals(c) || a.equals(d2) || bX.equals(c) || bX.equals(d2))
continue;
if (segIntersect(a, bX, c, d2)) {
ok = false;
break;
}
}
if (!ok)
break;
}
if (!ok)
continue;
Map<EdgeKey, EdgeVal> edges = new HashMap<>(baseEdges);
Pair<Point, Point>[] cwEdges = new Pair[] {
new Pair<>(p1, p4), new Pair<>(p4, p3), new Pair<>(p3, p2), new Pair<>(p2, p1)
};
for (Pair<Point, Point> e : cwEdges) {
EdgeKey ek = makeEdgeKey(e.first, e.second);
if (edges.containsKey(ek)) {
edges.remove(ek);
} else {
edges.put(ek, new EdgeVal(e.first, e.second));
}
}
List<List<Integer>> cycles = new ArrayList<>();
if (!extractCyclesFromEdges(edges, dirs, cycles))
continue;
if (cycles.isEmpty()) {
total = total.add(BigInteger.ONE);
continue;
}
BigInteger ways = BigInteger.ONE;
for (List<Integer> cycle : cycles) {
if (signedArea(cycle, dirs) < 0)
cycle = toCcw(cycle);
ways = ways.multiply(dfs(stepsKey(cycle)));
}
total = total.add(ways);
}
memo.put(key, total);
return total;
}
BigInteger solve(int a, int b) {
List<Integer> steps = new ArrayList<>(8 * (a + b));
int[] lens = { a, b, a, b, a, b, a, b };
for (int i = 0; i < 8; i++) {
for (int j = 0; j < lens[i]; j++) {
steps.add(i);
}
}
return dfs(stepsKey(steps));
}
}
public static String solve() {
Solver s = new Solver();
return s.solve(4, 2).toString();
}
public static void main(String[] args) {
System.out.println(solve());
}
}