# -*- coding: utf-8 -*-
"""
SELF-CONTAINED VERIFIER -- 2-adic valuation of the constant of the
Prouhet-Tarry-Escott problem, SYMMETRIC ideal solutions, ODD sizes.
New bound at n = 15.

David Niedbala Giraudin, ORCID 0009-0009-1526-1178.
Technical assistance: Claude (Anthropic).
Companion article [V7]: D. Niedbala Giraudin, The 2-adic valuation of the constant
for symmetric ideal Prouhet-Tarry-Escott solutions, hexagon:2611.00057.

Replays every numerical assertion from scratch and prints PASS or FAIL.
Pure Python 3: no dependency (no numpy, no scipy, no sympy).
No command-line argument. No intermediate file. Prints as it goes.
Exponentiation written with a double star throughout.

WHAT IS CHECKED
  B  the two-product identity      v_2(C) = v_2( prod(1+a_i) - prod(1-a_i) )
  C  the index identity            sum_i v_2(a_j + a_i) = v_2(C), for EVERY j
  D  a block of equal valuation is BIPARTITE modulo 4   (exhaustive)
  E  the block criterion                                 (exhaustive, m = 2..5)
  F  the rho forced by the index identity, including rho = 4 at n = 11 [5,1,1]
  G  NEGATIVE calibration : at n = 11 the two slots at t = 7 collapse
  H  POSITIVE calibration : the profiles of the REAL witnesses do not collapse
  I  THE RESULT : at n = 15 all four slots at t = 11 collapse
  J  t = 12 is empty, t = 13 holds a single slot, hence v_2(C) >= 14 at n = 15

SETTINGS -- change nothing but these
"""
BUDGET_SECONDS = 480      # total budget; heavy blocks are SKIPPED past it
FRONT_CAP = 60000        # ceiling on the lifting front (memory guard)
PREC = 150                # 2-adic working precision for the witnesses

import time
from collections import Counter

T0 = time.time()
PASS = [0]
FAIL = [0]
SKIP = [0]
INF = 10 ** 9


def say(s=""):
    print(s, flush=True)


def dit(nom, ok, detail=""):
    if ok:
        PASS[0] += 1
        say("   PASS  %-50s %s" % (nom, detail))
    else:
        FAIL[0] += 1
        say("   FAIL  %-50s %s" % (nom, detail))


def saute(nom, pourquoi):
    SKIP[0] += 1
    say("   SKIP  %-50s %s" % (nom, pourquoi))


def reste():
    return BUDGET_SECONDS - (time.time() - T0)


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


def binom(n, k):
    if k < 0 or k > n:
        return 0
    r = 1
    for i in range(k):
        r = r * (n - i) // (i + 1)
    return r


def polmul(A, B):
    C = [0] * (len(A) + len(B) - 1)
    for i, a in enumerate(A):
        if a:
            for j, b in enumerate(B):
                C[i + j] += a * b
    return C


MM = 1 << PREC


def sym(x):
    x %= MM
    return x if x < MM // 2 else x - MM


def sqrt2adic(d, P):
    """2-adic square root of d = 1 mod 8, modulo 2 to the power P."""
    s = 1
    for m in range(3, P):
        if (s * s - d) % (1 << (m + 1)):
            s += 1 << (m - 1)
    return s % (1 << P)


def hensel(coef, r0, P):
    """Newton refinement of a root of coef (decreasing degree) modulo 2**P."""
    M = 1 << (P + 80)

    def e(c, x):
        y = 0
        for a in c:
            y = (y * x + a) % M
        return y

    n = len(coef) - 1
    der = [coef[i] * (n - i) for i in range(n)]
    r = r0 % M
    for _ in range(120):
        f = e(coef, r)
        if f == 0:
            break
        fp = e(der, r)
        s = v2(fp)
        u = fp >> s
        rn = (r - ((f >> s) * pow(u, -1, M)) % M) % M
        if rn == r:
            break
        r = rn
    return r % (1 << P)


say("=" * 78)
say("VERIFIER -- PTE, 2-adic valuation, symmetric solutions, odd sizes")
say("=" * 78)
say("budget %d s; front ceiling %d; precision 2**%d"
    % (BUDGET_SECONDS, FRONT_CAP, PREC))
say()

