# -*- coding: utf-8 -*- """ SELF-CONTAINED VERIFIER -- 2-adic valuation of the constant of the Prouhet-Tarry-Escott problem, SYMMETRIC ideal solutions, EVEN sizes. New bound at n = 20. 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 A the six known even solutions, with k_0 and R = v_2(C) - 3 k_0 B Q = P + d, the two polynomials differ ONLY in the constant term C the internal identity P'(u_i) = Q'(u_i), as an identity over the integers D the internal lower bound, and its equality case E modulo 8 every odd number is its own inverse, hence S_i = k_0 xi_i - sum eta_j modulo 8 F EXHAUSTIVE enumeration modulo 8 of the congruences in the proof of Theorem R: no solution for k_0 odd nor for k_0 = 2 modulo 4, while k_0 = 4 survives, so the argument is not vacuous (the theorem itself needs k_0 >= 4) G no known solution has R = 0, so nothing above can contradict them H the bound engine replayed: the seven even bounds of [V7] are recovered EXACTLY, both branches, BEFORE any change I with the new constraint (R = 0 and k_0 >= 4 require k_0 = 0 modulo 4), the bound at n = 20 rises from 18 to 19 and no other bound moves J the known solutions stay compatible with the new bounds SETTINGS -- change nothing but these """ BUDGET_SECONDS = 14400 # total budget; heavy blocks are SKIPPED past it PROFILE_CAP = 8 # ceiling on the valuations in the case sigma = tau; in the case # of distinct sums the larger profile runs over all entries up to # s' + m + 3, which is exhaustive by the saturation lemma of [V7] import time from math import comb from itertools import combinations, combinations_with_replacement as cwr, product from functools import lru_cache from collections import Counter from fractions import Fraction 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 x = abs(int(x)) k = 0 while x % 2 == 0: x //= 2 k += 1 return k def psi(s, t): return 2 * min(s, t) + (3 if s == t else 0) def polmul(X, Y): Z = [0] * (len(X) + len(Y) - 1) for i, a in enumerate(X): if a: for j, b in enumerate(Y): Z[i + j] += a * b return Z say("=" * 78) say("VERIFIER -- PTE, 2-adic valuation, symmetric solutions, even sizes") say("=" * 78) say("budget %d s; profile ceiling %d (sigma = tau); distinct sums up to s' + m + 3" % (BUDGET_SECONDS, PROFILE_CAP)) say() # ============================================================ A say("-" * 78) say("A. THE SIX KNOWN EVEN SOLUTIONS") say("-" * 78) EVEN_SOLS = [ ("Rees-Smyth n=8", [5, 14, 23, 24], [2, 16, 21, 25], 8), ("BLP 2003 (a) n=10", [99, 100, 188, 301, 313], [71, 131, 180, 307, 308], 14), ("BLP 2003 (b) n=10", [103, 189, 366, 452, 515], [18, 245, 331, 471, 508], 11), ("KMS 1999 n=12", [22, 61, 86, 127, 140, 151], [35, 47, 94, 121, 146, 148], 12), ("Broadhurst 2007 n=12", [257, 891, 1109, 1618, 1896, 2058], [472, 639, 1294, 1514, 1947, 2037], 14), ("Choudhry-Wr. n=12", [107, 622, 700, 1075, 1138, 1511], [293, 413, 886, 953, 1180, 1510], 12), ] say(" solution n k_0 v_2(C) R sigma = tau ?") DATA = [] for nom, A, B, NU in EVEN_SOLS: m = len(A) n = 2 * m FA = [x for a in A for x in (a, -a)] FB = [x for b in B for x in (b, -b)] ideal = all(sum(x ** k for x in FA) == sum(x ** k for x in FB) for k in range(1, n)) U = [a * a for a in A] V = [b * b for b in B] pu, pv = 1, 1 for x in U: pu *= x for x in V: pv *= x C = pu - pv k0 = sum(1 for a in A if a % 2) R = NU - 3 * k0 S = tuple(sorted(v2(a) for a in A)) Tt = tuple(sorted(v2(b) for b in B)) DATA.append((nom, A, B, NU, U, V, C, k0, R, S, Tt)) say(" %-22s %4d %4d %6d %3d %s" % (nom, n, k0, NU, R, "oui" if S == Tt else "non")) dit("%s : ideal, v_2(C) = %d" % (nom, NU), ideal and v2(C) == NU) # ============================================================ B say() say("-" * 78) say("B. Q = P + d : THE TWO POLYNOMIALS DIFFER ONLY IN THE CONSTANT TERM") say("-" * 78) say(" The even power sums agreeing up to degree 2(m-1) means e_j(u) = e_j(v)") say(" for j = 1..m-1, so P and Q share every coefficient but the last.") for (nom, A, B, NU, U, V, C, k0, R, S, Tt) in DATA: m = len(U) P = [1] for x in U: P = polmul(P, [1, -x]) Q = [1] for x in V: Q = polmul(Q, [1, -x]) meme = all(P[i] == Q[i] for i in range(m)) d = 1 for j in range(m): d *= (U[0] - V[j]) pareil = True for i in range(m): p = 1 for j in range(m): p *= (U[i] - V[j]) if p != d: pareil = False dit("%s : coefficients 0..%d equal, d is the same for all i" % (nom, m - 1), meme and pareil and v2(d) == NU, "d = %s C, v_2(d) = %d" % ("+" if d == C else "-", v2(d))) # ============================================================ C say() say("-" * 78) say("C. THE INTERNAL IDENTITY P'(u_i) = Q'(u_i), OVER THE INTEGERS") say("-" * 78) say(" P and Q differ by a constant, so P' = Q'. At y = u_i this reads") say(" product_{k != i} (u_i - u_k) = sum_j product_{l != j} (u_i - v_l),") say(" an identity over the integers, equivalently") say(" product_{k != i} (u_i - u_k) = d . sum_j (u_i - v_j)**(-1).") for (nom, A, B, NU, U, V, C, k0, R, S, Tt) in DATA: m = len(U) ok = True for i in range(m): g = 1 for k in range(m): if k != i: g *= (U[i] - U[k]) s = 0 for j in range(m): p = 1 for l in range(m): if l != j: p *= (U[i] - V[l]) s += p if g != s: ok = False for j in range(m): g = 1 for l in range(m): if l != j: g *= (V[j] - V[l]) s = 0 for i in range(m): p = 1 for k in range(m): if k != i: p *= (V[j] - U[k]) s += p if g != s: ok = False dit("%s : the %d internal identities" % (nom, 2 * m), ok) # ============================================================ D say() say("-" * 78) say("D. THE INTERNAL LOWER BOUND AND ITS EQUALITY CASE") say("-" * 78) say(" sum_{k != i} v_2(a_i**2 - a_k**2) >= v_2(C) - max_j v_2(a_i**2 - b_j**2),") say(" with EQUALITY when that maximum is attained exactly once.") for (nom, A, B, NU, U, V, C, k0, R, S, Tt) in DATA: m = len(U) tout = True serres = 0 uniques = 0 for i in range(m): D = sum(v2(U[i] - U[k]) for k in range(m) if k != i) croise = [v2(U[i] - V[j]) for j in range(m)] mu = max(croise) if D < NU - mu: tout = False if D == NU - mu: serres += 1 if croise.count(mu) == 1: uniques += 1 if D != NU - mu: tout = False dit("%s : bound holds, tight exactly where the max is unique" % nom, tout, "%d tight rows, %d unique maxima, out of %d" % (serres, uniques, m)) # ============================================================ E say() say("-" * 78) say("E. MODULO 8 : EVERY ODD NUMBER IS ITS OWN INVERSE") say("-" * 78) dit("x * x = 1 modulo 8 for every odd x", all((x * x) % 8 == 1 for x in (1, 3, 5, 7))) say(" Hence, with S_i = sum_j (xi_i - eta_j)**(-1) over the k_0 unit indices,") say(" S_i = k_0 xi_i - sum_j eta_j modulo 8.") import random random.seed(31) ecarts = 0 for _ in range(20000): kk = random.choice((2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)) xi0 = 2 * random.randrange(0, 2 ** 9) eta = [2 * random.randrange(0, 2 ** 9) + 1 for _ in range(kk)] g = sum(pow((xi0 - e) % 8, -1, 8) for e in eta) % 8 if g != (kk * xi0 - sum(eta)) % 8: ecarts += 1 dit("the congruence holds on 20000 random draws", ecarts == 0, "discrepancies: %d" % ecarts) # ============================================================ F say() say("-" * 78) say("F. R = 0 : EXHAUSTIVE ENUMERATION MODULO 8") say("-" * 78) say(" If R = 0 then every xi_i - eta_j is odd, so the xi and the eta have") say(" OPPOSITE parities; and D_i = 3(k_0-1) + sum v(xi_i - xi_k) >= 4(k_0-1)") say(" >= 3 k_0 = v_2(C) as soon as k_0 >= 4, so v(sum_j (u_i-v_j)**(-1)) >= 0,") say(" which forces S_i = 0 modulo 8 on every row AND every column.") say(" Only the DIFFERENCES matter, so xi_1 is normalised.") say() say(" k_0 k_0 mod 4 configurations survivors verdict") PAIRS = (0, 2, 4, 6) ODDS = (1, 3, 5, 7) SURV = {} for k0 in (2, 3, 4, 5, 6): if reste() < 60: saute("k_0 = %d" % k0, "budget") continue total = 0 surv = 0 for parite in (0, 1): XS = PAIRS if parite == 0 else ODDS ES = ODDS if parite == 0 else PAIRS for xi in product(XS, repeat=k0 - 1): X = (XS[0],) + xi for eta in product(ES, repeat=k0): total += 1 bon = True for i in range(k0): if sum((X[i] - e) % 8 for e in eta) % 8: bon = False break if not bon: continue for e in eta: if sum((e - X[i]) % 8 for i in range(k0)) % 8: bon = False break if bon: surv += 1 SURV[k0] = surv say(" %6d %9d %14d %9d %s" % (k0, k0 % 4, total, surv, "IMPOSSIBLE" if surv == 0 else "possible")) dit("the congruences have no solution for k_0 = 3 and 5", SURV.get(3) == 0 and SURV.get(5) == 0) dit("the congruences have no solution for k_0 = 2 and 6 (k_0 = 2 modulo 4)", SURV.get(2) == 0 and SURV.get(6) == 0) dit("the congruences have solutions for k_0 = 4, so the argument is not vacuous", SURV.get(4, 0) > 0, "%d survivors" % SURV.get(4, 0)) # ============================================================ G say() say("-" * 78) say("G. NO KNOWN SOLUTION HAS R = 0") say("-" * 78) dit("every known solution has R >= 2", all(x[8] >= 2 for x in DATA), "R values: %s" % [x[8] for x in DATA]) # ============================================= the bound engine @lru_cache(maxsize=None) def node(a, b, cap): if a == 0 and b == 0: return frozenset() if b == 0: return frozenset({(0, None)}) if a == 0: return frozenset({(None, 0)}) base = set() for a1 in range(a + 1): for b1 in range(b + 1): a2, b2 = a - a1, b - b1 if (a1, b1) == (a, b) or (a2, b2) == (a, b): continue for (x1, y1) in node(a1, b1, cap): for (x2, y2) in node(a2, b2, cap): xs = [z for z in (x1, x2) if z is not None] ys = [z for z in (y1, y2) if z is not None] if len(xs) == 2 and xs[0] != xs[1]: continue if len(ys) == 2 and ys[0] != ys[1]: continue x = b + xs[0] if xs else None y = a + ys[0] if ys else None if (x is not None and x > cap) or (y is not None and y > cap): continue base.add((x, y)) out, fr = set(base), set(base) while fr: nx_ = set() for (x, y) in fr: nx = None if x is None else x + b ny = None if y is None else y + a if (nx is not None and nx > cap) or (ny is not None and ny > cap): continue if (nx, ny) not in out: nx_.add((nx, ny)) out |= nx_ fr = nx_ return frozenset(out) @lru_cache(maxsize=None) def R_realisable(k, cap=40): out = set() for a1 in range(k + 1): for b1 in range(k + 1): a2, b2 = k - a1, k - b1 L1 = list(node(a1, b1, cap)) or [(None, None)] L2 = list(node(a2, b2, cap)) or [(None, None)] for (x1, y1) in L1: for (x2, y2) in L2: if (a1 > 0) != (x1 is not None): continue if (a2 > 0) != (x2 is not None): continue if (b1 > 0) != (y1 is not None): continue if (b2 > 0) != (y2 is not None): continue xs = [z for z in (x1, x2) if z is not None] ys = [z for z in (y1, y2) if z is not None] if len(xs) == 2 and xs[0] != xs[1]: continue if len(ys) == 2 and ys[0] != ys[1]: continue if xs and ys and xs[0] == ys[0]: out.add(xs[0]) return frozenset(z for z in out if z <= cap) def index_bound(S, T): return max(max(sum(psi(s, t) for s in S) for t in T), max(sum(psi(s, t) for t in T) for s in S)) def ej_predicted(sig, j): m = len(sig) sums = sorted(sum(sig[i] for i in I) for I in combinations(range(m), j)) S = sums[0] N = sum(1 for x in sums if x == S) sup = [x for x in sums if x > S] dd = (sup[0] - S) if sup else 10 ** 9 if N % 8 and v2(N) < 2 * dd: return 2 * S + v2(N), True return 2 * S + min(3, 2 * dd), False def ej_fast(sig, j): """Same output as ej_predicted, computed without enumerating the j-subsets: the minimal weight is the sum of the j smallest entries, N is a binomial coefficient at the boundary value b, and the gap to the next attainable weight is the least increase obtained by one exchange.""" a = sorted(sig) S = sum(a[:j]) b = a[j-1] cnt_b = a.count(b) below = sum(1 for x in a if x < b) N = comb(cnt_b, j - below) cands = [] if j - below < cnt_b: lower = [x for x in a[:j] if x < b] if lower: cands.append(b - max(lower)) above = [x for x in a if x > b] if above: cands.append(min(above) - b) d = min(cands) if cands else 10**9 if N % 8 and v2(N) < 2*d: return 2*S + v2(N), True return 2*S + min(3, 2*d), False import random as _rnd _rnd.seed(12345) _bad = 0 for _ in range(20000): _m = _rnd.randint(2, 10) _sig = [_rnd.choice((0, 0, 0, 1, 1, 2, 3, 5, 9, 14)) for _ in range(_m)] for _j in range(1, _m): if ej_fast(_sig, _j) != ej_predicted(_sig, _j): _bad += 1 dit("fast predictor = direct predictor on 20000 random profiles", _bad == 0, "discrepancies: %d" % _bad) def feasible(a, b): for j in range(1, len(a)): x, dx = ej_fast(a, j) y, dy = ej_fast(b, j) if dx and dy and x != y: return False if dx and (not dy) and x < y: return False if dy and (not dx) and y < x: return False return True def lo_de(k0, neuf): # k_0 odd: R >= 1 by Lemma 3.4 of [V7] when sigma = tau, by parity when the sums differ. # new: R >= 1 when k_0 >= 4 and k_0 = 2 modulo 4 (Theorem R of the note). if neuf: return 1 if (k0 % 2 or (k0 % 4 == 2 and k0 >= 4)) else 0 return 1 if k0 % 2 else 0 def full_bound(S, T, neuf): k0 = Counter(S)[0] b = index_bound(S, T) b = max(b, 2 * sum(S) + 3) if S == T else max(b, 2 * min(sum(S), sum(T))) lo = lo_de(k0, neuf) cand = [3 * k0 + R for R in sorted(R_realisable(k0)) if R >= lo and 3 * k0 + R >= b] return min(cand) if cand else b def profiles(m, cap): return [P for P in cwr(range(cap + 1), m) if min(P) == 0] def branche_egale(n, neuf): m = n // 2 return min(full_bound(S, S, neuf) for S in profiles(m, PROFILE_CAP)) def branche_inegale(n, cible, neuf): m = n // 2 taus = [] for T in cwr(range((cible - 1) // 2 + 1), m): if min(T) != 0: continue st = sum(T) if 2 * st >= cible or st == 0: continue k0 = Counter(T)[0] R = 2 * st - 3 * k0 if R < lo_de(k0, neuf) or R not in R_realisable(k0): continue taus.append(T) bad = [] for T in taus: st = sum(T) k0 = Counter(T)[0] alpha = [w for w in range(st + m + 4) if sum(psi(w, t) for t in T) <= 2 * st] if 0 not in alpha: continue pos = [w for w in alpha if w >= 1] for P in cwr(pos, m - k0): S = tuple(sorted((0,) * k0 + P)) if Counter(S)[0] != k0 or sum(S) <= st: continue if 2 * st < index_bound(S, T): continue if feasible(S, T): bad.append((2 * st, S, T)) return min([x[0] for x in bad], default=None) # ============================================================ H say() say("-" * 78) say("H. THE BOUND ENGINE REPLAYED, BOTH BRANCHES, BEFORE ANY CHANGE") say("-" * 78) PUBLIE = {8: 8, 10: 11, 12: 12, 14: 12, 16: 16, 18: 18, 20: 18} # the bounds of [V7] say(" n published sigma = tau distinct sums engine") ANCIEN = {} for n in (8, 10, 12, 14, 16, 18, 20): if reste() < 90: saute("n = %d" % n, "budget") continue be = branche_egale(n, False) mi = branche_inegale(n, PUBLIE[n] + 2, False) b = min(be, mi) if mi is not None else be ANCIEN[n] = b say(" %4d %9d %11d %13s %6d" % (n, PUBLIE[n], be, mi if mi is not None else "-", b)) dit("n = %2d : the engine recovers the published %d" % (n, PUBLIE[n]), b == PUBLIE[n], "computed %d" % b) # ============================================================ I say() say("-" * 78) say("I. WITH THE NEW CONSTRAINT : R = 0 AND k_0 >= 4 REQUIRE k_0 = 0 MODULO 4") say("-" * 78) say(" n before sigma = tau distinct sums after gain") NOUVEAU = {} for n in (8, 10, 12, 14, 16, 18, 20): if n not in ANCIEN: continue if reste() < 60: saute("n = %d" % n, "budget") continue be = branche_egale(n, True) mi = branche_inegale(n, PUBLIE[n] + 2, True) b = min(be, mi) if mi is not None else be NOUVEAU[n] = b g = b - ANCIEN[n] say(" %4d %6d %11d %13s %5d %s" % (n, ANCIEN[n], be, mi if mi is not None else "-", b, ("+%d" % g) if g > 0 else "0")) for n in sorted(NOUVEAU): if n == 20: dit("n = 20 : the bound rises from 18 to 19", ANCIEN.get(20) == 18 and NOUVEAU.get(20) == 19, "before %s, after %s" % (ANCIEN.get(20), NOUVEAU.get(20))) else: dit("n = %2d : unchanged at %s" % (n, ANCIEN[n]), NOUVEAU[n] == ANCIEN[n], "after %d" % NOUVEAU[n]) # ============================================================ J say() say("-" * 78) say("J. THE KNOWN SOLUTIONS STAY COMPATIBLE") say("-" * 78) for (nom, A, B, NU, U, V, C, k0, R, S, Tt) in DATA: bn = full_bound(S, Tt, True) dit("%s : new bound %d <= v_2(C) = %d" % (nom, bn, NU), bn <= NU, "sigma %s, tau %s" % (list(S), list(Tt))) # ===================================================== 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, EVEN SIZES") say(" n Kleiman published this verifier") def vfact(n): s, p = 0, 2 while p <= n: s += n // p p *= 2 return s for n in (8, 10, 12, 14, 16, 18, 20): nv = NOUVEAU.get(n, None) say(" %3d %7d %9d %13s%s" % (n, vfact(n - 1), PUBLIE[n], nv if nv is not None else "-", " <== new" if (nv is not None and nv > PUBLIE[n]) else ""))