three-circles.sage

## FIRST SESSION: Find equation of circle through three points:
R.<x,y,x1,y1,x2,y2,x3,y3,u,v,w> = PolynomialRing(QQ,11)
eqn = x^2 + y^2 + u*x + v*y + w
eqn1 = eqn.subs({x:x1,y:y1})
eqn2 = eqn.subs({x:x2,y:y2})
eqn3 = eqn.subs({x:x3,y:y3})
M = Matrix(R,3,3,[[x1,x2,x3],[y1,y2,y3],[1,1,1]])
(R^3)((u,v,w)) * M + (R^3)((x1^2+y1^2, x2^2+y2^2, x3^2+y3^2)) == (R^3)((eqn1,eqn2,eqn3))
(ufrac,vfrac,wfrac) = - (R^3)((x1^2+y1^2, x2^2+y2^2, x3^2+y3^2)) * M.inverse()
## Results:
denom = -y1*x2 + x1*y2 + y1*x3 - y2*x3 - x1*y3 + x2*y3
unumer = y1*x2^2 - x1^2*y2 - y1^2*y2 + y1*y2^2 - y1*x3^2 + y2*x3^2 + x1^2*y3 + y1^2*y3 - x2^2*y3 - y2^2*y3 - y1*y3^2 + y2*y3^2
vnumer = x1^2*x2 + y1^2*x2 - x1*x2^2 - x1*y2^2 - x1^2*x3 - y1^2*x3 + x2^2*x3 + y2^2*x3 + x1*x3^2 - x2*x3^2 + x1*y3^2 - x2*y3^2
wnumer = -y1*x2^2*x3 + x1^2*y2*x3 + y1^2*y2*x3 - y1*y2^2*x3 + y1*x2*x3^2 - x1*y2*x3^2 - x1^2*x2*y3 - y1^2*x2*y3 + x1*x2^2*y3 + x1*y2^2*y3 + y1*x2*y3^2 - x1*y2*y3^2
ufrac = unumer/denom
vfrac = vnumer/denom
wfrac = wnumer/denom