# ======================================================= A. THE WITNESSES
say("-" * 78)
say("A. THE WITNESSES, LIFTED")
say("-" * 78)
QUINT = [1, 0, -4278125, 0, 3664306448884, -975527373189600]
TEMOINS = []
TEMOINS.append(("n=5  (9,5,1,-7,-8)", [9, 5, 1, -7, -8]))
TEMOINS.append(("n=7", [-50, -38, -13, -7, 24, 33, 51]))
TEMOINS.append(("Letac I  n=9", [-98, -82, -58, -34, 13, 16, 69, 75, 99]))
for nom, d, rac, fab in (
    ("CMSV A1  n=11", 26601, lambda s: (s - 163) // 2,
     lambda a: [-95, -68 - a, -52, -48, -13, -9, 30, 34, 61, 65, 95 + a]),
    ("CMSV A2  n=11", 400569, lambda s: (633 - s) // 2,
     lambda a: [-589, -546 + a, -363, -185, -170, 41, 87 - a, 234, 355, 548, 588])):
    al = rac(sqrt2adic(d, PREC + 10)) % MM
    TEMOINS.append((nom, [x % MM for x in fab(al)]))
TEMOINS.append(("witness n=13",
                [x % MM for x in [-1827, -989, -715, -405, 64, 902, 1330, 1640]]
                + [hensel(QUINT, r, PREC) for r in
                   (2723974027, 1263895509, 3276945658, 3386758734, 2233327960)]))

DONNEES = []
for nom, A in TEMOINS:
    A = [a % MM for a in A]
    n = len(A)
    M = (n - 3) // 2
    P = [sum(pow(a, 2 * j + 1, MM) for a in A) % MM for j in range(M + 1)]
    nuls = all(min(p, MM - p) < (1 << 40) for p in P)
    sg = [sym(a) for a in A]
    t = sum(v2(a) for a in sg if a % 2 == 0)
    w = sum(1 for a in sg if a % 2)
    prof = tuple(sorted((v2(a) for a in sg if a % 2 == 0), reverse=True))
    DONNEES.append((nom, A, n, M, t, w, prof))
    dit("%s : odd power sums vanish, v_2(C) = %d" % (nom, 1 + t), nuls,
        "n = %2d, w = %2d, profile %s" % (n, w, list(prof)))

# ============================================ B. THE TWO-PRODUCT IDENTITY
say()
say("-" * 78)
say("B. THE TWO-PRODUCT IDENTITY")
say("-" * 78)
say("   F(x) = x H(x**2) - c gives F(1) = prod(1-a_i) = H(1) - c and")
say("   prod(1+a_i) = H(1) + c, hence 2c = prod(1+a_i) - prod(1-a_i) and")
say("   v_2(C) = v_2(2c) = v_2( prod(1+a_i) - prod(1-a_i) ).")
for (nom, A, n, M, t, w, prof) in DONNEES:
    P1, P2 = 1, 1
    for a in A:
        P1 = (P1 * (1 + a)) % MM
        P2 = (P2 * (1 - a)) % MM
    dit("%s : v_2(A - B) = 1 + t = %d" % (nom, 1 + t),
        v2((P1 - P2) % MM) == 1 + t, "v_2(A) = %s, v_2(B) = %s"
        % ("inf" if v2(P1) >= INF else v2(P1), "inf" if v2(P2) >= INF else v2(P2)))

# ================================================== C. THE INDEX IDENTITY
say()
say("-" * 78)
say("C. THE INDEX IDENTITY, FOR EVERY INDEX")
say("-" * 78)
say("   Putting x_0 = a_j in prod(x_0+a_i) - prod(x_0-a_i) = 2c kills the second")
say("   product, so prod_i (a_j + a_i) = 2c and sum_i v_2(a_j + a_i) = v_2(C).")
for (nom, A, n, M, t, w, prof) in DONNEES:
    lignes = []
    for j in range(n):
        s = 0
        for i in range(n):
            d = (A[j] + A[i]) % MM
            s += v2(d) if min(d, MM - d) >= (1 << 40) or d == 0 else v2(d)
        lignes.append(s)
    dit("%s : all %d rows equal %d" % (nom, n, 1 + t),
        all(x == 1 + t for x in lignes),
        "rows = %s" % sorted(set(lignes)))

