nonic case

# 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, p8, p9, x, poly_x0, poly_x1, cops_x0, cops_x1, cops_x2, cops_x3, cops_x4, 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 + p8*x^8 + p9*x^9

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

# sum of three 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)) + (1/2) * ((-cops_x0+cops_x3)*(-cops_y0+cops_y4) - (-cops_y0+cops_y3)*(-cops_x0+cops_x4))

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),
    diff(lam*cops_area-poly_area, p8),
    diff(lam*cops_area-poly_area, p9)
)

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

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

# [1]