## SECOND SESSION: Find condition for three circles to go through the same point:
R.<x,y,ua,va,wa,ub,vb,wb,uc,vc,wc,ta,tb,tc> = PolynomialRing(QQ,14)
eqna = x^2 + y^2 + ua*x + va*y + wa
eqnb = x^2 + y^2 + ub*x + vb*y + wb
eqnc = x^2 + y^2 + uc*x + vc*y + wc
concurrence = R.ideal([eqna,eqnb,eqnc]).elimination_ideal([x,y]).gens()[0]
concurrence_hom = sum([c*mon*ta^(2-mon.degree(ua)-mon.degree(va)-mon.degree(wa))*tb^(2-mon.degree(ub)-mon.degree(vb)-mon.degree(wb))*tc^(2-mon.degree(uc)-mon.degree(vc)-mon.degree(wc)) for (c,mon) in list(concurrence)])
## Result:
concurrence = va*wa*ub*vb*uc - ua*wa*vb^2*uc - va^2*ub*wb*uc + ua*va*vb*wb*uc - va*wa*vb*uc^2 + wa*vb^2*uc^2 + va^2*wb*uc^2 - va*vb*wb*uc^2 - va*wa*ub^2*vc + ua*wa*ub*vb*vc + ua*va*ub*wb*vc - ua^2*vb*wb*vc + va*wa*ub*uc*vc + ua*wa*vb*uc*vc - 2*wa*ub*vb*uc*vc - 2*ua*va*wb*uc*vc + va*ub*wb*uc*vc + ua*vb*wb*uc*vc - ua*wa*ub*vc^2 + wa*ub^2*vc^2 + ua^2*wb*vc^2 - ua*ub*wb*vc^2 + va^2*ub^2*wc - 2*ua*va*ub*vb*wc + ua^2*vb^2*wc - va^2*ub*uc*wc + ua*va*vb*uc*wc + va*ub*vb*uc*wc - ua*vb^2*uc*wc + ua*va*ub*vc*wc - va*ub^2*vc*wc - ua^2*vb*vc*wc + ua*ub*vb*vc*wc + wa^2*ub^2 + wa^2*vb^2 - 2*ua*wa*ub*wb - 2*va*wa*vb*wb + ua^2*wb^2 + va^2*wb^2 - 2*wa^2*ub*uc + 2*ua*wa*wb*uc + 2*wa*ub*wb*uc - 2*ua*wb^2*uc + wa^2*uc^2 - 2*wa*wb*uc^2 + wb^2*uc^2 - 2*wa^2*vb*vc + 2*va*wa*wb*vc + 2*wa*vb*wb*vc - 2*va*wb^2*vc + wa^2*vc^2 - 2*wa*wb*vc^2 + wb^2*vc^2 + 2*ua*wa*ub*wc - 2*wa*ub^2*wc + 2*va*wa*vb*wc - 2*wa*vb^2*wc - 2*ua^2*wb*wc - 2*va^2*wb*wc + 2*ua*ub*wb*wc + 2*va*vb*wb*wc - 2*ua*wa*uc*wc + 2*wa*ub*uc*wc + 2*ua*wb*uc*wc - 2*ub*wb*uc*wc - 2*va*wa*vc*wc + 2*wa*vb*vc*wc + 2*va*wb*vc*wc - 2*vb*wb*vc*wc + ua^2*wc^2 + va^2*wc^2 - 2*ua*ub*wc^2 + ub^2*wc^2 - 2*va*vb*wc^2 + vb^2*wc^2
concurrence_hom = wa*vb^2*uc^2*ta - va*vb*wb*uc^2*ta - 2*wa*ub*vb*uc*vc*ta + va*ub*wb*uc*vc*ta + ua*vb*wb*uc*vc*ta + wa*ub^2*vc^2*ta - ua*ub*wb*vc^2*ta + va*ub*vb*uc*wc*ta - ua*vb^2*uc*wc*ta - va*ub^2*vc*wc*ta + ua*ub*vb*vc*wc*ta + wb^2*uc^2*ta^2 + wb^2*vc^2*ta^2 - 2*ub*wb*uc*wc*ta^2 - 2*vb*wb*vc*wc*ta^2 + ub^2*wc^2*ta^2 + vb^2*wc^2*ta^2 - va*wa*vb*uc^2*tb + va^2*wb*uc^2*tb + va*wa*ub*uc*vc*tb + ua*wa*vb*uc*vc*tb - 2*ua*va*wb*uc*vc*tb - ua*wa*ub*vc^2*tb + ua^2*wb*vc^2*tb - va^2*ub*uc*wc*tb + ua*va*vb*uc*wc*tb + ua*va*ub*vc*wc*tb - ua^2*vb*vc*wc*tb - 2*wa*wb*uc^2*ta*tb - 2*wa*wb*vc^2*ta*tb + 2*wa*ub*uc*wc*ta*tb + 2*ua*wb*uc*wc*ta*tb + 2*wa*vb*vc*wc*ta*tb + 2*va*wb*vc*wc*ta*tb - 2*ua*ub*wc^2*ta*tb - 2*va*vb*wc^2*ta*tb + wa^2*uc^2*tb^2 + wa^2*vc^2*tb^2 - 2*ua*wa*uc*wc*tb^2 - 2*va*wa*vc*wc*tb^2 + ua^2*wc^2*tb^2 + va^2*wc^2*tb^2 + va*wa*ub*vb*uc*tc - ua*wa*vb^2*uc*tc - va^2*ub*wb*uc*tc + ua*va*vb*wb*uc*tc - va*wa*ub^2*vc*tc + ua*wa*ub*vb*vc*tc + ua*va*ub*wb*vc*tc - ua^2*vb*wb*vc*tc + va^2*ub^2*wc*tc - 2*ua*va*ub*vb*wc*tc + ua^2*vb^2*wc*tc + 2*wa*ub*wb*uc*ta*tc - 2*ua*wb^2*uc*ta*tc + 2*wa*vb*wb*vc*ta*tc - 2*va*wb^2*vc*ta*tc - 2*wa*ub^2*wc*ta*tc - 2*wa*vb^2*wc*ta*tc + 2*ua*ub*wb*wc*ta*tc + 2*va*vb*wb*wc*ta*tc - 2*wa^2*ub*uc*tb*tc + 2*ua*wa*wb*uc*tb*tc - 2*wa^2*vb*vc*tb*tc + 2*va*wa*wb*vc*tb*tc + 2*ua*wa*ub*wc*tb*tc + 2*va*wa*vb*wc*tb*tc - 2*ua^2*wb*wc*tb*tc - 2*va^2*wb*wc*tb*tc + wa^2*ub^2*tc^2 + wa^2*vb^2*tc^2 - 2*ua*wa*ub*wb*tc^2 - 2*va*wa*vb*wb*tc^2 + ua^2*wb^2*tc^2 + va^2*wb^2*tc^2

