# -*- coding: utf-8 -*-
# =============================================================================
# VERIF_INERTES.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 9 (beyond 729: the law of inert primes) and Section "Verification".
# -----------------------------------------------------------------------------
# Statements checked, for every eigenform f of level 1 and weight k,
# kappa = k - 1, s_l = 1 + l**kappa:
#   (Thm R')  l == 2 (mod 3), j = v_3(l + 1): a_l**2 == s_l**2 (mod 3**min(2j+2, 10)).
#   (Thm R'') l == 1 (mod 27), 3 a cube modulo l: a_l == s_l (mod 3**7);
#             if k == 4 (mod 6) and l = x**2 + x y + 547 y**2: (mod 3**8),
#             and (mod 3**9) if moreover l == 1 (mod 81).
#   (Cor Y')  k not == 4 (mod 6), l == 2 (mod 3), v_3(l + 1) <= 4: v_3(a_l) = v_3(l + 1).
# RUN WITH PLAIN PYTHON 3. No library. No input needed.
# Running time: under a minute. Memory: a few MB. Progressive output.
#
# Part A: the identities of Lemma 16 in the dual-number representation
#         (x -> rho(x), y -> C rho(x) C), on random words (fixed seed),
#         modulo 3**14, and the valuations of the system of Lemma 17; witnesses: F can
#         be a unit if u != 0; if w = 0, beta/3 is not a function of mu/3 on
#         ker v (whereas it is if w >= 1).
# Part B: the law of inert primes and Corollary Y' on Hatada's lattice L_k,
#         12 <= k <= KMAX_B, l in PREMIERS_B (3-adic Newton polygons).
# Part C: exact values a_l (q-expansions) for the six weights of dimension 1
#         (k = 12, 16, 18, 20, 22, 26) and all primes l < LMAX_C:
#         Theorems R', R'' and Corollary Y'; witnesses 971 (k = 18) and 271 (k = 16).
# Part D: Pic of the order of conductor 27 (resp. 9) of Z[omega]: cyclic group
#         of order 9 (resp. 3), and consistency with x**2 + x y + 547 y**2 (resp. 61 y**2)
#         for all primes l == 1 (mod 3), l < LMAX_D.
# =============================================================================
import time
from fractions import Fraction
from math import comb, gcd
import random

