2607.19268v1 / tools/check_finite.c

all files

/*
 * Fast numerical companion to check_finite.py for the larger finite cases.
 *
 * This program is evidence only, not part of the Lean proof.  It accepts an
 * order and a graph6 file produced by plantri, then reports the three largest
 * adjacency spectral radii found by normalized power iteration.
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define MAXN 62
#define MAXLINE 2048

struct result {
    double radius;
    char code[MAXLINE];
    int degree[MAXN];
};

static int decode_graph6(const char *line, unsigned char adjacency[MAXN][MAXN])
{
    const unsigned char *cursor = (const unsigned char *)line;
    int n, upper, lower, bit_index = 0;
    memset(adjacency, 0, MAXN * MAXN);
    if (strncmp((const char *)cursor, ">>graph6<<", 10) == 0)
        cursor += 10;
    if (*cursor == '~')
        return -1;
    n = *cursor++ - 63;
    for (upper = 1; upper < n; ++upper) {
        for (lower = 0; lower < upper; ++lower) {
            int value = cursor[bit_index / 6] - 63;
            int bit = (value >> (5 - bit_index % 6)) & 1;
            if (bit)
                adjacency[lower][upper] = adjacency[upper][lower] = 1;
            ++bit_index;
        }
    }
    return n;
}

static double spectral_radius(
    int n, const unsigned char adjacency[MAXN][MAXN])
{
    double x[MAXN], y[MAXN], product[MAXN];
    double quotient = 0.0;
    int i, j, iteration;
    for (i = 0; i < n; ++i)
        x[i] = 1.0 / sqrt((double)n);
    for (iteration = 0; iteration < 500; ++iteration) {
        double norm = 0.0, next_quotient = 0.0;
        for (i = 0; i < n; ++i) {
            y[i] = 0.0;
            for (j = 0; j < n; ++j)
                y[i] += adjacency[i][j] * x[j];
            norm += y[i] * y[i];
        }
        norm = sqrt(norm);
        for (i = 0; i < n; ++i)
            x[i] = y[i] / norm;
        for (i = 0; i < n; ++i) {
            product[i] = 0.0;
            for (j = 0; j < n; ++j)
                product[i] += adjacency[i][j] * x[j];
            next_quotient += x[i] * product[i];
        }
        if (fabs(next_quotient - quotient) < 1e-14)
            return next_quotient;
        quotient = next_quotient;
    }
    return quotient;
}

static void sort_degrees(int n, int degree[MAXN])
{
    int i, j;
    for (i = 1; i < n; ++i) {
        int value = degree[i];
        for (j = i; j > 0 && degree[j - 1] < value; --j)
            degree[j] = degree[j - 1];
        degree[j] = value;
    }
}

static double target_radius(int n)
{
    unsigned char adjacency[MAXN][MAXN] = {{0}};
    int vertex;
    adjacency[0][1] = adjacency[1][0] = 1;
    for (vertex = 2; vertex < n; ++vertex) {
        adjacency[0][vertex] = adjacency[vertex][0] = 1;
        adjacency[1][vertex] = adjacency[vertex][1] = 1;
    }
    for (vertex = 2; vertex + 1 < n; ++vertex)
        adjacency[vertex][vertex + 1] =
            adjacency[vertex + 1][vertex] = 1;
    return spectral_radius(n, adjacency);
}

static void consider(
    int n, double radius, const char *code,
    const unsigned char adjacency[MAXN][MAXN], struct result best[3])
{
    int rank, i, j;
    struct result candidate;
    candidate.radius = radius;
    snprintf(candidate.code, sizeof(candidate.code), "%s", code);
    candidate.code[strcspn(candidate.code, "\r\n")] = '\0';
    for (i = 0; i < n; ++i) {
        candidate.degree[i] = 0;
        for (j = 0; j < n; ++j)
            candidate.degree[i] += adjacency[i][j];
    }
    sort_degrees(n, candidate.degree);
    for (rank = 0; rank < 3; ++rank) {
        if (candidate.radius > best[rank].radius) {
            for (i = 2; i > rank; --i)
                best[i] = best[i - 1];
            best[rank] = candidate;
            return;
        }
    }
}

int main(int argc, char **argv)
{
    FILE *input;
    char line[MAXLINE];
    unsigned char adjacency[MAXN][MAXN];
    struct result best[3] = {{-INFINITY, "", {0}},
                             {-INFINITY, "", {0}},
                             {-INFINITY, "", {0}}};
    long count = 0;
    int expected, n, rank, i;
    double target;
    if (argc != 3) {
        fprintf(stderr, "usage: %s ORDER GRAPH6_FILE\n", argv[0]);
        return 2;
    }
    expected = atoi(argv[1]);
    input = fopen(argv[2], "rb");
    if (!input) {
        perror(argv[2]);
        return 2;
    }
    while (fgets(line, sizeof(line), input)) {
        if (strncmp(line, ">>", 2) == 0)
            continue;
        n = decode_graph6(line, adjacency);
        if (n != expected) {
            fprintf(stderr, "expected order %d, got %d\n", expected, n);
            return 2;
        }
        consider(n, spectral_radius(n, adjacency), line, adjacency, best);
        ++count;
    }
    fclose(input);
    target = target_radius(expected);
    printf("n=%d count=%ld target=%.12f\n", expected, count, target);
    for (rank = 0; rank < 3; ++rank) {
        printf("  %d: radius=%.12f delta=%+.12e degrees=[",
               rank + 1, best[rank].radius, best[rank].radius - target);
        for (i = 0; i < expected; ++i)
            printf("%s%d", i ? ", " : "", best[rank].degree[i]);
        printf("] graph6=%s\n", best[rank].code);
    }
    return 0;
}