#!/usr/bin/env python3
"""Exact checks for the explicit branch and conformal-map certificates.

Standard library only. All checks remain enabled under python -O.
Finite rational samples test identities and examples, not general analytic theorems.
"""
from dataclasses import dataclass
from fractions import Fraction as F
from itertools import product, permutations
from pathlib import Path
import argparse
import json

@dataclass(frozen=True)
class C:
    r: F = F(0)
    i: F = F(0)
    def __post_init__(self):
        object.__setattr__(self,'r',F(self.r));object.__setattr__(self,'i',F(self.i))
    def __add__(self,z):
        z=cast(z);return C(self.r+z.r,self.i+z.i)
    __radd__=__add__
    def __neg__(self):return C(-self.r,-self.i)
    def __sub__(self,z):return self+-cast(z)
    def __rsub__(self,z):return cast(z)+-self
    def __mul__(self,z):
        z=cast(z);return C(self.r*z.r-self.i*z.i,self.r*z.i+self.i*z.r)
    __rmul__=__mul__
    def conjugate(self):return C(self.r,-self.i)
    def norm2(self):return self.r*self.r+self.i*self.i
    def __truediv__(self,z):
        z=cast(z);v=z.norm2()
        if not v:raise ZeroDivisionError('complex denominator is zero')
        n=self*z.conjugate();return C(n.r/v,n.i/v)
    def __rtruediv__(self,z):return cast(z)/self
    def __pow__(self,n):
        if not isinstance(n,int):raise TypeError('integer exponent required')
        if n<0:return (1/self)**(-n)
        out=C(1)
        for _ in range(n):out=out*self
        return out
    def __eq__(self,z):
        try:z=cast(z)
        except (TypeError,ValueError):return False
        return self.r==z.r and self.i==z.i
    def record(self):return {'real':str(self.r),'imag':str(self.i)}
def cast(z):return z if isinstance(z,C) else C(z)
ZERO=C();ONE=C(1);I=C(0,1)

@dataclass(frozen=True)
class Mobius:
    a:C;b:C;c:C;d:C
    def __post_init__(self):
        for k in ['a','b','c','d']:object.__setattr__(self,k,cast(getattr(self,k)))
        if self.det()==ZERO:raise ValueError('nonzero determinant required')
    def det(self):return self.a*self.d-self.b*self.c
    def __call__(self,z):
        # None denotes the single point at infinity.
        if z is None:return None if self.c==ZERO else self.a/self.c
        z=cast(z);den=self.c*z+self.d
        return None if den==ZERO else (self.a*z+self.b)/den
    def inverse(self):return Mobius(self.d,-self.b,-self.c,self.a)
    def compose(self,t):
        # self after t: its matrix is M_self M_t, in this order.
        return Mobius(self.a*t.a+self.b*t.c,self.a*t.b+self.b*t.d,
                      self.c*t.a+self.d*t.c,self.c*t.b+self.d*t.d)
    def derivative(self,z):return self.det()/(self.c*z+self.d)**2

def disk_shift(a):
    a=cast(a)
    if a.norm2()>=1:raise ValueError('center must be inside the disk')
    return Mobius(1,-a,-a.conjugate(),1)
def rho2(z,w):return ((z-w)/(1-w.conjugate()*z)).norm2()
counts={}
def check(group,condition):
    counts[group]=counts.get(group,0)+1
    if not condition:raise RuntimeError('failed '+group)
def expect_error(group,thunk):
    try:thunk()
    except (ValueError,ZeroDivisionError):check(group,True);return
    check(group,False)

