First 30 solutions of Pell's equation.

Question

$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?

Answer 1

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.

Answer 2

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$.

Answer 3

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)$.)

Answer 4

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?