Hi everyone!
Recently I started solving projecteuler, and while doing so I encountered two concepts about which I heard before, but I didn't really bother to learn them. I made a few notes to myself about how they work, and thought it could be useful for somebody else too. This blog focuses on Pythagorean triples and Pell's equations, which are recurrent concepts on projecteuler.
Great thanks to nor, Endagorion, Golovanov399 and Neodym for useful discussions about these topics.
Pythagorean triples
Pythagorean triple is a triple $$$a,b,c \in \mathbb Z$$$ such that
It is possible to parameterize them in a way that allows to find all triples such that $$$a,b,c \leq n$$$ in $$$O(\sqrt n)$$$.
Rational points on a unit circle
To do that, let's divide the equation by $$$c^2$$$ to get
This reduces the task of finding Pythagorean triples to finding rational points $$$(x, y)$$$ on a unit circle.
And there is a nice and simple way to find all such points. To do this, we should notice that a ray going from $$$(0, -1)$$$ into a rational point $$$(x, y)$$$ will intersect the line $$$y=0$$$ in a rational point $$$(\frac{p}{q}, 0)$$$. So, we can find all rational points $$$(x, y)$$$ by projecting the point $$$(\frac{p}{q}, 0)$$$ back onto unit circle for each $$$\frac{p}{q}$$$.
We consider a ray from $$$(0, -1)$$$ into $$$(\frac{p}{q}, 0)$$$ to find rational points on the unit circle
Note: The correct value for the numerator of $$$y$$$ in the image is $$$q^2 - p^2$$$
The direction vector of a ray from $$$(0, -1)$$$ to $$$(\frac{p}{q}, 0)$$$ is $$$(\frac{p}{q},1)$$$, so the equation of the ray is $$$py - xq = -p$$$. Let's solve the system
After expressing $$$y$$$ as a function of $$$x$$$ like
we put it back into the first equation and get the result:
As it turns out, the projection point on a unit circle is also always rational, so we found a bijection between rational points in $$$\mathbb R$$$ and rational points on the unit circle.
As a geometric inversion
It was also pointed to me that the bijection between the Pythagorean triples and the rational numbers can be perceived as an inversion transform on a complex plane. General formula for the inversion in the point $$$O$$$ with the radius $$$r$$$ is
where $$$\overline{z-O}$$$ is complex conjugate. Therefore,
In our particular case we use $$$O = -i$$$ and $$$r = \sqrt 2$$$. Then the image of a real number $$$t$$$ is
which rewrites as
Then it follows from geometric inversion properties that all such points lie on the unit circle. And for $$$t = \frac{p}{q}$$$ it is
Back to integers
This translates back into integer triples as
Note that for integers such parameterization is incomplete, as the fraction $$$x=\frac{a}{c}$$$ could as well correspond to $$$\frac{ka}{kc}$$$ for some integer $$$k$$$. Therefore, all possible triplets are obtained by multiplying vectors $$$(a,b,c)$$$ of the form defined above with all possible integer constants $$$k$$$.
Pell's equations
Another somewhat recurring theme is Pell's equations. These are the equations that look like
We're looking for positive integer solutions $$$(x, y)$$$ to it. To understand the solutions better, we can relax it to inequalities
The two inequalities above are equivalent to the equation when $$$(x, y)$$$ are restricted to integers. Now, by dividing with $$$y^2$$$ we turn it into
The only possible lattice points satisfying the inequalities above are actual solutions to $$$x^2 - Dy^2 = 1$$$. It means that the set of lattice points in the quater-plane $$$x, y \geq 0$$$, for which $$$x > y \sqrt D$$$, is a subset of the area surrounded by the curve $$$x^2 - Dy^2 = 1$$$. As this area is strictly convex (see the picture below), the subset can only intersect the boundary of the area in the convex hull of the subset. In other words, all solutions must lie on the convex hull of lattice points $$$x > y \sqrt D$$$ on $$$x, y \geq 0$$$.
Note that being on the convex hull of lattice points is a necessary, but generally not sufficient, condition for a solution.
Convergents (red) correspond to corners of the convex hull of lattice points
From the theory of continued fractions it is known that such points are exactly the convergents $$$\frac{x}{y}$$$ of $$$\sqrt D$$$. Knowing this, and trying all convergents in order we will eventually find at least one non-trivial solution, from which it will be possible to find all the others (see below).
Solutions as a subgroup of $$$2 \times 2$$$ matrices with determinant $$$1$$$
The text above shows that if solution exists, it must be a convergent of $$$\sqrt D$$$. But still, we would like to know which particular convergents are the solutions and what is their general structure. To better understand it, let's rewrite the equation as
That's a very nice representation because of the following facts:
- If an integer matrix has determinant $$$1$$$, its inverse is also an integer matrix with determinant $$$1$$$.
- If two matrices have determinant $$$1$$$, so does their product.
So, let's consider the inverse of the matrix under the determinant. It is
corresponding to the fact that if $$$(x, y)$$$ is an integer solution to the Pell's equation, so is $$$(x, -y)$$$. Now to the product:
In other words, if $$$(x_0, y_0)$$$ and $$$(x_1,y_1)$$$ are solutions, so is $$$(x_0 x_1 + D y_0 y_1, x_0 y_1 + y_0 x_1)$$$. This means that the full set of integer solutions to the Pell's equations in fact form a group with the group operation defined as above. Note that the point $$$(x, y) = (1, 0)$$$ corresponds to the identity element of such group. Note that for $$$(x_0, y_0) = (x_1, y_1)$$$ the result is $$$(x^2+Dy^2, 2xy)$$$.
Solutions as a subgroup of $$$\mathbb Z[\sqrt D]$$$
Another, perhaps simpler way to look on it is to consider numbers of kind $$$a + b \sqrt D$$$. For such numbers, it's natural to define multiplication
and the multiplicative norm $$$|a + b \sqrt D | = (a + b \sqrt D)(a - b \sqrt D) = a^2 - b^2 D$$$. Then it's also quite easy to see that the numbers satisfying $$$|a + b \sqrt D | = 1$$$ form a group, as the norm is multiplicative, that is $$$| AB | = |A | \cdot |B|$$$.
Periodicity of the continued fraction
Consider the continued fraction $$$\sqrt D = [a_0; a_1, a_2, \dots]$$$ and its complete quotients $$$s_k = [a_k; a_{k+1}, \dots]$$$. We can represent them as
where $$$x_k$$$ and $$$y_k$$$ are rational numbers. In fact, we can even prove that $$$x_k$$$ and $$$y_k$$$ are integers, and that
for any convergent $$$\frac{p_{k-1}}{q_{k-1}} = [a_0; a_1, \dots, a_{k-1}]$$$. The proof is a bit technically involved, thus I will place it under the spoiler.
The identity above means that $$$|D q_k^2 - p_k^2| = 1$$$ if and only if $$$y_{k+1} = 1$$$, that is $$$s_{k+1} = \sqrt D + x_{k+1}$$$. As it turns out, this is tightly related to the periodicity of the continued fraction of $$$\sqrt D$$$, and only ever happens at $$$k=-1$$$ or after we just completed another period. To better understand this, we should prove the following fact:
That is, $$$x + y \sqrt D \in \mathbb Q[\sqrt D]$$$ has a periodic fraction if and only if it's larger than $$$1$$$ and $$$-1 < x - y \sqrt D < 0$$$.
The result above allows us to prove that $$$\sqrt{D}$$$ is periodic starting with $$$a_1$$$, because $$$-1 < s_1^* < 0$$$ from the very beginning. It, in turn, means that when we get to $$$s_{k+1} = \sqrt{D} + x_{k+1}$$$ it must hold that $$$x_{k+1} = \lfloor \sqrt D \rfloor = a_0$$$, otherwise $$$s_{k+1}^*$$$ won't be in the interval $$$(-1, 0)$$$. But this immediately tells us that $$$s_{k+2} = s_1$$$.
Finding all solutions
In other words, the first non-trivial solution to $$$|x^2 - y^2 D|=1$$$ is given by the convergent $$$\frac{p_k}{q_k}$$$, where
and all the non-negative solutions are the convergents with indices $$$t(k+1)-1$$$ for integer $$$t \geq 0$$$. This observation also allows to find the explicit correspondence between a convergent $$$\frac{p_k}{q_k}$$$ and its matrix representation. It generally holds for convergents that
If $$$|p_k^2 - D q_k^2| = 1$$$, we can additionally derive that
It makes a lot of sense, as applying the transform defined by such matrix to any convergent $$$\small \frac{p_t}{q_t}$$$ will simply add another period in front of it, thus moving it into $$$\small \frac{p_{t+k+1}}{q_{t+k+1}}$$$. This gives explicit isomorphism between addition in convergent's indices and group operation in the solution group in terms of matrices. This also highlights that the solution group is cyclic, that is all solutions are obtained as powers of the smallest solution $$$\frac{p_k}{q_k}$$$, which can be found with continued fractions.
Summary
In other words, to find all solutions to the Pell's equation $$$x^2 - Dy^2=1$$$, one may find the smallest solution as the first convergent $$$\small \frac{p_k}{q_k}$$$ of $$$\sqrt D = [a_0; a_1, a_2, \dots]$$$ such that $$$p_k^2 - D q_k^2 = 1$$$. Then all the other solutions are obtained as
for all integer $$$t \geq 0$$$.
To find all primitive pythagorean triples up to some threshold, instead of using this parametrization and checking that $$$\mathrm{gcd}(p, q) = 1$$$, one can use this spell. It is not very difficult to formally prove why it works, but I have absolutely no idea how it works, how to find similar representations for other quadratic forms, why matrices should equal these and overall lack structural understanding.
The Wikipedia page mentions that one can also perceive it as a tree in pairs $$$(m, n)$$$ s.t.
where $$$m > n > 0$$$ are co-prime and have different parities. The suggested transforms from $$$\frac{m}{n}$$$ are:
They can be interpreted in terms of Calkin-Wilf tree, starting from $$$\frac{2}{1}$$$ as:
I suppose one could use Euclid-like inductive reasoning to prove that all such $$$\frac{m}{n} > 1$$$ are reached.
Note that the continued fractions encode the path to a fraction in the Calkin-Wilf tree, so the transforms are even more meaningful in terms of continued fractions. If $$$\frac{m}{n} = [a_0; a_1, \dots, a_n]$$$, then the transforms are
I might be wrong, but for finding primitive pythagorean triples we can use the parametrization, but p, q must be of different parities.
Yes, this is true. I didn't mention that I'm looking for primitive triples only though.
Very interesting and concise blog, many thanks! Concerning Pell's equation, there is another approach, which is essentially Dirichlet's unit theorem applied to $$$\mathbf{Q}(\sqrt{D})$$$. The proof can be significantly simplified for this field (as well as we need Minkowski's theorem only for $$$\mathbf{Z}^2$$$, which can be proved just with the pigeonhole principle) and doesn't need heavy theory. Though, I don't think it will be shorter than continued fractions.
I'm a simple man. I see adamant blog, I click it. :D
There's a more general way to find rational points on rational polynomial curves. In general, if you have a curve of degree $$$d$$$, and you pick $$$d - 1$$$ (possibly degenerate, i.e., tangents are allowed) rational points on the curve that are on a rational straight line, and intersect that line with the curve again, the resultant intersection point (if it exists) will also be rational. The proof is pretty straightforward — just use Vieta's relations. This method is also used to find rational points on conics and elliptic curves (in the elliptic curve case, it is related to point addition).
For Pell's equation, it seems that the treatment is quite concise (but lacking a few important details and structure), so I'll just list out how you'd go about proving and finding all solutions using the continued fraction of $$$\sqrt{d}$$$ itself. (Though I prefer the treatment using $$$x_0 - y_0 \sqrt{d}$$$ as the unit that generates all other units of the real quadratic field).
Seriously nice blog, understandable after first reading.