Pell's equation

It is not obvious at all that a hyperbola $x^2-dy^2=1$, which certainly has real points, must pass through a lattice point — this theorem guarantees it always does, for every non-square $d$, via a clever pigeonhole argument on rational approximations.

By Dirichlet's approximation theorem, for any integer $Q>0$ there exist integers $p,q$ with $1\le q\le Q$ and $\left|\sqrt d - \frac{p}{q}\right|<\frac{1}{qQ}$. Letting $Q\to\infty$ produces infinitely many pairs $(p,q)$ with $\left|\sqrt d-\frac{p}{q}\right|<\frac{1}{q^2}$.

For each such pair, $|p-q\sqrt d|<\frac{1}{q}$, so $|p^2-dq^2|=|p-q\sqrt d|\cdot|p+q\sqrt d|<\frac{1}{q}\left(2q\sqrt d+\frac{1}{q}\right)<2\sqrt d+1$. Thus $p^2-dq^2$ takes one of only finitely many integer values in $(-2\sqrt d-1,\,2\sqrt d+1)$, while there are infinitely many pairs $(p,q)$.

By the pigeonhole principle, some fixed nonzero integer $k$ in that finite range satisfies $p^2-dq^2=k$ for infinitely many pairs $(p,q)$. Among these infinitely many pairs, again by pigeonhole, infinitely many share the same residues $p\equiv p_0,\ q\equiv q_0\pmod{|k|}$.

Take two distinct such pairs $(p_1,q_1)\neq(p_2,q_2)$ with $p_1^2-dq_1^2=p_2^2-dq_2^2=k$ and matching residues mod $|k|$. Set $x+y\sqrt d = \frac{(p_1+q_1\sqrt d)(p_2-q_2\sqrt d)}{k}$; expanding shows $x=\frac{p_1p_2-dq_1q_2}{k}$ and $y=\frac{p_1q_2-p_2q_1}{k}$ are honest integers precisely because of the matching residues mod $|k|$, and multiplicativity of the norm $N(a+b\sqrt d)=a^2-db^2$ gives $x^2-dy^2=\frac{k\cdot k}{k^2}=1$. Since $(p_1,q_1)\neq(p_2,q_2)$ but they give the same ratio in the limit, one checks $y\neq 0$, and replacing $(x,y)$ by $(|x|,|y|)$ (still a solution, since only squares appear) gives a solution in positive integers.

It says the infinitely many solutions are not a mysterious scattered set but a completely predictable geometric-like sequence generated by repeatedly "multiplying" the smallest one — reducing an infinite search to finding just one number.

First check $(x_n,y_n)$ defined by $x_n+y_n\sqrt d=(x_1+y_1\sqrt d)^n$ is indeed a solution for every $n$: taking conjugates, $x_n-y_n\sqrt d=(x_1-y_1\sqrt d)^n$, so $x_n^2-dy_n^2=(x_n+y_n\sqrt d)(x_n-y_n\sqrt d)=\left[(x_1+y_1\sqrt d)(x_1-y_1\sqrt d)\right]^n=(x_1^2-dy_1^2)^n=1^n=1$.

Now suppose $(x,y)$ is any positive integer solution not of this form; since $x_n\to\infty$ as $n\to\infty$, there is a unique $n$ with $x_n+y_n\sqrt d \le x+y\sqrt d < x_{n+1}+y_{n+1}\sqrt d = (x_n+y_n\sqrt d)(x_1+y_1\sqrt d)$.

Divide through by $(x_n+y_n\sqrt d)$, i.e. multiply by its inverse $(x_n-y_n\sqrt d)$ (valid since $x_n^2-dy_n^2=1$): set $x'+y'\sqrt d = (x+y\sqrt d)(x_n-y_n\sqrt d)$. Then $1\le x'+y'\sqrt d < x_1+y_1\sqrt d$, and $x'^2-dy'^2=(x^2-dy^2)(x_n^2-dy_n^2)=1\cdot 1=1$, so $(x',y')$ is also a solution of Pell's equation.

A short computation using $x'+y'\sqrt d\ge 1$ and $x'^2-dy'^2=1$ shows $x'\ge 1$ and $y'\ge 0$ (a solution with $x'+y'\sqrt d\ge 1$ but $y'<0$ would force $x'>x_1$, contradicting $x'+y'\sqrt d<x_1+y_1\sqrt d$ combined with the norm equation). If $y'>0$, then $(x',y')$ is a positive solution with $x'+y'\sqrt d<x_1+y_1\sqrt d$, contradicting minimality of the fundamental solution. So $y'=0$, forcing $x'=1$, i.e. $x+y\sqrt d=x_n+y_n\sqrt d$, so $(x,y)=(x_n,y_n)$ after all — contradicting our assumption.

Hence every positive solution is exactly some $(x_n,y_n)$, proving the fundamental solution generates the entire solution set.