Genesect's Riddle

The continued fraction expansion of $\sqrt{649}$ is $\big[a_0,a_1,a_2,\ldots \big] \; = \; \big[25, \overline{2, 9, 1, 2, 3, 1, 1, 2, 1, 4, 1, 16, 6, 3, 4, 3, 6, 16, 1, 4, 1, 2, 1, 1, 3, 2, 1, 9, 2, 50} \big]$ In particular, we note that it is periodic of period $30$ . Since $30$ is even, it is standard continued fraction bookwork that positive integer solutions of the equation $x^2 - 649y^2 \; = \; 1$ are given by $x \; = \; p_{30u-1} \;, \; y \; = \; q_{30u-1} \hspace{2cm} u \in \mathbb{N}$ where $\frac{p_n}{q_n} \; = \; \big[a_0,a_1,\ldots,a_n\big]$ are the convergents of the continued fraction expansion. These are most readily determined by the recurrence relations $\begin{aligned} p_n & = a_n p_{n-1} + p_{n-2} & \hspace{1cm} p_0 = a_0 & \hspace{0.5cm} p_{-1} = 1 \\ q_n & = a_n q_{n-1} + q_{n-2} & \hspace{1cm} q_0 = 1 & \hspace{0.5cm} q_{-1} = 0 \end{aligned}$ Since we calculate $p_{29} \; = \; 1123593226162199 \hspace{2cm} q_{29} \; = \; \boxed{44104892095380}$ we have the answer.

We will generalize the problem to $x^2 - Dy^2 = 1$ , where $D$ is a positive integer and not a square number. (If $D$ is not positive, we have $-1 \le x \le 1$ ; just test them. If $D$ is a square number, we can factorize the left hand side into a product of two integers, and so the answer is again simple.)

If you know it, this is Pell's equation. There has been a lot of research about this; in fact, there is an algorithm where, given $D$ , it can generate all the solutions efficiently. The idea is that if $x^2 - Dy^2 = 1$ , then $x^2$ and $Dy^2$ are very close to each other, and so $\frac{x^2}{y^2} \approx D$ , or $\frac{x}{y} \approx \sqrt{D}$ . We can then use the continued fraction expansion of $\sqrt{D}$ to obtain the smallest solution, because we know approximations obtained from continued fractions are particularly accurate (e.g. $\frac{355}{113}$ approximation to $\pi$ is pretty accurate, because it's obtained from the continued fraction expansion).

For example, consider $D = 11$ . The continued fraction expansion is

$\sqrt{11} = 3 + \frac{1}{3 + \frac{1}{6 + \frac{1}{3 + \frac{1}{6 + \ldots}}}}$

Taking the approximations, we have:

$\begin{aligned} 3 &= \frac{3}{1} \\ 3 + \frac{1}{3} &= \frac{10}{3} \\ 3 + \frac{1}{3 + \frac{1}{6}} &= \frac{63}{19} \\ 3 + \frac{1}{3 + \frac{1}{6 + \frac{1}{3}}} &= \frac{199}{60} \\ \ldots \end{aligned}$

Thus, the candidate solutions are those approximations: $\frac{3}{1}, \frac{10}{3}, \frac{63}{19}, \frac{199}{60}, \ldots$ . By checking one by one, we find out that $(10, 3)$ is a solution to $x^2 - 11y^2 = 1$ . And later in the sequence, we find another solution $(199, 60)$ , and so on; continuing the sequence will give more solutions.

Moreover, computing these approximations is not too difficult, and it can be done way faster than brute-force search. First, we can generate the continued fraction expansion:

