# -*- coding: utf-8 -*-
# =============================================================================
# VERIF_A3.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 3
# (Lemma 3, Theorem F) and Section "Verification".
# -----------------------------------------------------------------------------
# Theorem F: for every eigenform f of level 1 and weight k >= 12,
#     a_3(f) == 0 (mod 9),  and  a_3(f) == 0 (mod 27) if 6 does not divide k,
# congruences in the ring of algebraic integers (a_3/9, resp. a_3/27, integral).
#
# RUN WITH PLAIN PYTHON 3. No library (no numpy, no scipy, no sympy).
# No input needed. Clean stop before the time budget is exhausted.
# Running time: well under a minute with the values below. Memory: a few MB.
#
# Part A: the valuation inequality in the proof of Lemma 3, for all even
#         w = k - 2, 10 <= w <= WMAX_A, and all odd i, 3 <= i <= w-1:
#           v3(c_i) >= v3(c_1) + 1   (and + 2 if w is not == 1 mod 3),
#           c_i = 3**i * C(w,i) * (1 + 2**(w-i)),  v3(c_1) = 2 + v3(w) + v3(w-1).
# Part B: Hatada's lattice L_k (odd period polynomials whose periods m_j are
#         integers), 12 <= k <= KMAX_B, 3-adic basis:
#   (i)   exact identity  sum_j C(w,j) 3**(w-j) (1 + 2**j) m_j = 0 ;
#   (ii)  m_1 == 0 mod 3, and mod 9 if 6 does not divide k ;
#   (iii) T_3(L_k) contained in 9 L_k, and in 27 L_k if 6 does not divide k ;
#   (iv)  when dim S_k = 1, the eigenvalue of T_3 is the known a_3.
# T_3 acts through the seven Manin matrices (ad - bc = 3, a > c >= 0, d > b >= 0).
# =============================================================================
import sys
import time
from fractions import Fraction
from math import comb, gcd

WMAX_A = int(3000)       # part A: even w up to WMAX_A
KMAX_B = int(100)        # part B: weights up to KMAX_B
BUDGET = float(480.0)    # seconds: clean stop beyond
PARTIE_B = True          # switch: False to run part A only

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

# ---------------------------------------------------------- substitutions ---
# Basis of Z[x,y]_w : x**j y**(w-j), j = 0..w.  The matrix (a b ; c d) acts by
# F(x,y) -> F(a x + c y, b x + d y).  Column j = coefficients of the image of x**j y**(w-j).
def image_monome(w, j, a, b, c, d):
    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
    return out

MAN3 = [(1, 0, 0, 3), (1, 1, 0, 3), (1, 2, 0, 3), (2, 1, 1, 2),
        (3, 0, 0, 1), (3, 0, 1, 1), (3, 0, 2, 1)]

