#!/usr/bin/env python3
"""Exact finite copula certificates, using only Python's standard library.
No simulations certify a limit. General proofs and the two infinite threshold
sequences are in the accompanying pages. --output writes only the named file.
"""
from fractions import Fraction as F
from itertools import product
from collections import Counter
from pathlib import Path
import argparse, json, math
COUNTS=Counter()
def check(ok, label):
    if not ok: raise RuntimeError(label)
    COUNTS[label]+=1

def clamp(x):return max(F(0),min(F(1),x))
def compositions(total,n):
    if n==1:yield (total,);return
    for k in range(total+1):
        for r in compositions(total-k,n-1):yield (k,)+r

def cumulative(p):
    s=F(0);out=[s]
    for a in p:s+=a;out.append(s)
    return out

def jitter_cdf(table,u,v,shared=False):
    rows=[sum(r)for r in table];cols=[sum(c)for c in zip(*table)]
    a=cumulative(rows);b=cumulative(cols);ans=F(0)
    for i,row in enumerate(table):
        for j,p in enumerate(row):
            if p:
                r=clamp((u-a[i])/rows[i]);s=clamp((v-b[j])/cols[j])
                ans+=p*(min(r,s)if shared else r*s)
    return ans

def rectangle(fun,a,b,c,d):return fun(b,d)-fun(a,d)-fun(b,c)+fun(a,c)
def clayton_one(u,v):return F(0)if u==v==0 else u*v/(u+v-u*v)
def independent(u,v):return u*v
def mixture(u,v):return (min(u,v)+max(u+v-1,0))/2

def sklar_checks():
    grid=[F(k,8)for k in range(9)];models=0
    for mass in compositions(4,6):
        table=[[F(mass[i*2+j],4)for j in range(2)]for i in range(3)]
        rows=[sum(r)for r in table];cols=[sum(c)for c in zip(*table)]
        a=cumulative(rows);b=cumulative(cols);models+=1
        for shared in (False,True):
            fun=lambda u,v:jitter_cdf(table,u,v,shared)
            for u in grid:
                check(fun(u,F(1))==u and fun(F(1),u)==u,'randomized_transform_uniform_margins')
            for i in range(4):
                for j in range(3):
                    old=sum((table[r][s]for r in range(i)for s in range(j)),F(0))
                    check(fun(a[i],b[j])==old,'sklar_discrete_grid_recovery')
            for i,j in product(range(8),repeat=2):
                check(rectangle(fun,grid[i],grid[i+1],grid[j],grid[j+1])>=0,'jitter_nonnegative_cell_mass')
        # Every positive-mass atom is recovered for interior jitter, even when other atoms have zero mass.
        for p,cum in [(rows,a),(cols,b)]:
            for i,w in enumerate(p):
                if w:
                    for r in (F(1,7),F(1,2),F(6,7)):
                        u=cum[i]+r*w
                        q=next(j for j in range(len(p)) if cum[j+1]>=u)
                        check(q==i,'generalized_inverse_recovers_atom')
    fair=[[F(1,4)]*2 for _ in range(2)]
    check(jitter_cdf(fair,F(1,4),F(1,4),False)==F(1,16),'fair_independent_jitter_value')
    check(jitter_cdf(fair,F(1,4),F(1,4),True)==F(1,8),'fair_shared_jitter_value')
    check(F(1,12)/4==F(1,48) and F(1,48)/F(1,12)==F(1,4),'shared_jitter_covariance')
    return {'tables':models,'independent_jitter_at_quarter':'1/16','shared_jitter_at_quarter':'1/8','shared_covariance':'1/48','shared_correlation':'1/4'}

# Circle arcs encode failures; conditional interval filling makes actual uniforms.
def circle_model(u):
    lengths=[1-x for x in u];starts=cumulative(lengths)[:-1]
    points={F(0),F(1)}
    for s,l in zip(starts,lengths):points.update((s%1,(s+l)%1))
    points=sorted(points);pieces=[]
    for a,b in zip(points,points[1:]):
        t=(a+b)/2
        failures=tuple((t-(s%1))%1<l for s,l in zip(starts,lengths))
        pieces.append((b-a,failures))
    return pieces

def circle_cdf(u,query,pieces):
    ans=F(0)
    for weight,fail in pieces:
        val=weight
        for cut,q,bad in zip(u,query,fail):
            if bad:
                if cut==1:raise RuntimeError('impossible all-failure branch')
                val*=clamp((q-cut)/(1-cut))
            else:
                if cut==0:raise RuntimeError('impossible no-failure branch')
                val*=clamp(q/cut)
        ans+=val
    return ans

