"""SymPy verification of the Weyl and Poisson lifts.

The identities JA = I and b_l^(ij) = 0 simultaneously verify the
defining relations for both constructions.
"""

from sympy import symbols, Matrix, Rational, expand, factor, diff, eye

x, y, z = symbols('x y z')
u = 1 + x*y

# Alpoge's map (posted 2026-07-20 UTC)
F1 = u**3 * z + y**2 * u * (4 + 3*x*y)
F2 = y + 3*x*u**2 * z + 3*x*y**2 * (4 + 3*x*y)
F3 = 2*x - 3*x**2*y - x**3*z
F = Matrix([F1, F2, F3])
X = Matrix([x, y, z])

# Step 1: Jacobian determinant
J = F.jacobian(X)
detJ = expand(J.det())
print("det J =", detJ)
assert detJ == -2

# Step 2: the 3-point collision (exact arithmetic)
pts = [(0, 0, Rational(-1,4)), (1, Rational(-3,2), Rational(13,2)), (-1, Rational(3,2), Rational(13,2))]
for p in pts:
    val = F.subs({x: p[0], y: p[1], z: p[2]})
    print("F", p, "=", list(val))
    assert list(val) == [Rational(-1,4), 0, 0]

# Step 3: inverse Jacobian A = J^{-1} = adj(J)/(-2), polynomial entries
A = J.adjugate().applyfunc(expand) * Rational(-1, 2)
# sanity: J*A == I
assert (J*A).applyfunc(expand) == eye(3)
assert (A*J).applyfunc(expand) == eye(3)
print("\nA = J^{-1} entries are polynomials, max total degree per entry:")
from sympy import Poly
for i in range(3):
    for j in range(3):
        print(f"  A[{i},{j}] degree {Poly(A[i,j], x, y, z).total_degree()}, {len(Poly(A[i,j], x, y, z).terms())} terms")

# Step 4: identities shared by the Weyl and Poisson lifts
# Relation [phi(d_i), phi(x_j)] = delta_ij  <=>  sum_k A[k,i] dF_j/dx_k = delta_ij (i.e. (J A)^T = I, done above)
# Relation [phi(d_i), phi(d_j)] = 0  <=> for all l:
#   sum_k ( A[k,i]*d(A[l,j])/dx_k - A[k,j]*d(A[l,i])/dx_k ) = 0
vars_ = [x, y, z]
ok = True
for i in range(3):
    for j in range(i+1, 3):
        for l in range(3):
            c = sum(A[k,i]*diff(A[l,j], vars_[k]) - A[k,j]*diff(A[l,i], vars_[k]) for k in range(3))
            c = expand(c)
            if c != 0:
                ok = False
                print("NONZERO commutator coefficient", i, j, l, c)
assert ok, "A nonzero commuting-vector-field coefficient was found"
print("\nWeyl relations verified:")
print("  [phi(d_i), phi(x_j)] = delta_ij")
print("  [phi(d_i), phi(d_j)] = 0")
print("\nPoisson relations verified from the same identities:")
print("  {Phi(xi_i), Phi(x_j)} = delta_ij")
print("  {Phi(xi_i), Phi(xi_j)} = 0")

# Step 5: dump the explicit endomorphism, factored where possible
print("\n=== Explicit endomorphism phi of A_3 ===")
for i, Fi in enumerate([F1, F2, F3], 1):
    print(f"phi(x{i}) =", factor(Fi))
print()
for i in range(3):
    terms = []
    for k in range(3):
        terms.append(f"({factor(A[k,i])}) * d{k+1}")
    print(f"phi(d{i+1}) =")
    for t in terms:
        print("    +", t)
    print()