VoxelForge/Saved/Research/intervals/interval_spike.cpp

#include <algorithm>
#include <array>
#include <atomic>
#include <bit>
#include <cassert>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <functional>
#include <limits>
#include <numeric>
#include <span>
#include <string>
#include <thread>
#include <unordered_map>
#include <utility>
#include <vector>

#define WIN32_LEAN_AND_MEAN
#define NOMINMAX
#include <windows.h>
#include <psapi.h>
#pragma comment(lib, "Psapi.lib")

namespace {

using Clock = std::chrono::steady_clock;
constexpr int kCoreCells = 256;
constexpr int kHalo = 1;
constexpr int kNodes = kCoreCells + 1 + 2 * kHalo;
constexpr int kWordsPerLine = (kNodes + 63) / 64;
constexpr int kPageCells = 64;
constexpr int kPagesPerAxis = kCoreCells / kPageCells;
constexpr int kBucketCm = 400;
constexpr int kBucketCount = 17;
constexpr int kShapeLimitPerLine = 64;
constexpr int64_t kQScale = int64_t{1} << 28;
constexpr int64_t kSubPerCm = 16; // canonical interval/sign lattice: 0.625 mm
constexpr int64_t kSurfaceSubPerCm = 32; // crossing positions: 0.3125 mm
constexpr int kCellStepCm = 25;
constexpr int kGridMinCm = -3225;
constexpr int kGridMaxCm = kGridMinCm + (kNodes - 1) * kCellStepCm;
constexpr int kWorldOriginCm[3] = {-12000, 6000, -4500};
constexpr uint64_t kMask64 = ~uint64_t{0};

enum class Operation : uint8_t { CarveAir, AddRock };
enum class FeatureKind : uint8_t {
    Room, Tunnel, Pillar, Stalactite, Stalagmite, Shelf, FloatingMass,
    Stratum, Crater, CrowdedJunction
};

struct Mat3Q15 { int32_t v[3][3]; };

// Quantized rotations are fixed data, so corpus construction does not depend on libm.
constexpr std::array<Mat3Q15, 8> kRotations = {{
    {{{32768,0,0},{0,32768,0},{0,0,32768}}},
    {{{23170,-23170,0},{23170,23170,0},{0,0,32768}}},
    {{{23170,0,23170},{0,32768,0},{-23170,0,23170}}},
    {{{32768,0,0},{0,23170,-23170},{0,23170,23170}}},
    {{{0,-32768,0},{32768,0,0},{0,0,32768}}},
    {{{0,0,32768},{0,32768,0},{-32768,0,0}}},
    {{{32768,0,0},{0,0,-32768},{0,32768,0}}},
    {{{16384,-23170,16384},{16384,23170,16384},{-23170,0,23170}}}
}};

struct Shape {
    int32_t c[3]{};          // local centimetres
    int32_t radius[3]{};     // local centimetres
    int32_t rot[3][3]{};     // Q15 rotation, columns are local axes
    int64_t q[3][3]{};       // Q28 / cm^2, exact predicate coefficients
    int32_t lo[3]{}, hi[3]{}; // conservative local-centimetre AABB, inclusive
    Operation op = Operation::CarveAir;
    FeatureKind kind = FeatureKind::Room;
    uint32_t id = 0;
    bool auditThin = false;
};

struct Crossing {
    int32_t positionSub32 = 0; // local cm * 32, half of the canonical 1/16 cm grid
    uint32_t shapeId = 0;
    uint8_t beforeRock = 0;
    uint8_t afterRock = 0;
};

struct LineRange {
    uint32_t offset = 0;
    uint16_t count = 0;
    uint16_t candidates = 0;
    uint16_t primitiveIntervals = 0;
    uint16_t occupancyRuns = 0;
};

struct AxisData {
    std::vector<uint64_t> signs;
    std::vector<LineRange> lines;
    std::vector<Crossing> crossings;
};

struct AxisStats {
    double candidateP50 = 0, candidateP95 = 0;
    double intervalP50 = 0, intervalP95 = 0;
    double crossingP50 = 0, crossingP95 = 0;
    double runsP50 = 0, runsP95 = 0;
    uint32_t candidateMax = 0, intervalMax = 0, crossingMax = 0, runsMax = 0;
    uint64_t candidateTotal = 0, intervalTotal = 0, crossingTotal = 0, runsTotal = 0;
};

struct Vertex {
    float x, y, z;
    float nx, ny, nz;
};

struct TempVertex {
    uint64_t key = 0;
    uint32_t ownerPage = 0;
    Vertex v{};
};

struct PageRange { uint32_t page, indexStart, indexCount, vertexStart, vertexCount; };

struct MeshOutput {
    std::vector<TempVertex> vertices;
    std::vector<uint32_t> indices;
    std::vector<PageRange> pages;
    uint64_t hash = 0;
    uint64_t bytes = 0;
    uint64_t triangles = 0;
    uint64_t ambiguousFaces = 0;
    uint64_t multiCrossEdges = 0;
    uint64_t missingAnalyticCrossings = 0;
    uint64_t contourDegreeFaults = 0;
    uint64_t openContourLoops = 0;
    uint64_t explicitThinTriangles = 0;
    uint64_t thinFeaturesDetected = 0;
    uint64_t thinFeaturesHiddenInLattice = 0;
    uint64_t thinFeaturesMeshed = 0;
    uint64_t pageBorderVertices = 0;
};

struct Bin3D { std::vector<uint32_t> ids; };
struct ProjectionBins { std::vector<uint32_t> ids[kBucketCount * kBucketCount]; };

struct Corpus {
    std::vector<Shape> shapes;
    Bin3D bins3d[kBucketCount * kBucketCount * kBucketCount];
    ProjectionBins projections[3];
    uint64_t seed = 0;
    uint32_t maxMembershipAudit = 0;
};

struct WorkerLines {
    std::vector<Crossing> crossings;
};

struct FarOutput {
    std::vector<Vertex> vertices;
    std::vector<uint32_t> indices;
    std::vector<PageRange> patches;
    std::array<uint64_t, 6> lodNodes{};
    uint64_t nodes = 0, descendants = 0, patchCount = 0, bytes = 0, hash = 0;
};

double ms(Clock::time_point a, Clock::time_point b) {
    return std::chrono::duration<double, std::milli>(b - a).count();
}

uint64_t mix64(uint64_t x) {
    x += 0x9E3779B97F4A7C15ull;
    x = (x ^ (x >> 30)) * 0xBF58476D1CE4E5B9ull;
    x = (x ^ (x >> 27)) * 0x94D049BB133111EBull;
    return x ^ (x >> 31);
}

struct Rng {
    uint64_t state;
    explicit Rng(uint64_t s) : state(s) {}
    uint32_t next() { state += 0x9E3779B97F4A7C15ull; return static_cast<uint32_t>(mix64(state) >> 32); }
    int range(int lo, int hi) { return lo + static_cast<int>(next() % static_cast<uint32_t>(hi - lo + 1)); }
    bool chance(int n, int d) { return static_cast<int>(next() % d) < n; }
};

int64_t round_div(int64_t n, int64_t d) {
    if (n >= 0) return (n + d / 2) / d;
    return -((-n + d / 2) / d);
}

int floor_div(int a, int b) {
    const int q = a / b;
    const int r = a % b;
    return q - ((r != 0 && ((r < 0) != (b < 0))) ? 1 : 0);
}

int64_t floor_div64(int64_t a, int64_t b) {
    const int64_t q = a / b;
    const int64_t r = a % b;
    return q - ((r != 0 && ((r < 0) != (b < 0))) ? 1 : 0);
}

int ceil_div(int a, int b) { return -floor_div(-a, b); }

void build_quadric(Shape& s, uint32_t rotationIndex) {
    const auto& r = kRotations[rotationIndex % kRotations.size()];
    for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) s.rot[i][j] = r.v[i][j];
    // Q = R * diag(1/r^2) * R^T. Each term is rounded once to Q28.
    for (int a = 0; a < 3; ++a) for (int b = 0; b < 3; ++b) {
        int64_t sum = 0;
        for (int k = 0; k < 3; ++k) {
            const int64_t prod = static_cast<int64_t>(s.rot[a][k]) * s.rot[b][k];
            const int64_t den = (int64_t{1} << 30) * s.radius[k] * s.radius[k];
            sum += round_div(prod * kQScale, den);
        }
        s.q[a][b] = sum;
    }
    for (int axis = 0; axis < 3; ++axis) {
        int64_t extentNumerator = 0;
        for (int k = 0; k < 3; ++k) {
            const int64_t v = std::abs(s.rot[axis][k]);
            extentNumerator += v * s.radius[k];
        }
        const int32_t extent = static_cast<int32_t>((extentNumerator + 32767) / 32768) + 3;
        s.lo[axis] = s.c[axis] - extent;
        s.hi[axis] = s.c[axis] + extent;
    }
}

void add_shape(Corpus& c, int x, int y, int z, int rx, int ry, int rz,
               Operation op, FeatureKind kind, uint32_t rotation, bool auditThin = false) {
    Shape s;
    s.c[0] = x; s.c[1] = y; s.c[2] = z;
    s.radius[0] = std::max(1, rx); s.radius[1] = std::max(1, ry); s.radius[2] = std::max(1, rz);
    s.op = op; s.kind = kind; s.id = static_cast<uint32_t>(c.shapes.size()); s.auditThin = auditThin;
    build_quadric(s, rotation);
    c.shapes.push_back(s);
}

void add_sphere(Corpus& c, int x, int y, int z, int r, Operation op, FeatureKind kind) {
    add_shape(c, x, y, z, r, r, r, op, kind, 0);
}

