# -*- coding: utf-8 -*-
# =============================================================================
# VERIF_HATADA2.py
# Author: David Niedbala Giraudin (ORCID 0009-0009-1526-1178)
# Technical assistance: Claude (Anthropic)
# Ancillary file to: "Hecke eigenvalues in level 1 modulo powers of 2 and 3:
# Hatada's question (1981) and the Coleman-Stein conjecture",
# Section 8 (a 2-adic Galois law) and Section "Verification".
# -----------------------------------------------------------------------------
# Statements checked, for every eigenform f of level 1 and weight k,
# s_p = 1 + p**(k-1), d_3 = a_3 - s_3:
#   (Lemma 13) T_5 - (1 + 5 Lam_5) == 0 (mod 512) on Hatada's lattice L_k, for
#             every even k: a_5 == s_5 (mod 512).  Certificate of 190 terms.
#   (Thm U)   p == 1 (mod 8): a_p == s_p (mod 2**10);  p == 5 (mod 8): (mod 2**8)
#             (part D checks 2**9, the depth given by Theorem X);
#             p == 3 (mod 8): a_p - s_p == d_3 (mod 2**8);
#             p == 7 (mod 16): a_p - s_p == 2 d_3 (mod 2**9); p == 15 (mod 16): (mod 2**9).
#   (Thm X)   p == -1 (mod 32): a_p == s_p (mod 2**10), hence Hatada's question for a = 10;
#             p == 5 (mod 8): a_p == s_p (mod 2**9).
#   (Cor Z)   v_2(a_p) = v_2(p + 1) if v_2(p + 1) <= 10: a_p = 0 implies p == -1 (mod 2**11)
#             (2**12 with part F, heavier, disabled by default).
# RUN WITH PLAIN PYTHON 3. No library. No input needed.
# Running time: under a minute. Memory: a few tens of MB. Progressive output.
#
# Part A: the certificate (criterion of Lemma 12, the analogue of Lemma 7 for
#         p = 2, w = 2t) modulo 512; witness: the target alone fails the criterion.
# Part B: direct check, action of the remainder on forms of degree w <= W_B_MAX,
#         with witnesses of odd degree.
# Part C: the Kummer symbols used in the proof (Euler criteria in F_9 and F_25)
#         for Q(sqrt 2) and Q(zeta_8).
# Part D: the five laws of Theorem U at the level of eigenvalues, on L_k,
#         12 <= k <= KMAX_D (2-adic Newton polygons), with witnesses.
# Part E: the finite group generated by c and t in GL_2((Z/2**n)[eps]) times the
#         units of Z/2**m (Lemma 15): exhaustive enumeration and check of the
#         guaranteed depths (moduli n = 16, m = 64, then n = 32, m = 512, for
#         Theorem X and Corollary Z).
# Part F: (optional, PARTIE_F = False) n = 64, m = 2048: sieve 2**12; about
#         500 000 elements, several hundred MB: do NOT run on a machine with
#         limited memory.
# =============================================================================
import sys
import time
from fractions import Fraction
from math import comb, gcd

W_B_MAX = int(40)
KMAX_D = int(40)
BUDGET = float(480.0)
PARTIE_F = False
T0 = time.time()

def dire(*a):
    print(*a, flush=True)

def v3(n):
    n = abs(int(n))
    if n == 0:
        return 10 ** 9
    v = 0
    while n % 3 == 0:
        n //= 3
        v += 1
    return v

def manin(n):
    out = []
    for a in range(1, n + 1):
        for d in range(1, n + 1):
            if a * d < n:
                continue
            for c in range(0, a):
                for b in range(0, d):
                    if a * d - b * c == n:
                        out.append((a, b, c, d))
    return out

def binom_entier(z, k):
    num = 1
    for t in range(k):
        num *= (z - t)
    f = 1
    for t in range(2, k + 1):
        f *= t
    return num // f

def image_monome(w, j, a, b, c, d, N=None):
    """coefficients of F(ax + cy, bx + dy) for F = x**j y**(w-j)"""
    A = [comb(j, s) * (a ** s) * (c ** (j - s)) for s in range(j + 1)]
    B = [comb(w - j, t) * (b ** t) * (d ** (w - j - t)) for t in range(w - j + 1)]
    out = [0] * (w + 1)
    for s, x in enumerate(A):
        if x:
            for t, y in enumerate(B):
                if y:
                    out[s + t] += x * y
    if N is not None:
        out = [x % N for x in out]
    return out

