"""D1258 C08.2. Exact Zorn Frobenius constraints and signed channel subsets."""

from itertools import product, combinations, permutations
from collections import Counter
import json
import platform
import sympy as S
if not __debug__:
    raise SystemExit("Do not use -O: assertions are required.")
def report(key, value):
    print(key + " = " + str(value))
def va(x,y): return tuple(a+b for a,b in zip(x,y))
def vn(x): return tuple(-a for a in x)
def vs(x,y): return va(x,vn(y))
def sc(k,x): return tuple(k*a for a in x)
def dot(x,y): return sum(a*b for a,b in zip(x,y))
def cross(x,y):
    return (x[1]*y[2]-x[2]*y[1],x[2]*y[0]-x[0]*y[2],x[0]*y[1]-x[1]*y[0])
def qm(p,r):
    a,b,c,d=p; e,f,g,h=r
    return (a*e-b*f-c*g-d*h,a*f+b*e+c*h-d*g,
            a*g-b*h+c*e+d*f,a*h+b*g-c*f+d*e)
def qb(p): return (p[0],-p[1],-p[2],-p[3])
O_NAMES=("1","e1","e2","e3","f1","f2","f3","l")
OB=[tuple(S.Integer(i==j) for i in range(8)) for j in range(8)]
OZ=(S.Integer(0),)*8
def obar(x): return (x[0],)+vn(x[1:])
def omul(x,y,epsilon=1):
    # Same basis as D1256: f_i=-e_i*l, not +e_i*l.
    p=x[:4]; q=(x[7],)+vn(x[4:7])
    r=y[:4]; s=(y[7],)+vn(y[4:7])
    first=va(qm(p,r),sc(epsilon,qm(qb(s),q)))
    second=va(qm(s,p),qm(q,qb(r)))
    return first+vn(second[1:])+(second[0],)
def onorm(x,epsilon=1):
    return sum(a*a for a in x[:4])-epsilon*sum(a*a for a in x[4:])
def to_zorn(x):
    return (x[0]+x[7],)+va(x[1:4],x[4:7])+vs(x[4:7],x[1:4])+(x[0]-x[7],)
def from_zorn(z):
    a=z[0]; u=z[1:4]; v=z[4:7]; b=z[7]
    return ((a+b)/2,)+sc(S.Rational(1,2),vs(u,v))+sc(S.Rational(1,2),va(u,v))+((a-b)/2,)
def zmul(z,w):
    a=z[0]; u=z[1:4]; v=z[4:7]; b=z[7]
    c=w[0]; U=w[1:4]; V=w[4:7]; d=w[7]
    return (a*c+dot(u,V),)+va(va(sc(a,U),sc(d,u)),cross(v,V))+vs(va(sc(c,v),sc(b,V)),cross(u,U))+(dot(v,U)+b*d,)
def zero_vector(v): return all(S.expand(a)==0 for a in v)
def inertia(M):
    """Exact rational symmetric congruence elimination; no eigenvalue tolerance."""
    A=S.Matrix(M); assert A==A.T
    positive=negative=null=0
    while A.rows:
        k=next((i for i in range(A.rows) if A[i,i]!=0),None)
        if k is None:
            pair=next(((i,j) for i in range(A.rows) for j in range(i+1,A.rows) if A[i,j]!=0),None)
            if pair is None:
                null+=A.rows; break
            i,j=pair
            P=S.eye(A.rows); P[j,i]=1
            A=P.T*A*P
            k=i
        inds=[k]+[i for i in range(A.rows) if i!=k]
        A=A.extract(inds,inds); d=A[0,0]
        assert d.is_positive or d.is_negative
        positive+=int(bool(d>0)); negative+=int(bool(d<0))
        v=A[1:,0]; A=A[1:,1:]-(v*v.T)/d
    return (positive,negative,null)

report("environment", {"python":platform.python_version(),"sympy":S.__version__})
report("data_basis","Native Zorn product from declared formula; no banked arrays.")
a,b=S.symbols("a b",real=True)
s=S.Matrix([1/S.sqrt(2),0,0,0,0,0,0,-1/S.sqrt(2)])
u=[S.eye(8)[:,i] for i in range(1,4)]
v=[S.eye(8)[:,i] for i in range(4,7)]
basis=[s]+u+v
C=S.Matrix.hstack(*basis)
ci=(C.T*C).inv()*C.T
def br(x,y):return S.Matrix(vs(zmul(tuple(x),tuple(y)),zmul(tuple(y),tuple(x))))
def co(z):
    q=S.simplify(ci*z); assert S.simplify(C*q-z)==S.zeros(8,1);return q
B=S.zeros(7);B[0,0]=b
for i in range(3): B[i+1,i+4]=B[i+4,i+1]=a
def pairing(x,y):return (x.T*B*y)[0]
std=[S.eye(7)[:,i] for i in range(7)]
bt=[[co(br(x,y)) for y in basis] for x in basis]
defects=[]
for i,j,k in product(range(7),repeat=3):
    d=S.simplify(pairing(bt[i][j],std[k])-pairing(std[i],bt[j][k]))
    if d!=0:defects.append((i,j,k,d))
constraint=S.groebner([t[3] for t in defects],a,b,extension=S.sqrt(2))
report("Frobenius.all_basis_triples",7**3)
report("Frobenius.nonzero_symbolic_defects",len(defects))
report("Frobenius.distinct_nonzero_defects",sorted(set(t[3] for t in defects),key=str))
report("Frobenius.constraint_ideal",list(constraint))
assert list(constraint)==[a-b]
lhs=pairing(bt[0][1],std[4]);rhs=pairing(std[0],bt[1][4])
report("Frobenius.witness_left_right",(lhs,rhs))
assert S.simplify(lhs-S.sqrt(2)*a)==0 and S.simplify(rhs-S.sqrt(2)*b)==0
trtotal=6*a+b;trtriplet=3*a;r=S.cancel(trtotal/trtriplet)
report("finite_trace.ratio",r)
report("finite_trace.normalized_ratio",S.cancel(r.subs(b,a)))
report("finite_trace.zero_weight","a=b=0 satisfies closure but the ratio is undefined.")
report("negative_control.a1_b2.witness",S.simplify((lhs-rhs).subs({a:1,b:2})))
# Actual algebraic imaginary metric, not an assumed all-positive count.
g=[onorm(e) for e in OB[1:]]
counts=Counter()
for select in product([0,1],repeat=7):
    counts[sum(q*t for q,t in zip(select,g))]+=1
report("imaginary.metric_diagonal",g)
report("imaginary.metric_inertia(+,-,0)",inertia(S.diag(*g)))
report("signed_subsets.count",sum(counts.values()))
report("signed_subsets.histogram",sorted(counts.items()))
report("signed_subsets.range",(min(counts),max(counts)))
report("full.signed_contraction",sum(g))
report("full.unweighted_multiplicity",len(g))
# Distinguish raw matrix trace from basis-invariant signature.
P=S.diag(2,1,1,1,1,1,1);G=S.diag(*g);Gp=P.T*G*P
report("nonorthonormal_control.raw_metric_matrix_trace",S.trace(Gp))
report("nonorthonormal_control.inertia(+,-,0)",inertia(Gp))
# Departing from the unit-weight hypothesis can reach seven.
report("outside_unit_weight_class.example",7*g[0])
report("physical_LSZ_residue","NOT COMPUTED; no spacetime action or LSZ limit supplied.")
report("registered_array_alignment","NOT VERIFIED; the declared Zorn object is tested directly.")