"""Exact-arithmetic sufficient criterion following the self-reviewed N=52 draft.

Geometric proof: docs/P56-ROW-CRITERION-2026-09-08.md.
An unresolved case is not a counterexample and is never accepted.
"""
from fractions import Fraction as F
from math import isqrt
from pathlib import Path
import argparse
import json

PL = F(3141592653589793238462643, 10**24)
PH = PL + F(1, 10**24)


def root_bounds(q):
    scale = 10**24
    k = isqrt(q.numerator * scale**2 // q.denominator)
    lo, hi = F(k, scale), F(k + 1, scale)
    assert lo**2 <= q <= hi**2
    return lo, hi


R3L, R3H = root_bounds(F(3))
_, R5H = root_bounds(F(5))


def ceil_fraction(q):
    return -(-q.numerator // q.denominator)


def sin_lower(x):
    assert 0 <= x <= 1
    return x - x**3 / 6


def cos_lower(x):
    assert 0 <= x <= 2
    return 1 - x**2 / 2 + x**4 / 24 - x**6 / 720


def cos_upper(x):
    assert 0 <= x <= 2
    return 1 - x**2 / 2 + x**4 / 24


def certify_pi():
    def atan_bounds(x):
        total, k = F(0), 0
        while True:
            total += (-1)**k * x**(2*k+1) / (2*k+1)
            k += 1
            term = x**(2*k+1)/(2*k+1)
            if term < F(1, 10**32):
                return (total, total+term) if k % 2 == 0 else (total-term, total)
    l5, h5 = atan_bounds(F(1,5))
    l239, h239 = atan_bounds(F(1,239))
    # Machin's identity and alternating arctangent series.
    assert PL < 16*l5-4*h239 < 16*h5-4*l239 < PH


def angle_checks(E, u, q):
    beta = 3*E/2
    A = PH/30+12*beta
    return dict(wall_cone=(5-R5H)/8-q/4,
        five_angles=2*PL-5*(PH/3+2*beta),
        seven_angles=7*(PL/3-beta)-2*PH,
        acute=PL/2-(PH/3+2*beta),
        wall_fan_one=3*PL/5-(PH/3+2*beta),
        wall_fan_five=5*(PL/3-beta)-7*PH/5,
        corner_fan_three=3*(PL/3-beta)-9*PH/10,
        mixed_modes=sum(sin_lower(2*PL/15-c*E) for c in [12,21])-u/2,
        same_sign_large=2*PL/3-16*beta-2*PH/5,
        adjacent_straight=sum(sin_lower(PL/12-c*E) for c in [3,12,21])-u/2,
        adjacent_zigzag=PL/10-20*beta,
        corner_color=PL/10-12*beta,
        corner_direction=2*PL/15-A,
        endpoint_monotonicity=1-q*A*A)


def representations(n):
    out = []
    for k in range(2, n // 2 + 1):
        if n % k == 0:
            m = n // k
            out.append(dict(family='A', k=k, m=m, counts=[m]*k, s=F(k-1, 2*m-1), d=F(2, 2*m-1)))
        for phase in [0, 1]:
            low_count = (k + phase) // 2
            if (n + low_count) % k:
                continue
            m = (n + low_count) // k
            # Only rows of at least two sites are relevant to the theorem.
            if m < 3:
                continue
            counts = [m - ((i+phase) % 2) for i in range(k)]
            assert sum(counts) == n
            out.append(dict(family=f'B{phase}', k=k, m=m, counts=counts, s=F(k-1, 2*m-2), d=F(1, m-1)))
    return out


def inspect(n):
    reps = representations(n)
    # Equality constructions are already handled by the universal lower bound.
    if any(3*r['s']**2 <= 1 and 3*(r['s'] + r['d']/3)**2 >= 1 for r in reps):
        return dict(n=n, status='universal-bound equality construction')
    compressed = [r for r in reps if 3*r['s']**2 > 1]
    if not compressed:
        return dict(n=n, status='no compressed candidate')
    best = min(compressed, key=lambda r: r['s'])
    S = best['s']
    q = 1 + S**2
    _, u = root_bounds(q)
    E = F(ceil_fraction((u - 2/R3H)*10**9), 10**9)
    if E > F(5481, 10**6):
        return dict(n=n, status='outside the tested local-angle range', candidateSquared=str(q))

    beta = 3*E/2
    kappa = 6*beta
    A = PH/30 + 12*beta
    alpha = 2*PL/15 - 8*beta
    B = u*(R3H/2 + beta)
    dpack, _ = root_bounds(F(n)*PL)
    _, dcover = root_bounds(F(n)*PH)
    D = 1/(dpack/2 - 1)  # strict delta upper bound
    d0 = 2/(u*dcover)     # strict delta lower bound
    wall_count = (1/(u*D)).numerator // (1/(u*D)).denominator + 1

    checks = dict(angle_checks(E, u, q),
        epsilon_cap=E-(u-2/R3H), h_lt_quarter=F(1,4)-u*D/2,
        wall_sites_at_least_four=F(wall_count-3))
    failed = [name for name, margin in checks.items() if margin <= 0]
    result = dict(n=n, candidateSquared=str(q), candidateApprox=float(q)**.5,
        candidate=dict(K=(2*best['m']-1 if best['family']=='A' else 2*best['m']-2),
            L=best['k']-1, phase=(1 if best['family']=='B1' else 0), rows=best['counts']),
        epsilonCap=str(E), deltaLower=str(d0), deltaUpper=str(D), wallSitesMinimum=wall_count,
        constants=dict(A=str(A), beta=str(beta), B=str(B), S=str(S)),
        localChecks={name: str(value) for name,value in checks.items()})
    if failed:
        return dict(result, status='local structural test unresolved', failed=failed)

    # A row cannot have t steps: before its t-th step it must still point right,
    # but these t steps would already have horizontal extent >1.
    total = F(0)
    barrier = None
    for t in range(1, n):
        theta = A+(t-1)*kappa
        if theta >= PL/2:
            break
        total += cos_lower(theta)
        if total*d0 > 1:
            barrier = dict(forbiddenSteps=t, angleUpper=str(theta),
                projectionLower=str(total*d0), projectionMarginApprox=float(total*d0-1))
            break
    if barrier:
        max_steps = barrier['forbiddenSteps']-1
        result['rowBarrier'] = barrier
    elif A+(n-2)*kappa < PL/2:
        max_steps = n-1
        result['rowBarrier'] = dict(finiteCardinality=True, maxSteps=max_steps)
    else:
        return dict(result, status='row projection test unresolved')

    outcomes = []
    for r in reps:
        counts, k = r['counts'], r['k']
        if k < wall_count or counts[0] < wall_count or counts[-1] < wall_count:
            continue
        steps = max(counts)-1
        aq = 1 + sum(F(1, counts[i]+counts[i+1]-1) for i in range(k-1))**2
        row = dict(family=r['family'], rows=counts, areaBoundSquared=str(aq))
        if steps > max_steps:
            row['result'] = 'excluded by row projection barrier'
        elif aq >= q:
            row['result'] = 'area bound reaches or exceeds candidate'
        else:
            theta = A+(steps-1)*kappa
            assert theta < PL/2
            # Both ends are controlled after the left-to-right structure is
            # proved. The i-th step is close to its nearer end, not necessarily
            # to the left end. Sum individual projection bounds.
            horizontal = sum(cos_lower(A+min(i, steps-1-i)*kappa) for i in range(steps))
            vertical = S+A+(k-1)*B*cos_upper(alpha)
            gap = horizontal-vertical
            row.update(exclusionMargin=str(gap), exclusionMarginApprox=float(gap))
            if gap > 0:
                row['result'] = 'excluded by incompatible spans'
            elif horizontal*d0 > 1:
                row['result'] = 'excluded by horizontal extent and covering area'
            elif vertical*D < 1:
                row['result'] = 'excluded by vertical extent and packing area'
            else:
                row['result'] = 'unresolved'
        outcomes.append(row)
    assert outcomes, 'A feasible supplied candidate cannot disappear from the row list.'
    result['representations'] = outcomes
    result['status'] = ('criterion passes (self-reviewed proof, not independent review)'
        if all(r['result'] != 'unresolved' for r in outcomes) else 'representation exclusion unresolved')
    return result


def run():
    parser = argparse.ArgumentParser()
    parser.add_argument('--max-n', type=int, default=400)
    args = parser.parse_args()
    assert 19 <= args.max_n <= 10000
    certify_pi()
    uniform = angle_checks(F(5481,10**6), F(1160182,10**6), F(1160182,10**6)**2)
    assert all(x > 0 for x in uniform.values())
    # An independent, deliberately simple count enumeration for regression.
    enumeration_test_max = min(args.max_n, 400)
    for n in range(19,enumeration_test_max+1):
        brute = set()
        for k in range(2,n//2+1):
            for m in range(2,n//k+2):
                patterns = [(m,)*k] + [tuple(m-((i+p)%2) for i in range(k)) for p in [0,1]]
                brute.update(p for p in patterns if min(p)>=2 and sum(p)==n)
        assert brute == {tuple(r['counts']) for r in representations(n)}, n
    results = [inspect(n) for n in range(19, args.max_n+1)]
    passed = [r for r in results if r['status'].startswith('criterion passes')]
    detailed = [r for r in results if 'localChecks' in r]
    counts = {s:sum(r['status']==s for r in results) for s in sorted(set(r['status'] for r in results))}
    output = dict(scope='Sufficient criterion for continuous optima; local research only. Failures are inconclusive.',
        selfTests=dict(piEnclosure='Machin identity + rational alternating series passed',
            integerEnumeration=f'Independent brute comparison at every N=19..{enumeration_test_max} passed'),
        uniformAngleChecks={name:str(value) for name,value in uniform.items()},
        maxN=args.max_n, statusCounts=counts, passedNs=[r['n'] for r in passed], detailed=detailed,
        allScreeningResults=[{k:r[k] for k in ['n','status','candidateSquared'] if k in r} for r in results])
    path = Path(__file__).resolve().parents[1]/f'docs/p56-row-criterion-to-{args.max_n}-2026-09-08.json'
    path.write_text(json.dumps(output, indent=2)+'\n', encoding='utf-8', newline='\n')
    print(json.dumps(dict(statusCounts=counts, passed=[dict(n=r['n'], M2=r['candidateSquared'],
        K=r['candidate']['K'], L=r['candidate']['L'], rowBarrier=r['rowBarrier']) for r in passed],
        unresolved=[dict(n=r['n'], status=r['status']) for r in detailed if not r['status'].startswith('criterion passes')]), indent=2))


if __name__ == '__main__':
    run()