Find point on rational circle for given angle

Question

I've a circle defined by 3 rationals points on the circle. I need to calculate another rational point on the circle given by an angle (a floating point number). The resulting angle does not need to be exact, but within certain limits (so the rational point should be close to the real point). It seems the default rational parametrisation of a circle $x=\frac{2t}{t^2+1}$ and $y=\frac{t^2-1}{t^2+1}$ is not helpful as it requires $t\to\infty$ for points close to the angle $\frac{\pi}{2}$.

How can I reliably find such rationals points on a rational circle?

Answer 1

Here is one possible approach. You can model the plane plus one extra point at infinity as the complex projective line $\mathbb{CP}^1$. So your three points become homogeneous coordinate vectors of the form $[x+iy:1]$. It is simple (see below for details) to find a projective transformation (i.e. Möbius transformation) which maps your circle to the real line in such a way that the three points you have correspond to $0$, $1$ and $\infty$. So as long as you can guarantee that the point you're looking for will be reasonably far away from one of your three known points, mapping that point to $\infty$ will avoid using infinite values for the point you're looking for. Rational points on the circle would simply correspond to rational numbers on the real number line.

You can then apply this transformation to a floating-point approximation of your desired point, and dehomogenize the result (i.e. divide first coordinate by second) to a real number (up to some imaginary noise due to numerics). That number is in fact the cross ratio of your desired point with respect to the three given points. In the special case where the three known points are west, east and north (or left, right and up) compass points, this cross ratio corresponds to the tangent half-angle parametrization you mentioned. So the cross-ratio parametrization may be seen as a generalization of the tangent half-angle parametrization, using the three given points instead of fixed compass points as the basis. You could approximate that floating-point number using rational numbers, e.g. using continued fractions to come up with a good approximation in terms of error compared to size of involved numbers. Then do the inverse Möbius transformation and you are back to your original circle. This may make your rational numbers more complicated, so I don't claim that a compact cross ratio approximation will lead to compact coordinates.

So how do you find the transformation? It's essentially a special case of this approach. Write your three points as complex numbers $A=x_A+i\,y_A$ and likewise $B$ and $C$. Now try to solve

$$\begin{pmatrix}A&B\\1&1\end{pmatrix}\cdot\begin{pmatrix}\lambda\\\mu\end{pmatrix}=\tau\begin{pmatrix}C\\1\end{pmatrix}$$

for some $\tau\in\mathbb C$. You can do this using the adjoint matrix, with $\tau$ equal to the determinant $(A-B)$:

$$\begin{pmatrix}\lambda\\\mu\end{pmatrix}= \begin{pmatrix}1&-B\\-1&A\end{pmatrix} \begin{pmatrix}C\\1\end{pmatrix}$$

Then scale the columns of the original matrix with these values to obtain

$$\begin{pmatrix}\lambda A&\mu B\\\lambda&\mu\end{pmatrix}$$

as the matrix which maps $\infty,0,1$ a.k.a. $[1:0],[0:1],[1:1]$ to $A,B,C$. Its inverse operation can be expressed using the adjoint again:

$$\begin{pmatrix}\mu&-\mu B\\-\lambda&\lambda A\end{pmatrix}$$

is the transformation from your circle to the real line, which maps $A,B,C$ to $\infty,0,1$ so make sure the point you are after is distinct from $A$ if you want to avoid infinity.

Example

To make all of this more accessible, let's have an exampe. Consider the circle through $(1,5)$, $(4,1)$ and $(5,4)$. Just for reference: the radius of the circle is $\frac5{26}\sqrt{170}\not\in\mathbb Q$, and its center is at $\left(\frac{69}{26},\frac{81}{26}\right)$.

So you have $A=1+5i, B=4+1i, C=5+4i$. Using these you compute

