heptic area

# sagemath script

# setting poly_x0 = -1, poly_x1 = 1 (without loss of generality) we will want to reduce the open logical formula "∃ cops. ∀ p. the areas are proportional with coefficient λ" to a single equation in only λ. it is therefore important that lam comes last in the following declaration:

R.<p0, p1, p2, p3, p4, p5, p6, p7, x, poly_x0, poly_x1, cops_x0, cops_x1, cops_x2, cops_x3, lam> = PolynomialRing(QQ, order='lex')

y = p0 + p1*x + p2*x^2 + p3*x^3 + p4*x^4 + p5*x^5 + p6*x^6 + p7*x^7

cops_y0 = y(x=cops_x0)
cops_y1 = y(x=cops_x1)
cops_y2 = y(x=cops_x2)
cops_y3 = y(x=cops_x3)

# sum of two triangles
cops_area = (1/2) * ((-cops_x0+cops_x1)*(-cops_y0+cops_y2) - (-cops_y0+cops_y1)*(-cops_x0+cops_x2)) + (1/2) * ((-cops_x0+cops_x2)*(-cops_y0+cops_y3) - (-cops_y0+cops_y2)*(-cops_x0+cops_x3))

poly_y0 = y(x=poly_x0)
poly_y1 = y(x=poly_x1)

# doing this the orthodox way
poly_area = -integral(y, x)(x=poly_x0) + integral(y, x)(x=poly_x1) - (1/2) * (-poly_x0 + poly_x1) * (poly_y0 + poly_y1)

I = ideal(
    # "∀ p. " is expressed in a mildly hacky way
    diff(lam*cops_area-poly_area, p0),
    diff(lam*cops_area-poly_area, p1),
    diff(lam*cops_area-poly_area, p2),
    diff(lam*cops_area-poly_area, p3),
    diff(lam*cops_area-poly_area, p4),
    diff(lam*cops_area-poly_area, p5),
    diff(lam*cops_area-poly_area, p6),
    diff(lam*cops_area-poly_area, p7)
)

# convenient coordinate system
J = I.subs(poly_x0=-1, poly_x1=1)

G = J.groebner_basis()
for i in G:
    print(str(i) + '\n')

"""
cops_x0 - 108/175*cops_x3^3*lam^2 + 35/18*cops_x3^3 + 8748/30625*cops_x3^2*lam^3 - 9/5*cops_x3^2*lam + 648/875*cops_x3*lam^2 - 11/6*cops_x3 - 26244/153125*lam^3 + 27/25*lam

cops_x1 + 216/175*cops_x3^3*lam^2 - 35/9*cops_x3^3 - 1296/875*cops_x3*lam^2 + 14/3*cops_x3

cops_x2 - 108/175*cops_x3^3*lam^2 + 35/18*cops_x3^3 - 8748/30625*cops_x3^2*lam^3 + 9/5*cops_x3^2*lam + 648/875*cops_x3*lam^2 - 11/6*cops_x3 + 26244/153125*lam^3 - 27/25*lam

cops_x3^4 - 6/5*cops_x3^2 + 27/175

lam^4 - 6125/972*lam^2 + 5359375/944784
"""