Locating three sets of collinear points

The operation that takes a point $P$ on the circle and find the other intersection point of the line $(AP)$ with the circle is a real automorphism of the circle.

The composition of $3$ such automorphisms is still a real automorphism. If you identify the circle with $\Bbb P^1(\Bbb R)$, it corresponds to a real homography $t \mapsto \frac {at+b}{ct+d}$.

In theory, given $3$ points on the circle and their images you can determine the parameters $a/d,b/d,c/d,\ldots \in \Bbb P^1(\Bbb R)$ and from then construct the line going through the two fixpoints (the parameters of the line are rational fractions in $a,b,c,d$ so they're a rational fraction in the coordinates of the three image points).

There probably even is a way to get a construction that doesn't depend on which $3$ points you choose on the circle, but at the moment I haven't checked on how to construct that line.

Given a point $A \in \Bbb P^2(\Bbb R)$, call $\sigma_A$ the involution of the circle obtained by taking a point $M$ and intersecting the line $(MA)$ with the circle.

Say an automorphism of the circle is direct if the output points turn in the same direction as the input point and indirect if it's not the case (this given by the sign of the determinant $ad-bc$).

The $\sigma_A$ are indirect involutions so they don't represent all the indirect automorphisms (a space of dimension $3$ while we have only $2$ dimensions by picking $A$). However, any direct automorphism can be given in infinitely many ways as a composition of two $\sigma_{A_i}$.

The fixpoints (possibly over $\Bbb C$) of $\sigma_A \circ \sigma_B$ are clearly the (possibly complex) intersections of $(AB)$ with the circle, so for any $C \in (AB)$ there is a unique point $D$ such that $\sigma_A \circ \sigma_B = \sigma_C \circ \sigma_D$. Constructing $D$ or $C$ from the other one is really easy, for example $C$ is the intersection of $(AB)$ with the line $(\sigma_D(M))( \sigma_A \circ \sigma_B (M))$, for any $M$ on the circle.

Next, we clearly have the simplification $\sigma_A \circ \sigma_A = id$.

So given $4$ points we can simplify $(\sigma_A \circ \sigma_B) \circ (\sigma_C \circ \sigma_D)$ by looking at the intersection of $(AB)$ and $(CD)$ (which is possibly at infinity) and moving both $B$ and $C$ to that point while changing $A$ and $D$ appropriately to $A'$ and $D'$. This gives us $(\sigma_A \circ \sigma_B) \circ (\sigma_C \circ \sigma_D) = \sigma_{A'} \circ \sigma_{D'}$

So now that we know how to represent direct automorphisms of the circle with couples of points and compute compositions, and because there is a direct relation between the fixpoints and the line formed by the points, all we have to do is to compute $\sigma_C \circ \sigma_B \circ \sigma_A \circ \sigma_C \circ \sigma_B \circ \sigma_A$ and simplify it to some $\sigma_X \circ \sigma_Y$. Then $(XY)$ intersects the circle at the two fixpoints $D_1$ and $D_2$ of $\sigma_C \circ \sigma_B \circ \sigma_A$.

Note that this construction should work if you replace the circle with any nondegenerate conic, so ellipses, parabolas, and hyperbolas too. You can also generalize this to an odd number $n$ of points instead of $3$, though the complexity of the construction increases linearly with $n$.

This is NOT an answer, but showing why there are two solutions.

Let $D$ be a point on the circle and let $E, F$ be the other intersections of lines $AD, BD$, resp., with the circle. When $D$ is rotated around the circle the points $E, F$ rotate around the circle in opposite direction. (The first image only shows $A, D$ and $E$.)

Hence also the line $EF$ rotates around the circle center in the opposite direction. For any point $C$ outside the circle the line $EF$ will run through $C$ twice during one full rotation: once while $E$ is closer to $C$, once while $F$ is closer to $C$.

[This is just a more step by step recipe based on mercio's solution for reference. It is not meant to take the bounty.]

Basic outline:

1. construct $B_1$ s.t. $\sigma_B \circ \sigma_A = \sigma_A \circ \sigma_{B_1}$
2. construct $C_1$ s.t. $\sigma_A \circ \sigma_C = \sigma_{C_1} \circ \sigma_{A}$
3. intersect line $BC$ with line $B_1C_1$ to get $S$
4. construct $X$ s.t. $\sigma_{C_1} \circ \sigma_{B_1} = \sigma_{S} \circ \sigma_{X}$
5. construct $Y$ s.t. $\sigma_{C} \circ \sigma_{B} = \sigma_{Y} \circ \sigma_{S}$
6. intersect line $XY$ with the circle to get $D_1$ and $D_2$, the two possible solutions for $D$.

After step 2 we have the points $B_1$ and $C_1$ s.t. $$\sigma_A \circ \sigma_C \circ \sigma_B \circ \sigma_A = \sigma_{C_1} \circ \sigma_A \circ \sigma_A \circ \sigma_{B_1} = \sigma_{C_1} \circ \sigma_{B_1} .$$ After step 5 we have the points $X$ and $Y$ s.t. $$ \sigma_{C} \circ \sigma_{B} \circ \sigma_{C_1} \circ \sigma_{B_1} = \sigma_{Y} \circ \sigma_{S} \circ \sigma_{S} \circ \sigma_{X} = \sigma_{Y} \circ \sigma_{X} .$$ Hence $$\sigma_C \circ \sigma_B \circ \sigma_A \circ \sigma_C \circ \sigma_B \circ \sigma_A = \sigma_{C} \circ \sigma_{B} \circ \sigma_{C_1} \circ \sigma_{B_1} = \sigma_{Y} \circ \sigma_{X} .$$ The fixed points of $\sigma_C \circ \sigma_B \circ \sigma_A$ are then also the fixed points of $\sigma_{Y} \circ \sigma_{X}$.

The details for construction steps 1, 2, 4, 5 are already given by mercio. I just give the details for step 1 for completeness:

1.1 choose any point $M$ on the circle and intersect line $AM$ with the circle to get the second intersection point $N = \sigma_A(M)$

1.2 intersect line $BN$ with the circle to get the second intersection point $P = \sigma_B(N) = \sigma_B(\sigma_A(M)) = (\sigma_B \circ \sigma_A)(M)$

1.3 intersect line $AP$ with the circle to get the second intersection point $Q = \sigma_A(P) = \sigma_A^{-1}(P)$

1.4 intersect line $QM$ with line $AB$ to get $B_1$, hence $\sigma_{B_1}(M)=Q=\sigma_A^{-1}(P)$ and $P = (\sigma_B \circ \sigma_A)(M) = (\sigma_{A} \circ \sigma_{B_1})(M)$.

[this is an addition to mercio's first answer: for my understanding there was still a proof missing, why the involved involutions on the circle like $\sigma_A$ (see mercio's answer or below for the precise definition) can be represented as a function of the form $t \mapsto \frac{at+b}{ct+d}$ in the projective line; so I want to share my results on this; I am not aware of any shorter argument]

A circle can be parametrized using stereographic projection. For our needs here it is easiest to consider the unit circle with center $M=(0,1)$:

The real parameter $t$ parametrizes the circle $x^2 + (y-1)^2 = 1$ by $$t \mapsto U(t) := \left(\frac{2t}{t^2+1},\frac{2}{t^2+1}\right) .$$ For a given point $(x,y)$ on the circle we get back $t$ by $t = \frac{x}{y}$. In this parametrization only the point $P$ is missing. The parametrization can be completed by going to homogenous coordinates $[t:1]$ and assigning $[1:0]$ to $P$.

Now let us consider the special situation when $A$ is on the $y$-axis. The involution $\sigma_A$ maps a point $T_0 = U(t_0)$ on the circle to the other intersection of the line $AT$ with the circle $T_1 = U(t_1)$.

Now how are $t_0$ and $t_1$ related? It turns out that $$t_1 = \frac{y_0-2}{y_0t_0} ,$$ where $y_0$ is the $y$-coordinate of $A$. (It takes some algebra to compute the second intersection of the circle $x^2 + (y-1)^2 = 1$ with the line through $A$ and $T_0$, whose equation is $$ y = \frac{\frac{2}{t_0^2+1}-y_0}{\frac{2t_0}{t_0^2+1}}x+y_0 .)$$ So in the special case of $A$ on the $y$-axis the involution $\sigma_A$ can indeed be represented in the form $\frac{0\cdot t+(y_0-2)}{y_0t+0}=\frac{at+b}{ct+d}$.

For the general case of the position of point $A$ it suffices to check that rotating a point on the circle corresponds to a transformation of parameter $t$ of the same form.

Now here the homogenous coordinates come in handy. From the above figure it can be seen that $$ U([t_1:h]^T) = \left(\begin{array}{cc} \cos\alpha & -\sin\alpha \\ \sin\alpha & \cos\alpha \end{array}\right) U([t_0:h]^T) = U\left( \left(\begin{array}{cc} \cos\frac{\alpha}{2} & -\sin\frac{\alpha}{2} \\ \sin\frac{\alpha}{2} & \cos\frac{\alpha}{2} \end{array}\right) [t_0:h]^T \right) ,$$ which shows that the transformation of $t_0$ to $t_1$ again is of the expected form. Since also concatentions of functions of the form $t \mapsto \frac{at+b}{ct+d}$ are of this form the - at least for me - missing details seem to be filled.