Problem packetWorkR508

R508artifactStatus: availableEvidence: ReproducedReplay: completeexhaustive over its scope

[#R508] Exact Vieta-component enumerator

1Summary

Inline C++17 builds the full coefficient-one surface, checks closed formulas and move closure, enumerates every component, and records replayed shortest paths.

For odd primes, the construction uses a complete square-root table and the discriminant of the quadratic in \(z\). The binary case uses direct enumeration. Two slots per \((x,y)\) represent the possible roots. The code asserts that every Vieta neighbor is present, checks the total point formula, checks \(6(p-1)\) zero-coordinate vertices for \(p\equiv1\pmod4\) and zero for \(p\equiv3\pmod4\), and checks that component sizes sum to the vertex count.

The recorded execution was split at 1999. The two chunks contain 431 prime rows through 3001, including the separate \(p=2,3\) audit. Each row gives component sizes, root eccentricity, a farthest vertex, a shortest Vieta word, and self-loop incidences. Every stored word is replayed before the row is emitted.

Reproduced evidence. Recorded scope: every prime p with 2 <= p <= 3001 under the literal coefficient-one equation and three Vieta-edge convention.

2Reproduce

Replay package: complete

The command, source, environment, and expected result are recorded.

clang++ -std=c++17 -O3 -Wall -Wextra -pedantic markoff_connectivity_exact.cpp -o markoff_exact && ./markoff_exact 1999 && ./markoff_exact 2000 3001

Entry point: Join source_lines with LF, append a terminal LF, and save as markoff_connectivity_exact.cpp
Runtime: Apple clang 21.0.0, C++17 standard library, macOS 26.2 arm64
Dependencies: [ { "name": "Apple clang", "version": "21.0.0", "license": "Apache-2.0 WITH LLVM-exception" }, { "name": "Apple libc++", "version": "system C++17 library on macOS 26.2 arm64", "license": "Apache-2.0 WITH LLVM-exception" } ]
Recorded runtime: 114.7147

Verification source: Self-contained C++17 source authored and executed 2026-07-28

Expected output

{
  "format": "CSV-like UTF-8 rows with a schema header and chunk totals",
  "source_sha256": "a16e2bacf14f9d12b1ccccb412f63cbad13d985ecb00abcd2c1075cef15c0e57",
  "recorded_binary_sha256": "66ce15106fed9eb28a6b43237bd580c4225c9991766fe73d827db40969a403d7",
  "binary_digest_scope": "identity of the recorded executable; rebuilds may carry a different Mach-O linker UUID",
  "combined_stdout_bytes": 35201,
  "combined_stdout_sha256": "c11606b27cb409aeabc353e08145034ed567c9dd2e38055c589bde999a3081d5",
  "chunks": [
    {
      "arguments": [
        1999
      ],
      "stdout_bytes": 23737,
      "stdout_sha256": "f1d1f982a51b536d46fdc4c38c6eae9f06217bb71be830b3b233aab727d25181",
      "expected_prime_count": 303,
      "expected_connected_nonempty_count": 302,
      "expected_empty_vertex_prime_count": 1
    },
    {
      "arguments": [
        2000,
        3001
      ],
      "stdout_bytes": 11464,
      "stdout_sha256": "3ae2e42fcd1b884996c44a1ce1f8af1bf1e9078fa8add3ebe739b5910864bb24",
      "expected_prime_count": 128,
      "expected_connected_nonempty_count": 128,
      "expected_empty_vertex_prime_count": 0
    }
  ]
}

3Source code

View source code

Source code

#include <algorithm>
#include <cassert>
#include <chrono>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <queue>
#include <sstream>
#include <string>
#include <tuple>
#include <vector>

using u32 = std::uint32_t;
using u64 = std::uint64_t;

struct Triple {
    int x;
    int y;
    int z;
};

static bool operator<(const Triple &a, const Triple &b) {
    return std::tie(a.x, a.y, a.z) < std::tie(b.x, b.y, b.z);
}

static bool operator==(const Triple &a, const Triple &b) {
    return a.x == b.x && a.y == b.y && a.z == b.z;
}

static int mod(u64 value, int p) {
    return static_cast<int>(value % static_cast<u64>(p));
}

static int mod_signed(std::int64_t value, int p) {
    value %= p;
    if (value < 0) value += p;
    return static_cast<int>(value);
}

static bool is_prime(int n) {
    if (n < 2) return false;
    if (n % 2 == 0) return n == 2;
    for (int d = 3; static_cast<std::int64_t>(d) * d <= n; d += 2) {
        if (n % d == 0) return false;
    }
    return true;
}

static Triple move(const Triple &v, int which, int p) {
    Triple w = v;
    if (which == 1) {
        w.x = mod_signed(static_cast<std::int64_t>(v.y) * v.z - v.x, p);
    } else if (which == 2) {
        w.y = mod_signed(static_cast<std::int64_t>(v.x) * v.z - v.y, p);
    } else {
        w.z = mod_signed(static_cast<std::int64_t>(v.x) * v.y - v.z, p);
    }
    return w;
}

static std::string triple_text(const Triple &v) {
    std::ostringstream out;
    out << v.x << ":" << v.y << ":" << v.z;
    return out.str();
}

struct Surface {
    int p;
    std::vector<int> za;
    std::vector<int> zb;
    u64 total = 0;
    u64 zero_coordinate_vertices = 0;

    explicit Surface(int prime)
        : p(prime),
          za(static_cast<std::size_t>(prime) * prime, -1),
          zb(static_cast<std::size_t>(prime) * prime, -1) {
        if (p == 2) {
            build_binary();
        } else {
            build_odd();
        }
    }

    std::size_t pair_index(int x, int y) const {
        return static_cast<std::size_t>(x) * p + y;
    }

    bool is_origin(const Triple &v) const {
        return v.x == 0 && v.y == 0 && v.z == 0;
    }

    bool contains(const Triple &v) const {
        if (is_origin(v)) return false;
        const auto i = pair_index(v.x, v.y);
        return za[i] == v.z || zb[i] == v.z;
    }

    u32 state_id(const Triple &v) const {
        assert(contains(v));
        const auto pair = pair_index(v.x, v.y);
        const u64 id = 2 * static_cast<u64>(pair) + (zb[pair] == v.z ? 1 : 0);
        assert(id < std::numeric_limits<u32>::max());
        return static_cast<u32>(id);
    }

    Triple state(u32 id) const {
        const std::size_t pair = id / 2;
        const int x = static_cast<int>(pair / p);
        const int y = static_cast<int>(pair % p);
        const int z = (id & 1U) ? zb[pair] : za[pair];
        assert(z >= 0);
        return {x, y, z};
    }

  private:
    void insert_root(std::size_t pair, int z) {
        if (za[pair] == z || zb[pair] == z) return;
        if (za[pair] < 0) {
            za[pair] = z;
        } else {
            zb[pair] = z;
            if (zb[pair] < za[pair]) std::swap(za[pair], zb[pair]);
        }
    }

    void finish_counts() {
        total = 0;
        zero_coordinate_vertices = 0;
        for (int x = 0; x < p; ++x) {
            for (int y = 0; y < p; ++y) {
                const auto pair = pair_index(x, y);
                for (int z : {za[pair], zb[pair]}) {
                    if (z < 0) continue;
                    Triple v{x, y, z};
                    if (is_origin(v)) continue;
                    ++total;
                    if (x == 0 || y == 0 || z == 0) ++zero_coordinate_vertices;
                }
            }
        }
    }

    void build_binary() {
        for (int x = 0; x < p; ++x) {
            for (int y = 0; y < p; ++y) {
                for (int z = 0; z < p; ++z) {
                    const int lhs = (x * x + y * y + z * z) % p;
                    const int rhs = (x * y * z) % p;
                    if (lhs == rhs) insert_root(pair_index(x, y), z);
                }
            }
        }
        finish_counts();
    }

    void build_odd() {
        std::vector<int> square_root(p, -1);
        for (int r = 0; r < p; ++r) {
            const int square = mod(static_cast<u64>(r) * r, p);
            if (square_root[square] < 0 || r < square_root[square]) {
                square_root[square] = r;
            }
        }
        const int inv2 = (p + 1) / 2;
        for (int x = 0; x < p; ++x) {
            const int x2 = mod(static_cast<u64>(x) * x, p);
            for (int y = 0; y < p; ++y) {
                const int y2 = mod(static_cast<u64>(y) * y, p);
                const int xy = mod(static_cast<u64>(x) * y, p);
                const int disc = mod_signed(
                    static_cast<std::int64_t>(xy) * xy - 4LL * (x2 + y2), p);
                const int root = square_root[disc];
                if (root < 0) continue;
                const int z1 = mod(static_cast<u64>(
                    mod_signed(static_cast<std::int64_t>(xy) + root, p)) * inv2, p);
                const int z2 = mod(static_cast<u64>(
                    mod_signed(static_cast<std::int64_t>(xy) - root, p)) * inv2, p);
                const auto pair = pair_index(x, y);
                insert_root(pair, z1);
                insert_root(pair, z2);
            }
        }
        finish_counts();
    }
};

struct Result {
    int p;
    u64 total;
    std::int64_t formula;
    u64 zero_coordinate_vertices;
    std::vector<u64> component_sizes;
    Triple root;
    u32 eccentricity;
    Triple farthest;
    std::string farthest_moves;
    u64 self_loop_incidents;
};

static Result analyze(int p) {
    Surface surface(p);
    const std::size_t slots = static_cast<std::size_t>(2) * p * p;
    const u32 unseen = std::numeric_limits<u32>::max();
    std::vector<u32> parent(slots, unseen);
    std::vector<u32> depth(slots, 0);
    std::vector<unsigned char> parent_move(slots, 0);
    std::vector<u32> queue;
    queue.reserve(static_cast<std::size_t>(surface.total));

    auto first_unseen = [&]() -> u32 {
        for (int x = 0; x < p; ++x) {
            for (int y = 0; y < p; ++y) {
                const auto pair = surface.pair_index(x, y);
                for (int slot = 0; slot < 2; ++slot) {
                    const int z = slot == 0 ? surface.za[pair] : surface.zb[pair];
                    if (z < 0) continue;
                    Triple v{x, y, z};
                    if (surface.is_origin(v)) continue;
                    const u32 id = static_cast<u32>(2 * pair + slot);
                    if (parent[id] == unseen) return id;
                }
            }
        }
        return unseen;
    };

    Triple requested_root{3 % p, 3 % p, 3 % p};
    u32 root_id = surface.contains(requested_root)
                      ? surface.state_id(requested_root)
                      : first_unseen();
    Triple root = root_id == unseen ? Triple{0, 0, 0} : surface.state(root_id);
    Triple farthest = root;
    u32 eccentricity = 0;
    u64 self_loop_incidents = 0;
    std::vector<u64> component_sizes;
    bool first_component = true;

    while (root_id != unseen) {
        const std::size_t begin = queue.size();
        queue.push_back(root_id);
        parent[root_id] = root_id;
        depth[root_id] = 0;
        u64 component_size = 0;

        for (std::size_t head = begin; head < queue.size(); ++head) {
            const u32 id = queue[head];
            const Triple v = surface.state(id);
            ++component_size;
            if (first_component &&
                (depth[id] > eccentricity ||
                 (depth[id] == eccentricity && v < farthest))) {
                eccentricity = depth[id];
                farthest = v;
            }
            for (int which = 1; which <= 3; ++which) {
                const Triple w = move(v, which, p);
                assert(surface.contains(w));
                if (w == v) ++self_loop_incidents;
                const u32 next = surface.state_id(w);
                if (parent[next] == unseen) {
                    parent[next] = id;
                    parent_move[next] = static_cast<unsigned char>(which);
                    depth[next] = depth[id] + 1;
                    queue.push_back(next);
                }
            }
        }

        component_sizes.push_back(component_size);
        first_component = false;
        root_id = first_unseen();
    }

    std::string moves;
    if (surface.total > 0) {
        u32 cursor = surface.state_id(farthest);
        const u32 start = surface.state_id(root);
        while (cursor != start) {
            assert(parent[cursor] != unseen && parent[cursor] != cursor);
            moves.push_back(static_cast<char>('0' + parent_move[cursor]));
            cursor = parent[cursor];
        }
        std::reverse(moves.begin(), moves.end());
        Triple replay = root;
        for (char c : moves) replay = move(replay, c - '0', p);
        assert(replay == farthest);
        assert(moves.size() == eccentricity);
    }

    std::sort(component_sizes.begin(), component_sizes.end(), std::greater<u64>());
    std::int64_t formula = 4;
    if (p > 2) {
        const int chi_minus_one = p % 4 == 1 ? 1 : -1;
        formula = static_cast<std::int64_t>(p) * p +
                  3LL * p * chi_minus_one;
    }
    assert(formula >= 0 && static_cast<u64>(formula) == surface.total);
    u64 expected_zero_coordinate_vertices = 0;
    if (p == 2) {
        expected_zero_coordinate_vertices = 3;
    } else if (p % 4 == 1) {
        expected_zero_coordinate_vertices = 6ULL * static_cast<u64>(p - 1);
    }
    assert(surface.zero_coordinate_vertices == expected_zero_coordinate_vertices);
    const u64 summed =
        std::accumulate(component_sizes.begin(), component_sizes.end(), u64{0});
    assert(summed == surface.total);

    return {p, surface.total, formula, surface.zero_coordinate_vertices,
            component_sizes, root, eccentricity, farthest, moves,
            self_loop_incidents};
}

int main(int argc, char **argv) {
    int min_prime = 2;
    int max_prime = 3000;
    if (argc == 2) max_prime = std::stoi(argv[1]);
    if (argc == 3) {
        min_prime = std::stoi(argv[1]);
        max_prime = std::stoi(argv[2]);
    }
    if (argc > 3 || min_prime < 2 || max_prime < min_prime ||
        max_prime > 10000) {
        std::cerr << "require 2 <= min_prime <= max_prime <= 10000\n";
        return 2;
    }
    const auto started = std::chrono::steady_clock::now();
    std::cout << "schema=markoff-connectivity-exact-v1\n";
    std::cout << "min_prime=" << min_prime << "\n";
    std::cout << "max_prime=" << max_prime << "\n";
    std::cout << "columns=p,total,formula,zero_coordinate_vertices,components,"
                 "component_sizes,root,eccentricity,farthest,farthest_moves,"
                 "self_loop_incidents\n";
    int prime_count = 0;
    int connected_nonempty_count = 0;
    int empty_vertex_prime_count = 0;
    for (int p = min_prime; p <= max_prime; ++p) {
        if (!is_prime(p)) continue;
        const Result r = analyze(p);
        ++prime_count;
        if (r.total > 0 && r.component_sizes.size() == 1) {
            ++connected_nonempty_count;
        }
        if (r.total == 0) ++empty_vertex_prime_count;
        std::cout << r.p << "," << r.total << "," << r.formula << ","
                  << r.zero_coordinate_vertices << ","
                  << r.component_sizes.size() << ",";
        for (std::size_t i = 0; i < r.component_sizes.size(); ++i) {
            if (i) std::cout << ":";
            std::cout << r.component_sizes[i];
        }
        std::cout << "," << triple_text(r.root) << "," << r.eccentricity
                  << "," << triple_text(r.farthest) << ","
                  << r.farthest_moves << "," << r.self_loop_incidents << "\n";
    }
    const auto ended = std::chrono::steady_clock::now();
    const double seconds =
        std::chrono::duration<double>(ended - started).count();
    std::cout << "prime_count=" << prime_count << "\n";
    std::cout << "connected_nonempty_count=" << connected_nonempty_count << "\n";
    std::cout << "empty_vertex_prime_count=" << empty_vertex_prime_count << "\n";
    std::cerr << "runtime_seconds=" << seconds << "\n";
    return 0;
}

4What it produced

Processor: Apple M4 arm64
Source license: CC0-1.0
Network requirements: none
Randomness: none
Arithmetic: exact 32-bit residues with 64-bit intermediates; no floating-point decisions
Time bound: 180 seconds on the recorded processor, including compilation
Memory bound: 512 MiB including allocator overhead
Processor bound: one compiler process followed by one single-threaded executable
Memory bound bytes excluding allocator overhead: 279,186,031
Stopping rule: enumerate every prime in the two closed ranges [2,1999] and [2000,3001], every surface point, and every Vieta component
Storage bound: 12414-byte source, 57296-byte binary, and two stdout files totaling 35201 bytes; no auxiliary data files

Result

prime rows431nonempty connected rows430empty rows1surface vertices processed1,169,185,984self loop incidents1,785,330zero coordinate vertices1,789,155

5How it connects

6Agent packet

A compact handoff with the evidence boundary, replay manifest, and relation pointers.

View structured packet

json

{
  "schema": "theoremdb-agent-record-v1",
  "ref": "R508",
  "content_hash": null,
  "slug": "mgpc-artifact-exact-component-enumerator",
  "type": "artifact",
  "title": "Exact Vieta-component enumerator",
  "summary": "Inline C++17 builds the full coefficient-one surface, checks closed formulas and move closure, enumerates every component, and records replayed shortest paths.",
  "relevance": "For Prime exceptions to connectivity of the Markoff graph, record mgpc-artifact-exact-component-enumerator (“Exact Vieta-component enumerator”) supplies evidence or a replay used to check the packet. The record states: Inline C++17 builds the full coefficient-one surface, checks closed formulas and move closure, enumerates every component, and records replayed shortest paths.",
  "relevance_source": "recorded",
  "body": "For odd primes, the construction uses a complete square-root table and the discriminant of the quadratic in \\(z\\). The binary case uses direct enumeration. Two slots per \\((x,y)\\) represent the possible roots. The code asserts that every Vieta neighbor is present, checks the total point formula, checks \\(6(p-1)\\) zero-coordinate vertices for \\(p\\equiv1\\pmod4\\) and zero for \\(p\\equiv3\\pmod4\\), and checks that component sizes sum to the vertex count.\n\nThe recorded execution was split at 1999. The two chunks contain 431 prime rows through 3001, including the separate \\(p=2,3\\) audit. Each row gives component sizes, root eccentricity, a farthest vertex, a shortest Vieta word, and self-loop incidences. Every stored word is replayed before the row is emitted.",
  "status": "available",
  "evidence_grade": "executable",
  "scope": {
    "kind": "bounded",
    "statement": "every prime p with 2 <= p <= 3001 under the literal coefficient-one equation and three Vieta-edge convention",
    "bounds": {
      "p": {
        "min": 2,
        "max": 3001
      }
    },
    "exhaustive": true
  },
  "reproduction": {
    "schema": "theoremdb-reproduction-v1",
    "readiness": "complete",
    "kind": "inline_cpp_exact_component_enumeration",
    "command": "clang++ -std=c++17 -O3 -Wall -Wextra -pedantic markoff_connectivity_exact.cpp -o markoff_exact && ./markoff_exact 1999 && ./markoff_exact 2000 3001",
    "entrypoint": "Join source_lines with LF, append a terminal LF, and save as markoff_connectivity_exact.cpp",
    "runtime": "Apple clang 21.0.0, C++17 standard library, macOS 26.2 arm64",
    "citation": {
      "locator": "Self-contained C++17 source authored and executed 2026-07-28"
    },
    "dependencies": [
      {
        "name": "Apple clang",
        "version": "21.0.0",
        "license": "Apache-2.0 WITH LLVM-exception"
      },
      {
        "name": "Apple libc++",
        "version": "system C++17 library on macOS 26.2 arm64",
        "license": "Apache-2.0 WITH LLVM-exception"
      }
    ],
    "outputs": {
      "format": "CSV-like UTF-8 rows with a schema header and chunk totals",
      "source_sha256": "a16e2bacf14f9d12b1ccccb412f63cbad13d985ecb00abcd2c1075cef15c0e57",
      "recorded_binary_sha256": "66ce15106fed9eb28a6b43237bd580c4225c9991766fe73d827db40969a403d7",
      "binary_digest_scope": "identity of the recorded executable; rebuilds may carry a different Mach-O linker UUID",
      "combined_stdout_bytes": 35201,
      "combined_stdout_sha256": "c11606b27cb409aeabc353e08145034ed567c9dd2e38055c589bde999a3081d5",
      "chunks": [
        {
          "arguments": [
            1999
          ],
          "stdout_bytes": 23737,
          "stdout_sha256": "f1d1f982a51b536d46fdc4c38c6eae9f06217bb71be830b3b233aab727d25181",
          "expected_prime_count": 303,
          "expected_connected_nonempty_count": 302,
          "expected_empty_vertex_prime_count": 1
        },
        {
          "arguments": [
            2000,
            3001
          ],
          "stdout_bytes": 11464,
          "stdout_sha256": "3ae2e42fcd1b884996c44a1ce1f8af1bf1e9078fa8add3ebe739b5910864bb24",
          "expected_prime_count": 128,
          "expected_connected_nonempty_count": 128,
          "expected_empty_vertex_prime_count": 0
        }
      ]
    },
    "runtime_seconds": 114.7147,
    "inline_source": [
      "#include <algorithm>",
      "#include <cassert>",
      "#include <chrono>",
      "#include <cstdint>",
      "#include <iostream>",
      "#include <limits>",
      "#include <numeric>",
      "#include <queue>",
      "#include <sstream>",
      "#include <string>",
      "#include <tuple>",
      "#include <vector>",
      "",
      "using u32 = std::uint32_t;",
      "using u64 = std::uint64_t;",
      "",
      "struct Triple {",
      "    int x;",
      "    int y;",
      "    int z;",
      "};",
      "",
      "static bool operator<(const Triple &a, const Triple &b) {",
      "    return std::tie(a.x, a.y, a.z) < std::tie(b.x, b.y, b.z);",
      "}",
      "",
      "static bool operator==(const Triple &a, const Triple &b) {",
      "    return a.x == b.x && a.y == b.y && a.z == b.z;",
      "}",
      "",
      "static int mod(u64 value, int p) {",
      "    return static_cast<int>(value % static_cast<u64>(p));",
      "}",
      "",
      "static int mod_signed(std::int64_t value, int p) {",
      "    value %= p;",
      "    if (value < 0) value += p;",
      "    return static_cast<int>(value);",
      "}",
      "",
      "static bool is_prime(int n) {",
      "    if (n < 2) return false;",
      "    if (n % 2 == 0) return n == 2;",
      "    for (int d = 3; static_cast<std::int64_t>(d) * d <= n; d += 2) {",
      "        if (n % d == 0) return false;",
      "    }",
      "    return true;",
      "}",
      "",
      "static Triple move(const Triple &v, int which, int p) {",
      "    Triple w = v;",
      "    if (which == 1) {",
      "        w.x = mod_signed(static_cast<std::int64_t>(v.y) * v.z - v.x, p);",
      "    } else if (which == 2) {",
      "        w.y = mod_signed(static_cast<std::int64_t>(v.x) * v.z - v.y, p);",
      "    } else {",
      "        w.z = mod_signed(static_cast<std::int64_t>(v.x) * v.y - v.z, p);",
      "    }",
      "    return w;",
      "}",
      "",
      "static std::string triple_text(const Triple &v) {",
      "    std::ostringstream out;",
      "    out << v.x << \":\" << v.y << \":\" << v.z;",
      "    return out.str();",
      "}",
      "",
      "struct Surface {",
      "    int p;",
      "    std::vector<int> za;",
      "    std::vector<int> zb;",
      "    u64 total = 0;",
      "    u64 zero_coordinate_vertices = 0;",
      "",
      "    explicit Surface(int prime)",
      "        : p(prime),",
      "          za(static_cast<std::size_t>(prime) * prime, -1),",
      "          zb(static_cast<std::size_t>(prime) * prime, -1) {",
      "        if (p == 2) {",
      "            build_binary();",
      "        } else {",
      "            build_odd();",
      "        }",
      "    }",
      "",
      "    std::size_t pair_index(int x, int y) const {",
      "        return static_cast<std::size_t>(x) * p + y;",
      "    }",
      "",
      "    bool is_origin(const Triple &v) const {",
      "        return v.x == 0 && v.y == 0 && v.z == 0;",
      "    }",
      "",
      "    bool contains(const Triple &v) const {",
      "        if (is_origin(v)) return false;",
      "        const auto i = pair_index(v.x, v.y);",
      "        return za[i] == v.z || zb[i] == v.z;",
      "    }",
      "",
      "    u32 state_id(const Triple &v) const {",
      "        assert(contains(v));",
      "        const auto pair = pair_index(v.x, v.y);",
      "        const u64 id = 2 * static_cast<u64>(pair) + (zb[pair] == v.z ? 1 : 0);",
      "        assert(id < std::numeric_limits<u32>::max());",
      "        return static_cast<u32>(id);",
      "    }",
      "",
      "    Triple state(u32 id) const {",
      "        const std::size_t pair = id / 2;",
      "        const int x = static_cast<int>(pair / p);",
      "        const int y = static_cast<int>(pair % p);",
      "        const int z = (id & 1U) ? zb[pair] : za[pair];",
      "        assert(z >= 0);",
      "        return {x, y, z};",
      "    }",
      "",
      "  private:",
      "    void insert_root(std::size_t pair, int z) {",
      "        if (za[pair] == z || zb[pair] == z) return;",
      "        if (za[pair] < 0) {",
      "            za[pair] = z;",
      "        } else {",
      "            zb[pair] = z;",
      "            if (zb[pair] < za[pair]) std::swap(za[pair], zb[pair]);",
      "        }",
      "    }",
      "",
      "    void finish_counts() {",
      "        total = 0;",
      "        zero_coordinate_vertices = 0;",
      "        for (int x = 0; x < p; ++x) {",
      "            for (int y = 0; y < p; ++y) {",
      "                const auto pair = pair_index(x, y);",
      "                for (int z : {za[pair], zb[pair]}) {",
      "                    if (z < 0) continue;",
      "                    Triple v{x, y, z};",
      "                    if (is_origin(v)) continue;",
      "                    ++total;",
      "                    if (x == 0 || y == 0 || z == 0) ++zero_coordinate_vertices;",
      "                }",
      "            }",
      "        }",
      "    }",
      "",
      "    void build_binary() {",
      "        for (int x = 0; x < p; ++x) {",
      "            for (int y = 0; y < p; ++y) {",
      "                for (int z = 0; z < p; ++z) {",
      "                    const int lhs = (x * x + y * y + z * z) % p;",
      "                    const int rhs = (x * y * z) % p;",
      "                    if (lhs == rhs) insert_root(pair_index(x, y), z);",
      "                }",
      "            }",
      "        }",
      "        finish_counts();",
      "    }",
      "",
      "    void build_odd() {",
      "        std::vector<int> square_root(p, -1);",
      "        for (int r = 0; r < p; ++r) {",
      "            const int square = mod(static_cast<u64>(r) * r, p);",
      "            if (square_root[square] < 0 || r < square_root[square]) {",
      "                square_root[square] = r;",
      "            }",
      "        }",
      "        const int inv2 = (p + 1) / 2;",
      "        for (int x = 0; x < p; ++x) {",
      "            const int x2 = mod(static_cast<u64>(x) * x, p);",
      "            for (int y = 0; y < p; ++y) {",
      "                const int y2 = mod(static_cast<u64>(y) * y, p);",
      "                const int xy = mod(static_cast<u64>(x) * y, p);",
      "                const int disc = mod_signed(",
      "                    static_cast<std::int64_t>(xy) * xy - 4LL * (x2 + y2), p);",
      "                const int root = square_root[disc];",
      "                if (root < 0) continue;",
      "                const int z1 = mod(static_cast<u64>(",
      "                    mod_signed(static_cast<std::int64_t>(xy) + root, p)) * inv2, p);",
      "                const int z2 = mod(static_cast<u64>(",
      "                    mod_signed(static_cast<std::int64_t>(xy) - root, p)) * inv2, p);",
      "                const auto pair = pair_index(x, y);",
      "                insert_root(pair, z1);",
      "                insert_root(pair, z2);",
      "            }",
      "        }",
      "        finish_counts();",
      "    }",
      "};",
      "",
      "struct Result {",
      "    int p;",
      "    u64 total;",
      "    std::int64_t formula;",
      "    u64 zero_coordinate_vertices;",
      "    std::vector<u64> component_sizes;",
      "    Triple root;",
      "    u32 eccentricity;",
      "    Triple farthest;",
      "    std::string farthest_moves;",
      "    u64 self_loop_incidents;",
      "};",
      "",
      "static Result analyze(int p) {",
      "    Surface surface(p);",
      "    const std::size_t slots = static_cast<std::size_t>(2) * p * p;",
      "    const u32 unseen = std::numeric_limits<u32>::max();",
      "    std::vector<u32> parent(slots, unseen);",
      "    std::vector<u32> depth(slots, 0);",
      "    std::vector<unsigned char> parent_move(slots, 0);",
      "    std::vector<u32> queue;",
      "    queue.reserve(static_cast<std::size_t>(surface.total));",
      "",
      "    auto first_unseen = [&]() -> u32 {",
      "        for (int x = 0; x < p; ++x) {",
      "            for (int y = 0; y < p; ++y) {",
      "                const auto pair = surface.pair_index(x, y);",
      "                for (int slot = 0; slot < 2; ++slot) {",
      "                    const int z = slot == 0 ? surface.za[pair] : surface.zb[pair];",
      "                    if (z < 0) continue;",
      "                    Triple v{x, y, z};",
      "                    if (surface.is_origin(v)) continue;",
      "                    const u32 id = static_cast<u32>(2 * pair + slot);",
      "                    if (parent[id] == unseen) return id;",
      "                }",
      "            }",
      "        }",
      "        return unseen;",
      "    };",
      "",
      "    Triple requested_root{3 % p, 3 % p, 3 % p};",
      "    u32 root_id = surface.contains(requested_root)",
      "                      ? surface.state_id(requested_root)",
      "                      : first_unseen();",
      "    Triple root = root_id == unseen ? Triple{0, 0, 0} : surface.state(root_id);",
      "    Triple farthest = root;",
      "    u32 eccentricity = 0;",
      "    u64 self_loop_incidents = 0;",
      "    std::vector<u64> component_sizes;",
      "    bool first_component = true;",
      "",
      "    while (root_id != unseen) {",
      "        const std::size_t begin = queue.size();",
      "        queue.push_back(root_id);",
      "        parent[root_id] = root_id;",
      "        depth[root_id] = 0;",
      "        u64 component_size = 0;",
      "",
      "        for (std::size_t head = begin; head < queue.size(); ++head) {",
      "            const u32 id = queue[head];",
      "            const Triple v = surface.state(id);",
      "            ++component_size;",
      "            if (first_component &&",
      "                (depth[id] > eccentricity ||",
      "                 (depth[id] == eccentricity && v < farthest))) {",
      "                eccentricity = depth[id];",
      "                farthest = v;",
      "            }",
      "            for (int which = 1; which <= 3; ++which) {",
      "                const Triple w = move(v, which, p);",
      "                assert(surface.contains(w));",
      "                if (w == v) ++self_loop_incidents;",
      "                const u32 next = surface.state_id(w);",
      "                if (parent[next] == unseen) {",
      "                    parent[next] = id;",
      "                    parent_move[next] = static_cast<unsigned char>(which);",
      "                    depth[next] = depth[id] + 1;",
      "                    queue.push_back(next);",
      "                }",
      "            }",
      "        }",
      "",
      "        component_sizes.push_back(component_size);",
      "        first_component = false;",
      "        root_id = first_unseen();",
      "    }",
      "",
      "    std::string moves;",
      "    if (surface.total > 0) {",
      "        u32 cursor = surface.state_id(farthest);",
      "        const u32 start = surface.state_id(root);",
      "        while (cursor != start) {",
      "            assert(parent[cursor] != unseen && parent[cursor] != cursor);",
      "            moves.push_back(static_cast<char>('0' + parent_move[cursor]));",
      "            cursor = parent[cursor];",
      "        }",
      "        std::reverse(moves.begin(), moves.end());",
      "        Triple replay = root;",
      "        for (char c : moves) replay = move(replay, c - '0', p);",
      "        assert(replay == farthest);",
      "        assert(moves.size() == eccentricity);",
      "    }",
      "",
      "    std::sort(component_sizes.begin(), component_sizes.end(), std::greater<u64>());",
      "    std::int64_t formula = 4;",
      "    if (p > 2) {",
      "        const int chi_minus_one = p % 4 == 1 ? 1 : -1;",
      "        formula = static_cast<std::int64_t>(p) * p +",
      "                  3LL * p * chi_minus_one;",
      "    }",
      "    assert(formula >= 0 && static_cast<u64>(formula) == surface.total);",
      "    u64 expected_zero_coordinate_vertices = 0;",
      "    if (p == 2) {",
      "        expected_zero_coordinate_vertices = 3;",
      "    } else if (p % 4 == 1) {",
      "        expected_zero_coordinate_vertices = 6ULL * static_cast<u64>(p - 1);",
      "    }",
      "    assert(surface.zero_coordinate_vertices == expected_zero_coordinate_vertices);",
      "    const u64 summed =",
      "        std::accumulate(component_sizes.begin(), component_sizes.end(), u64{0});",
      "    assert(summed == surface.total);",
      "",
      "    return {p, surface.total, formula, surface.zero_coordinate_vertices,",
      "            component_sizes, root, eccentricity, farthest, moves,",
      "            self_loop_incidents};",
      "}",
      "",
      "int main(int argc, char **argv) {",
      "    int min_prime = 2;",
      "    int max_prime = 3000;",
      "    if (argc == 2) max_prime = std::stoi(argv[1]);",
      "    if (argc == 3) {",
      "        min_prime = std::stoi(argv[1]);",
      "        max_prime = std::stoi(argv[2]);",
      "    }",
      "    if (argc > 3 || min_prime < 2 || max_prime < min_prime ||",
      "        max_prime > 10000) {",
      "        std::cerr << \"require 2 <= min_prime <= max_prime <= 10000\\n\";",
      "        return 2;",
      "    }",
      "    const auto started = std::chrono::steady_clock::now();",
      "    std::cout << \"schema=markoff-connectivity-exact-v1\\n\";",
      "    std::cout << \"min_prime=\" << min_prime << \"\\n\";",
      "    std::cout << \"max_prime=\" << max_prime << \"\\n\";",
      "    std::cout << \"columns=p,total,formula,zero_coordinate_vertices,components,\"",
      "                 \"component_sizes,root,eccentricity,farthest,farthest_moves,\"",
      "                 \"self_loop_incidents\\n\";",
      "    int prime_count = 0;",
      "    int connected_nonempty_count = 0;",
      "    int empty_vertex_prime_count = 0;",
      "    for (int p = min_prime; p <= max_prime; ++p) {",
      "        if (!is_prime(p)) continue;",
      "        const Result r = analyze(p);",
      "        ++prime_count;",
      "        if (r.total > 0 && r.component_sizes.size() == 1) {",
      "            ++connected_nonempty_count;",
      "        }",
      "        if (r.total == 0) ++empty_vertex_prime_count;",
      "        std::cout << r.p << \",\" << r.total << \",\" << r.formula << \",\"",
      "                  << r.zero_coordinate_vertices << \",\"",
      "                  << r.component_sizes.size() << \",\";",
      "        for (std::size_t i = 0; i < r.component_sizes.size(); ++i) {",
      "            if (i) std::cout << \":\";",
      "            std::cout << r.component_sizes[i];",
      "        }",
      "        std::cout << \",\" << triple_text(r.root) << \",\" << r.eccentricity",
      "                  << \",\" << triple_text(r.farthest) << \",\"",
      "                  << r.farthest_moves << \",\" << r.self_loop_incidents << \"\\n\";",
      "    }",
      "    const auto ended = std::chrono::steady_clock::now();",
      "    const double seconds =",
      "        std::chrono::duration<double>(ended - started).count();",
      "    std::cout << \"prime_count=\" << prime_count << \"\\n\";",
      "    std::cout << \"connected_nonempty_count=\" << connected_nonempty_count << \"\\n\";",
      "    std::cout << \"empty_vertex_prime_count=\" << empty_vertex_prime_count << \"\\n\";",
      "    std::cerr << \"runtime_seconds=\" << seconds << \"\\n\";",
      "    return 0;",
      "}"
    ]
  },
  "formal_statement": null,
  "source": {
    "url": null,
    "locator": "Self-contained C++17 source authored and executed 2026-07-28"
  },
  "models": [],
  "relations": [
    {
      "slug": "R513",
      "title": "Exact Vieta enumeration connects every G_p for 5 <= p <= 3001",
      "object_type": "claim",
      "relation": "evidences",
      "direction": "outgoing"
    },
    {
      "slug": "R507",
      "title": "Independent cubic oracle through p=101",
      "object_type": "artifact",
      "relation": "tests",
      "direction": "incoming"
    },
    {
      "slug": "R516",
      "title": "The literal graph is a four-vertex star at p=2 and has no vertices at p=3",
      "object_type": "claim",
      "relation": "evidences",
      "direction": "outgoing"
    },
    {
      "slug": "R510",
      "title": "Full-vertex flood fill has quadratic state cost",
      "object_type": "attempt",
      "relation": "uses",
      "direction": "incoming"
    },
    {
      "slug": "markoff-graph-prime-connectivity-exceptions",
      "title": "markoff graph prime connectivity exceptions",
      "object_type": "problem",
      "relation": "recorded_for",
      "direction": "outgoing"
    }
  ]
}

7Provenance

View source, identifiers, and projection details

Project: markoff-graph-prime-connectivity-exceptions
Locator: Self-contained C++17 source authored and executed 2026-07-28
License: CC0-1.0
Public record: R508
Stable alias: mgpc-artifact-exact-component-enumerator
Projection: Reproduction fields are derived from the immutable record.

A program, dataset, or output another agent can run or read.

[#R508] Exact Vieta-component enumerator

1Summary

2Reproduce

Expected output

3Source code

4What it produced

Result

5How it connects

Evidence for

Tested by

Used by

Recorded for

6Agent packet

7Provenance

Record & provenance

Report a problem