Most efficient checking algorithm for Pell's Equation

What is the most computationally efficient way to check, given $x,y,D$ that they satisfy Pell's equation (positive or negative) ($x^2-Dy^2=1$)? (Obviously the question is concerned with very large values of $x,y,D$.)

I know (I think) that it'll have to be checking mod $p$ but I just can't find the right balance between computation time required for factorizing $x,y$ or $D$ and brute-force calculations.

Update: What I said about mod $p$ had to do with the fact that I thought that with Pell especially there might be some computationally efficient way to get a distinguished set of primes (and then the question was how best to balance the 'finding' with the moding) - but as Carnahan indicates the case might not be more special than that of evaluating any binomial. I think I'm pretty much covered by what Carnahan says, though I would be interested to know how these bounds on the number of operations are obtained. Also does it really follow that there is no more efficient way to evaluate a 'close' possible triple (i.e. one that has passed a mod-small-prime or log test) than to use multiplication algorithms? Is that obvious?

• Start by comparing 2log(x) with log(D) + 2 log(y). Then consider multiplying the high order bits in comparison, and the low order bits. You can then try mod p computations if the first three checks show sanity. Gerhard "Ask Me About System Design" Paseman, 2011.04.01 – Gerhard Paseman Apr 1 '11 at 17:39
• This question warrants the question: what other operations on these large numbers are you doing, that computing $x^2-Dy^2$ is more expensive than? – Dror Speiser Apr 1 '11 at 20:29
• Computations mod $p$ don't need any factoring, so I don't see why you are worried about the computation time that it takes. Also, it would help if you explained what your end goal is: If you are just seeking solutions to the equation, then an algorithm using the $PGL_2$-structure is much better than a brute force check. If you have some other program that for some reason produces large triples $(x,y,D)$, then Paseman's suggestion in the previous comment is good for quickly removing bad candidates from consideration. Asymptotically, the Chinese remainder thm is about as good as brute force. – S. Carnahan Apr 2 '11 at 4:14
• The end goal is to simply check if a triple $(x,y,D)$ satisfies Pell's equation - it's not brute forcing $(x,y)$ given $D$. The log filter obviously gets rid of a bad triple very quickly. My problem is that if I go mod $p$ then I don't see a way to find a 'privileged' set of primes relevant to my triple (e.g. factors of $D$) without getting more expensive than simply checking mod $p$ for 'all' primes in ascending order. So I guess the general question is how to strike the most efficient balance between getting a privileged set and 'brute force' in the sense of moding for undistinguished primes – Chuck Apr 2 '11 at 9:22
• (cont.) That said, mod $p$ might not be so 'obviously' the most efficient route: I'm trying to get my head around Jagy's algorithm below (my computational savvy is embarassingly non-existent) - in any case my question is really not related to what I'm working on, it's more theoretical in nature, i.e. optimal way to do this, so I guess another question is how does one prove that this or that algorithm is optimal wrt checking if a triple satisfies Pell – Chuck Apr 2 '11 at 9:26

Based on the comments, it looks like this is not a question specific to Pell's equation, and that you just want to evaluate a single binomial with big inputs as quickly as possible.

If you check the equation directly using fast multiplication algorithms (e.g., Schönhage-Strassen), you can expect the calculation to require about $(\log z )\cdot (\log \log z)$ operations, where $z = \max \{ x, y, D \}$.

If you want a way to quickly find a negative answer, you can check relative sizes and leading digits, then try reduction modulo small integers. If you start with a random negative answer, you can expect it to be eliminated after at most a few divisions (i.e., about $\log z$ operations).

To find a positive answer using modular arithmetic, you can use the Chinese remainder theorem. To prove that the identity holds, it suffices to check it modulo $n$, for $n$ ranging over a collection of positive integers whose least common denominator is larger than $x^2$ and $Dy^2$. It is common to check modulo a large collection of small primes, and this will require about $\log z$ primes and $(\log z)^2$ operations. Another natural choice with a binary computer is Fermat numbers, of the form $2^{2^n}+1$, since the division-with-remainder can be optimized - this ends up looking a lot like direct calculation.

