"""test_626 -- MA open problem 7: the nine components of R reduce to a symmetric tensor W with tr(W P) = 0.

MA Proposition 8.1(5): R_ab^d (antisymmetric in a, b; nine components) satisfies R_ab^b = 0 and
sum_cyc(a,b,c) R_ab^d P_dc = 0. Put W^{cd} = (1/2) eps^{cab} R_ab^d, so R_ab^d = eps_{cab} W^{cd}.
  P1 [exact] the skew part of W is (up to a factor) the trace R_ab^b: eps_{ecd} W^{cd} = R_{de}^d. So the trace identity
     makes W symmetric (nine components to six).
  P2 [exact] eps^{abc} (cyclic sum) = 6 W^{cd} P_dc for every R and P, and the cyclic sum is totally antisymmetric, so the
     cyclic identity is the single scalar identity tr(W P) = 0 (six to five, the count of Theorem C).
Proposed by Gemini (12 November 2026); checked here, in Euclidean index position (orthonormal frame).
"""
import sys, itertools
import sympy as sp
FAIL = []
def check(name, ok, detail=""):
    ok = bool(ok)
    if not ok: FAIL.append(name)
    print(f"[{'PASS' if ok else 'FAIL'}] {name}" + (f"   ({detail})" if detail else ""), flush=True)
E = lambda *ix: sp.LeviCivita(*ix)
R = {}
for a, b, d in itertools.product(range(3), repeat=3):
    if a < b: R[(a, b, d)] = sp.Symbol(f"R{a}{b}{d}")
def Rf(a, b, d):
    if a == b: return 0
    return R[(a, b, d)] if a < b else -R[(b, a, d)]
W = sp.Matrix(3, 3, lambda c, d: sp.Rational(1, 2)*sum(E(c, a, b)*Rf(a, b, d) for a in range(3) for b in range(3)))
recon = all(sp.expand(sum(E(c, a, b)*W[c, d] for c in range(3)) - Rf(a, b, d)) == 0 for a, b, d in itertools.product(range(3), repeat=3))
skew = all(sp.expand(sum(E(e, c, d)*W[c, d] for c in range(3) for d in range(3)) - sum(Rf(d, e, d) for d in range(3))) == 0 for e in range(3))
check("P1 R_ab^d = eps_cab W^cd, and eps_ecd W^cd = R_de^d: the trace identity is exactly W's skew part", recon and skew)
P = sp.Matrix(3, 3, lambda i, j: sp.Symbol(f"P{min(i,j)}{max(i,j)}"))
cyc = lambda a, b, c: sum(Rf(a, b, d)*P[d, c] + Rf(b, c, d)*P[d, a] + Rf(c, a, d)*P[d, b] for d in range(3))
contr = sp.expand(sum(E(a, b, c)*cyc(a, b, c) for a in range(3) for b in range(3) for c in range(3)))
tot_anti = all(sp.expand(cyc(*p) - sp.combinatorics.Permutation(list(p)).signature()*cyc(0, 1, 2)) == 0 for p in itertools.permutations(range(3)))
check("P2 eps^abc (cyclic sum) = 6 tr(W P), and the cyclic sum is totally antisymmetric (one component)",
      sp.expand(contr - 6*sum(W[c, d]*P[d, c] for c in range(3) for d in range(3))) == 0 and tot_anti)
print("\nALL PASS" if not FAIL else f"\nFAILURES: {FAIL}")
sys.exit(1 if FAIL else 0)