p = 115792089237316195423570985008687907853269984665640564039457584007908834671663
Fp = GF(p)
E = EllipticCurve(Fp, [0, 7])

K.<x> = PolynomialRing(Fp, implementation="generic")
L.<y> = PolynomialRing(K, implementation="generic")

def call(f, P):
    Px, Py = P.xy()
    return f(x=Px, y=Py)

def slope_intercept(P0, P1):
    P0x, P0y = P0.xy()
    P1x, P1y = P1.xy()
    m = (P1y - P0y) / (P1x - P0x)
    c = P1y - m*P1x
    return m, c

inf = E(0, 1, 0)

P0 = E(2638891549212558194816434702774699814912586136468438548683319413291233982670, 10140456388517236016202114238224306675726576846632714430386916682749918607464)
P1 = E(17532914127565625088989484000002092349237591241902561765149751398626236254705, 88418181982451178952779569856811854926376929222674054022068577993668451054968)
P2 = -(P0 + P1)
assert P0 + P1 + P2 == inf

m, c = slope_intercept(P0, P1)
f = y - m*x - c
assert call(f, P0) == 0
assert call(f, P1) == 0
assert call(f, P2) == 0

A0 = E(71667150045698532747085079020221438975232032539323499361837608460679887058944, 112531462156938649599975140073878402074313288452223046891973447679626205036672)
A1 = E(24400044380857008437858416414134907169391225728426050882076550827717815954208, 31119844425475844378953714282286295361574062544261840901985213408230446258932)
A2 = -(A0 + A1)
assert A0 + A1 + A2 == inf

λ, μ = slope_intercept(A0, A1)
g = y - λ*x - μ
assert call(g, A0) == 0
assert call(g, A1) == 0
assert call(g, A2) == 0

s = call(f, A0) * call(f, A1) * call(f, A2)
t = call(g, P0) * call(g, P1) * call(g, P2)
assert s == -t
print(s)
print(-t)