Most efficient checking algorithm for Pell's Equation

Question

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?

Answer 1

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.

Answer 2

See my answer at:

Upper bound of period length of continued fraction representation of very composite number square root

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>

Answer 3

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.