# =================================== D. A BLOCK IS BIPARTITE MODULO FOUR
say()
say("-" * 78)
say("D. A BLOCK OF EQUAL VALUATION IS BIPARTITE MODULO FOUR  (exhaustive)")
say("-" * 78)
say("   For b, b' odd, (b+b') + (b-b') = 2b has valuation 1, so")
say("   min( v_2(b-b'), v_2(b+b') ) = 1, and more precisely")
say("   v_2(b+b') = 1 if and only if b = b' modulo 4.")
mauvais = 0
for b in range(1, 256, 2):
    for bp in range(1, 256, 2):
        if b == bp:
            continue
        s = (b + bp) % 256
        un = (v2(s) == 1) if s else False
        if un != (b % 4 == bp % 4):
            mauvais += 1
dit("over all odd pairs modulo 2**8", mauvais == 0, "discrepancies: %d" % mauvais)

# ========================================================= E. THE CRITERION
say()
say("-" * 78)
say("E. THE BLOCK CRITERION  (exhaustive, m = 2 to 5)")
say("-" * 78)
say("   Counting the cross rho from both sides gives p(R-p+1) = q(R-q+1),")
say("   hence p = q or R = m-1. Therefore")
say("      R realisable  <=>  R = m-1,  or  m even and R >= 3m/2 - 1.")
say("   In particular R = m is FORBIDDEN as soon as m >= 3.")


def bloc_realisable(m, R):
    if m == 1:
        return R == 0
    if R < m - 1:
        return False
    if R == m - 1:
        return True
    return m % 2 == 0 and R >= 3 * m // 2 - 1


def brute(m, K, Rmax):
    M_ = 1 << K
    out = set()

    def rec(bs):
        if len(bs) == m:
            rows = []
            for i in range(m):
                s = 0
                for j in range(m):
                    if i != j:
                        d = (bs[i] + bs[j]) % M_
                        s += v2(d) if d else K
                rows.append(s)
            if len(set(rows)) == 1 and rows[0] <= Rmax:
                out.add(rows[0])
            return
        for b in range(1, M_, 2):
            rec(bs + [b])
    rec([1])
    return out


for m, K, Rmax in ((2, 8, 9), (3, 7, 9), (4, 6, 11), (5, 5, 9)):
    at = sorted(brute(m, K, Rmax))
    lim = min(Rmax, K - 1)
    coh = all((R in at) == bloc_realisable(m, R) for R in range(0, lim + 1))
    dit("m = %d : brute force matches the criterion up to R = %d" % (m, lim), coh,
        "reached: %s" % at)

# ================================================ F. THE rho FORCED BY C
say()
say("-" * 78)
say("F. THE rho FORCED BY THE INDEX IDENTITY")
say("-" * 78)


def rhos(t, profil):
    cnt = Counter(profil)
    out = {}
    for f, mf in sorted(cnt.items()):
        autre = sum(mv * min(f, v) for v, mv in cnt.items() if v != f)
        out[f] = (mf, (1 + t) - (1 + f) - autre - (mf - 1) * f)
    return out


for n, t, pr in ((11, 7, (5, 1, 1)), (11, 7, (3, 2, 2)), (11, 10, (4, 2, 2, 1, 1)),
                 (13, 11, (7, 1, 1, 1, 1)), (13, 11, (5, 2, 2, 1, 1)),
                 (15, 11, (9, 1, 1)), (15, 11, (7, 2, 2)), (15, 11, (5, 3, 3))):
    r = rhos(t, pr)
    say("   n = %2d, t = %2d, %-20s %s"
        % (n, t, str(list(pr)),
           " ".join("f=%d: m=%d, R=%d" % (f, m, R) for f, (m, R) in r.items())))
dit("n = 11, profile [5,1,1] forces rho = 4", rhos(7, (5, 1, 1))[1][1] == 4)
dit("n = 11, profile [4,2,2,1,1] forces rho = 5 on the valuation-1 pair",
    rhos(10, (4, 2, 2, 1, 1))[1][1] == 5, "checked on 30 and 34: v_2(15+17) = %d" % v2(32))

