2607.19268v1 / triangulation_search.cpp
all files
#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <unordered_map>
#include <utility>
#include <vector>
// Independent finite check for the paper. Every maximal planar graph is
// generated from K4 by splitting a vertex along two positions of its rotation.
// Whitney uniqueness lets us canonically label a triangulation by trying every
// directed root edge and both orientations of its spherical rotation system.
struct Code {
std::uint64_t lo = 0, hi = 0;
bool operator==(const Code& o) const { return lo == o.lo && hi == o.hi; }
bool operator<(const Code& o) const {
return hi < o.hi || (hi == o.hi && lo < o.lo);
}
};
struct CodeHash {
std::size_t operator()(const Code& c) const {
std::uint64_t x = c.lo ^ (c.hi + 0x9e3779b97f4a7c15ULL +
(c.lo << 6) + (c.lo >> 2));
x ^= x >> 30; x *= 0xbf58476d1ce4e5b9ULL;
x ^= x >> 27; x *= 0x94d049bb133111ebULL;
return static_cast<std::size_t>(x ^ (x >> 31));
}
};
using Rotation = std::vector<std::vector<std::uint8_t>>;
static bool adjacent(const Rotation& r, int u, int v) {
const auto& ns = r[u];
return std::find(ns.begin(), ns.end(), v) != ns.end();
}
static Code rooted_code(const Rotation& r, int root, int first, int dir) {
const int n = static_cast<int>(r.size());
std::array<int, 14> label, vertex, parent;
label.fill(-1); vertex.fill(-1); parent.fill(-1);
label[root] = 0; vertex[0] = root; parent[root] = first;
label[first] = 1; vertex[1] = first; parent[first] = root;
int discovered = 2;
for (int q = 0; q < discovered; ++q) {
const int x = vertex[q];
const auto& ns = r[x];
int start = 0;
while (ns[start] != parent[x]) ++start;
for (int step = 0; step < static_cast<int>(ns.size()); ++step) {
int pos = start + dir * step;
pos %= static_cast<int>(ns.size());
if (pos < 0) pos += static_cast<int>(ns.size());
const int y = ns[pos];
if (label[y] < 0) {
label[y] = discovered;
vertex[discovered++] = y;
parent[y] = x;
}
}
}
Code out;
int bit = 0;
for (int i = 0; i < n; ++i) {
for (int j = i + 1; j < n; ++j, ++bit) {
if (adjacent(r, vertex[i], vertex[j])) {
if (bit < 64) out.lo |= 1ULL << bit;
else out.hi |= 1ULL << (bit - 64);
}
}
}
return out;
}
static Code canonical_code(const Rotation& r) {
bool have = false;
Code best;
for (int u = 0; u < static_cast<int>(r.size()); ++u) {
for (int v : r[u]) {
for (int dir : {-1, 1}) {
Code c = rooted_code(r, u, v, dir);
if (!have || c < best) { best = c; have = true; }
}
}
}
return best;
}
static void insert_after(std::vector<std::uint8_t>& xs, int oldv, int newv) {
auto it = std::find(xs.begin(), xs.end(), oldv);
xs.insert(it + 1, static_cast<std::uint8_t>(newv));
}
static void insert_before(std::vector<std::uint8_t>& xs, int oldv, int newv) {
auto it = std::find(xs.begin(), xs.end(), oldv);
xs.insert(it, static_cast<std::uint8_t>(newv));
}
static Rotation split_vertex(const Rotation& r, int v, int i, int j) {
const int n = static_cast<int>(r.size());
const auto old = r[v];
const int d = static_cast<int>(old.size());
const int a = old[i], b = old[j];
Rotation s = r;
s.push_back({});
std::vector<std::uint8_t> arc_a, arc_b;
for (int k = i; k <= j; ++k) arc_a.push_back(old[k]);
for (int k = j; k < d; ++k) arc_b.push_back(old[k]);
for (int k = 0; k <= i; ++k) arc_b.push_back(old[k]);
arc_a.push_back(static_cast<std::uint8_t>(n));
arc_b.push_back(static_cast<std::uint8_t>(v));
s[v] = std::move(arc_a);
s[n] = std::move(arc_b);
insert_after(s[a], v, n);
insert_before(s[b], v, n);
for (int k = j + 1; k < d; ++k) {
auto it = std::find(s[old[k]].begin(), s[old[k]].end(), v);
*it = static_cast<std::uint8_t>(n);
}
for (int k = 0; k < i; ++k) {
auto it = std::find(s[old[k]].begin(), s[old[k]].end(), v);
*it = static_cast<std::uint8_t>(n);
}
return s;
}
static double spectral_radius(const Rotation& r) {
const int n = static_cast<int>(r.size());
std::array<double, 14> x{}, y{};
for (int i = 0; i < n; ++i) x[i] = 1.0 / std::sqrt(n);
for (int iter = 0; iter < 120; ++iter) {
y.fill(0.0);
for (int i = 0; i < n; ++i)
for (int j : r[i]) y[i] += x[j];
double norm = 0.0;
for (int i = 0; i < n; ++i) norm += y[i] * y[i];
norm = std::sqrt(norm);
for (int i = 0; i < n; ++i) x[i] = y[i] / norm;
}
double rq = 0.0;
for (int i = 0; i < n; ++i)
for (int j : r[i]) rq += x[i] * x[j];
return rq;
}
static void print_graph(const Rotation& r) {
std::cout << " degrees=[";
std::vector<int> ds;
for (const auto& x : r) ds.push_back(static_cast<int>(x.size()));
std::sort(ds.begin(), ds.end(), std::greater<int>());
for (std::size_t i = 0; i < ds.size(); ++i)
std::cout << (i ? "," : "") << ds[i];
std::cout << "] edges=";
bool first = true;
for (int i = 0; i < static_cast<int>(r.size()); ++i)
for (int j : r[i]) if (i < j) {
if (!first) std::cout << ',';
first = false;
std::cout << i << '-' << j;
}
std::cout << '\n';
}
int main(int argc, char** argv) {
int target = 14;
if (argc > 1) target = std::stoi(argv[1]);
if (target < 4 || target > 14) return 2;
Rotation k4 = {{1,2,3}, {0,3,2}, {0,1,3}, {0,2,1}};
std::unordered_map<Code, Rotation, CodeHash> current;
current.emplace(canonical_code(k4), k4);
for (int n = 4; n <= target; ++n) {
double best_rho = -1.0;
const Rotation* best = nullptr;
for (const auto& kv : current) {
double rho = spectral_radius(kv.second);
if (rho > best_rho) { best_rho = rho; best = &kv.second; }
}
std::cout << "n=" << n << " triangulations=" << current.size()
<< " max_rho=" << best_rho;
if (best) print_graph(*best);
if (n == target) break;
auto start = std::chrono::steady_clock::now();
std::unordered_map<Code, Rotation, CodeHash> next;
next.reserve(current.size() * 7);
for (const auto& kv : current) {
const Rotation& r = kv.second;
for (int v = 0; v < n; ++v) {
int d = static_cast<int>(r[v].size());
for (int i = 0; i < d; ++i)
for (int j = i + 1; j < d; ++j) {
Rotation child = split_vertex(r, v, i, j);
Code code = canonical_code(child);
next.try_emplace(code, std::move(child));
}
}
}
auto stop = std::chrono::steady_clock::now();
std::cerr << "generated n=" << n + 1 << " in "
<< std::chrono::duration<double>(stop - start).count() << "s\n";
current = std::move(next);
}
}