# First 30 solutions of Pell's equation.

$$k^2 - 5z^2 = 1 (1)$$

The equation (1) is a Pell equation where $$D = 5$$ and the first solution $$\{k, z\} =\{9, 4\} = \{p, q\} (2)$$

$$k_n = \dfrac{(p + q\sqrt{D})^n + (p - q\sqrt{D})^n}{2} ; z_n = \dfrac{(p + q\sqrt{D})^n - (p - q\sqrt{D})^n}{2\sqrt{D}} (3)$$

Where n is index of a solution. Replace p and q in (2) with (3)

$$k_n = \dfrac{(9 + 4\sqrt{5})^n + (9 - 4\sqrt{5})^n}{2} ; z_n = \dfrac{(9 + 4\sqrt{5})^n - (9 - 4\sqrt{5})^n}{2\sqrt{5}} (4)$$

Equations (35,36) from Pell's equation on Wolfram

I wrote a small Python program to calculate 30 solutions. However, the first 12 solutions are correct then all others are incorrect.

import math
def pell(u,v,D,nMaxUpperBound=30):
lstResult = []
a = u + v*math.sqrt(D)
b = u - v*math.sqrt(D)
for i in range (1, nMaxUpperBound, 1) :
k = int(round((math.pow(a, i) + math.pow(b, i))/2))
z = int(round((math.pow(a, i) - math.pow(b, i))/(2 * math.sqrt(D))))
if k**2 - D*z**2 == 1:
lstResult.append((k,z))
else:
print("failed. i = ", i, "k,z => (%d,%d)" % (k,z))
return lstResult
lstResult = pell(9, 4, 5)
print(lstResult)


What should I do to improve the program or use a different approach to calculate all 30 solutions?

Floating point precision is not sufficient for the larger solutions. You should probably look at using a recursive definition for the $k_n$ and $z_n$. They should satisfy $a_{n+2}=18a_{n+1}-a_{n}$ (derived from the minimum polynomial of $9\pm4\sqrt{5}$). Try to find a tail recursive solution, and you shouldn't run into problems until you overflow the ints.

• Python converts from normal ints to arbitrary precision ints when required, so overflow is not a problem. One can keep going this way until running out of time or memory. Aug 13, 2017 at 23:21
• @RossMillikan thanks for the tip. Aug 13, 2017 at 23:23
• How to explain why minimum polynomial of $9+/-4\sqrt{5}$ equals 18? Aug 13, 2017 at 23:27
• Let $\alpha = 9+4\sqrt{5}$ then it's "conjugate" is $\tilde \alpha = 9-4\sqrt{5}$ They are roots of the equation with coefficient of the middle term equal to $-(\alpha + \tilde \alpha) = -18$ and constant term equal to $\alpha \cdot \tilde \alpha = 1$, so $x^2-18x+1$ This comes from multiplying out $(x-r_1)(x-r_2)$. Aug 13, 2017 at 23:34
• Or, Cayley Hamilton, where $18$ is the trace of $$\left( \begin{array}{cc} 9 & 20 \\ 4 & 9\end{array} \right)$$ Aug 13, 2017 at 23:36

One problem is that round-off error will gradually wreak havoc on any floating-point computation using $\sqrt{D}$. A surer way to generate solutions is to use Python's arbitrarily-long integers and do integer arithemetic, looking at the integer-coefficient recursion for the (integer!) coefficients of powers $(a+b\sqrt{D})^n$.

• I used Decimal Python package. Results are the same. It failed after 12th solution. Aug 13, 2017 at 22:54
• @AlexGawkins, but that does not overcome the general problem, namely, that things like $\sqrt{D}$ have no exact finite decimal expression. That is, any choice of finite decimal expansion will eventually have blow-up of error. There's just no numerical way to make computations with irrational numbers precise. Symbolic, yes. Or, often, the easier way is to convert the issues to (infinite-precision, but always finite, at every particular moment!) integer arithmetic. Aug 13, 2017 at 22:57

Any formula of the form $a_n=cx^{n}+dy^{n}$ can be computed recursively as:

$$a_{n+1}=(x+y)a_n - xya_{n-1}$$ if you know $a_0,a_1$.

In your case, you've got $x,y=9\pm 4\sqrt{5}$, and $x+y=18$, $xy=1$.

So you have $k_0=1,k_1=9,z_0=0,z_1=4$:

$$k_{n+1}=18k_n-k_{n-1}\\z_{n+1}=18z_n-z_{n-1}$$

If you only want to compute one pair, $(k_n,z_n)$ you can compute this in $O(\log n)$ time by computing:

$$\begin{pmatrix}0&1\\ -1&18 \end{pmatrix}^{n-1}\begin{pmatrix}1&0\\ 9&4\end{pmatrix}=\begin{pmatrix}k_{n-1}&z_{n-1}\\k_{n}&z_{n}\end{pmatrix}$$

