"""D1258 C02.1. Full generic cubic element certificate and operator controls.""" 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","Declared H3(CD-split) product; no banked arrays.") X=S.symbols("A0:27") XX=jp(X,X) ss=S.expand((tr(X)**2-tr(XX))/2) assert S.expand(ss-sig(X))==0 Q=sharp(X) R=vs(jp(X,Q),sc(jnorm(X),JI)) res=[S.Poly(S.expand(r),*X) for r in R] assert all(r.is_zero for r in res) report("generic.cubic_coordinate_polynomials_zero",len(res)) report("generic.norm_degree",S.Poly(S.expand(jnorm(X)),*X).total_degree()) report("generic.norm_monomials",len(S.Poly(S.expand(jnorm(X)),*X).terms())) # Differential identity: gradient N paired with trace metric equals adjoint. B=S.diag(*([1]*3+[2*onorm(e) for e in OB]*3)) grad=S.Matrix([S.diff(jnorm(X),x) for x in X]) gradres=grad-B*S.Matrix(Q) assert all(S.expand(t)==0 for t in gradres) report("generic.gradient_N_equals_trace_metric_sharp",True) report("trace_form.inertia(+,-,0)",inertia(B)) # Universal Jordan identity: [L_X,L_(X^2)]=0 coefficientwise for generic X. LX=L(X); LX2=L(tuple(S.expand(z) for z in XX)) comm=LX*LX2-LX2*LX assert all(S.expand(z)==0 for z in comm) report("Jordan_identity.generic_commutator_coordinate_polynomials_zero",len(comm)) # Exact linear-in-A trace-selfadjoint identity on every basis generator. for e in JB: Le=L(e) assert Le.T*B==B*Le report("trace_selfadjoint.basis_operator_checks",len(JB)) # A real multiplication matrix is trace-form self-adjoint, not Euclidean symmetric. D=tuple(S.Integer(i%5-2) for i in range(27)); M=L(D) assert M.T*B==B*M report("sample.L_dimension",M.shape) report("sample.L_trace_selfadjoint",True) report("sample.L_Euclidean_symmetric",M==M.T) # Diagonal counterexample to p_A(L_A)=0. diag=sc(1,JB[0]);diag=va(va(diag,sc(2,JB[1])),sc(3,JB[2])) Ld=L(diag); opres=Ld**3-tr(diag)*Ld**2+sig(diag)*Ld-jnorm(diag)*S.eye(27) assert not opres.is_zero_matrix report("diagonal123.Jordan_cubic_residual",brief(vs(jp(diag,sharp(diag)),sc(jnorm(diag),JI)))) report("diagonal123.L_spectrum",sorted(Ld.eigenvals().items(),key=lambda t:t[0])) report("diagonal123.operator_cubic_nonzero_entries",sum(x!=0 for x in opres)) report("diagonal123.operator_cubic_on_J12_scalar",opres[3,3]) # Hermitian split witness W12(l); l^2=1 but n(l)=-1. W=JB[10]; W2=jp(W,W) assert W2==vn(va(JB[0],JB[1])) assert va(va(W2,jp(JB[0],JB[0])),jp(JB[1],JB[1]))==JZ lam=S.symbols("lambda") p=lam**3-tr(W)*lam**2+sig(W)*lam-jnorm(W) report("split_witness.A_squared",brief(W2)) report("split_witness.nonzero_square_sum_zero",True) report("split_witness.element_characteristic",S.factor(p)) report("split_witness.element_roots",S.solve(p,lam)) LW=L(W) report("split_witness.L_characteristic",S.factor(LW.charpoly(lam).as_expr())) assert LW.T*B==B*LW report("split_witness.L_trace_selfadjoint",True) report("split_witness.L_Euclidean_symmetric",LW==LW.T) # Jordan identity sampled on explicitly deterministic integer pairs, not a universal proof. count=0 for k in range(6): A=tuple(S.Integer((i*i+3*k*i+k)%5-2) for i in range(27)) Y=tuple(S.Integer((3*i*i+k*i+2)%5-2) for i in range(27)) A2=jp(A,A) assert zero_vector(vs(jp(A2,jp(A,Y)),jp(A,jp(A2,Y)))) count+=1 report("Jordan_identity.sample_pairs_passed",count) report("Jordan_identity.universal_product_proof","Coefficientwise generic polynomial certificate, not inferred from samples.") report("banked_array_alignment","NOT VERIFIED")