# -*- coding: utf-8 -*- """ SELF-CONTAINED VERIFIER - 2-adic valuation of the constant of the Prouhet-Tarry-Escott problem, SYMMETRIC ideal solutions, even and odd sizes. David Niedbala Giraudin, ORCID 0009-0009-1526-1178. Technical assistance: Claude (Anthropic). Replays every numerical assertion from scratch and prints PASS or FAIL. Pure Python: no dependency (no numpy, no scipy, no sympy). No command-line argument. No intermediate file. Prints as it goes. SETTINGS - change nothing but these """ BUDGET_SECONDS = 7200 # total budget; blocks that are too heavy are SKIPPED PROFILE_CAP = 8 # ceiling for the profiles in the case sigma = tau (entries >= 9 # force nu_2(C) >= 21); in the case of distinct sums the larger # profile is enumerated with all entries up to s' + r + 3, which is # exhaustive by the saturation lemma (Lemma lem:sat of the paper) EVEN_SIZES = (6, 8, 10, 12, 14, 16, 18, 20) 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 T0 = time.time() PASS = 0 FAIL = 0 SKIP = 0 def say(txt=""): print(txt, flush=True) def verdict(name, ok, detail=""): global PASS, FAIL if ok: PASS += 1 say(" PASS %-58s %s" % (name, detail)) else: FAIL += 1 say(" FAIL %-58s %s" % (name, detail)) def skip_block(name, why): global SKIP SKIP += 1 say(" SKIP %-58s %s" % (name, why)) def remaining(): return BUDGET_SECONDS - (time.time() - T0) # ----------------------------------------------------------------- tools 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 esym(X): m = len(X) e = [1] + [0]*m for x in X: for j in range(m, 0, -1): e[j] += e[j-1]*x return e def sqrt2adic(a, prec): """2-adic square root of an a == 1 mod 8, modulo 2 to the power prec.""" assert a % 8 == 1 x = 1 for t in range(3, prec): if (x*x - a) % (1 << (t+1)): x += 1 << (t-1) assert (x*x - a) % (1 << (t+1)) == 0 return x % (1 << prec) def psi(s, t): return 2*min(s, t) + (3 if s == t else 0) # ------------------------------------------------------------------ data # half-families of the known symmetric ideal solutions (CMSV, arXiv:2304.11254) EVEN_SOLS = [ ("Rees-Smyth (4.7) 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. 2008 n=12", [107, 622, 700, 1075, 1138, 1511], [293, 413, 886, 953, 1180, 1510], 12), ] # odd sizes: B = -A. Letac 1942 (integral); CMSV (nine integers + a quadratic pair) ODD_SOLS = [ ("Letac 1942 (a) n=9 ", [-98, -82, -58, -34, 13, 16, 69, 75, 99], None, None, 9), ("Letac 1942 (b) n=9 ", [-174, -148, -132, -50, -8, 63, 119, 161, 169], None, None, 10), ("CMSV A1 n=11", [-95, -52, -48, -13, -9, 30, 34, 61, 65], 27, -6468, 11), ("CMSV A2 n=11", [-589, -363, -185, -170, 41, 234, 355, 548, 588], -459, -47472, 11), ] EVEN_BOUNDS = {6: 6, 8: 8, 10: 11, 12: 12, 14: 12, 16: 16, 18: 18, 20: 18} EXACT_AT = {8: 8, 10: 11, 12: 12} def kleiman(n): s, k = 0, 2 while k <= n-1: s += (n-1)//k k *= 2 return s say("="*76) say("VERIFIER - PTE, 2-adic valuation, symmetric solutions") say("="*76) say("budget %d s; profile ceiling %d (case sigma = tau); distinct sums up to s' + r + 3" % (BUDGET_SECONDS, PROFILE_CAP)) say() # ============================================================ PART A say("-"*76) say("A. THE KNOWN SOLUTIONS") say("-"*76) for lab, A, B, NU in EVEN_SOLS: m = len(A) FA = [x for a in A for x in (a, -a)] FB = [x for b in B for x in (b, -b)] n = 2*m ok_id = all(sum(x**k for x in FA) == sum(x**k for x in FB) for k in range(1, n)) pa, pb = 1, 1 for a in A: pa *= a*a for b in B: pb *= b*b C = pa - pb ok_nu = (v2(C) == NU) U = [a*a for a in A] V_ = [b*b for b in B] ok_red = all(sum(u**j for u in U) == sum(w**j for w in V_) for j in range(1, m)) S = sorted(v2(a) for a in A) T = sorted(v2(b) for b in B) verdict("%s: ideal, nu_2(C) = %d, reduction" % (lab, NU), ok_id and ok_nu and ok_red, "sigma %s, tau %s%s" % (S, T, "" if S == T else " (distinct profiles)")) say() for lab, ints, total, prod, NU in ODD_SOLS: n = len(ints) + (0 if prod is None else 2) if prod is None: ok_id = all(sum(x**k for x in ints) == sum((-x)**k for x in ints) for k in range(1, n)) p = 1 for a in ints: p *= a ok_nu = (v2(2*p) == NU) verdict("%s: ideal (B = -A), nu_2(C) = %d" % (lab, NU), ok_id and ok_nu) else: ok_s = (total == -sum(ints)) PREC = 44 MM = 1 << PREC disc = (total*total - 4*prod) % MM r = sqrt2adic(disc, PREC) r1 = ((total + r) % MM) >> 1 r2 = ((total - r) % MM) >> 1 p = 1 for a in ints: p *= a ok_nu = (v2(2*p*prod) == NU) verdict("%s: consistent sum, nu_2(C) = %d" % (lab, NU), ok_s and ok_nu, "valuations of the pair: %d and %d" % (v2(r1), v2(r2))) # ============================================================ PART B say() say("-"*76) say("B. THE TWO REGULAR SYSTEMS (row law)") say("-"*76) for lab, A, B, NU in EVEN_SOLS: Ao = [a for a in A if a % 2] Ae = [a for a in A if a % 2 == 0] Bo = [b for b in B if b % 2] Be = [b for b in B if b % 2 == 0] Mo = [[v2(a*a-b*b) for b in Bo] for a in Ao] Me = [[v2(a*a-b*b) for b in Be] for a in Ae] ro = (all(sum(r) == NU for r in Mo) and all(sum(Mo[i][k] for i in range(len(Ao))) == NU for k in range(len(Bo)))) re = (all(sum(r) == NU for r in Me) and all(sum(Me[i][k] for i in range(len(Ae))) == NU for k in range(len(Be)))) k0 = len(Ao) R = NU - 3*k0 xi = [(a*a-1)//8 for a in Ao] eta = [(b*b-1)//8 for b in Bo] Mm = [[v2(x-y) for y in eta] for x in xi] rm = (all(sum(r) == R for r in Mm) and all(sum(Mm[i][k] for i in range(k0)) == R for k in range(k0))) verdict("%s: odd system, even system, M = R" % lab, ro and re and rm, "k_0 = %d, R = %d" % (k0, R)) # ============================================================ PART C say() say("-"*76) say("C. THE INDEX BOUND AND THE PREDICTOR") say("-"*76) 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)) tight_count = 0 for lab, A, B, NU in EVEN_SOLS: S = [v2(a) for a in A] T = [v2(b) for b in B] b = index_bound(S, T) if b == NU: tight_count += 1 verdict("%s: bound %2d <= nu_2(C) = %2d" % (lab, b, NU), b <= NU, "TIGHT" if b == NU else "") verdict("the bound is tight on exactly four of the six", tight_count == 4, "tight: %d" % tight_count) 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] d = (sup[0]-S) if sup 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 mismatch = 0 for lab, A, B, NU in EVEN_SOLS: m = len(A) for X in (A, B): sig = tuple(sorted(v2(x) for x in X)) e = esym([x*x for x in X]) for j in range(1, m+1): pred, det = ej_predicted(sig, j) true_val = v2(e[j]) if det and pred != true_val: mismatch += 1 if (not det) and pred > true_val: mismatch += 1 verdict("predictor of nu_2(e_j) exact on the twelve half-solutions", mismatch == 0, "discrepancies: %d" % mismatch) def polygon(coeffs): pts = [(i, v2(c)) for i, c in enumerate(coeffs) if c != 0] h = [] for p in pts: while len(h) >= 2: (x1, y1), (x2, y2) = h[-2], h[-1] if (y2-y1)*(p[0]-x1) >= (p[1]-y1)*(x2-x1): h.pop() else: break h.append(p) sl = [] for (x1, y1), (x2, y2) in zip(h, h[1:]): sl += [(y1-y2)//(x2-x1)]*(x2-x1) return sorted(sl) okNP = True for lab, A, B, NU in EVEN_SOLS: m = len(A) for X in (A, B): e = esym([x*x for x in X]) f = [((-1)**(m-i))*e[m-i] for i in range(m+1)] if polygon(f) != sorted(2*v2(x) for x in X): okNP = False verdict("slopes of the Newton polygon = 2 sigma", okNP) # ============================================================ PART D say() say("-"*76) say("D. THE THEOREM FOR THE EVEN SIZES") say("-"*76) 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 verdict("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 @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 full_bound(S, T): 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 = 1 if k0 % 2 else 0 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] # --- calibration t1 = time.time() _ = [full_bound(P, P) for P in profiles(5, PROFILE_CAP)] dt5 = time.time() - t1 say(" calibration: m = 5 handled in %.2f s" % dt5) for n in EVEN_SIZES: m = n//2 target = EVEN_BOUNDS[n] nb = len(profiles(m, PROFILE_CAP)) predicted = dt5 * nb / max(1, len(profiles(5, PROFILE_CAP))) if predicted > remaining() * 0.45: skip_block("n = %d, case sigma = tau" % n, "predicted %.0f s, remaining %.0f s" % (predicted, remaining())) continue worst = None for S in profiles(m, PROFILE_CAP): b = full_bound(S, S) if worst is None or b < worst[0]: worst = (b, S) guard = 2*(PROFILE_CAP+1) + 3 verdict("n = %2d, case sigma = tau: minimum %2d" % (n, worst[0]), worst[0] >= target and guard >= target, "profile %s; guard %d" % (list(worst[1]), guard)) # --- case of distinct sums: sigma' is the profile of smaller sum s', nu_2(C) = 2 s'. # All pairs with 2 s' <= bound are collected (entries of sigma up to s' + r + 3). PAIRS = {} for n in EVEN_SIZES: m = n//2 target = EVEN_BOUNDS[n] if remaining() < 30: skip_block("n = %d, case of distinct sums" % n, "budget exhausted") continue t1 = time.time() found = [] ntau = 0 for T in cwr(range(target//2 + 1), m): if min(T) != 0: continue st = sum(T) if 2*st > target or st == 0: continue k0 = Counter(T)[0] R = 2*st - 3*k0 if R < (1 if k0 % 2 else 0) or R not in R_realisable(k0): continue ntau += 1 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): found.append((2*st, S, T)) PAIRS[n] = found below = [x for x in found if x[0] < target] verdict("n = %2d, distinct sums: no pair below %2d" % (n, target), not below, "%d profiles sigma' ; %.0f s%s" % (ntau, time.time() - t1, "" if not below else "; MINIMUM %d" % min(x[0] for x in below))) for n, val in sorted(EXACT_AT.items()): attained = [lab for lab, A, B, NU in EVEN_SOLS if 2*len(A) == n and NU == val] verdict("n = %2d: the bound %2d is attained" % (n, val), bool(attained), "by %s" % ", ".join(x.split()[0] for x in attained)) # ============================================================ PART E say() say("-"*76) say("E. THE ODD SIZES") say("-"*76) PREC = 44 MM = 1 << PREC for lab, ints, total, prod, NU in ODD_SOLS: if prod is None: el = list(ints) else: disc = (total*total - 4*prod) % MM r = sqrt2adic(disc, PREC) el = list(ints) + [((total + r) % MM) >> 1, ((total - r) % MM) >> 1] vs = sorted((v2(x) for x in el), reverse=True) V = sum(vs) ok_V = (1 + V == NU) simple_summit = (vs.count(vs[0]) == 1) c = Counter(vs[1:]) parity = all(mu % 2 == 0 for mu in c.values()) odds = [x for x in el if x % 2] k0 = len(odds) rows = [] for i, x in enumerate(odds): s = 0 for j, y in enumerate(odds): if i == j: continue d = (x + y) % MM s += v2(d) if d else 10**6 rows.append(s) regular = all(z == V for z in rows) verdict("%s: 1+V = nu_2(C), simple summit, parity, regularity" % lab, ok_V and simple_summit and parity and regular, "V = %d, k_0 = %d, profile %s" % (V, k0, vs)) # R = 1 impossible for k_0 = 3 found = None for xs in product(range(8), repeat=3): for ys in product(range(8), repeat=3): M = [] bad_ = False for x in xs: row = [] for y in ys: if x == y: bad_ = True break row.append(v2(x-y)) if bad_: break M.append(row) if bad_: continue if (all(sum(r) == 1 for r in M) and all(sum(M[i][k] for i in range(3)) == 1 for k in range(3))): found = (xs, ys) break if found: break verdict("R = 1 impossible for k_0 = 3 (brute force mod 8)", found is None) verdict("R = 1 absent from the computed set for k_0 = 3", 1 not in R_realisable(3), "realisable R: %s" % sorted(R_realisable(3))[:6]) verdict("R = 2 present for k_0 = 3 (consistent with BLP (b))", 2 in R_realisable(3)) # sieve on the odd profiles def admissible_odd(prof, V): c = Counter(prof) vals = sorted(c, reverse=True) if c[0] == 0: return False if c[vals[0]] != 1: return False for w in vals[1:]: if c[w] % 2: return False asc = sorted(c) for idx, w in enumerate(asc): if sum(c[w2]*(w2-w) for w2 in asc[idx+1:]) < c[w]-1: return False return V in R_odd(c[0]) @lru_cache(maxsize=None) def node_anti(a, b, cap=40): 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_anti(a1, b1, cap): for (x2, y2) in node_anti(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_odd(m0): if m0 < 2: return frozenset() s = {m0-1} for p in range(1, m0): for (x, y) in node_anti(p, m0-p): if x is not None and y is not None and x == y: s.add(m0-1+x) return frozenset(z for z in s if z <= 40) def enum_odd(n, V): out = [] def rec(pos, rest, lo, cur): if pos == n: if rest == 0 and admissible_odd(cur, V): out.append(tuple(sorted(cur, reverse=True))) return for w in range(lo, rest+1): if w*(n-pos) > rest: break rec(pos+1, rest-w, w, cur+[w]) rec(0, V, 0, []) return out CMSVprof = tuple(sorted((0,)*6 + (1, 1, 2, 2, 4), reverse=True)) EXPECTED = {(11, 10): 2, (13, 11): 4, (15, 11): 4} for (n, V), att in sorted(EXPECTED.items()): if remaining() < 20: skip_block("odd sieve n = %d" % n, "budget exhausted") continue L = enum_odd(n, V) verdict("n = %2d, V = %2d: %d admissible profiles" % (n, V, len(L)), len(L) == att, "expected %d%s" % (att, "; CMSV present" if CMSVprof in L else "")) # ============================================================ PART F say() say("-"*76) say("F. CLASSIFICATION OF THE EXTREMAL CONFIGURATIONS") say("-"*76) CLASS_EQUAL = {10: 3, 12: 4, 14: 1, 16: 0, 18: 14, 20: 5} CLASS_FORMS = {10: 0, 12: 3, 14: 1, 16: 1, 18: 6, 20: 2} LISTS = {10: [(0,0,0,1,1), (0,0,0,1,2), (0,0,0,1,3)], 12: [(0,0,0,0,1,1), (0,0,0,0,1,2), (0,0,0,0,1,3), (0,0,0,1,1,2)], 14: [(0,0,0,0,1,1,2)]} for n in [x for x in EVEN_SIZES if x in CLASS_EQUAL]: if remaining() < 20: skip_block("classification n = %d, equal profiles" % n, "budget exhausted") continue m = n//2 target = EVEN_BOUNDS[n] eg = sorted(S for S in profiles(m, PROFILE_CAP) if full_bound(S, S) == target) ok = (len(eg) == CLASS_EQUAL[n]) det = "expected %d" % CLASS_EQUAL[n] if n in LISTS: ok = ok and eg == sorted(LISTS[n]) det += "; list matches" if eg == sorted(LISTS[n]) else "; LIST DIFFERS %s" % eg verdict("n = %2d: %2d profile(s) sigma = tau at the bound" % (n, len(eg)), ok, det) for n in [x for x in EVEN_SIZES if x in CLASS_EQUAL]: if n not in PAIRS: skip_block("classification n = %d, forms" % n, "distinct-sums block skipped") continue m = n//2 target = EVEN_BOUNDS[n] fam = {} for (v, S, T) in PAIRS[n]: if v == target: fam.setdefault(T, []).append(S) # a form: sigma equals sigma' except in its largest entry, which takes every value from a # threshold up to the enumeration ceiling s' + r + 3 (hence every larger value as well) one_param = all(all(S[:-1] == T[:-1] for S in L) and sorted(S[-1] for S in L) == list(range(min(S[-1] for S in L), sum(T) + m + 4)) for T, L in fam.items()) verdict("n = %2d: %d form(s) with distinct profiles" % (n, len(fam)), len(fam) == CLASS_FORMS[n] and one_param, "expected %d; one-parameter forms with sigma[:-1] = sigma'[:-1]: %s" % (CLASS_FORMS[n], "yes" if one_param else "NO")) for T in sorted(fam): say(" sigma' = %-32s largest entry of sigma: every value >= %d" % (list(T), min(S[-1] for S in fam[T]))) # ========================================================== SUMMARY say() say("="*76) say("SUMMARY: %d PASS, %d FAIL, %d SKIPPED (%.1f s)" % (PASS, FAIL, SKIP, time.time()-T0)) say("="*76) say() say("TABLE OF BOUNDS, EVEN SIZES") say(" n Kleiman bound gain status") for n in EVEN_SIZES: b = EVEN_BOUNDS[n] k = kleiman(n) st = "EXACT" if n in EXACT_AT else "lower bound" say(" %2d %2d %2d +%-3d %s" % (n, k, b, b-k, st))