# ================================================ the lifting machine
def releve(n, t, profil, budget, bavard=False):
    """Bit by bit lifting of the odd parts of the even entries, the FIRST one
    normalised to 1. Constraints: the exact row sums of the rho forced by the
    index identity, and 2**(3m) dividing Delta**m E = sum c (c**2 - 1)**m.
    Returns (level, front size, peak front, reason)."""
    M = (n - 3) // 2
    u = len(profil)
    vmin = min(profil)
    LMAX = 3 * M - vmin + 3
    rh = rhos(t, profil)
    blocs = []
    for f, (mf, R) in rh.items():
        if R < mf - 1:
            return 0, 0, 0, "block f=%d impossible" % f
        blocs.append(([i for i, v in enumerate(profil) if v == f], R))
    libres = list(range(1, u))
    nl = len(libres)
    t1 = time.time()
    front = [tuple([1] * u)]
    pic = 1
    L = 1
    while L < LMAX:
        if time.time() - t1 > budget:
            return L, len(front), pic, "budget"
        nouveau = []
        MK = 1 << (L + 1)
        deuxL = 1 << L
        for base in front:
            for masque in range(1 << nl):
                bs = list(base)
                for q, i in enumerate(libres):
                    if (masque >> q) & 1:
                        bs[i] = bs[i] + deuxL
                bs = tuple(bs)
                bon = True
                for (idx, R) in blocs:
                    if len(idx) < 2:
                        continue
                    for i in idx:
                        connu = 0
                        inconnus = 0
                        for j in idx:
                            if j == i:
                                continue
                            s = (bs[i] + bs[j]) % MK
                            if s:
                                connu += v2(s)
                            else:
                                inconnus += 1
                        if connu + inconnus * (L + 1) > R:
                            bon = False
                            break
                        if inconnus == 0 and connu != R:
                            bon = False
                            break
                    if not bon:
                        break
                if not bon:
                    continue
                C = [(1 << profil[i]) * bs[i] for i in range(u)]
                prec = L + 1 + vmin
                for m in range(1, M + 1):
                    p = min(3 * m, prec)
                    mod = 1 << p
                    D = 0
                    for c in C:
                        D = (D + c * pow((c * c - 1) % mod, m, mod)) % mod
                    if D:
                        bon = False
                        break
                if bon:
                    nouveau.append(bs)
        front = nouveau
        L += 1
        pic = max(pic, len(front))
        if bavard:
            say("           level %2d : front %d" % (L, len(front)))
        if not front:
            return L, 0, pic, "EMPTY FRONT"
        if len(front) > FRONT_CAP:
            return L, len(front), pic, "front > %d" % FRONT_CAP
    return L, len(front), pic, "reached LMAX, front not empty"


# ===================================== G. NEGATIVE CALIBRATION, n = 11
say()
say("-" * 78)
say("G. NEGATIVE CALIBRATION : at n = 11 the two slots at t = 7 must collapse")
say("-" * 78)
for pr in ((5, 1, 1), (3, 2, 2)):
    if reste() < 60:
        saute("n = 11, %s" % str(list(pr)), "budget")
        continue
    L, taille, pic, motif = releve(11, 7, pr, budget=min(60.0, reste() - 20))
    dit("n = 11, t = 7, %s collapses" % str(list(pr)), taille == 0,
        "level %d, peak front %d [%s]" % (L, pic, motif))

# ===================================== H. POSITIVE CALIBRATION
say()
say("-" * 78)
say("H. POSITIVE CALIBRATION : the real witnesses must NOT collapse")
say("-" * 78)
for n, t, pr, nom in ((9, 8, (4, 1, 1, 1, 1), "Letac I"),
                      (11, 10, (4, 2, 2, 1, 1), "CMSV A1 / A2"),
                      (13, 16, (6, 3, 3, 1, 1, 1, 1), "witness n=13")):
    if reste() < 60:
        saute("%s" % nom, "budget")
        continue
    L, taille, pic, motif = releve(n, t, pr, budget=min(40.0, reste() - 20))
    dit("%s survives the lifting" % nom, taille > 0,
        "level %d, front %d [%s]" % (L, taille, motif))

