R.<x0,y0,x1,y1,x2,y2,x3,y3,r0,r1,r2,r3> = PolynomialRing(QQ, order='lex')

def tangent(Pi, Pj):
    xi, yi, ri = Pi
    xj, yj, rj = Pj
    return ideal((xi-xj)^2+(yi-yj)^2-(ri+rj)^2)

P = [(x0, y0, r0),
     (x1, y1, r1),
     (x2, y2, r2),
     (x3, y3, r3)]

I = (tangent(P[0], P[1]) +
     tangent(P[0], P[2]) +
     tangent(P[0], P[3]) +
     tangent(P[1], P[2]) +
     tangent(P[1], P[3]) +
     tangent(P[2], P[3]))

G = I.groebner_basis()

"""
sage: G
Polynomial Sequence with 65 Polynomials in 12 Variables
sage: G[-1]
r0^2*r1^2*r2^2 - 2*r0^2*r1^2*r2*r3 + r0^2*r1^2*r3^2 - 2*r0^2*r1*r2^2*r3 - 2*r0^2*r1*r2*r3^2 + r0^2*r2^2*r3^2 - 2*r0*r1^2*r2^2*r3 - 2*r0*r1^2*r2*r3^2 - 2*r0*r1*r2^2*r3^2 + r1^2*r2^2*r3^2
"""