"""D1258 C05.1. Exact chosen-vacuum identities; no physical readout."""

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)

def jmat(A):
    a,b,c=A[:3]; z=A[3:11]; y=A[11:19]; x=A[19:27]
    return [[sc(a,OB[0]),z,obar(y)],[obar(z),sc(b,OB[0]),x],[y,obar(x),sc(c,OB[0])]]
def jvec(M):
    assert all(zero_vector(M[i][i][1:]) for i in range(3))
    assert all(zero_vector(vs(M[i][j],obar(M[j][i]))) for i in range(3) for j in range(i+1,3))
    return tuple(M[i][i][0] for i in range(3))+tuple(M[0][1])+tuple(M[2][0])+tuple(M[1][2])
def mm(A,B):
    return [[va(va(omul(A[i][0],B[0][j]),omul(A[i][1],B[1][j])),omul(A[i][2],B[2][j]))
             for j in range(3)] for i in range(3)]
def jp(A,B):
    X=jmat(A); Y=jmat(B); XY=mm(X,Y); YX=mm(Y,X)
    return jvec([[sc(S.Rational(1,2),va(XY[i][j],YX[i][j])) for j in range(3)] for i in range(3)])
JB=[tuple(S.Integer(i==j) for i in range(27)) for j in range(27)]
JI=va(va(JB[0],JB[1]),JB[2]); JZ=(S.Integer(0),)*27
def tr(A): return sum(A[:3])
def sig(A):
    a,b,c=A[:3];z=A[3:11];y=A[11:19];x=A[19:27]
    return a*b+a*c+b*c-onorm(x)-onorm(y)-onorm(z)
def jnorm(A):
    a,b,c=A[:3];z=A[3:11];y=A[11:19];x=A[19:27]
    return a*b*c-a*onorm(x)-b*onorm(y)-c*onorm(z)+2*omul(omul(z,x),y)[0]
def sharp(A): return va(vs(jp(A,A),sc(tr(A),A)),sc(sig(A),JI))
def L(A): return S.Matrix.hstack(*[S.Matrix(jp(A,e)) for e in JB])
def brief(A): return {i:S.simplify(a) for i,a in enumerate(A) if S.simplify(a)!=0}

report("environment", {"python":platform.python_version(),"sympy":S.__version__})
report("data_basis","Positive root phi of t^2-t-1; declared H3(CD-split).")
phi=(1+S.sqrt(5))/2
simp=lambda z:S.simplify(S.expand(z))
J=(phi,S.Integer(1),1/phi)+(S.Integer(0),)*24
Js=tuple(simp(v) for v in sharp(J))
def B(x,y):return simp(tr(jp(x,y)))
G=S.Matrix([[B(J,J),B(J,Js)],[B(Js,J),B(Js,Js)]])
report("golden.polynomial_residual",simp(phi**2-phi-1))
report("golden.reciprocal_square_sum",simp(phi**2+phi**(-2)))
report("golden.determinant",simp(jnorm(J)))
report("golden.adjoint_diagonal",Js[:3])
report("golden.Gram",G.tolist())
report("golden.Gram_determinant",G.det())
report("golden.mismatch_cosine_squared",simp(G[0,1]**2/(G[0,0]*G[1,1])))
report("golden.mismatch_sine_squared",simp(G.det()/(G[0,0]*G[1,1])))
assert G==S.Matrix([[4,3],[3,4]])
assert G.det()/(G[0,0]*G[1,1])==S.Rational(7,16)
rank_frame=sum(jp(e,e)==e for e in JB[:3])
report("dimensions.octonion_Jordan_frame",(len(OB),len(JB),rank_frame))
report("dimensions.ratio",S.Rational(len(OB),rank_frame))
w=JB[3]
normblock=B(w,w)
gen_norm=dot((1,1,1),(1,1,1))
report("normalization.off_diagonal_unit_trace_square",normblock)
report("normalization.chosen_three_copy_vector_square",gen_norm)
report("normalization.product",simp(1/S.sqrt(normblock)/S.sqrt(gen_norm)))
report("normalization.generations","The three-copy readout is stipulated, not derived here.")
# Golden-only controls within unit-determinant reciprocal triples.
r=S.symbols("r",positive=True)
D=(r,1,1/r)+(S.Integer(0),)*24
Ds=tuple(simp(x) for x in sharp(D))
GD=S.Matrix([[B(D,D),B(D,Ds)],[B(Ds,D),B(Ds,Ds)]])
report("reciprocal_family.determinant",simp(jnorm(D)))
report("reciprocal_family.Gram_determinant",S.factor(GD.det()))
report("negative_control.r3over2.Gram_determinant",simp(GD.det().subs(r,S.Rational(3,2))))
report("negative_control.r3over2.mismatch_sine_squared",
       simp((GD.det()/(GD[0,0]*GD[1,1])).subs(r,S.Rational(3,2))))
# Exact unit determinant plus DET-7 countermodel to the extra equal-norm shorthand.
# u=((11+sqrt(105))/4)^(1/3)>0, spectrum (sqrt(u),sqrt(u),1/u).
u=S.symbols("u",positive=True)
P=2*u**6-11*u**3+2
Drep=(S.sqrt(u),S.sqrt(u),1/u)+(S.Integer(0),)*24
Dreps=tuple(simp(z) for z in sharp(Drep))
UU=B(Drep,Drep); VV=B(Dreps,Dreps)
assert simp(UU-(2*u+u**(-2)))==0 and simp(VV-(2/u+u*u))==0
assert simp(jnorm(Drep))==1 and B(Drep,Dreps)==3
q=(11+S.sqrt(105))/4
assert simp(2*q*q-11*q+2)==0 and q.is_positive
report("DET7_equal_norm_control.q_quadratic_residual",simp(2*q*q-11*q+2))
report("DET7_equal_norm_control.actual_determinant_and_pairing",(simp(jnorm(Drep)),B(Drep,Dreps)))
assert S.rem(S.together(UU*VV-16).as_numer_denom()[0],P,u)==0
norm4poly=S.together(UU-4).as_numer_denom()[0]
gcd=S.gcd(P,norm4poly)
assert gcd==1
report("DET7_equal_norm_control.positive_u_cubed",(11+S.sqrt(105))/4)
report("DET7_equal_norm_control.defining_polynomial",P)
report("DET7_equal_norm_control.Gram_product_mod_polynomial",
       S.rem(S.together(UU*VV-16).as_numer_denom()[0],P,u))
report("DET7_equal_norm_control.gcd_with_norm_equals4",gcd)
report("DET7_equal_norm_control.conclusion",
       "Unit determinant and Gram determinant 7 do not alone imply each squared norm is 4.")
report("registered_C1_coefficient_proof","NOT VERIFIED: source formula only; A558 matrices absent.")
report("registered_generator_uniqueness","NOT VERIFIED: A603/E110 construction absent.")
report("finite_trace_Frobenius_proof","NOT TESTED by this script; see companion C08.2, not dim/rank arithmetic.")
report("physical_attachment","NOT COMPUTED; no angle, mass or generation assignment inferred.")