# ===================================== I. THE RESULT, n = 15
say()
say("-" * 78)
say("I. THE RESULT : at n = 15 all four slots at t = 11 collapse")
say("-" * 78)
say("   The slot w = 8, profile [3,2,2,1,1,1,1] dies on the criterion alone:")
r = rhos(11, (3, 2, 2, 1, 1, 1, 1))
m4, R4 = r[1]
dit("n = 15, t = 11, [3,2,2,1,1,1,1] : m = %d, R = %d forbidden" % (m4, R4),
    not bloc_realisable(m4, R4), "R = m is forbidden from m = 3 on")
for pr in ((5, 3, 3), (7, 2, 2), (9, 1, 1)):
    if reste() < 70:
        saute("n = 15, %s" % str(list(pr)), "budget")
        continue
    L, taille, pic, motif = releve(15, 11, pr, budget=min(90.0, reste() - 20))
    dit("n = 15, t = 11, %s collapses" % str(list(pr)), taille == 0,
        "level %d, peak front %d [%s]" % (L, pic, motif))

# ===================================== J. THE LADDER AT n = 15
say()
say("-" * 78)
say("J. THE LADDER AT n = 15")
say("-" * 78)


def sval(n):
    s = 0
    for k in range(3, n - 1, 2):
        s = max(s, min(2 + v2(k - 1), k))
    return s


def partitions(tt, u, mx=None):
    if mx is None:
        mx = tt
    if u == 0:
        if tt == 0:
            yield ()
        return
    for cc in range(min(tt - (u - 1), mx), 0, -1):
        for r2 in partitions(tt - cc, u - 1, cc):
            yield (cc,) + r2


def rowbal(pr, t):
    for j in range(len(pr)):
        L = pr[j] + 1
        for i in range(len(pr)):
            if i != j:
                L += min(pr[i], pr[j]) if pr[i] != pr[j] else pr[j] + 1
        if L > 1 + t:
            return False
    return True


def creneaux(n, t, filtre=True):
    s = sval(n)
    mm = (n - 3) // 2
    out = []
    for w in range(2, n, 2):
        u = n - w
        if u < 1 or u % 2 == 0 or t < u:
            continue
        tau = t - (w - 1)
        if tau < 0:
            continue
        h = w // 2
        ok = (tau >= max(2, s - 1)) if w == 2 else \
             ((tau == 0 and w % 4 == 0) or (tau == h and h % 2 == 0)
              or tau >= h + (-(-h // 2)))
        if not ok:
            continue
        for pr in partitions(t, u):
            if pr.count(1) % 2:
                continue
            cm = min(pr)
            if pr.count(cm) % 2 and cm < 3 * mm:
                continue
            if not rowbal(pr, t):
                continue
            if filtre:
                if not bloc_realisable(w, t):
                    continue
                if any(not bloc_realisable(m, R) for f, (m, R) in rhos(t, pr).items()):
                    continue
            out.append((w, u, pr))
    return out


say("      t   v_2(C)   raw slots   after the criterion")
for t in range(11, 15):
    A = creneaux(15, t, False)
    B = creneaux(15, t, True)
    say("   %4d   %6d   %9d   %19d" % (t, t + 1, len(A), len(B)))
    for (w, u, pr) in A:
        say("            w = %2d, u = %d, profile %-26s %s"
            % (w, u, str(list(pr)), "kept" if (w, u, pr) in B else "killed by the criterion"))
dit("t = 12 is empty", not creneaux(15, 12, False), "raw count 0")
L13 = creneaux(15, 13, True)
dit("t = 13 holds a single slot", len(L13) == 1,
    "w = %d, profile %s" % (L13[0][0], list(L13[0][2])) if L13 else "")

# ===================================================== SUMMARY
say()
say("=" * 78)
say("SUMMARY:  %d PASS,  %d FAIL,  %d SKIPPED    (%.1f s)"
    % (PASS[0], FAIL[0], SKIP[0], time.time() - T0))
say("=" * 78)
say()
say("TABLE OF BOUNDS, ODD SIZES")
say("    n   Kleiman   previous   this verifier")


def vfact(n):
    s, p = 0, 2
    while p <= n:
        s += n // p
        p *= 2
    return s


for n, anc, nv in ((5, 4, 4), (7, 6, 6), (9, 9, 9), (11, 11, 11),
                   (13, 12, 12), (15, 12, 14)):
    say("  %3d   %7d   %8d   %13d%s"
        % (n, vfact(n - 1), anc, nv, "   <== new" if nv > anc else ""))