def action_directe(X, f, w, N):
    tot = [0] * (w + 1)
    for g, cg in X.items():
        a, b, c, d = g
        vp = [[1] + [0] * w]
        for m in range(1, w + 1):
            p = vp[-1]
            vp.append([(d * p[0]) % N] + [(b * p[e - 1] + d * p[e]) % N for e in range(1, w + 1)])
        Rr = [f[w] % N] + [0] * w
        for i in range(w - 1, -1, -1):
            V = vp[w - i]
            fi = f[i]
            Rr = [(c * Rr[0] + fi * V[0]) % N] + [(a * Rr[e - 1] + c * Rr[e] + fi * V[e]) % N for e in range(1, w + 1)]
        for e in range(w + 1):
            tot[e] = (tot[e] + cg * Rr[e]) % N
    return tot

def lignes_relations(w):
    rows = []
    for j in range(w + 1):
        imS = image_monome(w, j, 0, -1, 1, 0)
        imU = image_monome(w, j, 1, -1, 1, 0)
        imU2 = image_monome(w, j, 0, -1, 1, -1)
        imE = image_monome(w, j, -1, 0, 0, 1)
        rS = imS[:]
        rS[j] += 1
        rU = [imU[i] + imU2[i] for i in range(w + 1)]
        rU[j] += 1
        rE = imE[:]
        rE[j] += 1
        rows += [rS, rU, rE]
    return rows

