# VERIF_MELAS_M2.py  --  pure Python 3 (no numpy, no scipy). Run as is: no arguments.
# Launch with plain Python, NOT inside a Sage cell.
# Checks, for k = KMIN..KMAX_DIRECT and every c in F_{2**k}***:
#   (1) M2(k,c) computed directly from pattern counts of S_c(i) = Tr(a**i + c a**{-i})
#       equals 1 - 2**-k + 2**-k * W(c),  W(c) = sum_{h != 0, deg h < k} chi(kappa c / Nm(h)),
#       Nm(h) = h(a) h(1/a), for one constant kappa (found, then checked for ALL c);
#   (2) the representation law r(P) = (k - deg P) * (1*chi0)(P) for the norm map
#       h -> h(x)h(1/x) in F_2[s], s = x + 1/x, chi0(D) = (-1)**{D_1} if D_0 = 1, else 0;
#   (3) worst case over c of |M2(k,c) - 1|, for k up to KMAX_TABLE (via (1)).
# The character '**' is NOT used anywhere: xor is operator.xor.
import operator, time, sys
XOR = operator.xor
KMIN = int(4)
KMAX_DIRECT = int(10)      # direct pattern counting (slowest part)
KMAX_LAW = int(9)         # representation law
KMAX_TABLE = int(16)      # worst-case table
TIME_BUDGET = float(480.0)
PRIM = {3:0b1011,4:0b10011,5:0b100101,6:0b1000011,7:0b10000011,8:0b100011101,
        9:0b1000010001,10:0b10000001001,11:0b100000000101,12:0b1000001010011,
        13:0b10000000011011,14:0b100010001000011,15:0b1000000000000011,16:0b10001000000001011}
def say(*a):
    print(*a, flush=True)

# ---------------- arithmetic self-test ----------------
assert XOR(int(6), int(3)) == int(5)
assert (int(1) << int(10)) == int(1024)
assert int(7) // int(2) == int(3)
assert int(2) ** int(10) == int(1024)
say("self-test OK")
T0 = time.time()

def field(k):
    P = PRIM[k]; n = (int(1) << k) - int(1)
    ex = [0] * n; x = int(1)
    for i in range(n):
        ex[i] = x; x = x << 1
        if x >> k: x = XOR(x, P)
    assert x == 1, "polynomial not primitive"
    lg = [-1] * (1 << k)
    for i in range(n): lg[ex[i]] = i
    def mul(a, b):
        if a == 0 or b == 0: return 0
        return ex[(lg[a] + lg[b]) % n]
    def tr(v):
        t = v; s = v
        for _ in range(k - 1):
            s = mul(s, s); t = XOR(t, s)
        assert t in (0, 1)
        return t
    TR = [0] * (1 << k)
    basis = [tr(1 << j) for j in range(k)]
    for v in range(1 << k):
        t = 0
        for j in range(k):
            if (v >> j) & 1 and basis[j]: t = XOR(t, 1)
        TR[v] = t
    return n, ex, lg, TR, mul

def norm_logs(k, n, ex, lg):
    # log of Nm(h) = h(a) h(a**-1) for every h != 0 with deg h < k
    out = []
    for h in range(1, 1 << k):
        y1 = 0; y2 = 0
        for i in range(k):
            if (h >> i) & 1:
                y1 = XOR(y1, ex[i % n]); y2 = XOR(y2, ex[(-i) % n])
        out.append((lg[y1] + lg[y2]) % n)
    return out

def W_all(k, n, ex, TR, nl, shift):
    # W(c) for c = a**j, j = 0..n-1, with kappa = a**shift
    hist = [0] * n
    for e in nl: hist[e] += 1
    res = []
    for j in range(n):
        s = 0
        for e in range(n):
            if hist[e]:
                s += hist[e] * (1 - 2 * TR[ex[(j + shift - e) % n]])
        res.append(s)
    return res

