"""Reproduce a testing-versus-proof demonstration (Python standard library).
Source: Lander & Parkin (1966), DOI 10.1090/S0002-9904-1966-11654-3.
Bound 144 is deliberately informed by the published witness.
Observed original run: Python 3.12.14, macOS arm64.
Sequential 1,000,000: 0.148606417 s, no hits.
Random 1,000,000: 2.033453542 s, no hits.
Pair-sum 497,640 lookups + 10,296 builds: 0.082261459 s, one hit.
Times vary; these are single runs, not calibrated benchmarks.
"""
import itertools, math, random, time, sys, platform, json

N = 144
BUDGET = 1_000_000
p = [i**5 for i in range(N+1)]
roots = {p[e]: e for e in range(2, N+1)}
space = math.comb(N+3, 4)
results = []

t0 = time.perf_counter()
seq_hits = []
for seq_count, q in enumerate(itertools.islice(
        itertools.combinations_with_replacement(range(1, N+1), 4),
        BUDGET), 1):
    a,b,c,d = q
    e = roots.get(p[a]+p[b]+p[c]+p[d])
    if e is not None and d < e:
        seq_hits.append((*q,e))
seq_time = time.perf_counter()-t0
results.append(dict(method="sequential", candidates=seq_count,
                    seconds=seq_time, hits=seq_hits, last=q))

rng = random.Random(1966)
t0 = time.perf_counter()
random_hits = []
for _ in range(BUDGET):
    # Uniform 4-subsets of 0..146 map bijectively to sorted
    # quadruples with repetition in 1..144. Draws are with replacement.
    t = sorted(rng.sample(range(N+3),4))
    a,b,c,d = (t[i]+1-i for i in range(4))
    e = roots.get(p[a]+p[b]+p[c]+p[d])
    if e is not None and d < e:
        random_hits.append((a,b,c,d,e))
random_time = time.perf_counter()-t0
results.append(dict(method="random", candidates=BUDGET,
                    seconds=random_time, hits=random_hits, seed=1966))

t0 = time.perf_counter()
pairs = {}
pair_count = 0
for a in range(1,N):
    for b in range(a,N):
        # Keep every pair in case distinct pairs share a sum.
        pairs.setdefault(p[a]+p[b], []).append((a,b))
        pair_count += 1
build_time = time.perf_counter()-t0
lookups = matches = 0
structured_hits = []
for e in range(2,N+1):
    for c in range(1,e):
        for d in range(c,e):
            lookups += 1
            for a,b in pairs.get(p[e]-p[c]-p[d], ()):
                matches += 1
                if b <= c:
                    structured_hits.append((a,b,c,d,e))
structured_time = time.perf_counter()-t0
results.append(dict(method="pair-sum", candidates=lookups,
                    pair_builds=pair_count, raw_pair_matches=matches,
                    seconds=structured_time, build_seconds=build_time,
                    hits=structured_hits))
witness = (27,84,110,133,144)
assert sum(x**5 for x in witness[:4]) == witness[4]**5
assert witness in structured_hits
print(json.dumps(dict(N=N, quadruple_space=space, results=results,
                      powers=[x**5 for x in witness],
                      python=sys.version, platform=platform.platform()), indent=2))