\begin{align*} \lambda&=C-B=(5+4i)-(4+1i)=1+3i \\ \mu&=A-C=(1+5i)-(5+4i)=-4+1i \\ \lambda A&=(1+3i)(1+5i)=(1-3\cdot 5)+(3+5)i=-14+8i \\ \mu B&=(-4+1i)(4+1i)=(-4\cdot4-1)+(-4+4)i=-17+0i \end{align*}

Now suppose you want a point at $\varphi=123°\pm1''$. So we expect the point to lie close to $(1.2883,5.2183)$. Let's convert this to a real parameter:

\begin{multline*} \begin{pmatrix}\mu&-\mu B\\-\lambda&\lambda A\end{pmatrix}\cdot \begin{pmatrix}x+iy\\1\end{pmatrix}= \begin{pmatrix}-4+1i&17\\-1-3i&-14+8i\end{pmatrix}\cdot \begin{pmatrix}1.2883+5.2183i\\1\end{pmatrix} \\=\begin{pmatrix} (-4+1i)(1.2883+5.2183i)+17\\ (-1-3i)(1.2883+5.2183i)-14+8i \end{pmatrix} \\=\begin{pmatrix} (-4\cdot1.2883-1\cdot5.2183+17)+ (-4\cdot5.2183+1\cdot1.2883)i\\ (-1\cdot1.2883+3\cdot5.2183-14)+ (-1\cdot5.2183-3\cdot1.2883+8)i \end{pmatrix} \\=\begin{pmatrix} 6.6285-19.5849i\\ 0.3666-1.0832i \end{pmatrix} \end{multline*}

This is a homogeneous parameter vector. To turn this into a single real number, divide the first coordinate by the second.

\begin{multline*} \frac{6.6285-19.5849i}{0.3666-1.0832i} = \frac{(6.6285-19.5849i)(0.3666+1.0832i)}{(0.3666-1.0832i)(0.3666+1.0832i)} \\=\frac{(6.6285\cdot0.3666+19.5849\cdot1.0832)+ (6.6285\cdot1.0832-19.5849\cdot0.3666)i}{0.3666^2+1.0832^2} \\=\frac{23.64437178+0.00016686i}{1.3077178}\approx18.0806+0.0001i \end{multline*}

The imaginary part here is because our approximate position wasn't exactly on the circle; we simply ignore that. So our parameter $t$ for this point is approximately $18.0806$ according to the above computation. But that computation was pretty imprecise (to keep the numbers readable). In fact it was too imprecise to satisfy the tight error bounds I assumed. If we do the same computation with higher precision, and with the bounds of the desired interval, we obtain

$$18.084342<t<18.085554$$

as a sufficient condition (i.e. if $t$ lies in that interval, the error is guaranteed to be smaller than one arc second). We can also look at the continued fraction representation of the parameter you get for exactly $123°$, again computed with higher precision.

$$t\approx18.0849478026=18+\cfrac{1}{11+\cfrac{1}{1+\cfrac{1}{3+\cfrac{1}{2+\cfrac{1}{1+\cfrac{1}{1+\cdots}}}}}}$$

Trying one convergent (i.e. truncated version) after the other, one can verify that the first approximation to lie within the desired interval is

$$t=18+\cfrac{1}{11+\cfrac{1}{1+\cfrac{1}{3}}}=\frac{850}{47} \approx18.0851$$

Now we turn this back into a point on the circle. To keep the notation simple and avoid fractions, I'll write numerator and denominator of the computed number as two components of the homogeneous coordinate vector. Scalar multiples of homogeneous vectors describe the same point.

\begin{gather*} \begin{pmatrix}\lambda A&\mu B\\\lambda&\mu\end{pmatrix}\cdot 47\begin{pmatrix}t\\1\end{pmatrix}= \begin{pmatrix}-14+8i&-17\\1+3i&-4+1i\end{pmatrix}\cdot \begin{pmatrix}850\\47\end{pmatrix}= \begin{pmatrix}-12699+6800i\\662+2597i\end{pmatrix}\\[2ex] \frac{-12699+6800i}{662+2597i}= \frac{(-12699+6800i)(662-2597i)}{662^2+2597^2}= %\frac{9760822+37208903i}{7182653}= \frac{544286}{422509}+\frac{2204759}{422509}i \end{gather*}

