🍎🍌🍍.sage

R.<a,b,c> = PolynomialRing(QQ, 3)
R2.<x> = PolynomialRing(QQ)

pol = a / (b + c) + b / (c + a) + c / (a + b) - 4

pol *= pol.denominator()
F = R(pol)


def find_initial(F):
    for a0, b0 in cartesian_product([list(range(-20, 20))] * 2):
        if a0 + b0 == 0:
            continue
        for cs0 in pol(a0, b0, x).roots():
            for c0 in cs0:
                if any(t + c0 == 0 for t in [a0, b0]):
                    continue

                if pol(a0, b0, R(c0)) == 0:
                    P = [a0, b0, c0]
                    f = EllipticCurve_from_cubic(F, P)
                    E = f.codomain()
                    return E, f


def find_solutions(F):
    E, f = find_initial(F)
    p0 = E.gen(0)
    p = p0
    finv = f.inverse()
    sols = list()
    while True:
        q = finv(p)
        q.clear_denominators()
        if all(t > 0 for t in q):
            yield tuple(q)
        p = p + p0


for sol in find_solutions(F):
    for k,v in zip('🍎🍌🍍', sol):
        print(f'{k} = {v}')

    break