The reason this is an $O(\log n)$ operation is that you can apply exponentiation by squaring. to raise the matrix to the $n$th power.

But if you need all values $k_1,\dots,k_n$ and $z_1,\dots,z_n$, you'll of course require $O(n)$ time (you can't compute a list of $2n$ values in less than $O(n)$ time, after all.)

Another matrix approach, probably more direct, is to note that $$k_{n+1}+z_{n+1}\sqrt{5}=(k_n+z_n\sqrt5)(9+4\sqrt{5})=(9k_n+20z_n)+(4k_n+9z_n)\sqrt{5}$$

So we can write:

$$\begin{pmatrix}k_n\\z_n\end{pmatrix}=\begin{pmatrix}9&20\\4&9\end{pmatrix}^n\begin{pmatrix}1\\0\end{pmatrix}$$

Again, exponentiation by squaring allows $O(\log n)$ time calculation. (It might be closer to $O(\log^2 n)$ in both cases because of the number of bits in the integer multiplications, but still much better than $O(n)$.)

I also found recursive solution at Additional solutions from the fundamental solution $$\displaystyle x_{k}+y_{k}{\sqrt {D}}=(x_{1}+y_{1}{\sqrt {D}})^{k}(1)$$ that led them to $$\displaystyle x_{k+1}=x_{1}x_{k}+Dy_{1}y_{k} (2)$$ and $$\displaystyle y_{k+1}=x_{1}y_{k}+y_{1}x_{k} (3)$$

First of all let's proof (1)

To do this I shall use equation from my question

$$k_n = \dfrac{(p + q\sqrt{D})^n + (p - q\sqrt{D})^n}{2} ; z_n = \dfrac{(p + q\sqrt{D})^n - (p - q\sqrt{D})^n}{2\sqrt{D}} (3)$$

Equations (35,36) from Pell's equation on Wolfram or see (31,32) on the same page in Wolfram.

Calculate $$k_n+z_n\sqrt{D} = \dfrac{(p + q\sqrt{D})^n + (p - q\sqrt{D})^n}{2} + \dfrac{(p + q\sqrt{D})^n - (p - q\sqrt{D})^n}{2\sqrt{D}}\sqrt{D}= \dfrac{(p + q\sqrt{D})^n + (p - q\sqrt{D})^n}{2} + \dfrac{(p + q\sqrt{D})^n - (p - q\sqrt{D})^n}{2} = \dfrac{(p + q\sqrt{D})^n + (p - q\sqrt{D})^n + (p + q\sqrt{D})^n - (p - q\sqrt{D})^n}{2} = \dfrac{2(p + q\sqrt{D})^n}{2} = (p + q\sqrt{D})^n$$

What was required to prove !

Here's the proof:

From (1) for $$k+1$$ solutions

$$\displaystyle x_{k+1}+y_{k+1}{\sqrt {D}}=(x_{1}+y_{1}{\sqrt {D}})^{k+1}=(x_{1}+y_{1}{\sqrt {D}})^{k}*(x_{1}+y_{1}{\sqrt {D}}) (4)$$

First parentheses in (4) represents $$k$$-th solution , so (4) can be rewritten as

$$\displaystyle x_{k+1}+y_{k+1}{\sqrt {D}}=(x_{k}+y_{k}{\sqrt {D}})*(x_{1}+y_{1}{\sqrt {D}}) (5)$$.

Then open parentheses

$$x_{k+1}+y_{k+1}{\sqrt {D}} = x_kx_1+ x_ky_1\sqrt {D}+y_kx_1\sqrt {D}+y_ky_1D(6)$$

Combine similar terms

$$x_{k+1}+y_{k+1}{\sqrt {D}} = (x_kx_1+y_ky_1D) + \sqrt {D}(x_ky_1+y_kx_1) (7)$$

Now it is easy to see that

$$\displaystyle x_{k+1}=x_{1}x_{k}+Dy_{1}y_{k} (2)$$

and

$$\displaystyle y_{k+1}=x_{1}y_{k}+y_{1}x_{k} (3)$$

are true.

Let's find formula expressing $$x_{k+2}$$ by using equations (2) and (3)

from (2) $$x_{k+2}=x_{1}x_{k+1}+Dy_{1}y_{k+1} (8)$$

substitute $$y_{k+1}$$ from (3) in (8)

$$x_{k+2}=x_{1}x_{k+1}+Dy_{1}(x_{1}y_{k}+y_{1}x_{k}) (9)$$

simplify (9)

$$x_{k+2}=x_{1}x_{k+1}+Dy_{1}x_{1}y_{k}+Dy_{1}^2x_{k} (10)$$

$$x_{k+2}=x_{1}(x_{k+1}+Dy_{1}y_{k})+Dy_{1}^2x_{k} (11)$$

use (2) to simplify

$$x_{k+2}=x_{1}(x_{k+1}+x_{k+1}-x_{1}x_{k})+Dy_{1}^2x_{k} (12)$$

