Pasted

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
"""