def M2_direct(k, n, ex, TR):
    s = [TR[ex[i]] for i in range(n)]
    N = 1 << k
    out = []
    for j in range(n):                       # c = a**j ; S_c(i) = s[i] xor s[j - i]
        seq = [XOR(s[i], s[(j - i) % n]) for i in range(n)]
        cnt = [0] * N
        w = 0
        for t in range(k): w |= seq[t % n] << t
        for i in range(n):
            cnt[w] += 1
            w = (w >> 1) | (seq[(i + k) % n] << (k - 1))
        out.append(sum(c * (c - 1) for c in cnt))   # = 2**k * M2
    return out

say("== (1) M2 formula, all c ==")
for k in range(KMIN, KMAX_DIRECT + 1):
    t = time.time()
    n, ex, lg, TR, mul = field(k)
    nl = norm_logs(k, n, ex, lg)
    D = M2_direct(k, n, ex, TR)
    target = [d - (2 ** k - 1) for d in D]        # should equal W(c)
    found = None
    for sh in range(n):
        Wv = W_all(k, n, ex, TR, nl, sh) if k <= 6 else None
        if Wv is None: break
        if Wv == target: found = sh; break
    if found is None:
        # faster search: test a few c first
        hist = [0] * n
        for e in nl: hist[e] += 1
        def Wc(j, sh):
            return sum(hist[e] * (1 - 2 * TR[ex[(j + sh - e) % n]]) for e in range(n) if hist[e])
        for sh in range(n):
            if all(Wc(j, sh) == target[j] for j in range(4)):
                if all(Wc(j, sh) == target[j] for j in range(n)):
                    found = sh; break
    ok = found is not None
    say("k=%d  identity M2 = 1-2**-k+2**-k W(c) for all %d values of c: %s (kappa=a**%s)  %.1fs" %
        (k, n, "OK" if ok else "FAIL", found, time.time() - t))
    if not ok: sys.exit("FAILED")
    if time.time() - T0 > TIME_BUDGET: sys.exit("time budget reached")

say("== (2) representation law ==")
def pmul(a, b):
    r = 0
    while b:
        if b & 1: r = XOR(r, a)
        a = a << 1; b = b >> 1
    return r
def dg(a): return a.bit_length() - 1
def pdivmod(a, m):
    q = 0; dm = dg(m)
    while a and dg(a) >= dm:
        sh = dg(a) - dm; q |= 1 << sh; a = XOR(a, m << sh)
    return q, a
def normpoly(g):
    d = dg(g); co = [0] * (2 * d + 1)
    for i in range(d + 1):
        if (g >> i) & 1:
            for j in range(d + 1):
                if (g >> j) & 1: co[i - j + d] = XOR(co[i - j + d], 1)
    Dm = [0, 0b10]
    for j in range(1, d): Dm.append(XOR(pmul(0b10, Dm[j]), Dm[j - 1]))
    P = co[d]
    for j in range(1, d + 1):
        if co[d + j]: P = XOR(P, Dm[j])
    return P
def chi0(D):
    if not (D & 1): return 0
    return -1 if (D >> 1) & 1 else 1
for k in range(3, KMAX_LAW + 1):
    t = time.time()
    cnt = {}
    for g in range(1, 1 << k):
        P = normpoly(g); cnt[P] = cnt.get(P, 0) + 1
    ok = True
    for P in range(1, 1 << k):
        r0 = 0
        for D in range(1, 1 << (dg(P) + 1)):
            q, rem = pdivmod(P, D)
            if rem == 0: r0 += chi0(D)
        if cnt.get(P, 0) != (k - dg(P)) * r0: ok = False; break
    say("k=%d  r(P) = (k-deg P)(1*chi0)(P) for all %d P: %s  %.1fs" % (k, (1 << k) - 1, "OK" if ok else "FAIL", time.time() - t))
    if not ok: sys.exit("FAILED")

