# -*- coding: utf-8 -*-
"""
VERIF_PTE_2adic_sharp_10_12.py
David Niedbala Giraudin, ORCID 0009-0009-1526-1178 -- 2026-10-11
Technical assistance: Claude (Anthropic).

Statement checked
-----------------
There exist 2-adic ideal solutions of the Prouhet-Tarry-Escott problem
  - of size 10 with 2**9  exactly dividing the constant C,
  - of size 12 with 2**11 exactly dividing the constant C.
Hence the lower bounds 9 <= nu_2(C_10) and 11 <= nu_2(C_12) of
R. Bera and S. G. Nair (Acta Arith. 213 (2024), 75-95) are attained 2-adically.
The size-10 witness lies in their Case 2.3(ii), sub-case (I), and satisfies
every divisibility property they derive there.

A 2-adic ideal solution of size n is a pair of monic polynomials F, G of
degree n, with all n roots of each in Z_2, such that F - G is a nonzero constant C.

Certificate
-----------
  F10(z) = prod_{a in {0,4,8,16,17}} (z**2 + z - 2a)
  G10(z) = prod_{b in {1,2,10,14,18}} (z**2 + z - 2b)          F10 - G10 = 161280 = 2**9 * 315
  F12(z) = prod_{a in {0,4,5}} ((z**2+z)**2 + 2(z**2+z) - 8a)
  G12(z) = prod_{b in {1,2,6}} ((z**2+z)**2 + 2(z**2+z) - 8b)     F12 - G12 = 6144   = 2**11 * 3
Calibration at size 8:
  F8, G8 built in the same way on {0,4,7,11} / {1,2,9,10}:     F8 - G8 = -2880    = -2**6 * 45

Pure Python, no dependency, no input file, all parameters fixed below.
Runs in well under a second. Prints as it goes.
"""

import sys

PREC = 100          # working 2-adic precision for the Hensel certificates (bits)
MOD = 2 ** PREC

PASS = 0
FAIL = 0


def report(label, ok, detail=""):
    global PASS, FAIL
    if ok:
        PASS += 1
    else:
        FAIL += 1
    tag = "PASS" if ok else "FAIL"
    print("  [" + tag + "] " + label + ("   (" + detail + ")" if detail else ""), flush=True)


def v2(n):
    n = abs(int(n))
    if n == 0:
        return None
    k = 0
    while n % 2 == 0:
        n //= 2
        k += 1
    return k


def v2_mod(n):
    """valuation of a class mod 2**PREC (returns PREC for the zero class)"""
    n = int(n) % MOD
    if n == 0:
        return PREC
    return v2(n)


# ---------------------------------------------------------------- polynomials
# a polynomial is the list of its integer coefficients, constant term first

def pmul(a, b):
    r = [0] * (len(a) + len(b) - 1)
    for i in range(len(a)):
        for j in range(len(b)):
            r[i + j] += int(a[i]) * int(b[j])
    return r


def padd(a, b):
    n = max(len(a), len(b))
    a = list(a) + [0] * (n - len(a))
    b = list(b) + [0] * (n - len(b))
    return [int(x) + int(y) for x, y in zip(a, b)]


def pneg(a):
    return [-int(x) for x in a]


def ptrim(a):
    a = list(a)
    while len(a) > 1 and a[-1] == 0:
        a.pop()
    return a


def pprod(polys):
    r = [1]
    for p in polys:
        r = pmul(r, p)
    return r


def peval(p, z, mod=None):
    acc = 0
    for c in reversed(p):
        acc = acc * int(z) + int(c)
        if mod is not None:
            acc %= mod
    return acc


def pderiv(p):
    return [int(i) * int(p[i]) for i in range(1, len(p))]


