#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
verify_hgl4.py — HGL_4(F_4) のハミルトン閉路の**独立検証器**（第三者用）。

  依存は Python 3 の標準ライブラリだけ（numpy も scipy も要らない）。
  読むファイルは `hgl4_cycle.txt` ただ一つ。
  使い方:   python3 verify_hgl4.py [hgl4_cycle.txt]

検証するのは次の二つだけである。この二つで主張は尽きる。

  (i)  ファイルの 38,080 行は相異なり、可逆エルミート 4x4 行列の全体
       HGL_4(F_4) をちょうど尽くす。
  (ii) 巡回して隣り合う二行の差はすべて階数 1 である。

したがって、これらの行は HGL_4(F_4)（Orel 2015 の定義：頂点は可逆エルミート行列、
{A,B} が辺 ⟺ rk(A−B) = 1）のハミルトン閉路を与える。

## この検証器が使う判定法（他の実装と別経路にしてある）

  * F4 の乗法は「手で書いた 4x4 の乗積表」（体の公理 x^2 = x+1 から直接）。
    結合律・分配律・可換性・逆元の存在をすべて総当たりで確認してから使う。
  * 可逆性は掃き出しでも小行列式でもなく **ライプニッツ展開**
    （S_4 の 24 個の置換の総和。標数 2 なので符号は不要）。
  * 隣接判定は階数を測らない。差 D が
      S = { x x^* : x ∈ F4^4, x ≠ 0 }
    に属するかを集合所属で判定する。D = x x^* (x≠0) ならば D の列はすべて
    x の倍元だから rk(D) = 1。すなわち「D ∈ S」は rk(D)=1 の**十分条件**であり、
    辺の証明にはこれで足りる（逆向き「階数 1 のエルミート行列は必ず x x^* の形」は
    使わない）。S の各元が非零で全 2x2 小行列式 0 であることも総当たりで確かめる。
  * 差は符号のビット演算ではなく、行列の成分ごとの引き算（標数 2 なので XOR）。

## 出力

  各検査の結果と、閉路の三つの sha256（ファイル／保存順の符号列／
  回転・反転で正規化した符号列）。最後の一つは開始点と向きに依らないので、
  別に見つけた閉路が同じものかどうかを第三者が判定できる。