In summary, the advantage of checking modulo small primes is that it lets you quickly eliminate negative answers, and the disadvantage is that (if I'm not mistaken) it is roughly quadratically slower than direct calculation when you have a positive answer. You can choose your method depending on exactly what sort of calculation you plan to do.

• Thanks for this Scott - that pretty much answers my question. But isn't Paseman's log test more efficient that small-prime-moding to rule out negative answers? In which case (contrary to what I thought) there is absolutely no reason to stick with a mod $p$ method even with Pell...That seems counter-intuitive to me, but then again my 'computation' instincts are very....primitive – Chuck Apr 3 '11 at 14:30
• The log test is fastest, so it should be applied first. The small primes (or small Fermat numbers) are the next fastest, and are likely to catch the negative candidates that slip through the cracks of the log test, so they should be applied second. – S. Carnahan Apr 3 '11 at 16:17

See my answer at:

Instead of using continued fractions to solve Pell's equation in the first place, thus needing high accuracy real numbers repeatedly, and needing a means to check the solution at the end, one may do the whole process in "reduced" indefinite binary quadratic forms.

The additional information used, over the earlier answer posted, is that the little 2 by 2 matrices (in the answer I call them $R$) show how to update your $(x,y)$ pair. At the end of the "cycle", you have come back to the reduced form $(1, 2 a_0, a_0^2 - n).$ The product of all the $R$ matrices in order give an automorph of that form, in particular the left hand column are the first $(x,y)$ such that $x^2 + 2 a_0 x y + (a_0^2 - n) y^2 = 1.$ All you need now is to pull back, $( x + a_0 y)^2 - n y^2 = 1.$ Here $a_0^2 < n < (a_0 + 1)^2.$

Here is a sample from my home computer. I am assuming you have arbitrary precision integer arithmetic, I do not.

My own program is not quite what I described, I give a separate automorph for the original Pell form $x^2 - n y^2$, then the unit is just the left column of that. let me know if you would like the C++ program.

As you can see, the upper right entry in the Pell automorph is oversize for C++, I had to fill in 61 * 226153980 and other values by hand.

phoebus:~/Cplusplus> ./Pell
Input n for Pell
61

0  form   1 14 -12   delta  -1
1  form   -12 10 3   delta  4
2  form   3 14 -4   delta  -3
3  form   -4 10 9   delta  1
4  form   9 8 -5   delta  -2
5  form   -5 12 5   delta  2
6  form   5 8 -9   delta  -1
7  form   -9 10 4   delta  3
8  form   4 14 -3   delta  -4
9  form   -3 10 12   delta  1
10  form   12 14 -1   delta  -14
11  form   -1 14 12   delta  1
12  form   12 10 -3   delta  -4
13  form   -3 14 4   delta  3
14  form   4 10 -9   delta  -1
15  form   -9 8 5   delta  2
16  form   5 12 -5   delta  -2
17  form   -5 8 9   delta  1
18  form   9 10 -4   delta  -3
19  form   -4 14 3   delta  4
20  form   3 10 -12   delta  -1
21  form   -12 14 1   delta  14
22  form   1 14 -12

disc   244

Automorph, written on right of Gram matrix:
-183241189   -2713847760
-226153980   -3349396909

Pell automorph
-1766319049  -13795392780
-226153980   -1766319049

Pell unit
1766319049^2 - 61 * 226153980^2 = 1

=========================================
phoebus:~/Cplusplus>

It probably depends on the situation. If are convinced that indeed $x^2-Dy^2=1$, but you have to provide a proof, and have no special knowledge about $x,y,D$, then I'm not sure that anything would be faster than a direct calculation. You can check $\mod m$ for various $m$. There is no need for $m$ to be prime. If the numbers are provided in decimal then it would be easy to check $\mod 10^k+1$ or $10^k-1.$ If you know the basic solution $x_0^2-Dy_0^2=1$ then you could figure out the correct power $t$, if any with $(x_0+\sqrt{D}y_0)^t=x+\sqrt{D}y$ and do effective powering by squaring to see if it works. One could pick a base $b$ and and check $\mod b^i$ for $i=1,2,3,\cdots$ until $b^i$ is large enough. In the event that it is not true, you might find out quickly. If it is true, that might not be faster than direct calculation.