Can't find fundamental solution for $x^2 - 61y^2 = 1$ through continued fractions

Question

I'm trying to find the fundamental solution for the $n = 61$ case of Pell's equation $x^2 - ny^2 = 1$ through continued fractions. I know the lowest solution is $x = 1766319049$, $y = 226153980$, but the continued fraction expansion for $\sqrt{61}$ appears to miss the case where this pair of values for the nominator and denominator occurs. This upsets me because the Wikipedia page for Pell's equation explains that calculating the continued fraction for $\sqrt{n}$ should solve the problem. What am I doing wrong?

This is the code I'm using:

double r = sqrt(61);
unsigned long long h1 = 0;
unsigned long long h2 = 1;
unsigned long long k1 = 1;
unsigned long long k2 = 0;

for(int i=0;i<100;i++) {
    unsigned long long b = r;
    unsigned long long temp = b * h2 + h1;
    h1 = h2;
    h2 = temp;

    temp = b * k2 + k1;
    k1 = k2;
    k2 = temp;

    r = (double) 1 / (double) (r - b);
}

Answer 1

To answer your question (and to echo @LordSharktheUnknown's comment), what you're doing wrong is right here:

double r = sqrt(61);
....
r = (double) 1 / (double) (r - b);

That number, $r$, that you're computing, is not $\sqrt{61}$. It's an approximation to it using double-precision arithmetic. Playing around with a different language, perhaps using floats instead of doubles, I find that $61 - r^2$, which should be $0$, is about $7 \times 10^{-15}$, which is definitely not zero. So you're finding (sort of) the continued fraction for a number that is not $\sqrt{61}$; you should not be surprised that the results you get are of little use in saying anything about $\sqrt{61}$.

Answer 2

The period of the continued fraction expansion for $\sqrt{61}$ is $11$, and that convergent, $\frac{p_{11}}{q_{11}}=\frac{29718}{3805}$, yields $-1$. Thus the solutions to the Pell equation $x^2-61y^2=1$ are $\frac{p_{22}}{q_{22}}$, $\frac{p_{44}}{q_{44}}$, $\frac{p_{66}}{q_{66}}, \ldots$.

Answer 3

Here is the "back of the envelope" way to construct the CF of $\sqrt{61}$. $$a_0=\left\lfloor{\sqrt{61}}\right\rfloor=7$$ $$x_1=\frac1{\sqrt{61}-a_0}=\frac{\sqrt{61}+7}{12}$$ $$a_1=\left\lfloor{x_1}\right\rfloor=1$$ $$x_2=\frac1{x_1-a_1}=\frac{\sqrt{61}+5}3$$ $$a_2=\left\lfloor{x_2}\right\rfloor=4$$ $$x_3=\frac1{x_2-a_2}=\frac{\sqrt{61}+7}4$$ $$a_3=\left\lfloor{x_3}\right\rfloor=3$$ $$x_4=\frac1{x_3-a_3}=\frac{\sqrt{61}+5}9$$ $$a_4=\left\lfloor{x_4}\right\rfloor=1$$ $$x_5=\frac1{x_4-a_4}=\frac{\sqrt{61}+4}5$$ $$a_5=\left\lfloor{x_5}\right\rfloor=2$$ $$x_6=\frac1{x_5-a_5}=\frac{\sqrt{61}+6}5$$ $$a_6=\left\lfloor{x_6}\right\rfloor=2$$ $$x_7=\frac1{x_6-a_6}=\frac{\sqrt{61}+4}9$$ $$a_7=\left\lfloor{x_7}\right\rfloor=1$$ $$x_8=\frac1{x_7-a_7}=\frac{\sqrt{61}+5}4$$ $$a_8=\left\lfloor{x_8}\right\rfloor=3$$ $$x_9=\frac1{x_8-a_8}=\frac{\sqrt{61}+7}3$$ $$a_9=\left\lfloor{x_9}\right\rfloor=4$$ $$x_{10}=\frac1{x_9-a_9}=\frac{\sqrt{61}+5}{12}$$ $$a_{10}=\left\lfloor{x_{10}}\right\rfloor=1$$ $$x_{11}=\frac1{x_{10}-a_{10}}=\sqrt{61}+7$$ $$a_{11}=\left\lfloor{x_{11}}\right\rfloor=14$$ $$x_{12}=\frac1{x_{11}-a_{11}}=\frac{\sqrt{61}+7}{12}=x_1$$ and here the periodicity kicks in: $x_{n+11}=x_n$ and $a_{n+11}=a_n$ for $n\ge1$. So the CF is $$[7,1,4,3,1,2,2,1,3,4,1,14,1,4,3,1,2,2,1,3,4,1,14,\ldots].$$

Note I have only used two-digit integer arithmetic to calculate the CF. In general for a non-square positive integer $d$, the $x_n$ in the CF for $\sqrt d$ will have the form $$x_n=\frac{\sqrt d+b_n}{c_n}$$ where $b_n$ and $c_n$ are positive integers, and $c_n\mid(d-b_n^2)$.

Of course in this example, to compute the convergents one still needs integer arithmetic, but with more than two digits.

Answer 4

There are 53 integer bits in a double-precision floating point number.

$x^2 > 10^{18} \gt 2^{59}$ and this is more digits than can be exactly computed.

Answer 5

Expanding on rogeri's answer, once the minimum solution $\;(p_{22},q_{22})\;$ is computed, the other solutions may be alternatively (and perhaps more easily) computed via

Define $(P_k,Q_k)\;$ by $\;\left(p_{22} + q_{22}\sqrt{61}\right)^k\;$ is equal to $P_k + Q_k\sqrt{61}.\;\;$

Then $\;(P_2,Q_2) = (p_{44},q_{44}), \;(P_3,Q_3) = (p_{66},q_{66}), ...$