"""
import sys
import time
import hashlib
import itertools

t0 = time.time()
n = 4

# ---------------------------------------------------------------- F4 -------
# 元 0,1,2,3 を 0, 1, x, x+1 と読む。加法は 2 ビットの XOR。
# 乗法は体の公理から：1*a = a, x*x = x+1, x*(x+1) = 1, (x+1)*(x+1) = x
MUL = [
    [0, 0, 0, 0],
    [0, 1, 2, 3],
    [0, 2, 3, 1],
    [0, 3, 1, 2],
]
for a in range(4):
    for b in range(4):
        assert MUL[a][b] == MUL[b][a], "可換性"
        for c in range(4):
            assert MUL[MUL[a][b]][c] == MUL[a][MUL[b][c]], "結合律"
            assert MUL[a][b ^ c] == MUL[a][b] ^ MUL[a][c], "分配律"
for a in range(1, 4):
    assert any(MUL[a][b] == 1 for b in range(1, 4)), "逆元"

CONJ = [MUL[a][a] for a in range(4)]        # a ↦ a^2 は F4/F2 の非自明な自己同型
assert CONJ == [0, 1, 3, 2]
for a in range(4):
    assert CONJ[CONJ[a]] == a, "共役は対合"
assert [a for a in range(4) if CONJ[a] == a] == [0, 1], "不動体は F2"

PERMS = list(itertools.permutations(range(n)))
assert len(PERMS) == 24


def det_leibniz(A):
    """det A = Σ_{σ∈S_4} Π_i A[i][σ(i)]。標数 2 なので符号は不要。"""
    s = 0
    for p in PERMS:
        t = 1
        for i in range(n):
            t = MUL[t][A[i][p[i]]]
            if t == 0:
                break
        s ^= t
    return s


def is_hermitian(A):
    return all(A[j][i] == CONJ[A[i][j]] for i in range(n) for j in range(n))


# ------------------------------------------------- 符号化（.txt を読むための約束）
# 16 ビット: bit i (i=0..3) が対角 A[i][i] ∈ F2、
#            bit 4+2k, 4+2k+1 が PAIRS[k]=(i,j) の A[i][j] ∈ F4（下位ビットが先）。
# 2^16 = 65536 = 2^4 · 4^6 = |H_4(F_4)|。これは列挙とハッシュのための番号づけであって、
# 検証の中身ではない。
PAIRS = [(i, j) for i in range(n) for j in range(i + 1, n)]
NBITS = n + 2 * len(PAIRS)
assert NBITS == 16
NSTATE = 1 << NBITS


def decode(code):
    A = [[0] * n for _ in range(n)]
    for i in range(n):
        A[i][i] = (code >> i) & 1
    b = n
    for (i, j) in PAIRS:
        v = (code >> b) & 3
        b += 2
        A[i][j] = v
        A[j][i] = CONJ[v]
    return A


def encode(A):
    code = 0
    for i in range(n):
        assert A[i][i] in (0, 1), "対角は F2 でなければならない"
        code |= A[i][i] << i
    b = n
    for (i, j) in PAIRS:
        code |= A[i][j] << b
        b += 2
    return code


# ---------------------------------------------------------------- 入力 -----
path = sys.argv[1] if len(sys.argv) > 1 else 'hgl4_cycle.txt'
rows = []
with open(path) as f:
    raw = f.read()
for line in raw.splitlines():
    line = line.strip()
    if not line or line.startswith('#'):
        continue
    toks = line.replace(' ', '')
    assert len(toks) == 16, "一行は 16 個の F4 の元でなければならない: " + line
    A = [[int(toks[4 * i + j]) for j in range(n)] for i in range(n)]
    assert all(0 <= v <= 3 for r in A for v in r), "F4 の元は 0..3"
    rows.append(A)

print("=" * 78)
print("verify_hgl4.py — HGL_4(F_4) のハミルトン閉路の独立検証")
print("=" * 78)
print(f"[0] 入力 {path} : {len(rows)} 行")

# ------------------------------------------- [1] 各行がエルミートかつ可逆であること
nherm = sum(1 for A in rows if is_hermitian(A))
ninv = sum(1 for A in rows if det_leibniz(A) != 0)
print(f"[1] エルミートである行 = {nherm} / {len(rows)}")
print(f"[1] 可逆である行（ライプニッツ展開） = {ninv} / {len(rows)}")

# ------------------------------- [2] 頂点集合 HGL_4(F_4) を独立に作り、行の集合と比べる
V = []
for c in range(NSTATE):
    A = decode(c)
    if det_leibniz(A) != 0:
        V.append(c)
Vset = set(V)
print(f"[2] |H_4(F_4)| = {NSTATE}   (= 2^4 · 4^6)")
print(f"[2] |HGL_4(F_4)| = {len(V)}   "
      f"(Orel 2015 の公式 2^{{n(n-1)/2}} Π_j (2^j+(-1)^j) は n=4 で 64·1·5·7·17 = 38080)")
print(f"    ({time.time()-t0:.1f}s)")

codes = [encode(A) for A in rows]
for c, A in zip(codes, rows):
    assert decode(c) == A, "符号化と復号が逆写像でない"
ok_len = (len(codes) == len(V))
ok_dist = (len(set(codes)) == len(codes))
ok_set = (set(codes) == Vset)
print(f"[3] 行数 = |HGL_4(F_4)| : {ok_len}")
print(f"[3] 全行が相異なる      : {ok_dist}")
print(f"[3] 行の集合 = HGL_4(F_4) : {ok_set}")

# ------------------------------------------------------ [4] 階数 1 の集合 S
S = set()
for xs in itertools.product(range(4), repeat=n):
    if all(v == 0 for v in xs):
        continue
    M = tuple(tuple(MUL[xs[i]][CONJ[xs[j]]] for j in range(n)) for i in range(n))
    assert is_hermitian([list(r) for r in M])
    S.add(M)
print(f"[4] 外積 x x^* の相異なる値 |S| = {len(S)}   (= (4^4-1)/(4-1) = 85)")
bad = 0
for M in S:
    if all(v == 0 for r in M for v in r):
        bad += 1
    for i in range(n):
        for j in range(i + 1, n):
            for k in range(n):
                for l in range(k + 1, n):
                    if MUL[M[i][k]][M[j][l]] ^ MUL[M[i][l]][M[j][k]] != 0:
                        bad += 1
print(f"[4] S の元で「非零かつ全 2x2 小行列式 0」を破るもの = {bad} 個"
      "  ⇒ S の元の階数はちょうど 1")

# ------------------------------------------------------------ [5] 全辺の検査
nbad = 0
first_bad = None
used = set()
N = len(rows)
for idx in range(N):
    A = rows[idx]
    B = rows[(idx + 1) % N]
    D = tuple(tuple(A[i][j] ^ B[i][j] for j in range(n)) for i in range(n))
    if D not in S:
        nbad += 1
        if first_bad is None:
            first_bad = idx
    else:
        used.add(D)
print(f"[5] 閉路の辺の総数（末尾→先頭を含む） = {N}")
print(f"[5] 差 D = A − B が外積 x x^* の形でない辺 = {nbad} 本"
      + (f"（最初の破れは {first_bad} 行目）" if first_bad is not None else ""))
print(f"[5] 使われた階数 1 の元の種類 = {len(used)} / {len(S)}")

# --------------------------------------------------- [6] 次数（参考・版面との照合）
import random
random.seed(924)
degs = set()
for c in random.sample(V, 100):
    A = decode(c)
    d = 0
    for M in S:
        C = [[A[i][j] ^ M[i][j] for j in range(n)] for i in range(n)]
        if det_leibniz(C) != 0:
            d += 1
    degs.add(d)
print(f"[6] 無作為 100 頂点の次数の集合 = {sorted(degs)}"
      "   (Orel Prop.12: (2^n−(−1)^n)(2^{n−1}−(−1)^{n−1})/3 = 15·9/3 = 45)")

# ------------------------------------------------------------- [7] ハッシュ
h_file = hashlib.sha256(raw.encode()).hexdigest()
h_canon = hashlib.sha256(",".join(str(v) for v in codes).encode()).hexdigest()
m = min(range(N), key=lambda i: codes[i])
fwd = codes[m:] + codes[:m]
rev = [fwd[0]] + fwd[1:][::-1]
h_canon2 = hashlib.sha256(
    ",".join(str(v) for v in (fwd if fwd[1] < rev[1] else rev)).encode()).hexdigest()
print()
print(f"[7] sha256({path})                   = {h_file}")
print(f"[7] sha256(符号のコンマ区切り十進・保存順) = {h_canon}")
print(f"[7] sha256(回転・反転で正規化した符号列)   = {h_canon2}")

ok = (nherm == N and ninv == N and ok_len and ok_dist and ok_set
      and bad == 0 and nbad == 0)
print()
print("[判定] " + ("**合格**" if ok else "**不合格**")
      + f"  頂点 {len(V)} / 行 {N} / 相異なる {ok_dist} / 集合一致 {ok_set} / 悪い辺 {nbad}")
print(f"所要 {time.time()-t0:.1f} 秒")
sys.exit(0 if ok else 1)