def noyau_rationnel(rows, n):
    M = [[Fraction(x) for x in r] for r in rows if any(r)]
    piv_cols = []
    r = 0
    for c in range(n):
        p = None
        for i in range(r, len(M)):
            if M[i][c] != 0:
                p = i
                break
        if p is None:
            continue
        M[r], M[p] = M[p], M[r]
        inv = 1 / M[r][c]
        M[r] = [x * inv for x in M[r]]
        for i in range(len(M)):
            if i != r and M[i][c] != 0:
                f = M[i][c]
                M[i] = [M[i][k] - f * M[r][k] for k in range(n)]
        piv_cols.append(c)
        r += 1
        if r == len(M):
            break
    libres = [c for c in range(n) if c not in piv_cols]
    base = []
    for fcol in libres:
        v = [Fraction(0)] * n
        v[fcol] = Fraction(1)
        for i, pc in enumerate(piv_cols):
            v[pc] = -M[i][fcol]
        den = 1
        for x in v:
            den = den * x.denominator // gcd(den, x.denominator)
        vi = [int(x * den) for x in v]
        g = 0
        for x in vi:
            g = gcd(g, x)
        base.append([x // g for x in vi])
    return base

def dependance_mod3(base):
    d = len(base)
    n = len(base[0])
    M = [[base[i][k] % 3 for k in range(n)] + [1 if t == i else 0 for t in range(d)] for i in range(d)]
    r = 0
    for c in range(n):
        p = None
        for i in range(r, d):
            if M[i][c] % 3:
                p = i
                break
        if p is None:
            continue
        M[r], M[p] = M[p], M[r]
        inv = 1 if M[r][c] == 1 else 2
        M[r] = [(x * inv) % 3 for x in M[r]]
        for i in range(d):
            if i != r and M[i][c] % 3:
                f = M[i][c]
                M[i] = [(M[i][k] - f * M[r][k]) % 3 for k in range(n + d)]
        r += 1
    for i in range(r, d):
        if all(M[i][k] % 3 == 0 for k in range(n)):
            return M[i][n:]
    return None

def sature_en_3(base):
    base = [b[:] for b in base]
    while True:
        c = dependance_mod3(base)
        if c is None:
            return base
        i0 = next(i for i in range(len(c)) if c[i] % 3)
        v = [sum(c[i] * base[i][k] for i in range(len(base))) for k in range(len(base[0]))]
        base[i0] = [x // 3 for x in v]

def reseau_hatada(w):
    return sature_en_3(noyau_rationnel(lignes_relations(w), w + 1))

def matrice_Tn(w, n):
    """rows: m_l(T_n P) = sum_i Mt[l][i] m_i"""
    Mt = [[0] * (w + 1) for _ in range(w + 1)]
    for (a, b, c, d) in manin(n):
        for l in range(w + 1):
            im = image_monome(w, l, a, b, c, d)
            for i, x in enumerate(im):
                if x:
                    Mt[l][i] += x
    return Mt

def coordonnees(base, vecs):
    """expresses each vector of vecs (in the Z-module spanned by base) in coordinates on base"""
    d = len(base)
    n = len(base[0])
    # choose d columns on which the base is invertible (rational elimination)
    M = [[Fraction(base[i][k]) for k in range(n)] for i in range(d)]
    piv = []
    r = 0
    for c in range(n):
        p = None
        for i in range(r, d):
            if M[i][c] != 0:
                p = i
                break
        if p is None:
            continue
        M[r], M[p] = M[p], M[r]
        inv = 1 / M[r][c]
        M[r] = [x * inv for x in M[r]]
        for i in range(d):
            if i != r and M[i][c] != 0:
                f = M[i][c]
                M[i] = [M[i][k] - f * M[r][k] for k in range(n)]
        piv.append(c)
        r += 1
        if r == d:
            break
    # solve sum_i x_i base[i][piv] = vec[piv]
    A = [[Fraction(base[i][c]) for i in range(d)] for c in piv]   # d x d
    out = []
    for v in vecs:
        rhs = [Fraction(v[c]) for c in piv]
        Mx = [A[i][:] + [rhs[i]] for i in range(d)]
        for c in range(d):
            p = next(i for i in range(c, d) if Mx[i][c] != 0)
            Mx[c], Mx[p] = Mx[p], Mx[c]
            inv = 1 / Mx[c][c]
            Mx[c] = [x * inv for x in Mx[c]]
            for i in range(d):
                if i != c and Mx[i][c] != 0:
                    f = Mx[i][c]
                    Mx[i] = [Mx[i][k] - f * Mx[c][k] for k in range(d + 1)]
        x = [Mx[i][d] for i in range(d)]
        for xi in x:
            if xi.denominator % 3 == 0:
                raise ValueError('coordinate not 3-integral')
        # full check (the base is saturated at 3 only: coordinates in Z_(3))
        for k in range(n):
            if sum(x[i] * base[i][k] for i in range(d)) != v[k]:
                raise ValueError('vector outside the lattice')
        out.append(x)
    return out

def polcar(A):
    """characteristic polynomial (Faddeev-LeVerrier), coefficients c_0..c_d in Z_(3), monic"""
    d = len(A)
    Mk = [[Fraction(0)] * d for _ in range(d)]
    c = [Fraction(0)] * (d + 1)
    c[d] = Fraction(1)
    for k in range(1, d + 1):
        # M_k = A M_{k-1} + c_{d-k+1} I
        Mk = [[sum(A[i][t] * Mk[t][j] for t in range(d)) + (c[d - k + 1] if i == j else 0) for j in range(d)] for i in range(d)]
        tr = sum(sum(A[i][t] * Mk[t][i] for t in range(d)) for i in range(d))
        c[d - k] = -tr / k
    return c

def sigma_k(n, e):
    return sum(d ** e for d in range(1, n + 1) if n % d == 0)


# =============================================================================
# Ring of matrices modulo 2**n, ideal I, criterion of Lemma 12 (Lemma 7 for p = 2)
# =============================================================================
def mulmod(g, h, N):
    a, b, c, d = g
    e, f, p, q = h
    return ((a * e + b * p) % N, (a * f + b * q) % N, (c * e + d * p) % N, (c * f + d * q) % N)

def elt(g, N, c=1):
    return {tuple(int(t) % N for t in g): int(c) % N}

def somme(N, *termes):
    out = {}
    for c, X in termes:
        for g, x in X.items():
            out[g] = (out.get(g, 0) + c * x) % N
    return {g: x for g, x in out.items() if x}

def produit(X, Y, N):
    out = {}
    for g, x in X.items():
        for h, y in Y.items():
            k = mulmod(g, h, N)
            out[k] = (out.get(k, 0) + x * y) % N
    return {g: x for g, x in out.items() if x}

REPS2 = []
for _a in range(2):
    for _b in range(2):
        for _c in range(2):
            for _d in range(2):
                if (_a * _d - _b * _c) % 2:
                    REPS2.append((_a, _b, _c, _d))
CL2 = {g: i for i, g in enumerate(REPS2)}
KAPPA2 = [(3, 0, 0, 1), (1, 2, 0, 1), (1, 0, 2, 1), (1, 0, 0, 3)]

class Anneau2:
    def __init__(self, n):
        self.n = n
        self.N = N = 2 ** n
        self.un = elt((1, 0, 0, 1), N)
        U = (1, N - 1, 1, 0)
        self.gens = {'S': somme(N, (1, self.un), (1, elt((0, -1, 1, 0), N))),
                     'U': somme(N, (1, self.un), (1, elt(U, N)), (1, elt(mulmod(U, U, N), N))),
                     'E': somme(N, (1, self.un), (1, elt((-1, 0, 0, 1), N)))}
        self.bfac = [somme(N, (1, elt(k, N)), (-1, self.un)) for k in KAPPA2]
        self.repinv = []
        for (a, b, c, d) in REPS2:
            di = pow((a * d - b * c) % N, -1, N)
            self.repinv.append(((d * di) % N, (-b * di) % N, (-c * di) % N, (a * di) % N))

    def zeta(self, texte):
        acc = {}
        toks = ''.join(texte).split()
        for tok in toks:
            h = int(tok[1])
            al = [int(tok[2 + t]) for t in range(4)]
            c = int(tok[6:9])
            X = elt(REPS2[h], self.N)
            for t in range(4):
                for _ in range(al[t]):
                    X = produit(X, self.bfac[t], self.N)
            for g, x in produit(self.gens[tok[0]], X, self.N).items():
                acc[g] = (acc.get(g, 0) + c * x) % self.N
        return {g: x for g, x in acc.items() if x}, len(toks)

    def coord_kappa(self, k):
        p, q, r, s = k
        n, N = self.n, self.N
        out = {}
        for a in range(n):
            for b in range(n - a):
                for m in range(n - a - b):
                    for m2 in range(n - a - b - m):
                        val = (r ** a) * (q ** b) * ((p - 1) ** m) * ((s - 1) ** m2) * comb(a + m, a) * comb(b + m2, b)
                        cle = (b - a, a + m, b + m2)
                        out[cle] = (out.get(cle, 0) + val) % N
        return out

    def decompose(self, g):
        i = CL2[tuple(t % 2 for t in g)]
        k = mulmod(self.repinv[i], g, self.N)
        if not (k[0] % 2 == 1 and k[1] % 2 == 0 and k[2] % 2 == 0 and k[3] % 2 == 1):
            raise ValueError('decomposition impossible')
        return i, k

    def critere_nul(self, X):
        """Lemma 12 (Lemma 7 for p = 2): Mahler coefficients (in i and t, j = 2t - i, total degree < n)
        zero modulo 2**n; then X acts by 0 modulo 2**n on the forms of even degree w."""
        n, N = self.n, self.N
        blocs = {}
        cache = {}
        for g, x in X.items():
            h, k = self.decompose(g)
            if k not in cache:
                cache[k] = self.coord_kappa(k)
            for (dl, A, B), v in cache[k].items():
                bl = blocs.setdefault((h, dl), {})
                bl[(A, B)] = (bl.get((A, B), 0) + x * v) % N
        for bl in blocs.values():
            vals = [[sum(c * comb(i, A) * binom_entier(2 * t - i, B) for (A, B), c in bl.items())
                     for t in range(n)] for i in range(n)]
            for A in range(n):
                for B in range(n - A):
                    tot = 0
                    for u in range(A + 1):
                        for v in range(B + 1):
                            tot += (-1) ** (A - u + B - v) * comb(A, u) * comb(B, v) * vals[u][v]
                    if tot % N:
                        return False
        return True

def cible_T5(R):
    N = R.N
    T = somme(N, *[(1, elt(g, N)) for g in manin(5)])
    return somme(N, (1, T), (-1, R.un), (-5, elt((5, 0, 0, 5), N)))

CERT_T5 = (
    "S00000377 S00001020 S00002002 S00010165 S00011012 S00020006 S00021031 "
    "S00100071 S00110015 S00111009 S00120010 S00121004 S00130001 S00131001 "
    "S00200020 S00201021 S00210020 S00211002 S00220001 S00221001 S00300022 "
    "S00301002 S00310001 S00400001 S00500001 S01000053 S01001002 S01010031 "
    "S01011002 S01020003 S01100031 S01110002 S01200005 S02000001 S02100001 "
    "S10000349 S10001099 S10003001 S10010104 S10011045 S10020039 S10021031 "
    "S10100238 S10101002 S10110051 S10111020 S10120024 S10121006 S10130007 "
    "S10200023 S10201023 S10210007 S10211001 S10220003 S10300004 S10301001 "
    "S10310002 S10320001 S10400001 S11000008 S11001002 S11010011 S11100010 "
    "S11110003 S11300001 S12000001 S12100001 S20000004 S20001006 S20010010 "
    "S20011002 S20020002 S20100004 S20101002 S20200002 S30000356 S30001108 "
    "S30002002 S30010026 S30011019 S30012001 S30020035 S30021029 S30030005 "
    "S30031003 S30040003 S30100242 S30101060 S30102001 S30110002 S30111012 "
    "S30120001 S30121006 S30130006 S30131001 S30200022 S30201020 S30210020 "
    "S30211003 S30220007 S30221001 S30300017 S30301003 S30310002 S30400001 "
    "S30410001 S30500001 S31000042 S31010009 S31011004 S31020003 S31030001 "
    "S31100030 S31110001 S31111001 S31200006 S31300001 S32000001 S40000032 "
    "S40001012 S40010004 S40100014 S40110002 S40200002 S50000025 S50001009 "
    "S50010012 S50020001 S50100012 S50110001 U00000133 U00001005 U00002002 "
    "U00010208 U00011051 U00020057 U00021001 U00100038 U00101040 U00110045 "
    "U00111009 U00120027 U00121007 U00130007 U00200008 U00201024 U00210021 "
    "U00211003 U00220005 U00300027 U00310003 U00400002 U00401001 U00410001 "
    "U00500001 U01000032 U01001002 U01010029 U01100022 U01300001 U02000003 "
    "U10000084 U10001009 U10010011 U10011019 U10020025 U10021001 U10100033 "
    "U10101010 U10110026 U10111002 U10120002 U10121002 U10130003 U10200028 "
    "U10201001 U10210002 U10300007 U10301001 U10310001 U10400001 U11000003 "
    "U11100007 U11200003 U20000048 U20001008 U20010010 U20100006 U30001002 "
    "U30010004 "
)

# =============================================================================
# Part A: the certificate
# =============================================================================
def partie_A():
    dire('')
    dire('Part A: certificate for T_5 - (1 + 5 Lam_5) modulo 512 (criterion of Lemma 12, p = 2, w even)')
    t1 = time.time()
    R = Anneau2(9)
    zeta, nt = R.zeta(CERT_T5)
    Xi = cible_T5(R)
    reste = somme(R.N, (1, Xi), (-1, zeta))
    ok = R.critere_nul(reste)
    temoin = not R.critere_nul(Xi)
    dire('  %d terms, remainder of %d matrices; criterion: %s; witness (target alone fails): %s  (%.1f s)'
         % (nt, len(reste), 'VERIFIED' if ok else 'FAILED', 'yes' if temoin else 'NO', time.time() - t1))
    return ok and temoin, reste

# =============================================================================
# Part B: direct action of the remainder
# =============================================================================
def partie_B(reste):
    dire('')
    dire('Part B: direct action of the remainder on the forms of degree w <= %d (modulo 512)' % W_B_MAX)
    N = 512
    ok = True
    res = []
    for w in list(range(2, W_B_MAX + 1, 2)) + [3, 7, 11]:
        f = [(7 * i * i + 3 * i + 5) % N for i in range(w + 1)]
        g = [(i * i * i + 11 * i + 1) % N for i in range(w + 1)]
        nul = not any(action_directe(reste, f, w, N)) and not any(action_directe(reste, g, w, N))
        if w % 2 == 0:
            ok = ok and nul
            res.append('w=%d:%s' % (w, 'yes' if nul else 'NO'))
        else:
            res.append('w=%d(odd):%s' % (w, 'no' if not nul else 'YES?'))
            ok = ok and not nul
    dire('  ' + ' '.join(res))
    return ok

# =============================================================================
# Part C: Kummer symbols (Euler's criteria)
# =============================================================================
def fp2_mul(x, y, p, d):
    """F_p[s]/(s**2 - d) : (a + b s)(c + e s)"""
    return ((x[0] * y[0] + d * x[1] * y[1]) % p, (x[0] * y[1] + x[1] * y[0]) % p)

def fp2_pow(x, e, p, d):
    r = (1, 0)
    while e:
        if e & 1:
            r = fp2_mul(r, x, p, d)
        x = fp2_mul(x, x, p, d)
        e >>= 1
    return r

def symbole(x, p, d):
    """+1 if x is a square in F_{p**2} = F_p[s]/(s**2 - d), -1 otherwise (Euler)"""
    y = fp2_pow(x, (p * p - 1) // 2, p, d)
    return 1 if y == (1, 0) else (-1 if y == (p - 1, 0) else 0)

def racine8(p, d):
    for a in range(p):
        for b in range(p):
            z = (a, b)
            if fp2_pow(z, 8, p, d) == (1, 0) and fp2_pow(z, 4, p, d) != (1, 0):
                return z
    return None

def partie_C():
    dire('')
    dire('Part C: Kummer symbols (Frobenius at the primes above 3 and 5)')
    ok = True
    # Q(sqrt 2): F_9 = F_3[s]/(s**2 - 2), F_25 = F_5[s]/(s**2 - 2); s = sqrt 2
    v3_ = [symbole(x, 3, 2) for x in [(2, 0), (1, 1), (0, 1)]]      # -1, 1 + sqrt2, sqrt2
    v5_ = [symbole(x, 5, 2) for x in [(4, 0), (1, 1), (0, 1)]]
    ok = ok and v3_ == [1, -1, 1] and v5_ == [1, 1, -1]
    dire('  Q(sqrt2), Frob_3**2: (-1, 1+sqrt2, sqrt2) -> %s (expected [1, -1, 1], image (0,1,0))' % v3_)
    dire('  Q(sqrt2), Frob_5**2: (-1, 1+sqrt2, sqrt2) -> %s (expected [1, 1, -1], image (0,0,1))' % v5_)
    # Q(zeta_8): residue fields F_9 = F_3[s]/(s**2 + 1) and F_25 = F_5[s]/(s**2 - 2)
    for p, d, att in [(3, -1, [-1, -1]), (5, 2, [-1, 1])]:
        z = racine8(p, d)
        zi = fp2_pow(z, 7, p, d)
        r2 = ((z[0] + zi[0]) % p, (z[1] + zi[1]) % p)                # sqrt 2 = zeta + zeta**-1
        if fp2_mul(r2, r2, p, d) != (2 % p, 0):
            ok = False
        eps = ((1 + r2[0]) % p, r2[1])
        vals = [symbole(z, p, d), symbole(eps, p, d)]
        ok = ok and vals == att
        dire('  Q(zeta8), Frobenius above %d: (zeta8, 1+sqrt2) -> %s (expected %s)' % (p, vals, att))
    return ok

# =============================================================================
# Part D: the five laws on L_k (eigenvalues)
# =============================================================================
def v2(n):
    n = abs(int(n))
    if n == 0:
        return 10 ** 9
    v = 0
    while n % 2 == 0:
        n //= 2
        v += 1
    return v

def v2frac(x):
    x = Fraction(x)
    if x == 0:
        return 10 ** 9
    return v2(x.numerator) - v2(x.denominator)

def pentes2_ok(c, seuil):
    """all roots of the monic polynomial sum c_i x**i have 2-adic valuation >= seuil"""
    d = len(c) - 1
    for i in range(d):
        if c[i] != 0 and v2frac(c[i]) < seuil * (d - i):
            return False
    return True

def matrice_sur_reseau(base, w, n):
    Mt = matrice_Tn(w, n)
    images = [[sum(r[i] * m[i] for i in range(w + 1) if r[i] and m[i]) for r in Mt] for m in base]
    A = coordonnees(base, images)
    d = len(base)
    return [[A[j][i] for j in range(d)] for i in range(d)]

def mcomb(*termes):
    d = None
    for c, M in termes:
        if M is not None:
            d = len(M)
    out = [[Fraction(0)] * d for _ in range(d)]
    for c, M in termes:
        for i in range(d):
            for j in range(d):
                out[i][j] += c * (M[i][j] if M is not None else (1 if i == j else 0))
    return out

PREMIERS_D = [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]

def partie_D():
    dire('')
    dire('Part D: laws of Theorem U on L_k, 12 <= k <= %d, p in %s' % (KMAX_D, PREMIERS_D))
    echecs = 0
    temoins = 0
    for k in range(12, KMAX_D + 1, 2):
        if k == 14:
            continue
        t1 = time.time()
        w = k - 2
        base = reseau_hatada(w)
        s = lambda l: 1 + l ** (k - 1)
        A = {p: matrice_sur_reseau(base, w, p) for p in PREMIERS_D}
        D3 = mcomb((1, A[3]), (-s(3), None))
        ok = True
        for p in PREMIERS_D:
            X = mcomb((1, A[p]), (-s(p), None))
            if p % 8 == 1:
                ok = ok and pentes2_ok(polcar(X), 10)
            elif p % 8 == 5:
                ok = ok and pentes2_ok(polcar(X), 9)
            elif p % 8 == 3:
                ok = ok and pentes2_ok(polcar(mcomb((1, X), (-1, D3))), 8)
            elif p % 16 == 7:
                ok = ok and pentes2_ok(polcar(mcomb((1, X), (-2, D3))), 9)
                if not pentes2_ok(polcar(X), 9):
                    temoins += 1
            else:
                ok = ok and pentes2_ok(polcar(X), 9)
        if not ok:
            echecs += 1
        dire('  k=%3d d=%2d : %s  (%.1f s)' % (k, len(base), 'ok' if ok else 'FAILED', time.time() - t1))
        if time.time() - T0 > BUDGET:
            dire('  budget reached during k = %d' % k)
            break
    dire('  weights failing: %d (expected: 0); cases p == 7 (mod 16) where the term 2 d_3 is needed: %d (expected: > 0)' % (echecs, temoins))
    return echecs == 0 and temoins > 0


# =============================================================================
# Parts E and F: the universal family and the finite group
# =============================================================================
def mul_eps(g, h, N, M):
    (a1, b1, c1, d1, e1, f1, g1, h1, x1), (a2, b2, c2, d2, e2, f2, g2, h2, x2) = g, h
    return ((a1 * a2 + b1 * c2) % N, (a1 * b2 + b1 * d2) % N, (c1 * a2 + d1 * c2) % N, (c1 * b2 + d1 * d2) % N,
            (a1 * e2 + b1 * g2 + e1 * a2 + f1 * c2) % N, (a1 * f2 + b1 * h2 + e1 * b2 + f1 * d2) % N,
            (c1 * e2 + d1 * g2 + g1 * a2 + h1 * c2) % N, (c1 * f2 + d1 * h2 + g1 * b2 + h1 * d2) % N,
            (x1 * x2) % M)

def groupe_fini(kappa, direction, N, M):
    """image of <c, t>: c -> (diag(1,-1), 0, -1), t -> (A_t, B_t, 5), A_t = [[1,1],[0,D]],
    B_t = [[1,0],[D+1,1]] (direction x = 2**8 X) or [[1,0],[D-1,-1]] (direction y = 2**7 Y)"""
    D = pow(5, kappa, N * M)
    c = (1, 0, 0, N - 1, 0, 0, 0, 0, M - 1)
    if direction == 'X':
        Bt = (1, 0, (D + 1) % N, 1)
    else:
        Bt = (1, 0, (D - 1) % N, N - 1)
    t = (1, 1, 0, D % N) + Bt + (5 % M,)
    vus = {c, t}
    pile = [c, t]
    while pile:
        g = pile.pop()
        for s in (c, t):
            h = mul_eps(g, s, N, M)
            if h not in vus:
                vus.add(h)
                pile.append(h)
    return vus

def profondeurs(N, M):
    """guaranteed depth of Delta (exponent of 2) for each class of chi, over all classes of kappa"""
    prof = {}
    ordre = N // 4                      # order of 5 modulo N: D = 5**kappa depends only on kappa mod N/4
    for kappa in range(1, ordre, 2):
        for direction, echelle in (('X', 7), ('Y', 6)):
            H = groupe_fini(kappa, direction, N, M)
            for g in H:
                ch = g[8]
                if ch % 4 == 3:
                    cle = ('-', min(v2((ch + 1) % M), v2(M)))
                elif ch % 8 == 1:
                    cle = ('+', 3)
                else:
                    cle = ('5', 0)
                tr = (g[4] + g[7]) % N
                d = echelle + (v2(tr) if tr else v2(N))
                prof[cle] = min(prof.get(cle, 99), d)
    return prof

def partie_E():
    dire('')
    dire('Part E: finite group of the universal family (Lemma 15)')
    t1 = time.time()
    ok = True
    temoin = False
    for kappa in (1, 3):                # N = 16: D depends only on kappa mod 4
        for direction, mod in (('X', 8), ('Y', 16)):
            H = groupe_fini(kappa, direction, 16, 64)
            c1 = all((g[4] + g[7]) % mod == 0 for g in H if g[8] % 8 == 1)
            c2 = all((g[4] + g[7]) % mod == 0 for g in H if g[8] % 32 == 31)
            if any((g[4] + g[7]) % mod != 0 for g in H if g[8] % 32 == 15):
                temoin = True
            dire('  kappa == %d (mod 4), direction %s: |H| = %d; chi == 1 (mod 8): %s; chi == -1 (mod 32): %s'
                 % (kappa, direction, len(H), 'zero mod 2**10' if c1 else 'FAILED', 'zero mod 2**10' if c2 else 'FAILED'))
            ok = ok and c1 and c2
    dire('  witness: chi == 15 (mod 32) is not always zero mod 2**10: %s' % ('yes' if temoin else 'NO'))
    # class 5, with t = Frob_5 and x = a_5 - s_5 divisible by 2**9 (Lemma 13): 2**8 tau_X + 2**6 tau_Y
    c5 = True
    for kappa in range(1, 8, 2):
        for direction, mod in (('X', 2), ('Y', 8)):
            H = groupe_fini(kappa, direction, 32, 512)
            if not all((g[4] + g[7]) % mod == 0 for g in H if g[8] % 8 == 5):
                c5 = False
    dire('  class 5 (x divisible by 2**9): Delta zero mod 2**9: %s' % ('yes' if c5 else 'NO'))
    ok = ok and c5
    prof = profondeurs(32, 512)
    crible = True
    for cle in prof:
        if cle[0] == '-' and cle[1] >= 2 and prof[cle] <= min(cle[1], 10):
            crible = False
    # Theorem X: depth >= min(m + 5, 11) on the class v_2(1 + chi) = m, and >= 10 on class 1
    bornes = all(prof[cle] >= min(cle[1] + 5, 11) for cle in prof if cle[0] == '-' and cle[1] >= 2)
    bornes = bornes and prof[('+', 3)] >= 10
    dire('  Theorem X: depth >= min(m + 5, 11) for v_2(1 + chi) = m >= 2, and >= 10 for chi == 1 (mod 8): %s'
         % ('yes' if bornes else 'NO'))
    crible = crible and bornes
    dire('  guaranteed depths (n = 32, m = 512): ' + ' '.join('%s%d:%d' % (k[0], k[1], prof[k]) for k in sorted(prof)))
    dire('  (class "-m": v_2(1 + chi) = m, the last bin collecting the m largest values; "+3": chi == 1 (mod 8))')
    dire('  sieve: v_2(a_p) = v_2(p + 1) for v_2(p + 1) <= 10: %s  (%.1f s)' % ('yes' if crible else 'NO', time.time() - t1))
    return ok and temoin and crible

def partie_F():
    dire('')
    dire('Part F (heavy): n = 64, m = 2048')
    prof = profondeurs(64, 2048)
    dire('  ' + ' '.join('%s%d:%d' % (k[0], k[1], prof[k]) for k in sorted(prof)))
    ok = all(prof[k] > min(k[1], 11) for k in prof if k[0] == '-' and k[1] >= 2)
    dire('  sieve 2**12: %s' % ('yes' if ok else 'NO'))
    return ok

# =============================================================================
# Self-test and summary
# =============================================================================
def auto_test():
    if pow(2, 8) != 256 or comb(10, 3) != 120 or v2(1216) != 6 or len(REPS2) != 6:
        return False
    base = reseau_hatada(10)
    if len(base) != 1:
        return False
    m = base[0]
    A = coordonnees(base, [[sum(r[i] * m[i] for i in range(11)) for r in matrice_Tn(10, 3)]])
    if A[0][0] != 252:                          # tau(3)
        return False
    A5 = coordonnees(base, [[sum(r[i] * m[i] for i in range(11)) for r in matrice_Tn(10, 5)]])
    if A5[0][0] != 4830:                        # tau(5)
        return False
    return v2(4830 - 1 - 5 ** 11) >= 8 and len(manin(2)) == 4 and len(manin(5)) == 15

dire('=' * 76)
dire('VERIF_HATADA2: 2-adic Galois law in level 1')
dire('Budget %.0f s; Part B: w <= %d; Part D: k <= %d' % (BUDGET, W_B_MAX, KMAX_D))
dire('=' * 76)
if not auto_test():
    dire('SELF-TEST FAILED: stop (environment not as expected).')
    sys.exit(1)
dire('Self-test: OK (period of Delta, tau(3) = 252, tau(5) = 4830, v2(tau(5) - 1 - 5**11) >= 8).')
rA, reste = partie_A()
rB = partie_B(reste)
rC = partie_C()
rD = partie_D()
rE = partie_E()
rF = partie_F() if PARTIE_F else True
dire('')
dire('SUMMARY: Part A (certificate) %s, Part B (direct action) %s, Part C (Kummer) %s, Part D (laws, k <= %d) %s, Part E (finite group) %s%s'
     % ('OK' if rA else 'FAILED', 'OK' if rB else 'FAILED', 'OK' if rC else 'FAILED', KMAX_D, 'OK' if rD else 'FAILED', 'OK' if rE else 'FAILED', (', Part F ' + ('OK' if rF else 'FAILED')) if PARTIE_F else ''))
dire('Total time: %.1f s' % (time.time() - T0))