$$x_{k+2}=x_{1}(2x_{k+1}-x_{1}x_{k})+Dy_{1}^2x_{k} (13)$$

$$x_{k+2}=2x_{1}x_{k+1}-x_{1}^2x_{k}+Dy_{1}^2x_{k} (14)$$

$$x_{k+2}=2x_{1}x_{k+1}-x_{k}(x_{1}^2-Dy_{1}^2) (15)$$

You can see that in parentheses $$x_{1}^2-Dy_{1}^2$$ is original Pell's equation , and it equals to 1, so

$$x_{k+2}=2x_{1}x_{k+1}-x_{k} (16)$$

Formula (16) for $$k^2-5z^2=1$$ with the first solution $$\{k,z\}=\{9,4\}=\{x_1,y_1\}$$ gives us a recursive family of solutions $$x_{k+2}=2*9*x_{k+1}-x_{k}$$ or $$x_{k+2}=18x_{k+1}-x_{k}$$

Let's find formula expressing $$y_{k+2}$$ by using equations (2) and (3)

$$\text{from (3) } y_{k+2}=x_{1}y_{k+1}+y_{1}x_{k+1} (17)$$

substitute $$x_{k+1}$$ from (2)

$$y_{k+2}=x_{1}y_{k+1}+y_{1}(x_{1}x_{k}+Dy_{1}y_{k}) (17)$$

simplify

$$y_{k+2}=x_{1}y_{k+1}+y_{1}x_{1}x_{k}+Dy_{1}^2y_{k} (18)$$

get $$x_k$$ from (3) and apply it to (18)

$$x_k=\frac{y_{k+1} - x_1y_k}{y_1} (19)$$

$$y_{k+2}=x_{1}y_{k+1}+y_{1}x_{1}(\frac{y_{k+1} - x_1y_k}{y_1})+Dy_{1}^2y_{k} (20)$$

simplify (20)

$$y_{k+2}=x_{1}y_{k+1}+x_{1}y_{k+1} - x_1^2y_k+Dy_{1}^2y_{k} (21)$$

or

$$y_{k+2}=2x_{1}y_{k+1} - y_k(x_1^2-Dy_{1}^2) (22)$$

Again, you can see that in parentheses $$x_{1}^2-Dy_{1}^2$$ is original Pell's equation , and it equals to 1, so

$$y_{k+2}=2x_{1}y_{k+1} - y_k (23)$$

The updated program looks like this based on formulas (2) and (3)

def pell2(u,v,D,nMaxUpperBound=30):
lstResult = []
x1 = u
y1 = v
lstResult.append((x1,y1))
for i in range (1, nMaxUpperBound, 1) :
x2 = u*x1 + D*v*y1
y2 = u*y1 + v * x1
lstResult.append((x2,y2))
x1 = x2
y1 = y2
return lstResult
lstResult = pell2(9, 4, 5)
for j in range(0, len(lstResult)):
print(j+1, ' => ', lstResult[j])


Results:

1  =>  (9, 4)
2  =>  (161, 72)
3  =>  (2889, 1292)
4  =>  (51841, 23184)
5  =>  (930249, 416020)
6  =>  (16692641, 7465176)
7  =>  (299537289, 133957148)
8  =>  (5374978561, 2403763488)
9  =>  (96450076809, 43133785636)
10  =>  (1730726404001, 774004377960)
11  =>  (31056625195209, 13888945017644)
12  =>  (557288527109761, 249227005939632)
13  =>  (10000136862780489, 4472197161895732)
14  =>  (179445175002939041, 80250321908183544)
15  =>  (3220013013190122249, 1440033597185408060)
16  =>  (57780789062419261441, 25840354427429161536)
17  =>  (1036834190110356583689, 463686346096539499588)
18  =>  (18605234632923999244961, 8320513875310281831048)
19  =>  (333857389202521629825609, 149305563409488533459276)
20  =>  (5990827771012465337616001, 2679179627495483320435920)
21  =>  (107501042489021854447262409, 48075927731509211234387284)
22  =>  (1929027937031380914713107361, 862687519539670318898535192)
23  =>  (34615001824075834610388670089, 15480299423982556528939246172)
24  =>  (621141004896333642072282954241, 277782702112146347202007895904)
25  =>  (11145923086309929722690704506249, 4984608338594651693107202880100)
26  =>  (200005474548682401366360398158241, 89445167392591584128727643945896)
27  =>  (3588952618789973294871796462342089, 1605028404728053862623990388146028)
28  =>  (64401141663670836906325975923999361, 28801066117712377943103099342682608)
29  =>  (1155631597327285091018995770169646409, 516814161714094749113231797780140916)
30  =>  (20736967610227460801435597887129636001, 9273853844735993106095069260699853880)


Questions?

• Proof is simply by definition of the $k$-th power. Are you sure you can't get it by yourself? =) Aug 14, 2017 at 8:58