say("== (3) worst case over c of |M2(k,c)-1| ==")
for k in range(KMIN, KMAX_TABLE + 1):
    # predicted cost ~ 2**k * k (Walsh transform) ; skip if budget exceeded
    if time.time() - T0 > TIME_BUDGET: say("time budget reached, stop"); break
    t = time.time()
    n, ex, lg, TR, mul = field(k)
    N = 1 << k
    # h(v) = number of h with 1/Nm(h) = v ; W(c) = sum_v h(v) chi(kappa c v): Walsh over trace pairing
    inv_hist = [0] * N
    for e in norm_logs(k, n, ex, lg): inv_hist[ex[(-e) % n]] += 1
    # Walsh-Hadamard in the coordinate basis gives all values sum_v h(v)(-1)**{<y,v>};
    # every nonzero linear functional v -> Tr(c v) is one <y,.> with y != 0, so the max over y != 0
    # equals the max over c != 0 (kappa is a fixed nonzero constant).
    f = inv_hist[:]
    L = 1
    while L < N:
        for i in range(0, N, 2 * L):
            for j in range(i, i + L):
                a = f[j]; b = f[j + L]; f[j] = a + b; f[j + L] = a - b
        L *= 2
    Wmax = max(abs(x) for x in f[1:])
    dev = max(abs(x - 1) for x in f[1:])          # M2 - 1 = (W - 1) / 2**k
    say("k=%2d  max_c |M2-1| = %.6f   (max|W| = %d, max|W|/2**(k/2) = %.2f)  %.1fs" %
        (k, dev / float(N), Wmax, Wmax / float(2 ** (k / 2.0)), time.time() - t))

say("== (4) duality: columns (N_w(S_c))_c  versus fibres of Phi_delta(y) = y psi(y) + delta y ==")
for k in range(KMIN, int(8) + 1):
    t = time.time()
    n, ex, lg, TR, mul = field(k); N = 1 << k
    sq = [TR[ex[i]] for i in range(n)]
    table = [[0] * n for _ in range(N)]           # table[w][j] = N_w(S_c), c = a**j
    for j in range(n):
        seq = [XOR(sq[i], sq[(j - i) % n]) for i in range(n)]
        w = 0
        for tt in range(k): w |= seq[tt % n] << tt
        for i in range(n):
            table[w][j] += 1
            w = (w >> 1) | (seq[(i + k) % n] << (k - 1))
    cols = sorted(tuple(sorted(r)) for r in table)
    psi = [0] * N
    for g in range(N):
        y1 = 0; y2 = 0
        for i in range(k):
            if (g >> i) & 1:
                y1 = XOR(y1, ex[i % n]); y2 = XOR(y2, ex[(-i) % n])
        psi[y1] = y2
    fib = []
    for d in range(N):
        h = [0] * N
        for y in range(N):
            h[XOR(mul(y, psi[y]), mul(d, y))] += 1
        fib.append(tuple(sorted(h[1:])))
    ok = cols == sorted(fib)
    say("k=%d  multiset of columns == multiset of fibre vectors: %s  %.1fs" % (k, "OK" if ok else "FAIL", time.time() - t))
    if not ok: sys.exit("FAILED")

say("== (5) pair count Sigma_0 = #{(a,b) != 0 : a psi(b) = b psi(a)} <= 6 k**2 2**k ==")
for k in range(3, int(8) + 1):
    n, ex, lg, TR, mul = field(k); N = 1 << k
    psi = [0] * N
    for g in range(N):
        y1 = 0; y2 = 0
        for i in range(k):
            if (g >> i) & 1:
                y1 = XOR(y1, ex[i % n]); y2 = XOR(y2, ex[(-i) % n])
        psi[y1] = y2
    # count via ratio classes rho(a) = psi(a)/a
    cl = {}
    for a in range(1, N):
        r = (lg[psi[a]] - lg[a]) % n
        cl[r] = cl.get(r, 0) + 1
    S0 = sum(v * v for v in cl.values())
    say("k=%d  Sigma_0 = %d   bound 6k**2 2**k = %d   %s" % (k, S0, 6 * k * k * N, "OK" if S0 <= 6 * k * k * N else "FAIL"))
    if S0 > 6 * k * k * N: sys.exit("FAILED")
say("ALL CHECKS PASSED  total %.1fs" % (time.time() - T0))