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?
 A: 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$.
A: 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)$.)
A: 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.