So one rational point satisfying all the requirements would be

$$\left(\frac{544286}{422509},\frac{2204759}{422509}\right)$$

Answer 2

It seems the default rational parametrisation of a circle $x=\frac{2t}{t^2+1}$ and $y=\frac{t^2-1}{t^2+1}$ is not helpful as it requires $t\to\infty$ for points close to the angle $\frac{\pi}{2}$.

One elegant way of dealing with this problem is by using a homogeneous approach: introduce a second parameter $u$ like this:

$$ x = \frac{2tu}{t^2+u^2} \qquad y = \frac{t^2-u^2}{t^2+u^2} $$

The old $t$ corresponds to $\frac tu$ in the new parameters. The situation of $t\to\infty$ can now be described by $t=1,u=0$ which represents division by zero in a well-defined way.

Suppose you have a point $(x,y)$ on the unit circle, i.e. with $x^2+y^2=1$. This is what you'd get from having an angle, with $x=\cos\varphi$ and $y=\sin\varphi$. Then you can choose $t=x,u=1-y$. To verify this claim:

$$\frac{2x(1-y)}{x^2+(1-y)^2} = \frac{2x(1-y)}{(x^2+y^2)+1-2y} = \frac{2x(1-y)}{2-2y} = x \\ \frac{x^2-(1-y)^2}{x^2+(1-y)^2} = \frac{x^2-1+2y-y^2}{(x^2+y^2)+1-2y} = \frac{x^2-(x^2+y^2)+2y-y^2}{2-2y} = \frac{2y-2y^2}{2-2y} = y$$

As usual for homogeneous coordinates, multiples of the parameter vector $(t,u)$ describe the same point. So if $\frac tu\in\mathbb Q$ (rational) you can even choose $t,u\in\mathbb Z$ (integers). You might compute $t$ and $u$ from an angle $\varphi$ as described above, to obtain a first floating-point approximation. Then you can check which of these is larger, and try to approximate either $\frac tu$ or $\frac ut$ by a rational number, e.g. using continued fraction approximations. The numerator and denominator of the resulting fraction can then be used as parameters $t$ and $u$ for the formula above.

Personally I'd usually use $x=\frac{u^2-t^2}{u^2+t^2}$ and $y=\frac{2tu}{u^2+t^2}$, so that $t=0,u=1$ represents $0°$, $t=1,u=1$ is $90°$, $t=-1,u=1$ is $-90°$ and $t=1,u=0$ is $\pm180°$. That way you have $\frac tu=\tan\frac\varphi2$, which gives rise to the term tangent half-angle substitution. In that case you'd choose $t=y$ and $u=1+x$. It's just a minor change of coordinate system, though.

As noted in the comment, the above description works for the unit circle, and by a simple scaling also works for circles with rational radius. But circles defined by three rational points may well have an irrational radius. To deal with that fact, you could establish a different coordinate system. The center of the circle (which has rational coordinates) would be the origin of the new coordinate system, and the vector from there to one of the defining points could serve as the first basis vector. The second basis vector would be perpendicular to that.

So if the center of the circle is at $(x_0,y_0)$ and the first of the three defining points is at $(x_1,y_1)$ you could use a parametrization like

\begin{align*} x &= x_0 + \frac{2tu(x_1-x_0)-(t^2-u^2)(y_1-y_0)}{t^2+u^2} \\ y &= y_0 + \frac{2tu(y_1-y_0)+(t^2-u^2)(x_1-x_0)}{t^2+u^2} \end{align*}

When reading the parameters $t$ and $u$ off the angle $\varphi$, you'd first have to subtract the angle of the vector from center to first defining point, as that is your first coordinate axis in the new coordinate system.