KMAX_B = int(40)
PREMIERS_B = [5, 11, 17, 53, 71, 107]
LMAX_C = int(3000)
LMAX_D = int(20000)
N_MOTS = int(2000)
BUDGET = float(480.0)
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):
    """express 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 where 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 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


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

def pente_min(c):
    """smallest 3-adic valuation of the roots of the monic polynomial sum c_i x**i"""
    d = len(c) - 1
    m = Fraction(10 ** 9)
    for i in range(d):
        if c[i] != 0:
            m = min(m, Fraction(v3f(c[i]), d - i))
    return m

def pentes_egales(c, j):
    """all roots have valuation exactly j"""
    d = len(c) - 1
    return v3f(c[0]) == d * j and pente_min(c) >= j

# =============================================================================
# Self-test
# =============================================================================
def tau_table(N):
    sig = [0] * (N + 1)
    for d in range(1, N + 1):
        for m in range(d, N + 1, d):
            sig[m] += d
    t = [0] * (N + 1)
    t[1] = 1
    for n in range(2, N + 1):
        acc = 0
        for k in range(1, n):
            acc += sig[k] * t[n - k]
        t[n] = (-24 * acc) // (n - 1)
    return t

def auto_test():
    t = tau_table(12)
    ok = (t[2] == -24 and t[3] == 252 and t[5] == 4830 and t[7] == -16744 and t[11] == 534612)
    ok = ok and v3(54) == 3 and v3f(Fraction(9, 2)) == 2 and (3 ** 4) == 81 and int(7) * int(7) == 49
    ok = ok and pente_min([Fraction(27), Fraction(0), Fraction(1)]) == Fraction(3, 2)
    return ok

# =============================================================================
# Part A: dual-number representation (Lemma 16)
# =============================================================================
PA = 3 ** 14

def m2(X, Y):
    return ((X[0] * Y[0] + X[1] * Y[2]) % PA, (X[0] * Y[1] + X[1] * Y[3]) % PA,
            (X[2] * Y[0] + X[3] * Y[2]) % PA, (X[2] * Y[1] + X[3] * Y[3]) % PA)

def a2(X, Y):
    return tuple((p + q) % PA for p, q in zip(X, Y))

def mul_d(U, V):
    return (m2(U[0], V[0]), a2(m2(U[0], V[1]), m2(U[1], V[0])))

def conj_C(U):
    A, B = U
    return ((A[0], (-A[1]) % PA, (-A[2]) % PA, A[3]), (B[0], (-B[1]) % PA, (-B[2]) % PA, B[3]))

def puiss_d(U, e):
    R = ((1, 0, 0, 1), (0, 0, 0, 0))
    while e:
        if e & 1:
            R = mul_d(R, U)
        U = mul_d(U, U)
        e >>= 1
    return R

def inv(a):
    return pow(int(a) % PA, -1, PA)

def donnees(U):
    """(psi1, psi2, theta, beta, gamma, z, e, trZ) for rho(h) = rho0(h)(1 + eps Z(h))"""
    A, B = U
    p1, p2 = A[0], A[3]
    det = (A[0] * A[3] - A[1] * A[2]) % PA
    di = inv(det)
    Ai = ((A[3] * di) % PA, (-A[1] * di) % PA, (-A[2] * di) % PA, (A[0] * di) % PA)
    Z = m2(Ai, B)
    th = p1 * inv(p2) % PA
    beta = A[1] * inv(p2) % PA
    gam = Z[2] * inv(th) % PA
    z = Z[3]
    e = (B[0] + B[3]) % PA
    return p1, p2, th, beta, gam, z, e, (Z[0] + Z[3]) % PA

def partie_A():
    dire('')
    dire('Part A: identities of Lemma 16 (dual-number representation, modulo 3**14)')
    rnd = random.Random(20261022)
    tout = True
    for (kappa, s) in [(11, 45), (13, 18), (17, 0), (23, 9), (15, 12), (21, 3), (27, 21)]:
        t1 = time.time()
        w = v3(2 * s - kappa)
        p1 = pow(7, s, PA)
        p2 = pow(7, kappa - s, PA)
        th = p1 * inv(p2) % PA
        A = (p1, p2, 0, p2)
        rx = (A, (p1, 0, p1, 0))          # rho0(x) (1 + eps [[0,0],[theta(x),0]])
        ry = conj_C(rx)
        def mot():
            U = ((1, 0, 0, 1), (0, 0, 0, 0))
            nx = ny = 0
            for _ in range(rnd.randint(1, 4)):
                if rnd.random() < 0.5:
                    e = rnd.randrange(1, 3 ** 5)
                    U = mul_d(U, puiss_d(rx, e))
                    nx += e
                else:
                    e = rnd.randrange(1, 3 ** 5)
                    U = mul_d(U, puiss_d(ry, e))
                    ny += e
            return U, nx, ny
        def_b = def_c = def_d = def_e = def_tr = 0
        tem_d = tem_e = False
        assoc = {}
        for _ in range(N_MOTS):
            U, nx, ny = mot()
            q1, q2, t, be, ga, z, e, trZ = donnees(U)
            if trZ % PA or (e - q2 * ((1 - t) * z + t * be * ga)) % PA:
                def_tr += 1
            F = ((1 + t) * z - t * be * ga) % PA
            W = mul_d(conj_C(U), U)
            eW = (W[1][0] + W[1][3]) % PA
            if (eW - 2 * q2 * q2 * (1 - t) * F) % PA:
                def_c += 1
            m = nx - ny
            if (be - m) % 3 or (ga - m) % 3:
                def_b += 1
            if (nx + ny) % 3 == 0:
                if F % 3:
                    def_d += 1
            elif F % 3:
                tem_d = True
            if w >= 1:
                if (be - m) % 9 or (ga - m) % 9:
                    def_e += 1
            elif (be - m) % 9 or (ga - m) % 9:
                tem_e = True
            if m % 3 == 0:
                assoc.setdefault((m // 3) % 3, set()).add((be // 3) % 3)
        # on ker v, is beta/3 a function of mu/3?
        tem_f = any(len(x) > 1 for x in assoc.values())
        # system of Lemma 17 (lambda(x) normalized to 1, a unit)
        _, _, _, _, _, _, ex, _ = donnees(rx)
        xy = mul_d(rx, ry)
        q1, q2, _, _, _, _, exy, _ = donnees(xy)
        er_x = (p1 - p2) % PA
        er_xy = 2 * (q1 - q2) % PA
        Dt = (er_x * exy - ex * er_xy) % PA
        sys_ok = (v3(ex) == 0 and v3(er_x) == 1 + w and v3(er_xy) == 1 + w and v3(exy) >= 1 + w and v3(Dt) == 1 + w)
        ok = (def_tr == 0 and def_c == 0 and def_b == 0 and def_d == 0 and def_e == 0 and sys_ok and tem_d
              and (w >= 1 or (tem_e and tem_f)) and (w == 0 or not tem_f))
        tout = tout and ok
        dire('  kappa=%2d s=%2d w=%d : trace/det %d, (c) %d, beta=gamma=v (mod 3) %d, F (u=0) %d, (e) mod 9 %d failures ;'
             ' system %s ; witness F unit (u!=0) %s ; beta/3 function of mu/3 on ker v : %s  (%.1f s)'
             % (kappa, s, w, def_tr, def_c, def_b, def_d, def_e, 'ok' if sys_ok else 'FAILED',
                'yes' if tem_d else 'no', 'no' if tem_f else 'yes', time.time() - t1))
    return tout

# =============================================================================
# Part B: Hatada lattice
# =============================================================================
CACHE_MANIN = {}

def matrice_Tn_cache(w, n):
    if n not in CACHE_MANIN:
        CACHE_MANIN[n] = manin(n)
    Mt = [[0] * (w + 1) for _ in range(w + 1)]
    for (a, b, c, d) in CACHE_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 matrice_reseau(base, w, n):
    Mt = matrice_Tn_cache(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 carre(A):
    d = len(A)
    return [[sum(A[i][r] * A[r][j] for r in range(d)) for j in range(d)] for i in range(d)]

def partie_B():
    dire('')
    dire('Part B: law of inert primes and Corollary Y\' on L_k, 12 <= k <= %d, l in %s' % (KMAX_B, PREMIERS_B))
    echecs = 0
    exact = {l: False for l in PREMIERS_B}
    for k in range(12, KMAX_B + 1, 2):
        if k == 14:
            continue
        t1 = time.time()
        w = k - 2
        base = reseau_hatada(w)
        ok = True
        for l in PREMIERS_B:
            j = v3(l + 1)
            s = 1 + l ** (k - 1)
            A = matrice_reseau(base, w, l)
            X = mcomb((1, carre(A)), (-s * s, None))
            pm = pente_min(polcar(X))
            borne = min(2 * j + 2, 10)
            if pm < borne:
                ok = False
            if pm == borne:
                exact[l] = True
            if k % 6 != 4 and j <= 4:
                if not pentes_egales(polcar(A), j):
                    ok = False
        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('  failing weights: %d (expected: 0) ; bound min(2j+2, 10) attained for each l: %s' % (echecs, 'yes' if all(exact.values()) else 'no'))
    return echecs == 0 and all(exact.values())

# =============================================================================
# Part C: exact values in dimension 1
# =============================================================================
def rep(l, n):
    """l = x**2 + x y + n y**2 with x, y integers?"""
    y = 0
    while (4 * n - 1) * y * y <= 4 * l:
        D = y * y - 4 * (n * y * y - l)
        if D >= 0:
            r = int(D ** 0.5)
            while r * r > D:
                r -= 1
            while (r + 1) * (r + 1) <= D:
                r += 1
            if r * r == D and (r - y) % 2 == 0:
                return True
        y += 1
    return False

def premiers(N):
    crible = [True] * (N + 1)
    crible[0] = crible[1] = False
    for p in range(2, int(N ** 0.5) + 1):
        if crible[p]:
            for m in range(p * p, N + 1, p):
                crible[m] = False
    return [p for p in range(N + 1) if crible[p]]

def partie_C():
    dire('')
    dire('Part C: exact a_l for k = 12, 16, 18, 20, 22, 26 and all primes l < %d' % LMAX_C)
    t1 = time.time()
    N = LMAX_C
    t = tau_table(N)
    dire('  table of tau up to %d: %.1f s' % (N, time.time() - t1))
    def sigma_table(e):
        sg = [0] * (N + 1)
        for d in range(1, N + 1):
            de = d ** e
            for m in range(d, N + 1, d):
                sg[m] += de
        return sg
    eis = {0: None}
    for (m, cm) in [(4, 240), (6, -504), (8, 480), (10, -264), (14, -24)]:
        sg = sigma_table(m - 1)
        eis[m] = [1] + [cm * sg[n] for n in range(1, N + 1)]
    P = premiers(N - 1)
    ok = True
    compte = {'inertes': 0, 'decomposes': 0, 'k4_547': 0, 'corY': 0}
    temoins = {}
    for k in [12, 16, 18, 20, 22, 26]:
        e = eis[k - 12] if k > 12 else None
        for l in P:
            if l == 3:
                continue
            if e is None:
                a = t[l]
            else:
                a = sum(t[i] * e[l - i] for i in range(1, l + 1))
            s = 1 + l ** (k - 1)
            if l % 3 == 2:
                j = v3(l + 1)
                d = v3(a * a - s * s)
                compte['inertes'] += 1
                if d < min(2 * j + 2, 10):
                    ok = False
                    dire('  FAILED R\' : k=%d l=%d' % (k, l))
                if k % 6 != 4 and j <= 4:
                    compte['corY'] += 1
                    if v3(a) != j:
                        ok = False
                        dire('  FAILED Y\' : k=%d l=%d' % (k, l))
                temoins[(k, l)] = d
            elif l % 27 == 1 and pow(3, (l - 1) // 3, l) == 1:
                d = v3(a - s)
                compte['decomposes'] += 1
                if d < 7:
                    ok = False
                    dire('  FAILED R\'\' : k=%d l=%d' % (k, l))
                if k % 6 == 4 and rep(l, 547):
                    compte['k4_547'] += 1
                    if d < (9 if l % 81 == 1 else 8):
                        ok = False
                        dire('  FAILED R\'\' (547) : k=%d l=%d' % (k, l))
                temoins[(k, l)] = d
    t971 = (temoins.get((18, 971)), temoins.get((16, 971)))
    t271 = temoins.get((16, 271))
    tem = (t971 == (11, 12) and t271 == 7)
    dire('  cases checked: %d inert, %d split (of which %d with k == 4 (mod 6) and l = x**2+xy+547y**2), %d for Y\''
         % (compte['inertes'], compte['decomposes'], compte['k4_547'], compte['corY']))
    dire('  witnesses: v_3(a_971**2 - s_971**2) = %s for k = 18 (2j+2 = 12 fails) and %s for k = 16 ;'
         ' v_3(a_271 - s_271) = %s for k = 16 (271 not represented by the form 547)' % (t971[0], t971[1], t271))
    dire('  Part C summary: %s  (%.1f s)' % ('ok' if ok and tem else 'FAILED', time.time() - t1))
    return ok and tem

# =============================================================================
# Part D: Pic(Z + f Z[omega]) for f = 9 and 27
# =============================================================================
def partie_D():
    dire('')
    dire('Part D: (O_K / f)* / (Z/f)* O_K*, f = 9 and 27, and forms x**2 + xy + 61 y**2, x**2 + xy + 547 y**2')
    t1 = time.time()
    ok = True
    P = premiers(LMAX_D)
    for (f, n, h_att) in [(9, 61, 3), (27, 547, 9)]:
        def mu(p, q):
            # (a + b w)(c + d w), w**2 = -1 - w
            return ((p[0] * q[0] - p[1] * q[1]) % f, (p[0] * q[1] + p[1] * q[0] - p[1] * q[1]) % f)
        U = [(a, b) for a in range(f) for b in range(f) if (a * a - a * b + b * b) % 3]
        gens = [(m, 0) for m in range(1, f) if m % 3] + [(f - 1, 0), (0, 1)]
        H = {(1, 0)}
        pile = [(1, 0)]
        while pile:
            x = pile.pop()
            for g in gens:
                y = mu(x, g)
                if y not in H:
                    H.add(y)
                    pile.append(y)
        h = len(U) // len(H)
        ordre_max = 0
        for z in U:
            p = z
            o = 1
            while p not in H:
                p = mu(p, z)
                o += 1
            ordre_max = max(ordre_max, o)
        cyc = (h == h_att and ordre_max == h_att)
        incoh = 0
        for l in P:
            if l % 3 != 1:
                continue
            pi = None
            for b in range(1, int((4 * l / 3) ** 0.5) + 2):
                D = 4 * l - 3 * b * b
                if D < 0:
                    break
                r = int(D ** 0.5)
                while r * r > D:
                    r -= 1
                while (r + 1) * (r + 1) <= D:
                    r += 1
                if r * r == D and (r + b) % 2 == 0:
                    pi = ((r + b) // 2, b)     # a**2 - a b + b**2 = l
                    break
            triv = (pi[0] % f, pi[1] % f) in H
            if triv != rep(l, n):
                incoh += 1
            if f == 9 and triv != (pow(3, (l - 1) // 3, l) == 1):
                incoh += 1
        ok = ok and cyc and incoh == 0
        dire('  f = %2d : order %d, element of order %d (cyclic: %s) ; inconsistencies with x**2+xy+%d y**2%s : %d'
             % (f, h, ordre_max, 'yes' if cyc else 'no', n, ' and with the cube criterion' if f == 9 else '', incoh))
    dire('  (%.1f s)' % (time.time() - t1))
    return ok

# =============================================================================
if __name__ == '__main__':
    dire('=' * 76)
    dire('VERIF_INERTES: law of inert primes and split primes, level 1')
    dire('Budget %d s ; Part B: k <= %d ; Part C: l < %d ; Part D: l < %d' % (int(BUDGET), KMAX_B, LMAX_C, LMAX_D))
    dire('Expected memory: a few megabytes.')
    dire('=' * 76)
    if not auto_test():
        dire('Self-test: FAILED. Stop (non-conforming environment, e.g. Sage preparser).')
        raise SystemExit(1)
    dire('Self-test: OK (tau(2), tau(3), tau(5), tau(7), tau(11), valuations).')
    rA = partie_A()
    rD = partie_D()
    rC = partie_C()
    rB = partie_B()
    dire('')
    dire('SUMMARY: Part A (Lemma 16) %s, Part B (lattice) %s, Part C (dimension 1) %s, Part D (Pic) %s'
         % tuple('OK' if r else 'FAILED' for r in (rA, rB, rC, rD)))
    dire('Total time: %.1f s' % (time.time() - T0))