def lignes_relations(w):
    """A functional m (m_i = value on x**i y**(w-i)) is an odd period
    iff for every monomial F : m(F) + m(F.S) = 0, m(F) + m(F.U) + m(F.U2) = 0,
    m(F) + m(F.eps) = 0, with S = (0,-1;1,0), U = (1,-1;1,0), U2 = (0,-1;1,-1),
    eps = (-1,0;0,1)."""
    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):
    """Basis of the kernel over Q (primitive integer vectors)."""
    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):
    """Returns coefficients c (mod 3), not all zero, with sum c_i b_i == 0 mod 3, or None."""
    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]))]
        assert all(x % 3 == 0 for x in v)
        base[i0] = [x // 3 for x in v]

def reseau_hatada(k):
    w = k - 2
    base = noyau_rationnel(lignes_relations(w), w + 1)
    return sature_en_3(base)

def applique_T3(w, m):
    out = [0] * (w + 1)
    for (a, b, c, d) in MAN3:
        for l in range(w + 1):
            im = image_monome(w, l, a, b, c, d)
            out[l] += sum(im[i] * m[i] for i in range(w + 1) if im[i] and m[i])
    return out

def satisfait_relations(w, m, rows):
    return all(sum(r[i] * m[i] for i in range(w + 1)) == 0 for r in rows)

# ---------------------------------------------------------- self-test ---
A3_CONNUS = {12: 252, 16: -3348, 18: -4284, 20: 50652, 22: -128844, 26: -195804}

def auto_test():
    if v3(252) != 2 or v3(-3348) != 3 or comb(10, 3) != 120:
        return False
    if (1 + 2 ** 9) != 513 or v3(513) != 3:
        return False
    base = reseau_hatada(12)
    if len(base) != 1:
        return False
    m = base[0]
    impairs = [m[j] for j in range(1, 10, 2)]
    if impairs != [48, -25, 20, -25, 48] and impairs != [-48, 25, -20, 25, -48]:
        return False
    t = applique_T3(10, m)
    return all(t[j] == 252 * m[j] for j in range(11))

dire('=' * 76)
dire('VERIF_A3 : a_3 == 0 mod 9 (mod 27 if 6 does not divide k) in level 1')
dire('Budget %.0f s ; part A: w <= %d ; part B: k <= %d' % (BUDGET, WMAX_A, KMAX_B))
dire('=' * 76)
if not auto_test():
    dire('SELF-TEST FAILED: stop (non-conforming environment).')
    sys.exit(1)
dire('Self-test: OK (period of Delta (48,-25,20,-25,48), tau(3) = 252).')

# ------------------------------------------------------------- part A ---
dire('')
dire('Part A: valuation inequality, even w from 10 to %d' % WMAX_A)
MOD = 3 ** 40
v3_1p2 = {}
for n in range(1, WMAX_A + 1, 2):
    v3_1p2[n] = v3((1 + pow(2, n, MOD)) % MOD)

def somme_chiffres3(n):
    s = 0
    while n:
        s += n % 3
        n //= 3
    return s

S3 = [somme_chiffres3(n) for n in range(WMAX_A + 2)]

def v3_binom(n, i):          # Kummer: number of carries
    return (S3[i] + S3[n - i] - S3[n]) // 2

# cross-check of the computation of v3(C(n,i)) by Kummer on small cases
for n in range(1, 121):
    for i in range(n + 1):
        if v3_binom(n, i) != v3(comb(n, i)):
            dire('Kummer: inconsistency, stop.')
            sys.exit(1)

ecarts_A = 0
min_delta = {0: 10 ** 9, 1: 10 ** 9, 2: 10 ** 9}
for w in range(10, WMAX_A + 1, 2):
    vc1 = 1 + v3(w) + v3_1p2[w - 1]
    if vc1 != 2 + v3(w) + v3(w - 1):
        ecarts_A += 1
    seuil = 1 if w % 3 == 1 else 2
    for i in range(3, w, 2):
        vci = i + v3_binom(w, i) + v3_1p2[w - i]
        dl = vci - vc1
        if dl < min_delta[w % 3]:
            min_delta[w % 3] = dl
        if dl < seuil:
            ecarts_A += 1
    if time.time() - T0 > BUDGET / 3:
        dire('  budget of part A reached at w = %d' % w)
        break
dire('  discrepancies: %d ; min of v3(c_i) - v3(c_1) according to w mod 3 : %s' % (ecarts_A, min_delta))
dire('  (expected: 0 discrepancies; minimum 1 if w == 1 mod 3, at least 2 otherwise)')

# ------------------------------------------------------------- part B ---
if PARTIE_B:
    dire('')
    dire('Part B: Hatada lattice, weights 12 to %d' % KMAX_B)
    ecarts_B = 0
    for k in range(12, KMAX_B + 1, 2):
        if k == 14:
            continue
        t1 = time.time()
        w = k - 2
        rows = lignes_relations(w)
        base = reseau_hatada(k)
        d = len(base)
        e = 3 if k % 6 else 2           # 9 or 27
        ok_id = ok_m1 = ok_T3 = ok_stab = True
        vals = []
        for m in base:
            s = sum(comb(w, j) * 3 ** (w - j) * (1 + 2 ** j) * m[j] for j in range(1, w, 2))
            ok_id = ok_id and (s == 0)
            ok_m1 = ok_m1 and (m[1] % (3 ** (e - 1)) == 0)
            t = applique_T3(w, m)
            ok_T3 = ok_T3 and all(x % (3 ** e) == 0 for x in t)
            ok_stab = ok_stab and satisfait_relations(w, t, rows)
            if d == 1:
                l = next(j for j in range(w + 1) if m[j])
                vals.append(Fraction(t[l], m[l]))
        ok_val = True
        if d == 1 and k in A3_CONNUS:
            ok_val = (vals[0] == A3_CONNUS[k])
        bilan = ok_id and ok_m1 and ok_T3 and ok_stab and ok_val
        if not bilan:
            ecarts_B += 1
        dire('  k=%3d d=%2d  identity:%s  m_1 mod %d:%s  T3 in %d L:%s  stability:%s%s  (%.1f s)' % (
            k, d, 'ok' if ok_id else 'FAILED', 3 ** (e - 1), 'ok' if ok_m1 else 'FAILED',
            3 ** e, 'ok' if ok_T3 else 'FAILED', 'ok' if ok_stab else 'FAILED',
            ('  a_3=%s' % vals[0]) if d == 1 else '', time.time() - t1))
        if time.time() - T0 > BUDGET:
            dire('  budget reached after k = %d' % k)
            break
    dire('  failing weights: %d (expected: 0)' % ecarts_B)

dire('')
dire('Total time: %.1f s' % (time.time() - T0))