def box_delta(fun,lo,hi):
    out=F(0);d=len(lo)
    for bits in product((0,1),repeat=d):
        pt=tuple(hi[i]if bits[i]else lo[i]for i in range(d))
        out+=(-1)**(d-sum(bits))*fun(pt)
    return out

def frechet_checks():
    models=0
    for d in range(2,6):
        # Full ternary input grids include zero and one, disjoint and wrapped failures.
        for u in product((F(0),F(2,5),F(1)),repeat=d):
            pieces=circle_model(u);models+=1
            check(sum(w for w,b in pieces)==1,'circle_partition_mass')
            check(circle_cdf(u,u,pieces)==max(sum(u)-d+1,0),'circle_pointwise_lower_attainment')
            for i in range(d):
                for v in (F(1,7),F(1,2),F(6,7)):
                    q=[F(1)]*d;q[i]=v
                    check(circle_cdf(u,q,pieces)==v,'circle_filled_uniform_marginal')
    negatives={}
    for d in range(2,9):
        W=lambda u:max(sum(u)-len(u)+1,F(0))
        vol=box_delta(W,[F(1,2)]*d,[F(1)]*d)
        check(vol==1-F(d,2),'lower_function_half_cube_delta')
        if d>=3:check(vol<0,'high_dimensional_lower_rejected')
        negatives[str(d)]=str(vol)
    u=(F(3,5),F(7,10),F(4,5));pieces=circle_model(u)
    check(circle_cdf(u,u,pieces)==F(1,10),'capstone_three_threshold_attainment')
    mid=circle_cdf(u,[F(1,2)]*3,pieces)
    check(mid==F(25,672),'capstone_other_threshold_not_minimal')
    return {'constructed_models':models,'half_cube_masses':negatives,'target_mass':'1/10','same_model_center_mass':str(mid)}

def archimedean_checks():
    grid=[F(k,12)for k in range(13)];cells=0
    for name,fun in [('Clayton theta=1',clayton_one),('independent',independent),('countermonotone',lambda u,v:max(u+v-1,F(0)))]:
        for u in grid:check(fun(u,1)==fun(1,u)==u,'archimedean_margins')
        for i,j in product(range(12),repeat=2):
            check(rectangle(fun,grid[i],grid[i+1],grid[j],grid[j+1])>=0,'archimedean_exact_rectangle');cells+=1
    # Conditional inversion with exact squared rational innovations: sqrt(z)=s is exact.
    for iu,is_ in product(range(1,40),repeat=2):
        u=F(iu,40);s=F(is_,40);z=s*s;v=u*s/(1-(1-u)*s)
        check(0<v<1,'conditional_inverse_in_unit_interval')
        check(v*v/(u+v-u*v)**2==z,'conditional_inverse_equation')
    # Exact finite evaluations of derivative numerator identities, not symbolic proofs.
    for u,v in product((F(k,13)for k in range(1,13)),repeat=2):
        D=u+v-u*v
        check(v*D-u*v*(1-v)==v*v,'clayton_first_derivative_numerator')
        check(2*v*D-2*v*v*(1-u)==2*u*v,'clayton_density_numerator')
    h=lambda t:max(1-t*t,F(0))
    bad=h(0)-2*h(F(1,4))+h(F(1,2))
    check(bad==F(-1,8),'nonconvex_generator_rejected')
    # Finite convex inverse example h=max(1-t,0), verify midpoint and increment structure.
    H=lambda t:max(1-t,F(0))
    for s,a,b in product((F(k,4)for k in range(6)),repeat=3):
        check(H(s)-H(s+a)-H(s+b)+H(s+a+b)>=0,'convex_inverse_fixed_increment')
    check(clayton_one(F(1,2),F(1,4))==F(1,5),'capstone_continuous_joint_cdf')
    c=F(1,5);table=[c,F(1,2)-c,F(1,4)-c,1-F(1,2)-F(1,4)+c]
    check(table==[F(1,5),F(3,10),F(1,20),F(9,20)],'capstone_continuous_four_cells')
    check(F(1,2)*F(1,2)/(1-F(1,2)*F(1,2))==F(1,3),'capstone_conditional_sample')
    return {'rectangles':cells,'nonconvex_rectangle':str(bad),'four_cells':list(map(str,table)),'conditional_sample_v':'1/3'}

# Rational polygon clipping: independent geometric witness for modulo-one triple events.
def clip(poly,A,B,C):
    ans=[]
    if not poly:return ans
    for P,Q in zip(poly,poly[1:]+poly[:1]):
        p=A*P[0]+B*P[1]-C;q=A*Q[0]+B*Q[1]-C
        if p<=0:ans.append(P)
        if (p<0<q)or(q<0<p):
            t=p/(p-q);ans.append((P[0]+t*(Q[0]-P[0]),P[1]+t*(Q[1]-P[1])))
    return ans

