# -*- coding: utf-8 -*-
# =============================================================================
# VERIF_HATADA.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 2
# (Theorem A, Proposition B, Theorem C) and Section "Verification".
# -----------------------------------------------------------------------------
# Self-contained check of Hatada's question (Manuscripta Math. 34, 1981, p. 306)
#   "lambda_p == 1 + p (mod 2**a) if p == +-1 (mod 2**(a-1)) and 6 <= a <= 13 ?"
# for the six rational eigenforms of level 1 (weights 12, 16, 18, 20, 22, 26).
#
# RUN WITH PLAIN PYTHON 3 (numpy only; no scipy, no sympy). No input needed.
# Running time: a few seconds. Memory: a few MB.
# Computations modulo 2**23: all 2-adic valuations printed (< 23) are exact.
# Convolution sums stay below 2**60 (no int64 overflow).
# =============================================================================
import sys
import numpy as np

N = int(13313)
NB = int(23)
MOD = int(1) << NB

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

dire("Expected memory: about %d KB per series, about twenty series." % (int(8) * (N + 1) // 1024))

def eta(n):
    f = np.zeros(n + 1, dtype=np.int64)
    f[0] = 1
    k = 1
    while True:
        g1 = k * (3 * k - 1) // 2
        g2 = k * (3 * k + 1) // 2
        if g1 > n and g2 > n:
            break
        s = -1 if k % 2 else 1
        if g1 <= n:
            f[g1] += s
        if g2 <= n:
            f[g2] += s
        k += 1
    return f % MOD

def mul(a, b):
    return np.convolve(a, b)[:N + 1] % MOD

def sigma(e):
    s = np.zeros(N + 1, dtype=np.int64)
    for d in range(1, N + 1):
        s[d::d] += pow(d, int(e), MOD)
    return s % MOD

def v2(x):
    x = int(x) % MOD
    if x == 0:
        return NB
    v = 0
    while x % 2 == 0:
        x //= 2
        v += 1
    return v

def premier(n):
    if n < 2:
        return False
    d = 2
    while d * d <= n:
        if n % d == 0:
            return False
        d += 1
    return True

# --- series ---
dire("Computing Delta = q * prod(1-q**n)**24 ...")
e1 = eta(N); e2 = mul(e1, e1); e4 = mul(e2, e2); e8 = mul(e4, e4); e16 = mul(e8, e8)
F = mul(e16, e8)                              # F[m] = tau(m+1) mod 2**23
EIS = {4: (240, 3), 6: (-504, 5), 8: (480, 7), 10: (-264, 9), 14: (-24, 13)}
A = {12: F}
SIG = {12: sigma(11)}
for w, (c, e) in EIS.items():
    s = sigma(e)
    Ew = (int(c) * s) % MOD
    Ew[0] = 1
    A[12 + w] = mul(Ew, F)                    # Delta_{12+w} = E_w * Delta
    SIG[12 + w] = sigma(11 + w)
def coef(k, n):
    return int(A[k][n - 1])                   # a_n(Delta_k) mod 2**23

# --- self-test: refuses to continue if a control value is wrong ---
CONTROLE = {12: (-24, 252, -1472, 4830), 16: (216, -3348), 18: (-528, -4284),
            20: (456, 50652), 22: (-288, -128844), 26: (-48, -195804)}
ok = True
for k, vals in CONTROLE.items():
    for i, v in enumerate(vals):
        if coef(k, i + 2) != int(v) % MOD:
            ok = False
TAU_12289 = 13297575938818018265090          # exact value (PARI/GP: ramanujantau(12289))
if coef(12, 12289) != TAU_12289 % MOD:
    ok = False
if not ok:
    dire("SELF-TEST FAILED: computation stopped (results not reliable).")
    sys.exit(1)
dire("Self-test OK (known values + exact tau(12289) from PARI, residue mod 2**23).")
dire("")

# --- 1. class congruences: a_n == c_r * sigma_{k-1}(n) mod 2**A_r, n odd == r mod 8 ---
TABLE = {12: {1: (11, 1), 3: (13, 1217), 5: (12, 1537), 7: (14, 705)},
         16: {1: (10, 1), 3: (12, 609), 5: (11, 769), 7: (13, 353)},
         18: {1: (10, 1), 3: (12, 97), 5: (11, 769), 7: (13, 865)},
         20: {1: (12, 1), 3: (14, 3969), 5: (13, 3073), 7: (15, 2945)},
         22: {1: (12, 1), 3: (14, 897), 5: (13, 3073), 7: (15, 3969)},
         26: {1: (10, 1), 3: (12, 289), 5: (11, 257), 7: (13, 545)}}
NSTURM = int(2000)                            # Sturm bound for Gamma_0(64): 8k <= 208
dire("1. Class congruences, all odd n <= %d (Sturm bound 8k <= 208):" % NSTURM)
for k in sorted(TABLE):
    lig = []
    for r in (1, 3, 5, 7):
        Ar, cr = TABLE[k][r]
        m = int(1) << Ar
        bad = [n for n in range(r, NSTURM + 1, 8)
               if (int(cr) * int(SIG[k][n]) - coef(k, n)) % m != 0]
        lig.append("r=%d: 2**%d %s" % (r, Ar, "OK" if not bad else "FAILED"))
    dire("   Delta_%d : " % k + "  ".join(lig))
dire("")

# --- 2. the witnesses of the negative answer ---
dire("2. Witnesses (the question requires v2(a_p - 1 - p) >= a whenever p == +-1 mod 2**(a-1)):")
for p in (12289, 13313):
    m = max(v2(p - 1), v2(p + 1))
    dire("   p = %d  prime: %s   p-1 = 2**%d * %d" % (p, premier(p), v2(p - 1), (p - 1) >> v2(p - 1)))
    for k in sorted(A):
        v = v2(coef(k, p) - 1 - p)
        niv = [a for a in range(6, 14) if a - 1 <= m and v < a]
        dire("      Delta_%d : v2 = %2d   failing levels a: %s" % (k, v, niv if niv else "none"))
dire("")

# --- 3. scan: no violation for a <= 10 among all p <= N ---
dire("3. Scan of all odd primes p <= %d:" % N)
for k in sorted(A):
    pires = {}
    for p in range(3, N + 1, 2):
        if not premier(p):
            continue
        m = max(v2(p - 1), v2(p + 1))
        v = v2(coef(k, p) - 1 - p)
        for a in range(6, 14):
            if a - 1 <= m and v < a:
                pires.setdefault(a, p)
    dire("   Delta_%d : smallest witness for each level a: %s" % (k, pires if pires else "none"))
dire("")
dire("END.")