# ---------------------------------------------------------------- self-test
print("=" * 72, flush=True)
print("SELF-TEST (arithmetic of this interpreter)", flush=True)
print("=" * 72, flush=True)
st = (2 ** 10 == 1024 and (7 * 13) % 8 == 3 and v2(161280) == 9 and v2(6144) == 11
      and pmul([1, 1], [-1, 1]) == [-1, 0, 1] and 17 // 4 == 4 and (-3) % 8 == 5)
if not st:
    print("  self-test FAILED: this interpreter does not do integer arithmetic as expected.", flush=True)
    print("  Run this file with plain Python 3 (not inside a notebook cell).", flush=True)
    sys.exit(1)
print("  self-test OK", flush=True)


# ---------------------------------------------------------------- seeds
def is_ideal(X, Y):
    return all(sum(int(x) ** k for x in X) == sum(int(y) ** k for y in Y) for k in range(1, len(X)))


def const_of(X, Y):
    d = ptrim(padd(pprod([[-int(x), 1] for x in X]), pneg(pprod([[-int(y), 1] for y in Y]))))
    return d[0] if len(d) == 1 else None


SEEDS = {
    3: ([0, 4, 5], [1, 2, 6]),
    4: ([0, 4, 7, 11], [1, 2, 9, 10]),
    5: ([0, 4, 8, 16, 17], [1, 2, 10, 14, 18]),
}
print("=" * 72, flush=True)
print("INTEGER IDEAL SOLUTIONS USED AS SEEDS (classical)", flush=True)
print("=" * 72, flush=True)
for m in (3, 4, 5):
    X, Y = SEEDS[m]
    C = const_of(X, Y)
    report("size %d: %s =_%d %s is ideal" % (m, X, m - 1, Y), is_ideal(X, Y) and C is not None,
           "C = %s, nu_2(C) = %s" % (C, v2(C) if C else None))


# ---------------------------------------------------------------- the certificates
Z2Z = [0, 1, 1]                                   # z**2 + z
QUART = padd(pmul(Z2Z, Z2Z), [0, 2, 2])           # (z**2+z)**2 + 2(z**2+z)

def quad_factor(a):                               # z**2 + z - 2a
    return padd(Z2Z, [-2 * int(a)])

def quart_factor(a):                              # (z**2+z)**2 + 2(z**2+z) - 8a
    return padd(QUART, [-8 * int(a)])


def two_adic_roots_T(r):
    """the two solutions in Z/2**PREC of z**2 + z = 2r (r given mod 2**PREC), bit by bit"""
    out = []
    for z in (0, 1):
        for k in range(1, PREC):
            if (z * z + z - 2 * r) % (2 ** (k + 1)) != 0:
                z += 2 ** k
        if (z * z + z - 2 * r) % MOD != 0:
            return None
        out.append(z % MOD)
    return out


def hensel_certify(F, roots, label):
    """Each approximate root t satisfies v(F(t)) > 2 v(F'(t)); Hensel then gives a true root
    zeta in Z_2 with v(zeta - t) >= v(F(t)) - v(F'(t)). The n balls are pairwise disjoint,
    so F has n distinct roots in Z_2, i.e. F splits completely over Z_2."""
    n = len(F) - 1
    dF = pderiv(F)
    ok = len(roots) == n
    radius = []
    worst_sep = 0
    for t in roots:
        a = v2_mod(peval(F, t, MOD))
        b = v2_mod(peval(dF, t, MOD))
        if not (a > 2 * b):
            ok = False
        radius.append(a - b)
    for i in range(n):
        for j in range(i + 1, n):
            s = v2_mod(roots[i] - roots[j])
            worst_sep = max(worst_sep, s)
            if not (s < min(radius[i], radius[j])):
                ok = False
    report(label + ": %d Hensel-certified, pairwise distinct roots in Z_2" % n, ok,
           "max v(t_i - t_j) = %d, min certified radius = %d" % (worst_sep, min(radius) if radius else -1))
    return ok


def build(size_label, Xseed, Yseed, factor, expected_C, expected_v):
    print("=" * 72, flush=True)
    print(size_label, flush=True)
    print("=" * 72, flush=True)
    F = pprod([factor(a) for a in Xseed])
    G = pprod([factor(b) for b in Yseed])
    n = len(F) - 1
    report("F and G are monic of degree %d with integer coefficients" % n,
           F[-1] == 1 and G[-1] == 1 and len(G) == len(F))
    D = ptrim(padd(F, pneg(G)))
    is_const = len(D) == 1
    report("F - G is the constant %s" % (D[0] if is_const else D), is_const and D[0] == expected_C)
    report("nu_2(F - G) = %d" % expected_v, is_const and v2(D[0]) == expected_v,
           "%d = 2**%d * %d" % (D[0], v2(D[0]), D[0] // 2 ** v2(D[0])) if is_const else "")
    return F, G, (D[0] if is_const else None)


def roots_level1(seed):
    out = []
    for a in seed:
        out += two_adic_roots_T(int(a) % MOD)
    return out


def roots_level2(seed):
    out = []
    for a in seed:
        for r in two_adic_roots_T(int(a) % MOD):     # u = 2r is a root of u*u + 2u - 8a, with u = z*z + z
            out += two_adic_roots_T(r)
    return out


# -------- size 8 (calibration: Filaseta-Markovich lower bound 6)
F8, G8, C8 = build("SIZE 8 (calibration: Filaseta-Markovich 2017 prove nu_2(C_8) >= 6)",
                   SEEDS[4][0], SEEDS[4][1], quad_factor, -2880, 6)
hensel_certify(F8, roots_level1(SEEDS[4][0]), "F8")
hensel_certify(G8, roots_level1(SEEDS[4][1]), "G8")

# -------- size 10
F10, G10, C10 = build("SIZE 10 (Bera-Nair: 9 <= nu_2(C_10) <= 11)",
                      SEEDS[5][0], SEEDS[5][1], quad_factor, 161280, 9)
RF10 = roots_level1(SEEDS[5][0])
RG10 = roots_level1(SEEDS[5][1])
hensel_certify(F10, RF10, "F10")
hensel_certify(G10, RG10, "G10")

# -------- size 12
F12, G12, C12 = build("SIZE 12 (Bera-Nair: 11 <= nu_2(C_12) <= 12)",
                      SEEDS[3][0], SEEDS[3][1], quart_factor, 6144, 11)
RF12 = roots_level2(SEEDS[3][0])
RG12 = roots_level2(SEEDS[3][1])
hensel_certify(F12, RF12, "F12")
hensel_certify(G12, RG12, "G12")


# ---------------------------------------------------------------- position in Bera-Nair's case analysis
print("=" * 72, flush=True)
print("POSITION OF THE SIZE-10 WITNESS IN BERA-NAIR'S CASE 2.3(ii)", flush=True)
print("=" * 72, flush=True)
# their X is our G10 side (valuation profile 0,0,0,0,0,1,2,2,2,2), their Y our F10 side
T0 = 2 ** 7                                  # translation avoiding the root 0; touches no valuation < 7
X = [(t + T0) % MOD for t in RG10]
Y = [(t + T0) % MOD for t in RF10]
px = sorted(v2_mod(t) for t in X)
py = sorted(v2_mod(t) for t in Y)
report("profile of X: five odd, 2 || x6, 4 || x7..x10", px == [0, 0, 0, 0, 0, 1, 2, 2, 2, 2], str(px))
report("profile of Y: five odd, 2 || y6, 8 | y7..y10", py[:6] == [0, 0, 0, 0, 0, 1] and min(py[6:]) >= 3, str(py))
x6 = [t for t in X if v2_mod(t) == 1][0]
y6 = [t for t in Y if v2_mod(t) == 1][0]
report("x6 = y6 mod 32, exactly (their consequence of 2**9 || C)", v2_mod(x6 - y6) == 5,
       "v_2(x6 - y6) = %d" % v2_mod(x6 - y6))
def xprime(S):
    return [(t // 4) % 8 for t in S if v2_mod(t) == 2]
n1 = sum(1 for c in xprime(X) if c % 4 == 1)
report("an odd number of x'_j = x_j/4 are 1 mod 4 (their sub-case (ii))", n1 % 2 == 1, "%d of them" % n1)
if n1 == 3:                                  # their normalisation: translate by 8
    X = [(t + 8) % MOD for t in X]
    Y = [(t + 8) % MOD for t in Y]
one = [c for c in xprime(X) if c % 4 == 1]
if one and one[0] % 8 == 5:                  # their normalisation: translate by 16
    X = [(t + 16) % MOD for t in X]
    Y = [(t + 16) % MOD for t in Y]
cls = sorted(xprime(X))
c8 = sum(1 for t in Y if v2_mod(t) == 3)
others = sorted(c for c in cls if c != 1)
report("after their normalisation: classes of x'_j mod 8 = %s, and %d y_j with 8 || y_j" % (cls, c8),
       one is not None and len(one) == 1)
report("this is their sub-case (I): {3,3,7} or {7,7,7}, with exactly three y_j such that 8 || y_j",
       others in ([3, 3, 7], [7, 7, 7]) and c8 == 3)

# every quantity their sub-case (I), and (17) of sub-case (ii)(b), derives for such a solution
def coeffs(roots):
    """coefficients, constant term first, of the monic product over the given roots, mod 2**PREC"""
    p = [1]
    for r in roots:
        q = [0] * (len(p) + 1)
        for i, c in enumerate(p):
            q[i + 1] = (q[i + 1] + c) % MOD
            q[i] = (q[i] - c * r) % MOD
        p = q
    return p

xs = sorted(X, key=v2_mod)
ys = sorted(Y, key=v2_mod)
u = coeffs(xs[6:])          # their u_j, on x_7 .. x_10
v = coeffs(ys[6:])          # their v_j, on y_7 .. y_10
a = coeffs(X)               # their a_j
b = coeffs(Y)               # their b_j
report("sub-case (I): nu_2(u_1) = 7", v2_mod(u[1]) == 7, "= %d" % v2_mod(u[1]))
report("sub-case (I): nu_2(u_2) >= 7", v2_mod(u[2]) >= 7, "= %d" % v2_mod(u[2]))
report("sub-case (I): nu_2(a_2) = 7", v2_mod(a[2]) == 7, "= %d" % v2_mod(a[2]))
report("(17): nu_2(v_0) >= 13, nu_2(v_1) = 9, nu_2(v_2) = 6, hence nu_2(b_2) = 7",
       v2_mod(v[0]) >= 13 and v2_mod(v[1]) == 9 and v2_mod(v[2]) == 6 and v2_mod(b[2]) == 7,
       "v_0:%d v_1:%d v_2:%d b_2:%d" % (v2_mod(v[0]), v2_mod(v[1]), v2_mod(v[2]), v2_mod(b[2])))


# ---------------------------------------------------------------- verdict
print("=" * 72, flush=True)
print("TOTAL: %d PASS, %d FAIL" % (PASS, FAIL), flush=True)
if FAIL == 0:
    print("All checks passed. F10, G10 (resp. F12, G12) split completely over Z_2 with distinct", flush=True)
    print("roots and differ by a constant of 2-adic valuation 9 (resp. 11).", flush=True)
print("=" * 72, flush=True)