def area(poly):
    if len(poly)<3:return F(0)
    return abs(sum((p[0]*q[1]-p[1]*q[0]for p,q in zip(poly,poly[1:]+poly[:1])),F(0)))/2

def modulo_cdf(u,v,w):
    rect=[(F(0),F(0)),(u,F(0)),(u,v),(F(0),v)];out=F(0)
    for k in (0,1):
        poly=clip(rect,-1,-1,-k);poly=clip(poly,1,1,k+w);out+=area(poly)
    return out

def tail_checks():
    fair=[[F(1,4)]*2 for _ in range(2)]
    for n in range(3,90):
        u=F(1,n)
        check(mixture(u,u)/u==F(1,2),'mixture_lower_tail_ratio')
        v=1-u
        check((1-2*v+mixture(v,v))/(1-v)==F(1,2),'mixture_upper_tail_ratio')
        check(jitter_cdf(fair,u,u,True)/u==F(1,2),'shared_jitter_lower_tail_ratio')
        check((1-2*v+jitter_cdf(fair,v,v,True))/(1-v)==F(1,2),'shared_jitter_upper_tail_ratio')
        check(clayton_one(u,u)/u==1/(2-u),'clayton_lower_diagonal_ratio')
        check((1-2*v+clayton_one(v,v))/(1-v)==2*(1-v)/(2-v),'clayton_upper_diagonal_ratio')
    check((1-2*F(9,10)+mixture(F(9,10),F(9,10)))==F(1,20),'zero_correlation_joint_upper_tail')
    check(F(1,2)*F(1,3)+F(1,2)*F(1,6)==F(1,4),'zero_correlation_moment')
    # Exact infinite-construction subsequences evaluated only at their symbolic locations.
    seq=[]
    for k in range(0,10,2):
        a=F(1,2**(2**k));b=a*a;mid=(a+b)/2
        # Current reflected interval contributes length max(0,2*mid-a-b)=0;
        # all lower intervals together have mass b, by the preserved-block proof.
        cmid=b+max(F(0),2*mid-a-b)
        check(cmid/mid==2*a/(1+a),'nonexistent_tail_reflected_midpoint_ratio')
        check(a/a==1,'nonexistent_tail_endpoint_ratio')
        seq.append({'k':k,'endpoint_ratio':'1','midpoint_ratio':str(cmid/mid)})
    grid=[F(k,10)for k in range(11)]
    for u,v in product(grid,repeat=2):
        check(modulo_cdf(u,v,F(1))==u*v,'modulo_pair_UV_independent')
        check(modulo_cdf(u,F(1),v)==u*v,'modulo_pair_UW_independent')
        check(modulo_cdf(F(1),u,v)==u*v,'modulo_pair_VW_independent')
    for n in range(2,65):
        e=F(1,n);val=modulo_cdf(e,e,e)
        check(val==e*e/2,'modulo_triple_triangle_mass')
        check(val/(e*e)==F(1,2),'modulo_triple_conditional_probability')
    # A complete coarse 3D CDF cell partition independently confirms nonnegative masses summing to one.
    grid=[F(k,4)for k in range(5)];total=F(0)
    for idx in product(range(4),repeat=3):
        lo=[grid[k]for k in idx];hi=[grid[k+1]for k in idx]
        vol=box_delta(lambda x:modulo_cdf(*x),lo,hi)
        check(vol>=0,'modulo_nonnegative_3D_cells');total+=vol
    check(total==1,'modulo_3D_mass_one')
    return {'reflection_subsequences':seq,'pairwise_independent_triple_at_one_tenth':'1/200','fully_independent_triple_at_one_tenth':'1/1000','triple_conditional':'1/2'}

def run():
    COUNTS.clear()
    a=sklar_checks();b=frechet_checks();c=archimedean_checks();d=tail_checks()
    return {'status':'PASS','checks':sum(COUNTS.values()),'categories':dict(sorted(COUNTS.items())),'sklar':a,'frechet':b,'archimedean':c,'tail':d,'scope':'Exact finite algebra, finite-support jitter integration, circle partitions and polygon clipping. Grid checks do not prove universal convexity or existence of a tail limit; those are supplied by the accompanying proofs.'}
if __name__=='__main__':
    ap=argparse.ArgumentParser();ap.add_argument('--output');args=ap.parse_args();result=json.dumps(run(),ensure_ascii=False,indent=2)+'\n'
    if args.output:Path(args.output).write_text(result)
    else:print(result,end='')