## THIRD SESSION: Combine these:
R.<x1,y1,x2,y2,x3,y3,xa1,ya1,xa2,ya2,xa3,ya3,xb1,yb1,xb2,yb2,xb3,yb3,xc1,yc1,xc2,yc2,xc3,yc3,ua,va,wa,ub,vb,wb,uc,vc,wc,ta,tb,tc> = PolynomialRing(QQ,36)
denom = -y1*x2 + x1*y2 + y1*x3 - y2*x3 - x1*y3 + x2*y3
unumer = y1*x2^2 - x1^2*y2 - y1^2*y2 + y1*y2^2 - y1*x3^2 + y2*x3^2 + x1^2*y3 + y1^2*y3 - x2^2*y3 - y2^2*y3 - y1*y3^2 + y2*y3^2
vnumer = x1^2*x2 + y1^2*x2 - x1*x2^2 - x1*y2^2 - x1^2*x3 - y1^2*x3 + x2^2*x3 + y2^2*x3 + x1*x3^2 - x2*x3^2 + x1*y3^2 - x2*y3^2
wnumer = -y1*x2^2*x3 + x1^2*y2*x3 + y1^2*y2*x3 - y1*y2^2*x3 + y1*x2*x3^2 - x1*y2*x3^2 - x1^2*x2*y3 - y1^2*x2*y3 + x1*x2^2*y3 + x1*y2^2*y3 + y1*x2*y3^2 - x1*y2*y3^2
dicta = {x1:xa1,y1:ya1,x2:xa2,y2:ya2,x3:xa3,y3:ya3}
dictb = {x1:xb1,y1:yb1,x2:xb2,y2:yb2,x3:xb3,y3:yb3}
dictc = {x1:xc1,y1:yc1,x2:xc2,y2:yc2,x3:xc3,y3:yc3}
denoma = denom.subs(dicta) ; unumera = unumer.subs(dicta) ; vnumera = vnumer.subs(dicta) ; wnumera = wnumer.subs(dicta)
denomb = denom.subs(dictb) ; unumerb = unumer.subs(dictb) ; vnumerb = vnumer.subs(dictb) ; wnumerb = wnumer.subs(dictb)
denomc = denom.subs(dictc) ; unumerc = unumer.subs(dictc) ; vnumerc = vnumer.subs(dictc) ; wnumerc = wnumer.subs(dictc)
# concurrence_hom = ... ## copy line above
fulleqn = concurrence_hom.subs({ua:unumera,va:vnumera,wa:wnumera,ta:denoma, ub:unumerb,vb:vnumerb,wb:wnumerb,tb:denomb, uc:unumerc,vc:vnumerc,wc:wnumerc,tc:denomc})
## This gives a polynomial of degree 18 in 18 variables having 27873486 terms