def run():
    counts.clear()
    # Exact polynomial roots and integer periods for the specified factorization.
    for r,i in product(range(-3,4),repeat=2):
        z=C(F(r,2),F(i,2))
        if z in [C(1),C(-1)]:continue
        f=(z-1)**2*(z+1)**4;h=(z-1)*(z+1)**2
        check('factorized_square_root',h*h==f and (-h)**2==f)
        hp=(z+1)**2+2*(z-1)*(z+1)
        check('logarithmic_derivative',2*hp/h==2/(z-1)+4/(z+1))
    for wp,wm in product(range(-5,6),repeat=2):
        period=2*wp+4*wm
        check('integer_period_samples',period%2==0 and -period==2*(-wp)+4*(-wm))
    check('single_loop_obstructions',2!=0 and 2%3!=0 and 4%3!=0)
    check('changed_exponent_obstruction',3%2!=0)
    check('basepoint_branch',(C(2)-1)*(C(2)+1)**2==C(9))

    mats=[Mobius(1,0,0,1),Mobius(1,I,0,1),Mobius(0,1,1,0),
          Mobius(1,-I,1,I),Mobius(2,1,I,1),disk_shift(C(F(1,3))),
          disk_shift(C(F(1,4),F(-1,3)))]
    points=[None,ZERO,ONE,-ONE,I,-I,C(2,3),C(F(1,2),F(-1,4))]
    for m,z in product(mats,points):
        check('sphere_inverse',m.inverse()(m(z))==z)
        scaled=Mobius(I*m.a,I*m.b,I*m.c,I*m.d)
        check('scalar_matrix_equivalence',scaled(z)==m(z))
    for s,t,z in product(mats,mats,points):
        check('matrix_composition_order',s.compose(t)(z)==s(t(z)))
    for a,b,c in permutations([ZERO,ONE,I,-ONE,None],3):
        if c is None:m=Mobius(1,-a,0,b-a)
        elif b is None:m=Mobius(1,-a,1,-c)
        elif a is None:m=Mobius(0,b-c,1,-c)
        else:m=Mobius(b-c,-a*(b-c),b-a,-c*(b-a))
        check('three_point_calibration',m(a)==ZERO and m(b)==ONE and m(c) is None)
    expect_error('invalid_mobius',lambda:Mobius(1,1,1,1))
    expect_error('boundary_disk_center',lambda:disk_shift(ONE))

    cayley=Mobius(1,-I,1,I)
    quadrant=Mobius(1,-2*I,1,2*I)
    right=Mobius(1,-1,1,1)
    for u in range(-6,7):
      for v in range(1,7):
        z=C(F(u,3),F(v,3));w=cayley(z)
        check('upper_half_plane_identity',1-w.norm2()==4*z.i/(z+I).norm2())
        check('upper_half_plane_inverse',cayley.inverse()(w)==z and w.norm2()<1)
        z=C(F(v,3),F(u,3));w=right(z)
        check('right_half_plane_identity',1-w.norm2()==4*z.r/(z+1).norm2())
        check('slit_square_inverse',right.inverse()(w)**2==z*z and right.inverse()(w).r>0)
    for u,v in product(range(1,9),repeat=2):
        z=C(F(u,4),F(v,4));square=z*z;w=quadrant(square)
        check('quadrant_target',square.i>0 and w.norm2()<1)
        check('quadrant_rational_inverse',quadrant.inverse()(w)==square)
        check('quadrant_disk_identity',1-w.norm2()==8*square.i/(square+2*I).norm2())
    a=C(1,1);deriv=quadrant.derivative(a*a)*2*a
    check('quadrant_basepoint',quadrant(a*a)==ZERO)
    check('normalized_derivative',deriv==C(F(1,2),F(-1,2)) and C(1,1)*deriv==ONE)
    check('conformal_radius_squared',1/deriv.norm2()==2)
    # Multiplication by (1+i)/sqrt(2) has norm one; the exact product above is 1.
    check('slit_base_derivative',right.derivative(ONE)/2==C(F(1,4)))
    check('transport_second_input',quadrant(6*I)==C(F(1,2)))

    for q in [F(k,3) for k in range(-12,13)]:
        u=C((1-q*q)/(1+q*q),2*q/(1+q*q));z=1+u
        check('rational_unit_circle',u.norm2()==1)
        if z!=ZERO:check('circle_through_pole_to_line',(1/z).r==F(1,2))

    grid=[C(F(r,3),F(i,3)) for r,i in product(range(-2,3),repeat=2) if r*r+i*i<9]
    centers=[ZERO,C(F(1,3)),C(F(1,4),F(1,4))]
    scales=[ZERO,C(F(1,2)),C(F(2,3),F(1,3)),ONE,I]
    for a,b,c in product(centers,centers,scales):
        f=disk_shift(b).inverse().compose(Mobius(c,0,0,1)).compose(disk_shift(a)) if c!=ZERO else None
        value=(lambda z:f(z)) if f else (lambda z:b)
        for z in grid:
            fz=value(z);df=f.derivative(z) if f else ZERO
            check('constructed_disk_target',fz.norm2()<1)
            left=df.norm2()*(1-z.norm2())**2;rightbound=(1-fz.norm2())**2
            check('derivative_pick_bound',left<=rightbound)
            check('derivative_rigidity', (left==rightbound)==(c.norm2()==1))
            for w in grid:
                lhs=rho2(fz,value(w));rhs=rho2(z,w)
                check('two_point_pick_bound',lhs<=rhs)
                if z!=w:check('two_point_rigidity',(lhs==rhs)==(c.norm2()==1))
    for z in grid:
        check('square_map_strict_derivative',4*z.norm2()*(1-z.norm2())**2<(1-z.norm2()**2)**2)

    cases=[];b0=C(F(1,3));shift=disk_shift(b0)
    for b,expected in [(F(2,3),'feasible'),(F(5,7),'unique_automorphism'),(F(3,4),'impossible')]:
        beta=shift(C(b));val=beta.norm2();out={'b':str(b),'normalized_value':str(beta.r),'classification':expected}
        actual='impossible' if val>F(1,4) else 'unique_automorphism' if val==F(1,4) else 'feasible'
        check('capstone_interpolation',actual==expected)
        check('feasible_output_disk',(b-F(9,35))**2<=F(16,35)**2 if actual!='impossible' else (b-F(9,35))**2>F(16,35)**2)
        if actual!='impossible':
            c=2*beta;f=shift.inverse().compose(Mobius(c,0,0,1))
            check('capstone_interpolation_values',f(ZERO)==b0 and f(C(F(1,2)))==C(b))
            check('transport_derivative_squared',(f.derivative(ZERO)*deriv).norm2()==f.derivative(ZERO).norm2()/2)
            out.update(scale=str(c.r),disk_derivative=str(f.derivative(ZERO).r),quadrant_derivative_norm_squared=str((f.derivative(ZERO)*deriv).norm2()))
        cases.append(out)
    for n in range(1,40):
        s=F(n,40)
        check('missing_value_strict_gain',(1+s)**2/(4*s)>1)
        check('gain_excess_identity',(1+s)**2/(4*s)-1==(1-s)**2/(4*s))
    return {'status':'PASS','checks':sum(counts.values()),'groups':dict(counts),
            'integer_periods':{'around_plus_one':2,'around_minus_one':4,'logarithm':False,'root_orders':[1,2],'positive_square_root_at_2':9},
            'explicit_maps':{'quadrant_basepoint':[1,1],'unnormalized_derivative':deriv.record(),'normalized_derivative':'1/sqrt(2)','quadrant_conformal_radius':'sqrt(2)','slit_conformal_radius':'4'},
            'interpolation':cases,'arithmetic':'exact Fraction pairs; no float comparison',
            'limitations':['The loop integer formula and global branch classification are proved in prose, not inferred from sampled paths.','Inverse checks concern explicit maps and their stated branches; no general domain mapping algorithm is claimed.','Finite Pick checks do not prove the theorem for all holomorphic maps.','Infinity is one sphere point, not a signed or directional infinity.']}

def main():
    ap=argparse.ArgumentParser();ap.add_argument('--output',type=Path);args=ap.parse_args()
    result=run();s=json.dumps(result,ensure_ascii=False,indent=2)+'\n'
    if args.output:args.output.parent.mkdir(parents=True,exist_ok=True);args.output.write_text(s)
    print(s,end='')
if __name__=='__main__':main()