1. Start with $P = 0, Q = 1$ . We're going to find the next term in the continued fraction expansion of $\frac{\sqrt{D} + P}{Q}$ .
2. Find the largest integer $a$ not exceeding $\frac{\sqrt{D} + P}{Q}$ and append it to the continued fraction expansion so far. The rest of the number is $\frac{\sqrt{D} + P - aQ}{Q}$ , which we must now take the inverse to continue with the next term of the expansion.
3. Rationalize $\frac{Q}{\sqrt{D} + P - aQ} = \frac{Q}{D - (P - aQ)^2} \cdot (\sqrt{D} - P + aQ)$ . The first factor (the fraction) magically always cancels into a unit fraction, so our number is $\frac{\sqrt{D} - P + aQ}{k}$ for some integer $k$ .
4. Set $P \gets -P + aQ, Q \gets k$ .
5. Go to step 2.

This will give a never-ending stream of the expansion. How can we use it to compute the approximation? Easy. Each time we obtain a term, we try cutting it off there and finding the current approximation:

1. Suppose the continued fraction expansion of $\sqrt{D}$ is $[a_0; a_1, a_2, a_3, \ldots, a_n]$ so far (we got $n+1$ terms of the expansion).
2. Start $i = n, P = a_n, Q = 1$ . This is the innermost fraction of the expansion.
3. Decrement $i$ . This takes us one level back; now we need to simplify the fraction $a_i + \frac{1}{P/Q}$ .
4. Set $P \gets a_i P + Q, Q \gets P$ , since that's what our fraction evaluates to.
5. If $i > 0$ , go to step 3. If $i = 0$ , we now have our candidate $(P, Q)$ .

Finally, it's trivial to check whether our candidate is indeed a solution or not. This might cause problems in languages that don't support arbitrary-precision integers, though. You can either move to a language that supports it, like Python or Java, or you can rig in some way to support larger integers. Fortunately, for this problem, it's possible to do with a small array in C++ to represent a big number; details follow.

Either way, we have found a solution. How can we guarantee that's the minimum? That's yet more theory. The intuitive reason is precisely the same reason why we're looking at continued fraction expansions: approximations obtained in any other way are weaker and so the difference between $\frac{x^2}{y^2}$ and $D$ will become bigger, which means $x^2 - Dy^2$ will also become bigger, far from 1.

Anyway, for $D = 649$ the smallest solution is $(1123593226162199, \boxed{44104892095380})$ .

---

This is a way to use the fact that C++ supports 64-bit integers to make a data structure that is effectively 120-bit integers. Any such 120-bit integer will in fact be an array of four 64-bit integers a[3..0] , which will represent the number $2^{90} a[3] + 2^{60} a[2] + 2^{30} a[1] + 2^0 a[0]$ . Whenever we do something to this array, we will immediately perform necessary carries so all a[i] remain 30-bit. (We don't use 32 bits because there's a subtlety later that we need to handle.)

Now, why will this be enough? Straightforward checking requires us to multiply $x$ by itself, as well as 649 with $y^2$ . Since we're told that $y$ has 14 digits, this means $y$ has at most 47 bits. Since $x \approx y \sqrt{649}$ , $x$ has at most 53 bits. This fits inside 60 bits, so if we represent $x$ in our array method, we will have x[3], x[2] zero. Likewise, y[3], y[2] will also be zero.

To compute $x^2$ and store it into x2 , we simply use schoolbook multiplication: x2[2] = x[1] * x[1], x2[1] = x[1] * x[0] + x[0] * x[1], x2[0] = x[0] * x[0] . (The reason we use 30 bits per element is clear here: if we use the full 32 bits, x2[1] can overflow.) This is schoolbook multiplication except each digit ranges up to $2^{30}$ instead of up to 9. Remember that after we do this multiplication, we need to perform carries: if x2[0] is larger than $2^{30}$ , we take the quotient when dividing by $2^{30}$ and add it to x2[1] , so x2[0] can simply store the remainder.

We compute $649y^2$ similarly. One method that is easy to prove is to compute $649y$ times $y$ ; we can find out that $649y$ fits in 57 bits, so it's not a problem. We obtain y2 that stores this result.

Finally, to check the computation, we can simply do schoolbook subtraction to compute x2 - y2 . Just compute the difference element-by-element, then perform any necessary carries. If the result is 1 (that is, element 0 has value 1, and the rest have value 0), then we have found our answer.