In the Math Stack Exchange post titled Six points lie on a circle, Oai Thanh Đào asked for a proof of a certain cocircularity given some others. As the setup is invariant under Möbius transformations, I restricted myself to a scenario where three of the points are at fixed positions.
Original problem statement
Here is what Oai Thanh Đào originally asked:
I am looking for a proof of a problem as follows:
Let six points lie on a circle. Define be any circle through points . Let we take modulo 6. Let lie on a circle . Let . Show that also lie on the circle and six points lie on a circle.
(Question is licensed cc-by-sa 3.0 as per Stack Exchange license conditions)
My answer
Here is the answer I wrote to the question above:
Choice of coordinates
I like coordinates. The incidence relations you are interested in are invariant under Möbius transformations, so a natural setup for this problem would be the one-point compactification of the plane. Think if you want to, but I’ll not use complex numbers for coordinates in this answer. Since your setup is invariant under Möbius transformations, you can choose specific coordinates of three points without loss of generality. So I’ll choose to be the point at infinity, and fixed coordinates on the axis for two other points. Thus the “circle” through all the is the axis in my setup, and I can use parameters for the coordinates of the remaining points.
The interactive illustration above shows this specific version of your configuration. You can move the points around and see the configuration change accordingly. I alternated between fixed and variable points in an attempt to make things more symmetric, in the hope of making individual expressions not too complicated.
Next we can pick an arbitrary point anywhere in the plane. This results in the
computed as the second point of intersection between the corresponding circles.
First condition for cocircularity of the Bi
Choosing will define the circle through . On the other hand this allows us to construct as the intersection of and . We can and check whether that point lies on the same circle. In general it does not. The condition for to be cocircular can be factored into the following independent conditions:
The first three of these conditions can be seen as non-degeneracy constraints. If then , which is a degenerate situation. The second condition expresses the situation where lies on a certain circle around , and characterizes the situations where and would coincide. So if we can assume that all points in the configuration are distinct, then we can concentrate on the last condition.
It describes a circle, so we now know that the point must lie on a certain circle if is to be cocircular with . Actually this circle already passes through and . So knowing all the and already fixes the circle on which all the must lie.
Second condition for cocircularity of the Bi
We can play the same game for and instead of and . What’s the conditions for being cocircular if we choose ?
The last conditon is the same as above, so if the cocircularity condition is satisfied for and all points are distinct, and if is chosen from that same circle, too, then the cocircularity condition must be satisfied for as well, and all six are cocircular.
Conditions for cocircularity of the Ci
A similar game can be played for the as well. Their coordinates are considerably more complicated, though, and so are the conditions for cocircularity of four of these. But assuming the to be cocircular (i.e. assuming the respective last condition to hold for and ), one can demonstrate that all the points must be cocircular, too. In particular, if the condition for holds, then one can show that and are cocircular, and if the condition for holds, then and are cocircular. If both hold, then all points must be cocircular. You could start by defining the circle via , then add due to , add due to and add due to . All again assuming distinct points, of course.
My computation
This is the Sage code which I used to compute the answer above.
First define a polynomial ring for all operations:
PR.<a,b,c,d,e,x1,y1,x2,y2> = ZZ[]
Declare the matrix of the bilinear form of our Möbius geometry:
L = matrix([[2, 0, 0, 0], [0, 2, 0, 0], [0, 0, 0, -1], [0, 0, -1, 0]]); L
[ 2 0 0 0]
[ 0 2 0 0]
[ 0 0 0 -1]
[ 0 0 -1 0]
Declare the point at infinity:
Pinf = vector([0, 0, 1, 0])
Define some utility functions:
def simpl(v, simplify=True):
"""Simplify homogeneous coordinate vectors by canceling common factors."""
if not simplify:
return v
if v.is_zero():
return v
g = gcd(v.list())
if g.degree() > 0:
print("Canceling {}".format(g.factor()))
return v.parent()(v / g)
def hom2circ(x, y, z):
"""Turn a point with given homogeneous coordinates into a circle of radius zero."""
res = vector([x*z, y*z, x^2+y^2, z^2])
assert not (res*L*res)
return res
def circThrough(a, b, c, simplify=True):
"""If a,b,c are points construct circle through them.
More generally construct circle perpendicular to a,b,c."""
d, e, f, g = (matrix([a, b, c])*L).minors(3)
res = simpl(vector([-g, f, -e, d]), simplify)
assert not (a*L*res or b*L*res or c*L*res)
return res
def otherPOI(a, b, p, simplify=True):
"""Other point of intersection"""
assert not (a*L*p or b*L*p)
# we want aLq = bLq = qLq = 0 but q ≠ p
h = circThrough(a, b, Pinf, simplify) # line connecting centers
q = simpl((h*L*h)*p - 2*(p*L*h)*h, simplify)
assert not (a*L*q or b*L*q or q*L*q)
return q
def cocircular(a, b, c, d):
"""Check whether 4 points are cocircular. Returns zero if that is the case.
More generally whether there exists a circle perpendicular to the 4 given ones."""
return matrix([a, b, c, d]).det()
Choose coordinates on the real axis for all the , with three points fixed and three points controlled by a single parameter. All of this is without loss of generality.
A1 = hom2circ(1, 0, 0)
A2 = hom2circ(a, 0, 1)
A3 = hom2circ(0, 0, 1)
A4 = hom2circ(b, 0, 1)
A5 = hom2circ(1, 0, 1)
A6 = hom2circ(c, 0, 1)
Now pick two generic coordinates for and anywhere in the plane.
B1 = hom2circ(d, e, 1)
B2 = hom2circ(x1, y1, 1)
This defines four circles, where I write Cij for .
C12 = circThrough(A1, A2, B1)
C45 = circThrough(A4, A5, B1)
C23 = circThrough(A2, A3, B2)
C56 = circThrough(A5, A6, B2)
Canceling 2 * (b - 1)
Canceling 2 * a
Canceling 2 * (c - 1)
These circles in turn define two more points, as extra points of intersection.
B4 = otherPOI(C12, C45, B1)
B5 = otherPOI(C23, C56, B2)
Canceling 2^2 * e
Canceling 2^2 * y1
To give an example of how such a point of intersection would actually be formulated, we have a closer look at .
B4[:2]/B4[-1]
((a^2*b - a^2*d - a*b*d + a*d^2 + a*e^2 + a^2 - a*b - a*d + b*d)/(a^2 - 2*a*d + d^2 + e^2),
(a^2*e - a*b*e - a*e + b*e)/(a^2 - 2*a*d + d^2 + e^2))
(a*(b-d)+(a-b))*(a-d)+a*e^2 - B4[0]
0
(a-b)*(a-1)*e - B4[1]
0
(a-d)^2+e^2 - B4[-1]
0
So here we know that can be written as
Next, investigate the condition that are cocircular.
cond1 = cocircular(B1, B2, B4, B5) # whole condition
cond1.factor()
(a - 1)
* (a*b - 2*a*d + d^2 + e^2 + a - b)
* (-a*x1^2 + c*x1^2 - a*y1^2 + c*y1^2 + a*c - 2*c*x1 + x1^2 + y1^2)
* (a*b*e*x1 - b*c*e*x1 - a*e*x1^2 + c*e*x1^2 - a*b*d*y1 + b*c*d*y1 + a*d^2*y1 - c*d^2*y1 + a*e^2*y1 - c*e^2*y1 - a*e*y1^2 + c*e*y1^2 - a*c*e + b*c*e + a*e*x1 - b*e*x1 + a*c*y1 - b*c*y1 - a*d*y1 + b*d*y1)
Verify that the second and third factors describe situations where pairs of coincide.
gcd(matrix([B1, B4]).minors(2)) == cond1.factor()[1][0] # B1 and B4 coincide
True
gcd(matrix([B2, B5]).minors(2)) == cond1.factor()[2][0] # B2 and B5 coincide
True
So now we can concentrate on the last factor, which is an actual requirement for .
cond2 = -cond1.factor()[-1][0] # last factor
for _ in cond1.factor()[-1][0].polynomial(x1):
print _.polynomial(_.parent(y1))
(-a*e + c*e)*y1^2 + (-a*b*d + b*c*d + a*d^2 - c*d^2 + a*e^2 - c*e^2 + a*c - b*c - a*d + b*d)*y1 - a*c*e + b*c*e
a*b*e - b*c*e + a*e - b*e
-a*e + c*e
I manually simplified these terms, and used the following to check my own version:
(1
* (a-1)
* ((d-a)^2+e^2-(1-a)*(b-a))
* ((a-c-1)*(x1^2+y1^2) + 2*c*x1 - a*c)
* (((a-c)*e)*(x1^2+y1^2)
+ (-a*b + b*c - a + b)*e*x1
+ ((a-c)*((b-d)*d - e^2) - (a-b)*(c-d))*y1
+ (a-b)*c*e)) == cond1
True
Check that and lie on the circle described by the last factor.
cond2(x1=B1[0]/B1[-1], y1=B1[1]/B1[-1])
0
cond2(x1=B4[0]/B4[-1], y1=B4[1]/B4[-1])
0
We might also want to express this circle as a circle vector.
CB = vector([
(-a*b+b*c-a+b)*e,
(a-c)*((b-d)*d-e^2)-(a-b)*(c-d),
-2*(a-b)*c*e,
-2*(a-c)*e
])
Make sure the vector describes the same circle as the condition.
CB*L*B2 == 2*cond2
True
Using this circle formulation we have another way to express and lying on the same circle:
[CB*L*B1, CB*L*B4]
[0, 0]
Next, pick as another generic point in the plane and do the same steps as above:
B3 = hom2circ(x2, y2, 1)
C34 = circThrough(A3, A4, B3)
C61 = circThrough(A6, A1, B3)
B6 = otherPOI(C34, C61, B3)
Canceling 2 * b
Canceling 2^2 * y2
cond3 = cocircular(B1, B3, B4, B6)
cond3.factor()
(-1)
* (b*c - 2*c*x2 + x2^2 + y2^2)
* (a*b - 2*a*d + d^2 + e^2 + a - b)
* (a*b*e*x2 - b*c*e*x2 - a*e*x2^2 + c*e*x2^2 - a*b*d*y2 + b*c*d*y2 + a*d^2*y2 - c*d^2*y2 + a*e^2*y2 - c*e^2*y2 - a*e*y2^2 + c*e*y2^2 - a*c*e + b*c*e + a*e*x2 - b*e*x2 + a*c*y2 - b*c*y2 - a*d*y2 + b*d*y2)
cond4 = -cond3.factor()[-1][0]
cond4(x2=x1, y2=y1) == cond2
True
The above shows that the last condition is the same for being cocircular and for being cocircular. Again I reformulated the whole condition, and want to make sure I didn’t have a sign wrong.
cond3 == (1
* ((x2-c)^2+y2^2-(c-b)*c)
* ((d-a)^2+e^2-(1-a)*(b-a))
* cond2(x1=x2, y1=y2))
True
Now on to the points which result from yet another bunch of circle-circle intersections.
C2 = otherPOI(C12, C23, A2)
C3 = otherPOI(C23, C34, A3)
C4 = otherPOI(C34, C45, A4)
C5 = otherPOI(C45, C56, A5)
C6 = otherPOI(C56, C61, A6)
Unfortunately can’t be computed like this, since is the point at infinity.
C1 = otherPOI(C61, C12, A1); C1
Canceling (-1) * 2^2 * (-c*e + e*x2 + a*y2 - d*y2)
(0, 0, 0, 0)
Instead we can do a line-line intersection like this:
C1 = simpl(hom2circ(*(diagonal_matrix([1,1,-2])*C61[:3].cross_product(C12[:3]))))
[C1*L*C61, C1*L*C12]
[0, 0]
Let’s check for cocirularities again. Doing one example first, to know what kind of expressions we’re looking at.
cocircular(C1, C2, C3, C4).factor()
y2
* e
* (a*b*e*x2 - b*c*e*x2 - a*e*x2^2 + c*e*x2^2 - a*b*d*y2 + b*c*d*y2 + a*d^2*y2 - c*d^2*y2 + a*e^2*y2 - c*e^2*y2 - a*e*y2^2 + c*e*y2^2 - a*c*e + b*c*e + a*e*x2 - b*e*x2 + a*c*y2 - b*c*y2 - a*d*y2 + b*d*y2)
* (-a*c*e^2*x1 + c*e^2*x1^2 + a^2*c*e*y1 - a*c*d*e*y1 + c*e^2*y1^2 + a*e^2*x1*x2 - e^2*x1^2*x2 - a^2*e*y1*x2 + a*d*e*y1*x2 - e^2*y1^2*x2 + a^2*e*x1*y2 - a*d*e*x1*y2 - a*e*x1^2*y2 + d*e*x1^2*y2 - a^2*c*y1*y2 + 2*a*c*d*y1*y2 - c*d^2*y1*y2 + a*e^2*y1*y2 - c*e^2*y1*y2 - a*e*y1^2*y2 + d*e*y1^2*y2)
* (a*b*e*x1*y1*x2 - b*e*x1^2*y1*x2 - a^2*b*y1^2*x2 + a*b*d*y1^2*x2 - b*e*y1^3*x2 - a*e*x1*y1*x2^2 + e*x1^2*y1*x2^2 + a^2*y1^2*x2^2 - a*d*y1^2*x2^2 + e*y1^3*x2^2 - a^2*e*x1^2*y2 + 2*a*e*x1^3*y2 - e*x1^4*y2 + a^2*b*x1*y1*y2 - a*b*d*x1*y1*y2 - a*b*x1^2*y1*y2 + b*d*x1^2*y1*y2 - a^2*e*y1^2*y2 + a*b*e*y1^2*y2 + 2*a*e*x1*y1^2*y2 - 2*e*x1^2*y1^2*y2 - a*b*y1^3*y2 + b*d*y1^3*y2 - e*y1^4*y2 - a*e*x1*y1*y2^2 + e*x1^2*y1*y2^2 + a^2*y1^2*y2^2 - a*d*y1^2*y2^2 + e*y1^3*y2^2)
* (a*b^2*e*y1*x2^2 - 2*a*b*e*y1*x2^3 + a*e*y1*x2^4 - a*b^2*e*x1*x2*y2 + b^2*e*x1^2*x2*y2 - a*b^2*d*y1*x2*y2 + a*b*d^2*y1*x2*y2 + a*b*e^2*y1*x2*y2 + b^2*e*y1^2*x2*y2 + a*b*e*x1*x2^2*y2 - b*e*x1^2*x2^2*y2 + a*b*d*y1*x2^2*y2 - a*d^2*y1*x2^2*y2 - a*e^2*y1*x2^2*y2 - b*e*y1^2*x2^2*y2 + a*b^2*d*x1*y2^2 - a*b*d^2*x1*y2^2 - a*b*e^2*x1*y2^2 - b^2*d*x1^2*y2^2 + b*d^2*x1^2*y2^2 + b*e^2*x1^2*y2^2 - b^2*d*y1^2*y2^2 + b*d^2*y1^2*y2^2 + b*e^2*y1^2*y2^2 - 2*a*b*e*y1*x2*y2^2 + 2*a*e*y1*x2^2*y2^2 + a*b*e*x1*y2^3 - b*e*x1^2*y2^3 + a*b*d*y1*y2^3 - a*d^2*y1*y2^3 - a*e^2*y1*y2^3 - b*e*y1^2*y2^3 + a*e*y1*y2^4 - b^2*e*y1*x2^2 + 2*b*e*y1*x2^3 - e*y1*x2^4 + a*b*e*x1*x2*y2 - b*e*x1^2*x2*y2 + a*b^2*y1*x2*y2 - a*b*d*y1*x2*y2 - b*e*y1^2*x2*y2 - a*e*x1*x2^2*y2 + e*x1^2*x2^2*y2 - a*b*y1*x2^2*y2 + a*d*y1*x2^2*y2 + e*y1^2*x2^2*y2 - a*b^2*x1*y2^2 + a*b*d*x1*y2^2 + b^2*x1^2*y2^2 - b*d*x1^2*y2^2 + a*b*e*y1*y2^2 - b^2*e*y1*y2^2 + b^2*y1^2*y2^2 - b*d*y1^2*y2^2 + 2*b*e*y1*x2*y2^2 - 2*e*y1*x2^2*y2^2 - a*e*x1*y2^3 + e*x1^2*y2^3 - a*b*y1*y2^3 + a*d*y1*y2^3 + e*y1^2*y2^3 - e*y1*y2^4)
Now let’s compute all the cocircularity conditions together. Since a single cocircularity condition is the value of a determinant, all of them can be computed as the list of all such determinants, or minors:
cond5 = matrix([C1, C2, C3, C4, C5, C6]).minors(4) # Warning, this takes some time!
The minors are returned in lexicographical order, so we can compute the quadruple of circles for each of them in the following way:
conditionLabels = sorted(Combinations(["C1", "C2", "C3", "C4", "C5", "C6"], 4))
conditionLabels
[['C1', 'C2', 'C3', 'C4'],
['C1', 'C2', 'C3', 'C5'],
['C1', 'C2', 'C3', 'C6'],
['C1', 'C2', 'C4', 'C5'],
['C1', 'C2', 'C4', 'C6'],
['C1', 'C2', 'C5', 'C6'],
['C1', 'C3', 'C4', 'C5'],
['C1', 'C3', 'C4', 'C6'],
['C1', 'C3', 'C5', 'C6'],
['C1', 'C4', 'C5', 'C6'],
['C2', 'C3', 'C4', 'C5'],
['C2', 'C3', 'C4', 'C6'],
['C2', 'C3', 'C5', 'C6'],
['C2', 'C4', 'C5', 'C6'],
['C3', 'C4', 'C5', 'C6']]
Now let’s find out which cocircularities are implied by each of our conditions.
[_[1] for _ in zip(cond5, conditionLabels) if cond2.divides(_[0])]
[['C1', 'C2', 'C5', 'C6'], ['C2', 'C3', 'C4', 'C5']]
[_[1] for _ in zip(cond5, conditionLabels) if cond4.divides(_[0])]
[['C1', 'C2', 'C3', 'C4'], ['C1', 'C4', 'C5', 'C6']]
All other cocircularities between any four of these can be deduced from these.
Coordinates
To make it easier to verify individual steps, I’ll just include coordinates for all the relevant objects.
C12
(e, a - d, 2*a*e, 0)
C23
(a*y1, -a*x1 + x1^2 + y1^2, 0, 2*y1)
-C34
(b*y2, -b*x2 + x2^2 + y2^2, 0, 2*y2)
C45
(b*e + e, -b*d + d^2 + e^2 + b - d, 2*b*e, 2*e)
-C56
(c*y1 + y1, -c*x1 + x1^2 + y1^2 + c - x1, 2*c*y1, 2*y1)
-C61
(y2, c - x2, 2*c*y2, 0)
B1.column()
[ d]
[ e]
[d^2 + e^2]
[ 1]
B2.column()
[ x1]
[ y1]
[x1^2 + y1^2]
[ 1]
B3.column()
[ x2]
[ y2]
[x2^2 + y2^2]
[ 1]
B4.column()
[ a^2*b - a^2*d - a*b*d + a*d^2 + a*e^2 + a^2 - a*b - a*d + b*d]
[ a^2*e - a*b*e - a*e + b*e]
[a^2*b^2 - 2*a^2*b*d + a^2*d^2 + a^2*e^2 + 2*a^2*b - 2*a*b^2 - 2*a^2*d + 2*a*b*d + a^2 - 2*a*b + b^2]
[ a^2 - 2*a*d + d^2 + e^2]
B5.column()
[ a^2*c*x1 - a*c^2*x1 - a*c*x1^2 + c^2*x1^2 - a*c*y1^2 + c^2*y1^2 + a*c^2 - a*c*x1 - c^2*x1 + c*x1^2 + c*y1^2]
[ a^2*c*y1 - a*c^2*y1 - a*c*y1 + c^2*y1]
[ a^2*c^2 - 2*a*c^2*x1 + c^2*x1^2 + c^2*y1^2]
[a^2*x1^2 - 2*a*c*x1^2 + c^2*x1^2 + a^2*y1^2 - 2*a*c*y1^2 + c^2*y1^2 + 2*a*c*x1 - 2*c^2*x1 - 2*a*x1^2 + 2*c*x1^2 - 2*a*y1^2 + 2*c*y1^2 + c^2 - 2*c*x1 + x1^2 + y1^2]
B6.column()
[ b*c^2 - b*c*x2 - c^2*x2 + c*x2^2 + c*y2^2]
[ -b*c*y2 + c^2*y2]
[b^2*c^2 - 2*b*c^2*x2 + c^2*x2^2 + c^2*y2^2]
[ c^2 - 2*c*x2 + x2^2 + y2^2]
C1.column()
[ a*c^2*e^2 - 2*a*c*e^2*x2 + a*e^2*x2^2 - a^2*c*e*y2 - a*c^2*e*y2 + a*c*d*e*y2 + c^2*d*e*y2 + a^2*e*x2*y2 + a*c*e*x2*y2 - a*d*e*x2*y2 - c*d*e*x2*y2 + a^2*c*y2^2 - 2*a*c*d*y2^2 + c*d^2*y2^2]
[ -a*c*e^2*y2 + c^2*e^2*y2 + a*e^2*x2*y2 - c*e^2*x2*y2 + a^2*e*y2^2 - a*c*e*y2^2 - a*d*e*y2^2 + c*d*e*y2^2]
[a^2*c^2*e^2 - 2*a^2*c*e^2*x2 + a^2*e^2*x2^2 - 2*a^2*c^2*e*y2 + 2*a*c^2*d*e*y2 + 2*a^2*c*e*x2*y2 - 2*a*c*d*e*x2*y2 + a^2*c^2*y2^2 - 2*a*c^2*d*y2^2 + c^2*d^2*y2^2 + a^2*e^2*y2^2 - 2*a*c*e^2*y2^2 + c^2*e^2*y2^2]
[ c^2*e^2 - 2*c*e^2*x2 + e^2*x2^2 - 2*a*c*e*y2 + 2*c*d*e*y2 + 2*a*e*x2*y2 - 2*d*e*x2*y2 + a^2*y2^2 - 2*a*d*y2^2 + d^2*y2^2]
C2.column()
[ a^2*e*x1*y1 - a*d*e*x1*y1 - a*e*x1^2*y1 + d*e*x1^2*y1 + a*e^2*y1^2 - a*e*y1^3 + d*e*y1^3]
[ -a*e^2*x1*y1 + e^2*x1^2*y1 + a^2*e*y1^2 - a*d*e*y1^2 + e^2*y1^3]
[a^2*e^2*x1^2 - 2*a*e^2*x1^3 + e^2*x1^4 + a^2*e^2*y1^2 - 2*a*e^2*x1*y1^2 + 2*e^2*x1^2*y1^2 + e^2*y1^4]
[ a^2*y1^2 - 2*a*d*y1^2 + d^2*y1^2 + e^2*y1^2]
C3.column()
[a*b^2*y1^2*x2^2 - 2*a*b*y1^2*x2^3 + a*y1^2*x2^4 - a^2*b*x1*y1*x2*y2 - a*b^2*x1*y1*x2*y2 + a*b*x1^2*y1*x2*y2 + b^2*x1^2*y1*x2*y2 + a*b*y1^3*x2*y2 + b^2*y1^3*x2*y2 + a^2*x1*y1*x2^2*y2 + a*b*x1*y1*x2^2*y2 - a*x1^2*y1*x2^2*y2 - b*x1^2*y1*x2^2*y2 - a*y1^3*x2^2*y2 - b*y1^3*x2^2*y2 + a^2*b*x1^2*y2^2 - 2*a*b*x1^3*y2^2 + b*x1^4*y2^2 - 2*a*b*x1*y1^2*y2^2 + 2*b*x1^2*y1^2*y2^2 + b*y1^4*y2^2 - 2*a*b*y1^2*x2*y2^2 + 2*a*y1^2*x2^2*y2^2 + a^2*x1*y1*y2^3 + a*b*x1*y1*y2^3 - a*x1^2*y1*y2^3 - b*x1^2*y1*y2^3 - a*y1^3*y2^3 - b*y1^3*y2^3 + a*y1^2*y2^4]
[ -a^2*b*y1^2*x2*y2 + a*b^2*y1^2*x2*y2 + a^2*y1^2*x2^2*y2 - a*b*y1^2*x2^2*y2 + a^2*b*x1*y1*y2^2 - a*b^2*x1*y1*y2^2 - a*b*x1^2*y1*y2^2 + b^2*x1^2*y1*y2^2 - a*b*y1^3*y2^2 + b^2*y1^3*y2^2 + a^2*y1^2*y2^3 - a*b*y1^2*y2^3]
[ a^2*b^2*y1^2*x2^2 - 2*a^2*b*y1^2*x2^3 + a^2*y1^2*x2^4 - 2*a^2*b^2*x1*y1*x2*y2 + 2*a*b^2*x1^2*y1*x2*y2 + 2*a*b^2*y1^3*x2*y2 + 2*a^2*b*x1*y1*x2^2*y2 - 2*a*b*x1^2*y1*x2^2*y2 - 2*a*b*y1^3*x2^2*y2 + a^2*b^2*x1^2*y2^2 - 2*a*b^2*x1^3*y2^2 + b^2*x1^4*y2^2 - 2*a*b^2*x1*y1^2*y2^2 + 2*b^2*x1^2*y1^2*y2^2 + b^2*y1^4*y2^2 - 2*a^2*b*y1^2*x2*y2^2 + 2*a^2*y1^2*x2^2*y2^2 + 2*a^2*b*x1*y1*y2^3 - 2*a*b*x1^2*y1*y2^3 - 2*a*b*y1^3*y2^3 + a^2*y1^2*y2^4]
[ b^2*y1^2*x2^2 - 2*b*y1^2*x2^3 + y1^2*x2^4 - 2*a*b*x1*y1*x2*y2 + 2*b*x1^2*y1*x2*y2 + 2*b*y1^3*x2*y2 + 2*a*x1*y1*x2^2*y2 - 2*x1^2*y1*x2^2*y2 - 2*y1^3*x2^2*y2 + a^2*x1^2*y2^2 - 2*a*x1^3*y2^2 + x1^4*y2^2 + a^2*y1^2*y2^2 - 2*a*b*y1^2*y2^2 + b^2*y1^2*y2^2 - 2*a*x1*y1^2*y2^2 + 2*x1^2*y1^2*y2^2 + y1^4*y2^2 - 2*b*y1^2*x2*y2^2 + 2*y1^2*x2^2*y2^2 + 2*a*x1*y1*y2^3 - 2*x1^2*y1*y2^3 - 2*y1^3*y2^3 + y1^2*y2^4]
C4.column()
[ b^2*e^2*x2^2 - 2*b*e^2*x2^3 + e^2*x2^4 - b^2*d*e*x2*y2 + b*d^2*e*x2*y2 + b*e^3*x2*y2 + b*d*e*x2^2*y2 - d^2*e*x2^2*y2 - e^3*x2^2*y2 - 2*b*e^2*x2*y2^2 + 2*e^2*x2^2*y2^2 + b*d*e*y2^3 - d^2*e*y2^3 - e^3*y2^3 + e^2*y2^4 + b^2*e*x2*y2 - b*d*e*x2*y2 - b*e*x2^2*y2 + d*e*x2^2*y2 + b*e^2*y2^2 - b*e*y2^3 + d*e*y2^3]
[ b^2*e^2*x2*y2 - b*e^2*x2^2*y2 - b^2*d*e*y2^2 + b*d^2*e*y2^2 + b*e^3*y2^2 - b*e^2*y2^3 - b*e^2*x2*y2 + e^2*x2^2*y2 + b^2*e*y2^2 - b*d*e*y2^2 + e^2*y2^3]
[ b^2*e^2*x2^2 - 2*b*e^2*x2^3 + e^2*x2^4 + b^2*e^2*y2^2 - 2*b*e^2*x2*y2^2 + 2*e^2*x2^2*y2^2 + e^2*y2^4]
[b^2*e^2*x2^2 - 2*b*e^2*x2^3 + e^2*x2^4 - 2*b^2*d*e*x2*y2 + 2*b*d^2*e*x2*y2 + 2*b*e^3*x2*y2 + 2*b*d*e*x2^2*y2 - 2*d^2*e*x2^2*y2 - 2*e^3*x2^2*y2 + b^2*d^2*y2^2 - 2*b*d^3*y2^2 + d^4*y2^2 - 2*b*d*e^2*y2^2 + 2*d^2*e^2*y2^2 + e^4*y2^2 - 2*b*e^2*x2*y2^2 + 2*e^2*x2^2*y2^2 + 2*b*d*e*y2^3 - 2*d^2*e*y2^3 - 2*e^3*y2^3 + e^2*y2^4 + 2*b^2*e*x2*y2 - 2*b*d*e*x2*y2 - 2*b*e*x2^2*y2 + 2*d*e*x2^2*y2 - 2*b^2*d*y2^2 + 4*b*d^2*y2^2 - 2*d^3*y2^2 + 2*b*e^2*y2^2 - 2*d*e^2*y2^2 - 2*b*e*y2^3 + 2*d*e*y2^3 + b^2*y2^2 - 2*b*d*y2^2 + d^2*y2^2 + e^2*y2^2]
C5.column()
[b*c^2*e^2*x1^2 - 2*b*c*e^2*x1^3 + b*e^2*x1^4 - b^2*c*d*e*x1*y1 - b*c^2*d*e*x1*y1 + b*c*d^2*e*x1*y1 + c^2*d^2*e*x1*y1 + b*c*e^3*x1*y1 + c^2*e^3*x1*y1 + b^2*d*e*x1^2*y1 + b*c*d*e*x1^2*y1 - b*d^2*e*x1^2*y1 - c*d^2*e*x1^2*y1 - b*e^3*x1^2*y1 - c*e^3*x1^2*y1 + b^2*c*d^2*y1^2 - 2*b*c*d^3*y1^2 + c*d^4*y1^2 - 2*b*c*d*e^2*y1^2 + 2*c*d^2*e^2*y1^2 + c*e^4*y1^2 - 2*b*c*e^2*x1*y1^2 + 2*b*e^2*x1^2*y1^2 + b^2*d*e*y1^3 + b*c*d*e*y1^3 - b*d^2*e*y1^3 - c*d^2*e*y1^3 - b*e^3*y1^3 - c*e^3*y1^3 + b*e^2*y1^4 - 2*b*c^2*e^2*x1 + 4*b*c*e^2*x1^2 - 2*b*e^2*x1^3 + b^2*c*d*e*y1 + b*c^2*d*e*y1 - b*c*d^2*e*y1 - c^2*d^2*e*y1 - b*c*e^3*y1 - c^2*e^3*y1 + b^2*c*e*x1*y1 + b*c^2*e*x1*y1 - b^2*d*e*x1*y1 - 2*b*c*d*e*x1*y1 - c^2*d*e*x1*y1 + b*d^2*e*x1*y1 + c*d^2*e*x1*y1 + b*e^3*x1*y1 + c*e^3*x1*y1 - b^2*e*x1^2*y1 - b*c*e*x1^2*y1 + b*d*e*x1^2*y1 + c*d*e*x1^2*y1 - 2*b^2*c*d*y1^2 + 4*b*c*d^2*y1^2 - 2*c*d^3*y1^2 + b^2*e^2*y1^2 + 2*b*c*e^2*y1^2 + c^2*e^2*y1^2 - 2*c*d*e^2*y1^2 - 2*b*e^2*x1*y1^2 - b^2*e*y1^3 - b*c*e*y1^3 + b*d*e*y1^3 + c*d*e*y1^3 + b*c^2*e^2 - 2*b*c*e^2*x1 + b*e^2*x1^2 - b^2*c*e*y1 - b*c^2*e*y1 + b*c*d*e*y1 + c^2*d*e*y1 + b^2*e*x1*y1 + b*c*e*x1*y1 - b*d*e*x1*y1 - c*d*e*x1*y1 + b^2*c*y1^2 - 2*b*c*d*y1^2 + c*d^2*y1^2]
[ -b^2*c*e^2*x1*y1 + b*c^2*e^2*x1*y1 + b^2*e^2*x1^2*y1 - b*c*e^2*x1^2*y1 + b^2*c*d*e*y1^2 - b*c^2*d*e*y1^2 - b*c*d^2*e*y1^2 + c^2*d^2*e*y1^2 - b*c*e^3*y1^2 + c^2*e^3*y1^2 + b^2*e^2*y1^3 - b*c*e^2*y1^3 + b^2*c*e^2*y1 - b*c^2*e^2*y1 - b^2*e^2*x1*y1 + 2*b*c*e^2*x1*y1 - c^2*e^2*x1*y1 - b*e^2*x1^2*y1 + c*e^2*x1^2*y1 - b^2*c*e*y1^2 + b*c^2*e*y1^2 - b^2*d*e*y1^2 + 2*b*c*d*e*y1^2 - c^2*d*e*y1^2 + b*d^2*e*y1^2 - c*d^2*e*y1^2 + b*e^3*y1^2 - c*e^3*y1^2 - b*e^2*y1^3 + c*e^2*y1^3 - b*c*e^2*y1 + c^2*e^2*y1 + b*e^2*x1*y1 - c*e^2*x1*y1 + b^2*e*y1^2 - b*c*e*y1^2 - b*d*e*y1^2 + c*d*e*y1^2]
[ b^2*c^2*e^2*x1^2 - 2*b^2*c*e^2*x1^3 + b^2*e^2*x1^4 - 2*b^2*c^2*d*e*x1*y1 + 2*b*c^2*d^2*e*x1*y1 + 2*b*c^2*e^3*x1*y1 + 2*b^2*c*d*e*x1^2*y1 - 2*b*c*d^2*e*x1^2*y1 - 2*b*c*e^3*x1^2*y1 + b^2*c^2*d^2*y1^2 - 2*b*c^2*d^3*y1^2 + c^2*d^4*y1^2 - 2*b*c^2*d*e^2*y1^2 + 2*c^2*d^2*e^2*y1^2 + c^2*e^4*y1^2 - 2*b^2*c*e^2*x1*y1^2 + 2*b^2*e^2*x1^2*y1^2 + 2*b^2*c*d*e*y1^3 - 2*b*c*d^2*e*y1^3 - 2*b*c*e^3*y1^3 + b^2*e^2*y1^4 - 2*b^2*c^2*e^2*x1 + 4*b^2*c*e^2*x1^2 - 2*b^2*e^2*x1^3 + 2*b^2*c^2*d*e*y1 - 2*b*c^2*d^2*e*y1 - 2*b*c^2*e^3*y1 + 2*b^2*c^2*e*x1*y1 - 2*b^2*c*d*e*x1*y1 - 2*b*c^2*d*e*x1*y1 + 2*b*c*d^2*e*x1*y1 + 2*b*c*e^3*x1*y1 - 2*b^2*c*e*x1^2*y1 + 2*b*c*d*e*x1^2*y1 - 2*b^2*c^2*d*y1^2 + 4*b*c^2*d^2*y1^2 - 2*c^2*d^3*y1^2 + 2*b^2*c*e^2*y1^2 + 2*b*c^2*e^2*y1^2 - 2*c^2*d*e^2*y1^2 - 2*b^2*e^2*x1*y1^2 - 2*b^2*c*e*y1^3 + 2*b*c*d*e*y1^3 + b^2*c^2*e^2 - 2*b^2*c*e^2*x1 + b^2*e^2*x1^2 - 2*b^2*c^2*e*y1 + 2*b*c^2*d*e*y1 + 2*b^2*c*e*x1*y1 - 2*b*c*d*e*x1*y1 + b^2*c^2*y1^2 - 2*b*c^2*d*y1^2 + c^2*d^2*y1^2 + b^2*e^2*y1^2 - 2*b*c*e^2*y1^2 + c^2*e^2*y1^2]
[ c^2*e^2*x1^2 - 2*c*e^2*x1^3 + e^2*x1^4 - 2*b*c*d*e*x1*y1 + 2*c*d^2*e*x1*y1 + 2*c*e^3*x1*y1 + 2*b*d*e*x1^2*y1 - 2*d^2*e*x1^2*y1 - 2*e^3*x1^2*y1 + b^2*d^2*y1^2 - 2*b*d^3*y1^2 + d^4*y1^2 + b^2*e^2*y1^2 - 2*b*c*e^2*y1^2 + c^2*e^2*y1^2 - 2*b*d*e^2*y1^2 + 2*d^2*e^2*y1^2 + e^4*y1^2 - 2*c*e^2*x1*y1^2 + 2*e^2*x1^2*y1^2 + 2*b*d*e*y1^3 - 2*d^2*e*y1^3 - 2*e^3*y1^3 + e^2*y1^4 - 2*c^2*e^2*x1 + 4*c*e^2*x1^2 - 2*e^2*x1^3 + 2*b*c*d*e*y1 - 2*c*d^2*e*y1 - 2*c*e^3*y1 + 2*b*c*e*x1*y1 - 2*b*d*e*x1*y1 - 2*c*d*e*x1*y1 + 2*d^2*e*x1*y1 + 2*e^3*x1*y1 - 2*b*e*x1^2*y1 + 2*d*e*x1^2*y1 - 2*b^2*d*y1^2 + 4*b*d^2*y1^2 - 2*d^3*y1^2 + 2*b*e^2*y1^2 + 2*c*e^2*y1^2 - 2*d*e^2*y1^2 - 2*e^2*x1*y1^2 - 2*b*e*y1^3 + 2*d*e*y1^3 + c^2*e^2 - 2*c*e^2*x1 + e^2*x1^2 - 2*b*c*e*y1 + 2*c*d*e*y1 + 2*b*e*x1*y1 - 2*d*e*x1*y1 + b^2*y1^2 - 2*b*d*y1^2 + d^2*y1^2]
C6.column()
[ c^2*x1*y1*y2 - c*x1^2*y1*y2 - c*y1^3*y2 - c*x1*y1*x2*y2 + x1^2*y1*x2*y2 + y1^3*x2*y2 + c*y1^2*y2^2 + c^2*y1^2 - 2*c*y1^2*x2 + y1^2*x2^2 - c^2*y1*y2 + c*x1*y1*y2 + c*y1*x2*y2 - x1*y1*x2*y2]
[ c^2*y1^2*y2 - c*y1^2*x2*y2 - c*x1*y1*y2^2 + x1^2*y1*y2^2 + y1^3*y2^2 - c*y1^2*y2 + y1^2*x2*y2 + c*y1*y2^2 - x1*y1*y2^2]
[c^2*x1^2*y2^2 - 2*c*x1^3*y2^2 + x1^4*y2^2 + c^2*y1^2*y2^2 - 2*c*x1*y1^2*y2^2 + 2*x1^2*y1^2*y2^2 + y1^4*y2^2 + 2*c^2*x1*y1*y2 - 2*c*x1^2*y1*y2 - 2*c*y1^3*y2 - 2*c*x1*y1*x2*y2 + 2*x1^2*y1*x2*y2 + 2*y1^3*x2*y2 - 2*c^2*x1*y2^2 + 4*c*x1^2*y2^2 - 2*x1^3*y2^2 + 2*c*y1^2*y2^2 - 2*x1*y1^2*y2^2 + c^2*y1^2 - 2*c*y1^2*x2 + y1^2*x2^2 - 2*c^2*y1*y2 + 2*c*x1*y1*y2 + 2*c*y1*x2*y2 - 2*x1*y1*x2*y2 + c^2*y2^2 - 2*c*x1*y2^2 + x1^2*y2^2]
[ c^2*y1^2 - 2*c*y1^2*x2 + y1^2*x2^2 + y1^2*y2^2]
Another figure
A user called Futurologist wrote a different answer, using a lot less brute-force computation. Here is an interactive version of the figure for that proof. The crucial new elements are Pascal’s line (black) and the orange circles with their common radical axis (blue). The point is the center of the black circle through all the , while the point is the center of the blue circle through all the .
You can drag both of these around, and also drag the and . The radius of the black circle can be modified, too. So there are 13 real degrees of freedom in this construction. My computation above used 9 variable parameters, but it also fixed the positions of three points (amounting to 6 real degrees of freedom) while at the same time imposing two extra conditions (removing two real degrees of freedom).