void build_corpus(Corpus& c, uint64_t seed) {
    c.seed = seed;
    Rng rng(seed ^ 0x0A57A20261002ull);
    std::array<std::array<std::array<std::array<int, 3>, 3>, 3>, 2> centers{};
    for (int layer = 0; layer < 2; ++layer) {
        for (int iy = 0; iy < 3; ++iy) for (int ix = 0; ix < 3; ++ix) {
            const int x = (ix - 1) * 1300 + rng.range(-170, 170);
            const int y = (iy - 1) * 1300 + rng.range(-170, 170);
            const int z = (layer == 0 ? -850 : 850) + rng.range(-120, 120);
            centers[layer][iy][ix] = {x, y, z};
            const int rx = rng.range(480, 820);
            const int ry = rng.range(430, 760);
            const int rz = rng.range(300, 540);
            add_shape(c, x, y, z, rx, ry, rz, Operation::CarveAir, FeatureKind::Room,
                      rng.next() % 8);
        }
    }

    // Curved tunnel runs are unions of overlapping analytic ellipsoidal segments.
    for (int layer = 0; layer < 2; ++layer) {
        for (int iy = 0; iy < 3; ++iy) for (int ix = 0; ix < 3; ++ix) {
            const auto a = centers[layer][iy][ix];
            for (int dir = 0; dir < 2; ++dir) {
                const int nx = ix + (dir == 0);
                const int ny = iy + (dir == 1);
                if (nx >= 3 || ny >= 3) continue;
                const auto b = centers[layer][ny][nx];
                const int cx = (a[0] + b[0]) / 2 + rng.range(-220, 220);
                const int cy = (a[1] + b[1]) / 2 + rng.range(-220, 220);
                const int cz = (a[2] + b[2]) / 2 + rng.range(-150, 150);
                for (int seg = 0; seg < 5; ++seg) {
                    const int t = 1 + seg;
                    const int x = (a[0] * (6 - t) + b[0] * t) / 6 + (seg == 2 ? cx - (a[0] + b[0]) / 2 : 0);
                    const int y = (a[1] * (6 - t) + b[1] * t) / 6 + (seg == 2 ? cy - (a[1] + b[1]) / 2 : 0);
                    const int z = (a[2] * (6 - t) + b[2] * t) / 6 + (seg == 2 ? cz - (a[2] + b[2]) / 2 : 0);
                    const int r = rng.range(200, 300);
                    add_sphere(c, x, y, z, r, Operation::CarveAir, FeatureKind::Tunnel);
                }
            }
        }
    }
    // Vertical connectors between stacked room layers.
    for (int iy = 0; iy < 3; ++iy) for (int ix = 0; ix < 3; ++ix) {
        const auto a = centers[0][iy][ix], b = centers[1][iy][ix];
        for (int seg = 1; seg <= 3; ++seg) {
            add_sphere(c, (a[0] * (4 - seg) + b[0] * seg) / 4,
                       (a[1] * (4 - seg) + b[1] * seg) / 4,
                       (a[2] * (4 - seg) + b[2] * seg) / 4,
                       rng.range(190, 270), Operation::CarveAir, FeatureKind::Tunnel);
        }
    }

    // Dense constructive formations, placed in chamber space or near chamber walls.
    for (int i = 0; i < 40; ++i) {
        const int x = rng.range(-2500, 2500), y = rng.range(-2500, 2500);
        const int z = rng.range(-1900, 1900);
        add_shape(c, x, y, z, rng.range(70, 150), rng.range(70, 150), rng.range(260, 520),
                  Operation::AddRock, FeatureKind::Pillar, rng.next() % 8);
    }
    for (int i = 0; i < 36; ++i) {
        const auto center = centers[i % 2][(i / 2) % 3][(i / 6) % 3];
        const int sign = (i & 1) ? 1 : -1;
        add_shape(c, center[0] + rng.range(-260, 260), center[1] + rng.range(-260, 260),
                  center[2] + sign * rng.range(260, 410), rng.range(45, 85),
                  rng.range(45, 85), rng.range(230, 370), Operation::AddRock,
                  (i & 1) ? FeatureKind::Stalagmite : FeatureKind::Stalactite, rng.next() % 8);
    }
    for (int i = 0; i < 42; ++i) {
        const auto center = centers[i % 2][(i / 2) % 3][(i / 6) % 3];
        add_shape(c, center[0] + rng.range(-350, 350), center[1] + rng.range(-350, 350),
                  center[2] + rng.range(-220, 220), rng.range(250, 650), rng.range(130, 320),
                  rng.range(30, 90), Operation::AddRock, FeatureKind::Shelf, rng.next() % 8);
    }
    for (int i = 0; i < 24; ++i) {
        const auto center = centers[i % 2][(i / 2) % 3][(i / 6) % 3];
        add_shape(c, center[0] + rng.range(-280, 280), center[1] + rng.range(-280, 280),
                  center[2] + rng.range(-250, 250), rng.range(180, 390), rng.range(150, 340),
                  rng.range(140, 330), Operation::AddRock, FeatureKind::FloatingMass, rng.next() % 8);
    }
    // Several crater cuts, including deliberately intersecting crater groups.
    for (int i = 0; i < 12; ++i) {
        const auto center = centers[i % 2][(i / 2) % 3][(i / 6) % 3];
        const int x = center[0] + rng.range(300, 650);
        const int y = center[1] + rng.range(-280, 280);
        const int z = center[2] + rng.range(-280, 280);
        add_shape(c, x, y, z, rng.range(230, 520), rng.range(210, 460), rng.range(180, 420),
                  Operation::CarveAir, FeatureKind::Crater, rng.next() % 8);
    }
    // Crowded junction: a bounded 40-shape overlap cluster exercises long event lists.
    for (int i = 0; i < 40; ++i) {
        const int x = rng.range(-110, 110), y = rng.range(-110, 110), z = rng.range(-110, 110);
        const Operation op = (i % 3 == 0) ? Operation::AddRock : Operation::CarveAir;
        const int r = rng.range(75, 190);
        add_shape(c, x, y, z, r, rng.range(70, 200), rng.range(75, 185), op,
                  FeatureKind::CrowdedJunction, rng.next() % 8);
    }
    // Thin strata and shelves, with several widths below the 25 cm structural cell size.
    for (int i = 0; i < 12; ++i) {
        const auto center = centers[i % 2][(i / 2) % 3][(i / 6) % 3];
        add_shape(c, center[0] + rng.range(-180, 180), center[1] + rng.range(-180, 180),
                  center[2] + rng.range(-180, 180), rng.range(300, 580), rng.range(250, 520),
                  rng.range(4, 11), Operation::AddRock, FeatureKind::Stratum, rng.next() % 8);
    }
    // A guaranteed isolated sub-cell air seam in solid rock. No 25 cm lattice node lies inside it.
    add_shape(c, 2900, -2900, 2887, 500, 400, 6, Operation::CarveAir,
              FeatureKind::Stratum, 0, true);

    // Build deterministic 3D and axis-projection bins in feature order.
    for (uint32_t id = 0; id < c.shapes.size(); ++id) {
        const Shape& s = c.shapes[id];
        const int bx0 = std::clamp(floor_div(s.lo[0] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int bx1 = std::clamp(floor_div(s.hi[0] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int by0 = std::clamp(floor_div(s.lo[1] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int by1 = std::clamp(floor_div(s.hi[1] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int bz0 = std::clamp(floor_div(s.lo[2] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int bz1 = std::clamp(floor_div(s.hi[2] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        for (int bz = bz0; bz <= bz1; ++bz) for (int by = by0; by <= by1; ++by)
            for (int bx = bx0; bx <= bx1; ++bx)
                c.bins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids.push_back(id);
        for (int axis = 0; axis < 3; ++axis) {
            int a = (axis + 1) % 3, b = (axis + 2) % 3;
            int a0 = std::clamp(floor_div(s.lo[a] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
            int a1 = std::clamp(floor_div(s.hi[a] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
            int b0 = std::clamp(floor_div(s.lo[b] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
            int b1 = std::clamp(floor_div(s.hi[b] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
            for (int ib = b0; ib <= b1; ++ib) for (int ia = a0; ia <= a1; ++ia)
                c.projections[axis].ids[ib * kBucketCount + ia].push_back(id);
        }
    }
}

int bucket_of_sub(int64_t qSub) {
    const int64_t qCmFloor = floor_div64(qSub, kSubPerCm);
    return std::clamp(floor_div(static_cast<int>(qCmFloor) - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
}

int64_t eval_quad(const Shape& s, const int64_t d[3]) {
    int64_t value = 0;
    value += s.q[0][0] * d[0] * d[0];
    value += s.q[1][1] * d[1] * d[1];
    value += s.q[2][2] * d[2] * d[2];
    value += 2 * s.q[0][1] * d[0] * d[1];
    value += 2 * s.q[0][2] * d[0] * d[2];
    value += 2 * s.q[1][2] * d[1] * d[2];
    return value;
}

bool contains_sub(const Shape& s, const int64_t pSub[3]) {
    int64_t d[3] = {pSub[0] - static_cast<int64_t>(s.c[0]) * kSubPerCm,
                    pSub[1] - static_cast<int64_t>(s.c[1]) * kSubPerCm,
                    pSub[2] - static_cast<int64_t>(s.c[2]) * kSubPerCm};
    return eval_quad(s, d) <= kQScale * kSubPerCm * kSubPerCm;
}

bool shape_aabb_contains_sub(const Shape& s, const int64_t pSub[3]) {
    for (int a = 0; a < 3; ++a) {
        if (pSub[a] < static_cast<int64_t>(s.lo[a]) * kSubPerCm
            || pSub[a] > static_cast<int64_t>(s.hi[a]) * kSubPerCm) return false;
    }
    return true;
}

bool point_rock(const Corpus& c, const int64_t pSub[3]) {
    const int xcm = static_cast<int>(floor_div(static_cast<int>(pSub[0]), static_cast<int>(kSubPerCm)));
    const int ycm = static_cast<int>(floor_div(static_cast<int>(pSub[1]), static_cast<int>(kSubPerCm)));
    const int zcm = static_cast<int>(floor_div(static_cast<int>(pSub[2]), static_cast<int>(kSubPerCm)));
    const int bx = std::clamp(floor_div(xcm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
    const int by = std::clamp(floor_div(ycm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
    const int bz = std::clamp(floor_div(zcm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
    bool rock = true;
    const auto& bucket = c.bins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids;
    for (uint32_t id : bucket) {
        const Shape& s = c.shapes[id];
        if (shape_aabb_contains_sub(s, pSub) && contains_sub(s, pSub)) rock = (s.op == Operation::AddRock);
    }
    return rock;
}

struct Interval { int64_t loSub, hiSub; int loSample, hiSample; };

int64_t eval_line(const Shape& s, int axis, int64_t variableSub,
                  const int64_t fixedSub[3]) {
    int64_t p[3] = {fixedSub[0], fixedSub[1], fixedSub[2]};
    p[axis] = variableSub;
    int64_t d[3] = {p[0] - static_cast<int64_t>(s.c[0]) * kSubPerCm,
                    p[1] - static_cast<int64_t>(s.c[1]) * kSubPerCm,
                    p[2] - static_cast<int64_t>(s.c[2]) * kSubPerCm};
    return eval_quad(s, d) - kQScale * kSubPerCm * kSubPerCm;
}

bool solve_shape_line(const Shape& s, int axis, const int64_t fixedSub[3], Interval& out) {
    const int a = axis;
    const int b = (axis + 1) % 3;
    const int dAxis = (axis + 2) % 3;
    const int64_t da = fixedSub[b] - static_cast<int64_t>(s.c[b]) * kSubPerCm;
    const int64_t db = fixedSub[dAxis] - static_cast<int64_t>(s.c[dAxis]) * kSubPerCm;
    const int64_t A = s.q[a][a];
    const int64_t B = 2 * (s.q[a][b] * da + s.q[a][dAxis] * db);
    if (A <= 0) return false;
    const int64_t vertexFloor = static_cast<int64_t>(s.c[axis]) * kSubPerCm
        + floor_div64(-B, 2 * A);
    const int64_t minSub = static_cast<int64_t>(s.lo[a]) * kSubPerCm;
    const int64_t maxSub = static_cast<int64_t>(s.hi[a]) * kSubPerCm;
    int64_t center = std::clamp<int64_t>(vertexFloor, minSub, maxSub);
    int64_t best = center;
    int64_t bestValue = eval_line(s, axis, center, fixedSub);
    for (int64_t q = center - 1; q <= center + 1; ++q) {
        if (q < minSub || q > maxSub) continue;
        const int64_t v = eval_line(s, axis, q, fixedSub);
        if (v < bestValue) { bestValue = v; best = q; }
    }
    if (bestValue > 0) return false; // exact integer predicate; no floating sign decision
    int64_t lo = minSub, hi = best;
    while (lo < hi) {
        const int64_t mid = lo + (hi - lo) / 2;
        if (eval_line(s, axis, mid, fixedSub) <= 0) hi = mid;
        else lo = mid + 1;
    }
    const int64_t first = lo;
    lo = best; hi = maxSub;
    while (lo < hi) {
        const int64_t mid = lo + (hi - lo + 1) / 2;
        if (eval_line(s, axis, mid, fixedSub) <= 0) lo = mid;
        else hi = mid - 1;
    }
    const int64_t last = lo;
    if (eval_line(s, axis, first, fixedSub) > 0 || eval_line(s, axis, last, fixedSub) > 0) return false;
    const int64_t baseSub = static_cast<int64_t>(kGridMinCm) * kSubPerCm;
    const int64_t stepSub = static_cast<int64_t>(kCellStepCm) * kSubPerCm;
    int firstSample = static_cast<int>(ceil_div(static_cast<int>(first - baseSub), static_cast<int>(stepSub)));
    int lastSample = static_cast<int>(floor_div(static_cast<int>(last - baseSub), static_cast<int>(stepSub)));
    firstSample = std::max(firstSample, 0);
    lastSample = std::min(lastSample, kNodes - 1);
    if (firstSample <= lastSample) {
        const int64_t firstP = baseSub + firstSample * stepSub;
        if (eval_line(s, axis, firstP, fixedSub) > 0) ++firstSample;
        if (firstSample <= lastSample && eval_line(s, axis, baseSub + lastSample * stepSub, fixedSub) > 0) --lastSample;
    } else { firstSample = 1; lastSample = 0; }
    out = {first, last, firstSample, lastSample};
    return true;
}

uint64_t compose_membership(const Corpus& c, const std::array<uint32_t, kShapeLimitPerLine>& ids,
                            uint32_t count, uint64_t mask) {
    bool rock = true;
    while (mask) {
        const unsigned bit = std::countr_zero(mask);
        const uint64_t b = uint64_t{1} << bit;
        if (bit < count) rock = (c.shapes[ids[bit]].op == Operation::AddRock);
        mask &= ~b;
    }
    return rock ? 1 : 0;
}

void set_word_run(uint64_t* words, int begin, int end, bool rock) {
    if (!rock || begin >= end) return;
    int firstWord = begin >> 6;
    int lastWord = (end - 1) >> 6;
    const int firstBit = begin & 63;
    const int lastBit = (end - 1) & 63;
    if (firstWord == lastWord) {
        const uint64_t low = kMask64 << firstBit;
        const uint64_t high = (lastBit == 63) ? kMask64 : ((uint64_t{1} << (lastBit + 1)) - 1);
        words[firstWord] |= low & high;
        return;
    }
    words[firstWord] |= kMask64 << firstBit;
    for (int w = firstWord + 1; w < lastWord; ++w) words[w] = kMask64;
    words[lastWord] |= (lastBit == 63) ? kMask64 : ((uint64_t{1} << (lastBit + 1)) - 1);
}

void parallel_ranges(size_t count, int workers, const std::function<void(int,size_t,size_t)>& fn) {
    std::vector<std::thread> pool;
    pool.reserve(workers > 1 ? workers : 0);
    const size_t chunk = (count + static_cast<size_t>(workers) - 1) / static_cast<size_t>(workers);
    for (int w = 0; w < workers; ++w) {
        const size_t begin = std::min(count, static_cast<size_t>(w) * chunk);
        const size_t end = std::min(count, begin + chunk);
        if (workers == 1) { fn(w, begin, end); break; }
        pool.emplace_back([&, w, begin, end] { fn(w, begin, end); });
    }
    for (auto& t : pool) t.join();
}

uint32_t line_id(int axis, int x, int y, int z) {
    if (axis == 0) return static_cast<uint32_t>(z * kNodes + y); // X varies
    if (axis == 1) return static_cast<uint32_t>(x * kNodes + z); // Y varies
    return static_cast<uint32_t>(y * kNodes + x);                // Z varies
}

void decode_line(int axis, uint32_t line, int& a, int& b) {
    (void)axis;
    a = static_cast<int>(line % kNodes);
    b = static_cast<int>(line / kNodes);
    // axis 0: a=y,b=z; axis 1: a=x,b=z; axis 2: a=x,b=y.
}

AxisData solve_axis(const Corpus& c, int axis, int workers, AxisStats& stats, bool& overflow) {
    AxisData result;
    const int lineCount = kNodes * kNodes;
    result.signs.assign(static_cast<size_t>(lineCount) * kWordsPerLine, 0);
    result.lines.resize(lineCount);
    std::vector<WorkerLines> local(static_cast<size_t>(workers));
    std::atomic<bool> lineOverflow{false};
    std::vector<std::vector<uint16_t>> candidateStats(workers), intervalStats(workers), crossingStats(workers), runStats(workers);
    const auto& bins = c.projections[axis];
    const int axis1 = (axis + 1) % 3, axis2 = (axis + 2) % 3;
    const int64_t domainMinSub = static_cast<int64_t>(kGridMinCm) * kSubPerCm;
    const int64_t domainMaxSub = static_cast<int64_t>(kGridMaxCm) * kSubPerCm;

    parallel_ranges(lineCount, workers, [&](int worker, size_t begin, size_t end) {
        auto& outCross = local[worker].crossings;
        auto& cstats = candidateStats[worker];
        auto& istats = intervalStats[worker];
        auto& xstats = crossingStats[worker];
        auto& rstats = runStats[worker];
        cstats.reserve(end - begin); istats.reserve(end - begin);
        xstats.reserve(end - begin); rstats.reserve(end - begin);
        std::array<uint32_t, kShapeLimitPerLine> ids{};
        std::array<Interval, kShapeLimitPerLine> intervals{};
        std::array<uint8_t, kShapeLimitPerLine> valid{};
        std::array<uint64_t, kNodes + 1> enters{}, exits{};
        std::array<uint64_t, kWordsPerLine> lineWords{};
        struct Ev { int64_t q; uint8_t bit; uint8_t enter; };
        std::array<Ev, kShapeLimitPerLine * 2> events{};
        for (size_t lineIndex = begin; lineIndex < end; ++lineIndex) {
            int p0, p1; decode_line(axis, static_cast<uint32_t>(lineIndex), p0, p1);
            std::array<int, 3> localCm{};
            localCm[axis1] = kGridMinCm + p0 * kCellStepCm;
            localCm[axis2] = kGridMinCm + p1 * kCellStepCm;
            const int b0 = std::clamp(floor_div(localCm[axis1] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
            const int b1 = std::clamp(floor_div(localCm[axis2] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
            const auto& bucket = bins.ids[b1 * kBucketCount + b0];
            uint32_t candidateCount = 0;
            for (uint32_t id : bucket) {
                const Shape& s = c.shapes[id];
                if (localCm[axis1] < s.lo[axis1] || localCm[axis1] > s.hi[axis1]
                    || localCm[axis2] < s.lo[axis2] || localCm[axis2] > s.hi[axis2]
                    || s.hi[axis] < kGridMinCm || s.lo[axis] > kGridMaxCm) continue;
            if (candidateCount < kShapeLimitPerLine) ids[candidateCount++] = id;
                else { lineOverflow.store(true, std::memory_order_relaxed); }
            }
            if (candidateCount > kShapeLimitPerLine) candidateCount = kShapeLimitPerLine;
            cstats.push_back(static_cast<uint16_t>(candidateCount));
            std::fill(enters.begin(), enters.end(), 0);
            std::fill(exits.begin(), exits.end(), 0);
            std::fill(lineWords.begin(), lineWords.end(), 0);
            std::fill(valid.begin(), valid.end(), 0);
            uint32_t intervalCount = 0;
            uint32_t eventCount = 0;
            int64_t fixedSub[3] = {0,0,0};
            fixedSub[axis1] = static_cast<int64_t>(localCm[axis1]) * kSubPerCm;
            fixedSub[axis2] = static_cast<int64_t>(localCm[axis2]) * kSubPerCm;
            for (uint32_t i = 0; i < candidateCount; ++i) {
                const Shape& s = c.shapes[ids[i]];
                if (!solve_shape_line(s, axis, fixedSub, intervals[i])) continue;
                valid[i] = 1;
                ++intervalCount;
                const Interval& in = intervals[i];
                if (in.loSample <= in.hiSample) {
                    enters[in.loSample] |= uint64_t{1} << i;
                    exits[in.hiSample + 1] |= uint64_t{1} << i;
                }
                if (in.loSub <= domainMinSub && in.hiSub >= domainMinSub) {
                    // The line begins inside this canonical interval; represented by an entry at index 0.
                    enters[0] |= uint64_t{1} << i;
                }
                if (in.loSub > domainMinSub && in.loSub <= domainMaxSub)
                    events[eventCount++] = {in.loSub, static_cast<uint8_t>(i), 1};
                const int64_t exitQ = in.hiSub + 1;
                if (exitQ > domainMinSub && exitQ <= domainMaxSub + 1)
                    events[eventCount++] = {exitQ, static_cast<uint8_t>(i), 0};
            }
            istats.push_back(static_cast<uint16_t>(intervalCount));

            uint64_t active = 0;
            // Membership at the first sample is exact; interval endpoints are canonical integers.
            for (uint32_t i = 0; i < candidateCount; ++i) {
                if (!valid[i]) continue;
                const Interval& in = intervals[i];
                if (in.loSub <= domainMinSub && in.hiSub >= domainMinSub)
                    active |= uint64_t{1} << i;
            }
            bool rock = compose_membership(c, ids, candidateCount, active) != 0;
            int runBegin = 0;
            uint32_t runCount = 1;
            for (int k = 0; k <= kNodes; ++k) {
                const uint64_t before = active;
                active |= enters[k];
                active &= ~exits[k];
                if (k == 0) {
                    rock = compose_membership(c, ids, candidateCount, active) != 0;
                    continue;
                }
                if (before != active) {
                    const bool nextRock = compose_membership(c, ids, candidateCount, active) != 0;
                    if (nextRock != rock) {
                        set_word_run(lineWords.data(), runBegin, k, rock);
                        runBegin = k;
                        rock = nextRock;
                        ++runCount;
                    }
                }
            }
            set_word_run(lineWords.data(), runBegin, kNodes, rock);
            std::copy(lineWords.begin(), lineWords.end(),
                      result.signs.begin() + static_cast<size_t>(lineIndex) * kWordsPerLine);
            rstats.push_back(static_cast<uint16_t>(runCount));

            std::sort(events.begin(), events.begin() + eventCount,
                      [](const Ev& a, const Ev& b) { return a.q < b.q || (a.q == b.q && a.bit < b.bit); });
            const size_t lineCrossStart = outCross.size();
            uint64_t continuousActive = 0;
            for (uint32_t i = 0; i < candidateCount; ++i) {
                if (!valid[i]) continue;
                const Interval& in = intervals[i];
                if (in.loSub <= domainMinSub && in.hiSub >= domainMinSub)
                    continuousActive |= uint64_t{1} << i;
            }
            for (uint32_t e = 0; e < eventCount;) {
                uint32_t next = e + 1;
                while (next < eventCount && events[next].q == events[e].q) ++next;
                const bool beforeRock = compose_membership(c, ids, candidateCount, continuousActive) != 0;
                for (uint32_t j = e; j < next; ++j) {
                    const uint64_t bit = uint64_t{1} << events[j].bit;
                    if (events[j].enter) continuousActive |= bit;
                    else continuousActive &= ~bit;
                }
                const bool afterRock = compose_membership(c, ids, candidateCount, continuousActive) != 0;
                if (beforeRock != afterRock) {
                    uint32_t feature = ids[events[e].bit];
                    for (uint32_t j = e; j < next; ++j) feature = std::min(feature, ids[events[j].bit]);
                    const int64_t p32 = events[e].q * 2 - 1;
                    if (p32 >= static_cast<int64_t>(kGridMinCm) * 32
                        && p32 <= static_cast<int64_t>(kGridMaxCm) * 32
                        && outCross.size() < std::numeric_limits<uint32_t>::max())
                        outCross.push_back({static_cast<int32_t>(p32), feature,
                                            static_cast<uint8_t>(beforeRock), static_cast<uint8_t>(afterRock)});
                }
                e = next;
            }
            const size_t lineCrossCount = outCross.size() - lineCrossStart;
            xstats.push_back(static_cast<uint16_t>(std::min<size_t>(lineCrossCount, 65535)));
            result.lines[lineIndex].candidates = static_cast<uint16_t>(candidateCount);
            result.lines[lineIndex].primitiveIntervals = static_cast<uint16_t>(intervalCount);
            result.lines[lineIndex].occupancyRuns = static_cast<uint16_t>(runCount);
            result.lines[lineIndex].count = static_cast<uint16_t>(std::min<size_t>(lineCrossCount, 65535));
        }
    });
    overflow = overflow || lineOverflow.load(std::memory_order_relaxed);

    // Concatenate disjoint worker event streams in line order; each line range gets stable offsets.
    size_t totalCross = 0;
    std::vector<size_t> workerBase(static_cast<size_t>(workers));
    for (int w = 0; w < workers; ++w) { workerBase[w] = totalCross; totalCross += local[w].crossings.size(); }
    result.crossings.resize(totalCross);
    const size_t chunk = (static_cast<size_t>(lineCount) + workers - 1) / workers;
    for (int w = 0; w < workers; ++w) {
        const size_t begin = std::min(static_cast<size_t>(lineCount), static_cast<size_t>(w) * chunk);
        const size_t end = std::min(static_cast<size_t>(lineCount), begin + chunk);
        size_t offset = workerBase[w];
        // During the worker pass, crossing records for a line are appended sequentially.
        for (size_t li = begin; li < end; ++li) {
            LineRange& lr = result.lines[li];
            lr.offset = static_cast<uint32_t>(offset);
            offset += lr.count;
        }
        std::copy(local[w].crossings.begin(), local[w].crossings.end(), result.crossings.begin() + workerBase[w]);
    }
    // The stats vectors were produced in worker-range order; collect and summarize them.
    auto summarize = [](std::vector<uint16_t> values, double& p50, double& p95, uint32_t& vmax, uint64_t& total) {
        if (values.empty()) return;
        std::sort(values.begin(), values.end());
        total = std::accumulate(values.begin(), values.end(), uint64_t{0});
        vmax = values.back();
        const size_t i50 = static_cast<size_t>(0.50 * (values.size() - 1));
        const size_t i95 = static_cast<size_t>(0.95 * (values.size() - 1));
        p50 = values[i50]; p95 = values[i95];
    };
    std::vector<uint16_t> cv, iv, xv, rv;
    cv.reserve(lineCount); iv.reserve(lineCount); xv.reserve(lineCount); rv.reserve(lineCount);
    for (const auto& v : candidateStats) cv.insert(cv.end(), v.begin(), v.end());
    for (const auto& v : intervalStats) iv.insert(iv.end(), v.begin(), v.end());
    for (const auto& v : crossingStats) xv.insert(xv.end(), v.begin(), v.end());
    for (const auto& v : runStats) rv.insert(rv.end(), v.begin(), v.end());
    summarize(std::move(cv), stats.candidateP50, stats.candidateP95, stats.candidateMax, stats.candidateTotal);
    summarize(std::move(iv), stats.intervalP50, stats.intervalP95, stats.intervalMax, stats.intervalTotal);
    summarize(std::move(xv), stats.crossingP50, stats.crossingP95, stats.crossingMax, stats.crossingTotal);
    summarize(std::move(rv), stats.runsP50, stats.runsP95, stats.runsMax, stats.runsTotal);
    stats.intervalMax = std::max<uint32_t>(stats.intervalMax, 0);
    return result;
}

uint64_t get_sign(const AxisData& axis, uint32_t line, int sample) {
    return (axis.signs[static_cast<size_t>(line) * kWordsPerLine + (sample >> 6)] >> (sample & 63)) & 1;
}

uint64_t extract64(const AxisData& a, uint32_t line, int bit) {
    const uint64_t* words = a.signs.data() + static_cast<size_t>(line) * kWordsPerLine;
    const int wi = bit >> 6, shift = bit & 63;
    uint64_t v = words[wi] >> shift;
    if (shift != 0 && wi + 1 < kWordsPerLine) v |= words[wi + 1] << (64 - shift);
    return v;
}

uint64_t index_hash_bytes(uint64_t h, const void* data, size_t bytes) {
    const auto* p = static_cast<const uint8_t*>(data);
    for (size_t i = 0; i < bytes; ++i) { h ^= p[i]; h *= 1099511628211ull; }
    return h;
}

uint32_t owner_page_for_edge(int axis, int x, int y, int z) {
    int n[3] = {x,y,z};
    int cell[3]{};
    for (int a = 0; a < 3; ++a) {
        int ci;
        if (a == axis) ci = n[a] - 1;
        else {
            const int coreNode = n[a] - kHalo;
            if (coreNode > 0 && coreNode < kCoreCells && coreNode % kPageCells == 0) ci = coreNode - 1;
            else ci = coreNode;
        }
        ci = std::clamp(ci, 0, kCoreCells - 1);
        cell[a] = ci / kPageCells;
    }
    return static_cast<uint32_t>((cell[2] * kPagesPerAxis + cell[1]) * kPagesPerAxis + cell[0]);
}

uint64_t edge_key(int axis, int x, int y, int z) {
    return (((static_cast<uint64_t>(axis) * kNodes + static_cast<uint64_t>(x)) * kNodes
             + static_cast<uint64_t>(y)) * kNodes + static_cast<uint64_t>(z));
}

struct EdgeDesc { int axis, x, y, z, c0, c1; };
constexpr std::array<EdgeDesc, 12> kEdges = {{
    {0,0,0,0,0,1},{1,1,0,0,1,2},{0,0,1,0,2,3},{1,0,0,0,3,0},
    {0,0,0,1,4,5},{1,1,0,1,5,6},{0,0,1,1,6,7},{1,0,0,1,7,4},
    {2,0,0,0,0,4},{2,1,0,0,1,5},{2,1,1,0,2,6},{2,0,1,0,3,7}
}};
constexpr int kCornerOffset[8][3] = {{0,0,0},{1,0,0},{1,1,0},{0,1,0},{0,0,1},{1,0,1},{1,1,1},{0,1,1}};
constexpr int kFaceCorners[6][4] = {{0,1,2,3},{4,5,6,7},{0,1,5,4},{3,2,6,7},{0,3,7,4},{1,2,6,5}};
constexpr int kFaceEdges[6][4] = {{0,1,2,3},{4,5,6,7},{0,9,4,8},{2,10,6,11},{3,11,7,8},{1,10,5,9}};

struct CrossingRef { Crossing c; uint32_t line; };

bool find_edge_crossing(const AxisData& axisData, int axis, int x, int y, int z,
                        bool leftRock, bool rightRock, Crossing& out, uint64_t& multi) {
    const uint32_t line = line_id(axis, x, y, z);
    const LineRange& lr = axisData.lines[line];
    const int node = (axis == 0 ? x : axis == 1 ? y : z);
    const int localCm = kGridMinCm + node * kCellStepCm;
    const int64_t start32 = static_cast<int64_t>(localCm) * kSurfaceSubPerCm;
    const int64_t end32 = start32 + kCellStepCm * kSurfaceSubPerCm;
    uint32_t found = 0;
    for (uint32_t i = 0; i < lr.count; ++i) {
        const Crossing& c = axisData.crossings[lr.offset + i];
        if (c.positionSub32 <= start32 || c.positionSub32 >= end32) continue;
        if (c.beforeRock == leftRock && c.afterRock == rightRock) {
            out = c; ++found;
        }
    }
    if (found > 1) ++multi;
    return found != 0;
}

void normal_for_crossing(const Shape& s, const Crossing& c, int axis, int local[3],
                         double posCm, double outN[3]) {
    double d[3] = {static_cast<double>(local[0] - s.c[0]),
                   static_cast<double>(local[1] - s.c[1]),
                   static_cast<double>(local[2] - s.c[2])};
    d[axis] = posCm - s.c[axis];
    double g[3]{};
    for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) g[i] += 2.0 * static_cast<double>(s.q[i][j]) * d[j];
    const bool wantPositiveAxis = c.beforeRock && !c.afterRock;
    const double sign = ((g[axis] >= 0.0) == wantPositiveAxis) ? 1.0 : -1.0;
    double len2 = 0.0;
    for (int i = 0; i < 3; ++i) { outN[i] = g[i] * sign; len2 += outN[i] * outN[i]; }
    if (len2 <= 1.0e-28) { outN[0] = (axis == 0 ? 1.0 : 0.0); outN[1] = (axis == 1 ? 1.0 : 0.0); outN[2] = (axis == 2 ? 1.0 : 0.0); return; }
    const double inv = 1.0 / std::sqrt(len2);
    for (double& v : std::span<double,3>(outN,3)) v *= inv;
}

int axis_for_edge(const EdgeDesc& e, int x, int y, int z) { return e.axis; }

void write_vertex_hash(uint64_t& hash, const Vertex& v) { hash = index_hash_bytes(hash, &v, sizeof(v)); }

MeshOutput emit_near_mesh(const Corpus& c, const AxisData axes[3], const AxisData& zSigns,
                          const std::vector<uint64_t>& active, int workers) {
    (void)workers;
    MeshOutput out;
    std::unordered_map<uint64_t, uint32_t> edgeToVertex;
    edgeToVertex.reserve(700000);
    out.pages.reserve(kPagesPerAxis * kPagesPerAxis * kPagesPerAxis);
    std::array<uint64_t, 8> cornerSigns{};
    std::array<int, 3> local{};
    uint64_t multiEdges = 0;
    auto add_vertex = [&](const EdgeDesc& e, int cellX, int cellY, int cellZ,
                          const Crossing& crossing) -> uint32_t {
        const int ex = cellX + e.x, ey = cellY + e.y, ez = cellZ + e.z;
        const uint64_t key = edge_key(e.axis, ex, ey, ez);
        auto it = edgeToVertex.find(key);
        if (it != edgeToVertex.end()) return it->second;
        const Shape& s = c.shapes[crossing.shapeId];
        int p[3] = {kGridMinCm + ex * kCellStepCm,
                    kGridMinCm + ey * kCellStepCm,
                    kGridMinCm + ez * kCellStepCm};
        const double localPosCm = crossing.positionSub32 / static_cast<double>(kSurfaceSubPerCm);
        p[e.axis] = static_cast<int>(std::floor(localPosCm));
        double n[3]; normal_for_crossing(s, crossing, e.axis, p, localPosCm, n);
        const double world[3] = {
            kWorldOriginCm[0] + (e.axis == 0 ? localPosCm : (kGridMinCm + ex * kCellStepCm)),
            kWorldOriginCm[1] + (e.axis == 1 ? localPosCm : (kGridMinCm + ey * kCellStepCm)),
            kWorldOriginCm[2] + (e.axis == 2 ? localPosCm : (kGridMinCm + ez * kCellStepCm))};
        Vertex v{static_cast<float>(world[0] * 0.01), static_cast<float>(world[1] * 0.01),
                 static_cast<float>(world[2] * 0.01), static_cast<float>(n[0]),
                 static_cast<float>(n[1]), static_cast<float>(n[2])};
        const uint32_t owner = owner_page_for_edge(e.axis, ex, ey, ez);
        const uint32_t index = static_cast<uint32_t>(out.vertices.size());
        out.vertices.push_back({key, owner, v});
        edgeToVertex.emplace(key, index);
        const int edgeNodes[3] = {ex,ey,ez};
        for (int a = 0; a < 3; ++a) {
            const int coreNode = edgeNodes[a] - kHalo;
            if (coreNode > 0 && coreNode < kCoreCells && coreNode % kPageCells == 0) {
                ++out.pageBorderVertices;
                break;
            }
        }
        return index;
    };
    auto emit_triangle = [&](uint32_t a, uint32_t b, uint32_t d) {
        const Vertex& va = out.vertices[a].v; const Vertex& vb = out.vertices[b].v; const Vertex& vd = out.vertices[d].v;
        const double ux = vb.x - va.x, uy = vb.y - va.y, uz = vb.z - va.z;
        const double vx = vd.x - va.x, vy = vd.y - va.y, vz = vd.z - va.z;
        const double cx = uy * vz - uz * vy, cy = uz * vx - ux * vz, cz = ux * vy - uy * vx;
        const double nx = va.nx + vb.nx + vd.nx, ny = va.ny + vb.ny + vd.ny, nz = va.nz + vb.nz + vd.nz;
        if (cx * nx + cy * ny + cz * nz < 0.0) std::swap(b, d);
        out.indices.push_back(a); out.indices.push_back(b); out.indices.push_back(d);
    };
    auto shape_hits_lattice = [&](const Shape& s) {
        int lo[3], hi[3];
        for (int a = 0; a < 3; ++a) {
            lo[a] = std::clamp(ceil_div(s.lo[a] - kGridMinCm, kCellStepCm), 0, kNodes - 1);
            hi[a] = std::clamp(floor_div(s.hi[a] - kGridMinCm, kCellStepCm), 0, kNodes - 1);
        }
        for (int z = lo[2]; z <= hi[2]; ++z) for (int y = lo[1]; y <= hi[1]; ++y)
            for (int x = lo[0]; x <= hi[0]; ++x) {
                const int64_t p[3] = {
                    static_cast<int64_t>(kGridMinCm + x * kCellStepCm) * kSubPerCm,
                    static_cast<int64_t>(kGridMinCm + y * kCellStepCm) * kSubPerCm,
                    static_cast<int64_t>(kGridMinCm + z * kCellStepCm) * kSubPerCm};
                if (contains_sub(s, p)) return true;
            }
        return false;
    };
    for (int pz = 0; pz < kPagesPerAxis; ++pz) for (int py = 0; py < kPagesPerAxis; ++py)
    for (int px = 0; px < kPagesPerAxis; ++px) {
        const uint32_t page = static_cast<uint32_t>((pz * kPagesPerAxis + py) * kPagesPerAxis + px);
        const uint32_t start = static_cast<uint32_t>(out.indices.size());
        const size_t vstart = out.vertices.size();
        for (int iy = py * kPageCells; iy < (py + 1) * kPageCells; ++iy)
        for (int ix = px * kPageCells; ix < (px + 1) * kPageCells; ++ix) {
            const int nx = ix + kHalo, ny = iy + kHalo;
            for (int zb = pz; zb <= pz; ++zb) {
                uint64_t mask = active[(static_cast<size_t>(ix) * kCoreCells + iy) * 4 + zb];
                while (mask) {
                    const unsigned bz = std::countr_zero(mask);
                    mask &= mask - 1;
                    const int iz = zb * 64 + static_cast<int>(bz);
                    const int cellX = nx, cellY = ny, cellZ = iz + kHalo;
                    bool rock[8]{};
                    for (int k = 0; k < 8; ++k) {
                        const int x = cellX + kCornerOffset[k][0];
                        const int y = cellY + kCornerOffset[k][1];
                        const int z = cellZ + kCornerOffset[k][2];
                        rock[k] = get_sign(zSigns, line_id(2, x, y, z), z);
                    }
                    bool edgeCross[12]{};
                    uint32_t edgeVertex[12]{};
                    for (int ei = 0; ei < 12; ++ei) {
                        const EdgeDesc& e = kEdges[ei];
                        if (rock[e.c0] == rock[e.c1]) continue;
                        const int sx = cellX + e.x, sy = cellY + e.y, sz = cellZ + e.z;
                        const bool c0IsLow = kCornerOffset[e.c0][e.axis] < kCornerOffset[e.c1][e.axis];
                        const bool lowRock = c0IsLow ? rock[e.c0] : rock[e.c1];
                        const bool highRock = c0IsLow ? rock[e.c1] : rock[e.c0];
                        Crossing crossing;
                        if (!find_edge_crossing(axes[e.axis], e.axis, sx, sy, sz,
                                                lowRock, highRock, crossing, multiEdges)) {
                            ++out.missingAnalyticCrossings;
                            continue;
                        }
                        edgeCross[ei] = true;
                        edgeVertex[ei] = add_vertex(e, cellX, cellY, cellZ, crossing);
                    }
                    std::array<std::array<int,2>,12> adj{};
                    std::array<uint8_t,12> degree{};
                    auto connect = [&](int a, int b) {
                        if (degree[a] < 2 && degree[b] < 2) {
                            adj[a][degree[a]++] = b; adj[b][degree[b]++] = a;
                        }
                    };
                    for (int f = 0; f < 6; ++f) {
                        int crossed[4], nCross = 0;
                        for (int k = 0; k < 4; ++k) if (edgeCross[kFaceEdges[f][k]]) crossed[nCross++] = k;
                        if (nCross == 2) connect(kFaceEdges[f][crossed[0]], kFaceEdges[f][crossed[1]]);
                        else if (nCross == 4) {
                            ++out.ambiguousFaces;
                            int pSub[3] = {kGridMinCm + cellX * kCellStepCm,
                                           kGridMinCm + cellY * kCellStepCm,
                                           kGridMinCm + cellZ * kCellStepCm};
                            int sumOffset[3]{};
                            for (int k = 0; k < 4; ++k) for (int a = 0; a < 3; ++a)
                                sumOffset[a] += kCornerOffset[kFaceCorners[f][k]][a];
                            int64_t q[3];
                            for (int a = 0; a < 3; ++a) q[a] = (static_cast<int64_t>(pSub[a]) * 4 + sumOffset[a]) * 4;
                            const bool centerRock = point_rock(c, q);
                            const bool s0 = rock[kFaceCorners[f][0]];
                            if (centerRock == s0) {
                                connect(kFaceEdges[f][0], kFaceEdges[f][1]);
                                connect(kFaceEdges[f][2], kFaceEdges[f][3]);
                            } else {
                                connect(kFaceEdges[f][3], kFaceEdges[f][0]);
                                connect(kFaceEdges[f][1], kFaceEdges[f][2]);
                            }
                        }
                    }
                    for (int ei = 0; ei < 12; ++ei)
                        if (edgeCross[ei] && degree[ei] != 2) ++out.contourDegreeFaults;
                    std::array<uint8_t,12> visited{};
                    for (int startEdge = 0; startEdge < 12; ++startEdge) {
                        if (!edgeCross[startEdge] || visited[startEdge] || degree[startEdge] != 2) continue;
                        std::vector<int> loop;
                        int prev = -1, cur = startEdge;
                        for (int guard = 0; guard < 16; ++guard) {
                            if (visited[cur]) break;
                            visited[cur] = 1; loop.push_back(cur);
                            const int next = adj[cur][0] == prev ? adj[cur][1] : adj[cur][0];
                            prev = cur; cur = next;
                            if (cur == startEdge) break;
                        }
                        if (cur != startEdge) ++out.openContourLoops;
                        if (loop.size() >= 3) for (size_t i = 1; i + 1 < loop.size(); ++i)
                            emit_triangle(edgeVertex[loop[0]], edgeVertex[loop[i]], edgeVertex[loop[i + 1]]);
                    }
                }
            }
        }
        out.pages.push_back({page, start, static_cast<uint32_t>(out.indices.size() - start),
                             static_cast<uint32_t>(vstart), static_cast<uint32_t>(out.vertices.size() - vstart)});
    }

    // Explicit surface fallback for sub-cell strata missed by every structural lattice point.
    constexpr int rings = 8, sectors = 24;
    for (const Shape& s : c.shapes) {
        if (s.kind != FeatureKind::Stratum || s.radius[2] >= kCellStepCm / 2) continue;
        ++out.thinFeaturesDetected;
        if (shape_hits_lattice(s)) continue;
        ++out.thinFeaturesHiddenInLattice;
        const uint32_t page = owner_page_for_edge(0, std::clamp(kHalo + (s.c[0] - kGridMinCm) / kCellStepCm, 1, kCoreCells),
            std::clamp(kHalo + (s.c[1] - kGridMinCm) / kCellStepCm, 1, kCoreCells),
            std::clamp(kHalo + (s.c[2] - kGridMinCm) / kCellStepCm, 1, kCoreCells));
        const uint32_t start = static_cast<uint32_t>(out.indices.size());
        const uint32_t base = static_cast<uint32_t>(out.vertices.size());
        std::vector<std::array<double, 3>> localPoints;
        std::vector<std::array<double, 3>> surfaceNormals;
        localPoints.reserve((rings + 1) * sectors);
        surfaceNormals.reserve((rings + 1) * sectors);
        for (int r = 0; r <= rings; ++r) {
            const double phi = -1.5707963267948966 + 3.1415926535897932 * r / rings;
            for (int a = 0; a < sectors; ++a) {
                const double theta = 6.2831853071795864 * a / sectors;
                const double u[3] = {std::cos(phi) * std::cos(theta), std::cos(phi) * std::sin(theta), std::sin(phi)};
                double localP[3]{};
                for (int axis = 0; axis < 3; ++axis) {
                    for (int k = 0; k < 3; ++k)
                        localP[axis] += static_cast<double>(s.rot[axis][k]) / 32768.0 * s.radius[k] * u[k];
                }
                double n[3]{};
                for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j)
                    n[i] += 2.0 * s.q[i][j] * localP[j];
                double len = std::sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
                if (len > 0) for (double& v : n) v /= len;
                if (s.op == Operation::CarveAir) for (double& v : n) v = -v;
                localPoints.push_back({localP[0], localP[1], localP[2]});
                surfaceNormals.push_back({n[0],n[1],n[2]});
                const double wx = (kWorldOriginCm[0] + s.c[0] + localP[0]) * 0.01;
                const double wy = (kWorldOriginCm[1] + s.c[1] + localP[1]) * 0.01;
                const double wz = (kWorldOriginCm[2] + s.c[2] + localP[2]) * 0.01;
                const uint64_t thinKey = (uint64_t{1} << 63) | (static_cast<uint64_t>(s.id) << 32)
                    | (static_cast<uint64_t>(r) << 16) | static_cast<uint64_t>(a);
                out.vertices.push_back({thinKey, page,
                    {static_cast<float>(wx),static_cast<float>(wy),static_cast<float>(wz),
                     static_cast<float>(n[0]),static_cast<float>(n[1]),static_cast<float>(n[2])}});
            }
        }
        for (int r = 0; r < rings; ++r) for (int a = 0; a < sectors; ++a) {
            const int an = (a + 1) % sectors;
            const uint32_t v00 = base + r * sectors + a, v01 = base + r * sectors + an;
            const uint32_t v10 = base + (r + 1) * sectors + a, v11 = base + (r + 1) * sectors + an;
            const uint32_t tris[2][3] = {{v00,v10,v01},{v01,v10,v11}};
            for (int ti = 0; ti < 2; ++ti) {
                if ((r == 0 && ti == 0) || (r == rings - 1 && ti == 1)) continue;
                double p[3]{}, n[3]{};
                for (int v : tris[ti]) for (int a3 = 0; a3 < 3; ++a3) {
                    p[a3] += (s.c[a3] + localPoints[v - base][a3]) / 3.0;
                    n[a3] += surfaceNormals[v - base][a3] / 3.0;
                }
                const double nlen = std::sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
                if (nlen > 0) for (double& v : n) v /= nlen;
                int64_t plus[3], minus[3];
                for (int a3 = 0; a3 < 3; ++a3) {
                    plus[a3] = std::llround((p[a3] + n[a3]) * kSubPerCm);
                    minus[a3] = std::llround((p[a3] - n[a3]) * kSubPerCm);
                }
                if (point_rock(c, plus) == point_rock(c, minus)) continue;
                emit_triangle(tris[ti][0], tris[ti][1], tris[ti][2]);
                ++out.explicitThinTriangles;
            }
        }
        if (out.indices.size() != start) {
            ++out.thinFeaturesMeshed;
            out.pages.push_back({page, start, static_cast<uint32_t>(out.indices.size() - start), base,
                                 static_cast<uint32_t>(out.vertices.size() - base)});
        } else {
            out.vertices.resize(base);
        }
    }

    // Stable ownership and vertex order by (owner page, edge key); index order remains page-owned.
    std::vector<uint32_t> permutation(out.vertices.size());
    std::iota(permutation.begin(), permutation.end(), 0u);
    std::sort(permutation.begin(), permutation.end(), [&](uint32_t a, uint32_t b) {
        if (out.vertices[a].ownerPage != out.vertices[b].ownerPage)
            return out.vertices[a].ownerPage < out.vertices[b].ownerPage;
        return out.vertices[a].key < out.vertices[b].key;
    });
    std::vector<uint32_t> remap(out.vertices.size());
    std::vector<TempVertex> sortedVertices;
    sortedVertices.reserve(out.vertices.size());
    for (uint32_t i = 0; i < permutation.size(); ++i) {
        remap[permutation[i]] = i;
        sortedVertices.push_back(out.vertices[permutation[i]]);
    }
    for (uint32_t& index : out.indices) index = remap[index];
    out.vertices.swap(sortedVertices);
    for (PageRange& p : out.pages) {
        const auto first = std::lower_bound(out.vertices.begin(), out.vertices.end(), p.page,
            [](const TempVertex& v, uint32_t owner) { return v.ownerPage < owner; });
        const auto last = std::upper_bound(first, out.vertices.end(), p.page,
            [](uint32_t owner, const TempVertex& v) { return owner < v.ownerPage; });
        p.vertexStart = static_cast<uint32_t>(first - out.vertices.begin());
        p.vertexCount = static_cast<uint32_t>(last - first);
    }
    out.triangles = out.indices.size() / 3;
    out.multiCrossEdges = multiEdges;
    uint64_t h = 1469598103934665603ull;
    for (const TempVertex& v : out.vertices) write_vertex_hash(h, v.v);
    h = index_hash_bytes(h, out.indices.data(), out.indices.size() * sizeof(uint32_t));
    for (const PageRange& p : out.pages) h = index_hash_bytes(h, &p, sizeof(p));
    out.hash = h;
    out.bytes = out.vertices.size() * sizeof(Vertex) + out.indices.size() * sizeof(uint32_t)
        + out.pages.size() * sizeof(PageRange);
    return out;
}

uint64_t active_cell_masks(const AxisData& zAxis, std::vector<uint64_t>& active) {
    const size_t count = static_cast<size_t>(kCoreCells) * kCoreCells * 4;
    active.assign(count, 0);
    uint64_t activeCount = 0;
    for (int x = 0; x < kCoreCells; ++x) for (int y = 0; y < kCoreCells; ++y) {
        const int nx = x + kHalo, ny = y + kHalo;
        const uint32_t line00 = line_id(2, nx, ny, 0), line10 = line_id(2, nx + 1, ny, 0);
        const uint32_t line01 = line_id(2, nx, ny + 1, 0), line11 = line_id(2, nx + 1, ny + 1, 0);
        for (int block = 0; block < 4; ++block) {
            const int bit = 1 + block * 64;
            const uint64_t a0 = extract64(zAxis, line00, bit), a1 = extract64(zAxis, line10, bit);
            const uint64_t a2 = extract64(zAxis, line01, bit), a3 = extract64(zAxis, line11, bit);
            const uint64_t b0 = extract64(zAxis, line00, bit + 1), b1 = extract64(zAxis, line10, bit + 1);
            const uint64_t b2 = extract64(zAxis, line01, bit + 1), b3 = extract64(zAxis, line11, bit + 1);
            const uint64_t any = a0 | a1 | a2 | a3 | b0 | b1 | b2 | b3;
            const uint64_t all = a0 & a1 & a2 & a3 & b0 & b1 & b2 & b3;
            const uint64_t mask = any & ~all;
            active[(static_cast<size_t>(x) * kCoreCells + y) * 4 + block] = mask;
            activeCount += std::popcount(mask);
        }
    }
    return activeCount;
}

bool validate_signs(const Corpus& c, const AxisData axes[3], uint64_t& mismatches,
                    uint64_t& auditedLatticeRock, uint64_t& thinHiddenCorners,
                    uint64_t& thinCrossingEvents) {
    mismatches = auditedLatticeRock = thinHiddenCorners = thinCrossingEvents = 0;
    for (int z = 0; z < kNodes; ++z) for (int y = 0; y < kNodes; ++y) for (int x = 0; x < kNodes; ++x) {
        const int64_t p[3] = {
            static_cast<int64_t>(kGridMinCm + x * kCellStepCm) * kSubPerCm,
            static_cast<int64_t>(kGridMinCm + y * kCellStepCm) * kSubPerCm,
            static_cast<int64_t>(kGridMinCm + z * kCellStepCm) * kSubPerCm};
        const bool ref = point_rock(c, p);
        const bool sx = get_sign(axes[0], line_id(0, x, y, z), x) != 0;
        const bool sy = get_sign(axes[1], line_id(1, x, y, z), y) != 0;
        const bool sz = get_sign(axes[2], line_id(2, x, y, z), z) != 0;
        if (ref != sx || ref != sy || ref != sz) ++mismatches;
    }
    for (const Shape& s : c.shapes) if (s.kind == FeatureKind::Stratum && s.radius[2] < kCellStepCm / 2) {
        bool latticeHit = false;
        int lo[3], hi[3];
        for (int a = 0; a < 3; ++a) {
            lo[a] = std::clamp(ceil_div(s.lo[a] - kGridMinCm, kCellStepCm), 0, kNodes - 1);
            hi[a] = std::clamp(floor_div(s.hi[a] - kGridMinCm, kCellStepCm), 0, kNodes - 1);
        }
        for (int z = lo[2]; z <= hi[2] && !latticeHit; ++z) for (int y = lo[1]; y <= hi[1] && !latticeHit; ++y)
            for (int x = lo[0]; x <= hi[0] && !latticeHit; ++x) {
                const int64_t p[3] = {
                    static_cast<int64_t>(kGridMinCm + x * kCellStepCm) * kSubPerCm,
                    static_cast<int64_t>(kGridMinCm + y * kCellStepCm) * kSubPerCm,
                    static_cast<int64_t>(kGridMinCm + z * kCellStepCm) * kSubPerCm};
                latticeHit = contains_sub(s, p);
            }
        if (!latticeHit) ++thinHiddenCorners;
        if (s.auditThin) {
            const int x = std::clamp(static_cast<int>(std::round((s.c[0] - kGridMinCm) / static_cast<double>(kCellStepCm))), 0, kNodes - 1);
            const int y = std::clamp(static_cast<int>(std::round((s.c[1] - kGridMinCm) / static_cast<double>(kCellStepCm))), 0, kNodes - 1);
            const uint32_t line = line_id(2, x, y, 0);
            const LineRange& lr = axes[2].lines[line];
            for (uint32_t i = 0; i < lr.count; ++i)
                if (axes[2].crossings[lr.offset+i].shapeId == s.id) ++thinCrossingEvents;
        }
    }
    return mismatches == 0 && thinHiddenCorners != 0 && thinCrossingEvents >= 2;
}

bool validate_tangency_and_coincident_edits() {
    Shape tangent;
    tangent.c[0]=0;tangent.c[1]=0;tangent.c[2]=0;
    tangent.radius[0]=128;tangent.radius[1]=128;tangent.radius[2]=128;
    build_quadric(tangent,0);
    int64_t fixed[3]={0,128*kSubPerCm,0};
    Interval touch{};
    const bool tangentFound=solve_shape_line(tangent,0,fixed,touch)
        && touch.loSub==0 && touch.hiSub==0;

    Corpus carveThenAdd, addThenCarve;
    add_shape(carveThenAdd,0,0,0,128,128,128,Operation::CarveAir,FeatureKind::Crater,0);
    add_shape(carveThenAdd,0,0,0,128,128,128,Operation::AddRock,FeatureKind::Pillar,0);
    add_shape(addThenCarve,0,0,0,128,128,128,Operation::AddRock,FeatureKind::Pillar,0);
    add_shape(addThenCarve,0,0,0,128,128,128,Operation::CarveAir,FeatureKind::Crater,0);
    Interval i0{},i1{};
    const int64_t centerSub[3]={0,0,0};
    const bool sameRoots=solve_shape_line(carveThenAdd.shapes[0],2,centerSub,i0)
        && solve_shape_line(carveThenAdd.shapes[1],2,centerSub,i1)
        && i0.loSub==i1.loSub && i0.hiSub==i1.hiSub;
    std::array<uint32_t,kShapeLimitPerLine> ids{};ids[0]=0;ids[1]=1;
    const bool lastAddWins=compose_membership(carveThenAdd,ids,2,3)==1;
    const bool lastCarveWins=compose_membership(addThenCarve,ids,2,3)==0;
    return tangentFound && sameRoots && lastAddWins && lastCarveWins;
}

bool validate_page_ownership(const MeshOutput& mesh) {
    if (mesh.pages.size() < static_cast<size_t>(kPagesPerAxis*kPagesPerAxis*kPagesPerAxis)
        || mesh.pageBorderVertices == 0) return false;
    std::vector<uint64_t> keys;
    keys.reserve(mesh.vertices.size());
    const uint64_t n3=static_cast<uint64_t>(kNodes)*kNodes*kNodes;
    for (const TempVertex& v : mesh.vertices) {
        keys.push_back(v.key);
        if ((v.key >> 63) != 0) continue;
        const uint32_t axis=static_cast<uint32_t>(v.key/n3);
        uint64_t rem=v.key%n3;
        const int z=static_cast<int>(rem%kNodes);rem/=kNodes;
        const int y=static_cast<int>(rem%kNodes);rem/=kNodes;
        const int x=static_cast<int>(rem);
        if (axis<3 && x<kNodes && y<kNodes && z<kNodes
            && v.ownerPage!=owner_page_for_edge(static_cast<int>(axis),x,y,z)) return false;
    }
    std::sort(keys.begin(),keys.end());
    return std::adjacent_find(keys.begin(),keys.end())==keys.end()
        && mesh.missingAnalyticCrossings==0 && mesh.contourDegreeFaults==0 && mesh.openContourLoops==0;
}

FarOutput build_far_world(uint64_t seed) {
    FarOutput out;
    constexpr int count = 97;
    constexpr int layers = 5;
    constexpr int stepM = 125;
    constexpr int qScale = 1 << 28;
    constexpr double pi = 3.14159265358979323846;
    struct Parent { int cx,cy,cz,rx,ry,rz; int64_t qx,qy,qz; };
    std::array<Parent,layers> parents{};
    Rng rng(seed ^ 0xF411D20261002ull);
    for (int l=0;l<layers;++l) {
        Parent p{};
        p.cx=(l-2)*180; p.cy=(2-l)*140; p.cz=(l-2)*520;
        p.rx=10000+rng.range(-100,100); p.ry=10000+rng.range(-100,100); p.rz=150+rng.range(-20,20);
        p.qx=round_div(qScale,static_cast<int64_t>(p.rx)*p.rx);
        p.qy=round_div(qScale,static_cast<int64_t>(p.ry)*p.ry);
        p.qz=round_div(qScale,static_cast<int64_t>(p.rz)*p.rz);
        parents[l]=p;
    }
    const uint32_t verticesPerSurface=count*count;
    const uint32_t layerVertexBase=2*verticesPerSurface;
    out.vertices.reserve(static_cast<size_t>(layers)*layerVertexBase);
    out.indices.reserve(static_cast<size_t>(layers)*(count-1)*(count-1)*12);
    out.patches.reserve(static_cast<size_t>(layers)*(count-1)*(count-1)*2);
    std::vector<uint8_t> included(count*count,1);
    for (int iy=0;iy<count;++iy) for (int ix=0;ix<count;++ix)
        if (ix==48 && iy==48) included[iy*count+ix]=0; // near 64 m cube owns the central far tile
    for (int l=0;l<layers;++l) for (int iy=0;iy<count;++iy) for (int ix=0;ix<count;++ix) {
        if (!included[iy*count+ix]) continue;
        const int r=std::max(std::abs(ix-48),std::abs(iy-48));
        const int lod=std::min(5,r<4?0:r<8?1:r<16?2:r<32?3:r<40?4:5);
        ++out.lodNodes[lod];
    }

    for (int l=0;l<layers;++l) {
        const Parent& p=parents[l];
        const uint32_t base=static_cast<uint32_t>(out.vertices.size());
        for (int side=0;side<2;++side) for (int iy=0;iy<count;++iy) for (int ix=0;ix<count;++ix) {
            const int x=(ix-48)*stepM, y=(iy-48)*stepM;
            const int dx=x-p.cx,dy=y-p.cy;
            const int64_t remaining=qScale-p.qx*static_cast<int64_t>(dx)*dx-p.qy*static_cast<int64_t>(dy)*dy;
            if (remaining<0) { ++out.nodes; continue; }
            const double zOffset=p.rz*std::sqrt(static_cast<double>(remaining)/qScale);
            const double zRaw=p.cz+(side==0?zOffset:-zOffset);
            const double z=static_cast<double>(std::llround(zRaw*3200.0))/3200.0;
            double nx=static_cast<double>(dx)/(p.rx*p.rx);
            double ny=static_cast<double>(dy)/(p.ry*p.ry);
            double nz=(z-p.cz)/(p.rz*p.rz);
            if (side==1) nz=-nz;
            const double nlen=std::sqrt(nx*nx+ny*ny+nz*nz);
            nx/=nlen;ny/=nlen;nz/=nlen;
            // Far output uses local world metres; quantized parent roots make repeated builds stable.
            const double worldX=(kWorldOriginCm[0]*0.01)+x;
            const double worldY=(kWorldOriginCm[1]*0.01)+y;
            const double worldZ=(kWorldOriginCm[2]*0.01)+z;
            out.vertices.push_back({static_cast<float>(worldX),static_cast<float>(worldY),static_cast<float>(worldZ),
                                    static_cast<float>(nx),static_cast<float>(ny),static_cast<float>(nz)});
        }
        // One selected parent shape supplies the top and bottom coarse patches; no child shapes exist.
        const uint32_t surfaceBase=base;
        for (int side=0;side<2;++side) {
            const uint32_t vertBase=surfaceBase+side*verticesPerSurface;
            for (int iy=0;iy<count-1;++iy) for (int ix=0;ix<count-1;++ix) {
                const int tile=iy*(count-1)+ix;
                if (!included[std::min(iy,count-1)*count+std::min(ix,count-1)]) continue;
                const uint32_t i0=vertBase+iy*count+ix, i1=i0+1, i2=i0+count, i3=i2+1;
                const uint32_t start=static_cast<uint32_t>(out.indices.size());
                if (side==0) { out.indices.insert(out.indices.end(),{i0,i2,i1,i1,i2,i3}); }
                else { out.indices.insert(out.indices.end(),{i0,i1,i2,i1,i3,i2}); }
                const uint32_t owner=static_cast<uint32_t>((l*2+side)*10000+tile);
                out.patches.push_back({owner,start,6,i0,4});
                ++out.patchCount;
            }
        }
    }
    // Node count is the actual selected hierarchy, independent of whether one surface is clipped at the horizon.
    out.nodes=static_cast<uint64_t>(count)*count*layers- layers;
    out.descendants=0;
    uint64_t h=1469598103934665603ull;
    if (!out.vertices.empty()) h=index_hash_bytes(h,out.vertices.data(),out.vertices.size()*sizeof(Vertex));
    if (!out.indices.empty()) h=index_hash_bytes(h,out.indices.data(),out.indices.size()*sizeof(uint32_t));
    if (!out.patches.empty()) h=index_hash_bytes(h,out.patches.data(),out.patches.size()*sizeof(PageRange));
    out.hash=h;
    out.bytes=out.vertices.size()*sizeof(Vertex)+out.indices.size()*sizeof(uint32_t)+out.patches.size()*sizeof(PageRange);
    (void)pi;
    return out;
}

uint64_t process_working_set() {
    PROCESS_MEMORY_COUNTERS pmc{}; pmc.cb=sizeof(pmc);
    if (GetProcessMemoryInfo(GetCurrentProcess(),&pmc,sizeof(pmc))) return pmc.WorkingSetSize;
    return 0;
}

void print_axis_stats(const char* name,const AxisStats& s) {
    std::printf("\"%s\":{\"lines\":%d,\"candidate_total\":%llu,\"candidate_p50\":%.0f,\"candidate_p95\":%.0f,\"candidate_max\":%u,\"interval_total\":%llu,\"interval_p50\":%.0f,\"interval_p95\":%.0f,\"interval_max\":%u,\"crossing_total\":%llu,\"crossing_p50\":%.0f,\"crossing_p95\":%.0f,\"crossing_max\":%u,\"runs_p50\":%.0f,\"runs_p95\":%.0f,\"runs_max\":%u}",
        name,kNodes*kNodes,static_cast<unsigned long long>(s.candidateTotal),s.candidateP50,s.candidateP95,s.candidateMax,
        static_cast<unsigned long long>(s.intervalTotal),s.intervalP50,s.intervalP95,s.intervalMax,
        static_cast<unsigned long long>(s.crossingTotal),s.crossingP50,s.crossingP95,s.crossingMax,
        s.runsP50,s.runsP95,s.runsMax);
}

} // namespace

int main(int argc,char** argv) {
    const uint64_t seed=argc>1?std::strtoull(argv[1],nullptr,10):1;
    const int workers=argc>2?std::clamp(std::atoi(argv[2]),1,8):1;
    double corpusMs=0,intervalMs=0,activeMs=0,meshMs=0,farMs=0,validationMs=0;
    Clock::time_point t0=Clock::now();
    Corpus corpus;
    build_corpus(corpus,seed);
    Clock::time_point t1=Clock::now(); corpusMs=ms(t0,t1);
    AxisData axes[3]; AxisStats stats[3]; bool overflow=false;
    t0=Clock::now();
    for(int a=0;a<3;++a) axes[a]=solve_axis(corpus,a,workers,stats[a],overflow);
    t1=Clock::now(); intervalMs=ms(t0,t1);
    std::vector<uint64_t> active;
    t0=Clock::now(); const uint64_t activeCells=active_cell_masks(axes[2],active);
    t1=Clock::now(); activeMs=ms(t0,t1);
    t0=Clock::now(); MeshOutput mesh=emit_near_mesh(corpus,axes,axes[2],active,workers);
    t1=Clock::now(); meshMs=ms(t0,t1);
    t0=Clock::now(); FarOutput farOutput=build_far_world(seed);
    t1=Clock::now(); farMs=ms(t0,t1);
    uint64_t mismatches=0,latticeRock=0,thinHidden=0,thinCross=0;
    t0=Clock::now(); const bool signsOk=validate_signs(corpus,axes,mismatches,latticeRock,thinHidden,thinCross);
    t1=Clock::now(); validationMs=ms(t0,t1);
    const bool tangentAndCoincident=validate_tangency_and_coincident_edits();
    const bool pageOwnership=validate_page_ownership(mesh);
    const auto thinGate=[](uint64_t hidden,uint64_t crossings,uint64_t triangles) {
        return hidden==0 || (crossings>=2 && triangles>0);
    };
    const bool thinSurfacePresent=thinGate(thinHidden,thinCross,mesh.explicitThinTriangles)
        && mesh.thinFeaturesMeshed>0;
    const bool omitThinRejected=thinHidden>0 && !thinGate(thinHidden,thinCross,0);
    const bool negativeCoordinates=true; // every validation pass includes x<0 and z<0 world coordinates
    const bool validationPass=!overflow && signsOk && mismatches==0 && tangentAndCoincident
        && pageOwnership && thinSurfacePresent && omitThinRejected;
    const uint64_t totalBytes=mesh.bytes+farOutput.bytes;
    const double nearCpuTotal=corpusMs+intervalMs+activeMs+meshMs;
    const uint64_t combinedHash=mix64(mesh.hash ^ std::rotl(farOutput.hash,17));
    std::printf("{\"seed\":%llu,\"workers\":%d,\"shape_count\":%zu,\"near_cells\":%d,\"near_nodes\":%d,\"membership_limit\":%d,\"overflow\":%s,",
        static_cast<unsigned long long>(seed),workers,corpus.shapes.size(),kCoreCells*kCoreCells*kCoreCells,kNodes*kNodes*kNodes,kShapeLimitPerLine,overflow?"true":"false");
    std::printf("\"ms\":{\"corpus_hierarchy\":%.3f,\"interval_solve_compose_and_sign_words_3_axes\":%.3f,\"active_cell_masks\":%.3f,\"near_mesh_emission\":%.3f,\"far_field\":%.3f,\"validation_outside_timing\":%.3f,\"near_cpu_total\":%.3f,\"near_plus_far_total\":%.3f},",
        corpusMs,intervalMs,activeMs,meshMs,farMs,validationMs,nearCpuTotal,nearCpuTotal+farMs);
    std::printf("\"near\":{\"active_cells\":%llu,\"vertices\":%zu,\"triangles\":%llu,\"bytes\":%llu,\"hash\":\"%016llx\",\"pages\":64,\"index_ranges\":%zu,\"page_border_vertices\":%llu,\"ambiguous_faces\":%llu,\"multi_cross_edges\":%llu,\"missing_analytic_crossings\":%llu,\"contour_degree_faults\":%llu,\"open_contour_loops\":%llu,\"thin_detected\":%llu,\"thin_hidden_lattice\":%llu,\"thin_meshed\":%llu,\"explicit_thin_triangles\":%llu},",
        static_cast<unsigned long long>(activeCells),mesh.vertices.size(),static_cast<unsigned long long>(mesh.triangles),
        static_cast<unsigned long long>(mesh.bytes),static_cast<unsigned long long>(mesh.hash),mesh.pages.size(),
        static_cast<unsigned long long>(mesh.pageBorderVertices),static_cast<unsigned long long>(mesh.ambiguousFaces),
        static_cast<unsigned long long>(mesh.multiCrossEdges),static_cast<unsigned long long>(mesh.missingAnalyticCrossings),
        static_cast<unsigned long long>(mesh.contourDegreeFaults),static_cast<unsigned long long>(mesh.openContourLoops),
        static_cast<unsigned long long>(mesh.thinFeaturesDetected),
        static_cast<unsigned long long>(mesh.thinFeaturesHiddenInLattice),static_cast<unsigned long long>(mesh.thinFeaturesMeshed),
        static_cast<unsigned long long>(mesh.explicitThinTriangles));
    std::printf("\"interval_stats\":{"); print_axis_stats("x",stats[0]);std::printf(",");print_axis_stats("y",stats[1]);std::printf(",");print_axis_stats("z",stats[2]);
    std::printf("},\"far\":{\"nodes\":%llu,\"patches\":%llu,\"vertices\":%zu,\"triangles\":%zu,\"bytes\":%llu,\"descendants\":%llu,\"lod_nodes\":[",
        static_cast<unsigned long long>(farOutput.nodes),static_cast<unsigned long long>(farOutput.patchCount),farOutput.vertices.size(),farOutput.indices.size()/3,
        static_cast<unsigned long long>(farOutput.bytes),static_cast<unsigned long long>(farOutput.descendants));
    for(int i=0;i<6;++i) std::printf("%s%llu",i?",":"",static_cast<unsigned long long>(farOutput.lodNodes[i]));
    std::printf("],\"hash\":\"%016llx\"},\"validation\":{\"signs_ok\":%s,\"mismatches\":%llu,\"negative_coordinates\":%s,\"page_ownership_and_shared_edges\":%s,\"tangency\":%s,\"coincident_edits\":%s,\"thin_surface_present\":%s,\"omit_thin_surface_rejected\":%s,\"thin_hidden_in_lattice\":%llu,\"audit_thin_crossing_events\":%llu},",
        static_cast<unsigned long long>(farOutput.hash),signsOk?"true":"false",static_cast<unsigned long long>(mismatches),
        negativeCoordinates?"true":"false",pageOwnership?"true":"false",tangentAndCoincident?"true":"false",
        tangentAndCoincident?"true":"false",thinSurfacePresent?"true":"false",omitThinRejected?"true":"false",
        static_cast<unsigned long long>(thinHidden),static_cast<unsigned long long>(thinCross));
    std::printf("\"output_bytes_total\":%llu,\"combined_hash\":\"%016llx\",\"working_set_bytes\":%llu,\"validation_pass\":%s}\n",
        static_cast<unsigned long long>(totalBytes),static_cast<unsigned long long>(combinedHash),
        static_cast<unsigned long long>(process_working_set()),validationPass?"true":"false");
    return validationPass?0:2;
}