Exercise 7.8.9 (Least solution of $x^2 - 29 y^2 = \pm 1$)

Proof. We generalize the method of Problem 8 to give a program in pure Python to obtain the least solution of $x^2-d{y}^2=\pm 1$.

from math import sqrt

def pell_fermat(d):
    """input  : d > 1 ,not a perfect square
        output: least solution of x^2 - d y^2 = +/- 1
    """
    (h, k) = (0, 1)
    (H, K) = (1, 0)
    Vd = int(sqrt(d))
    m = Vd
    q = d - m * m
    a = Vd
    (h, k, H, K) = (H, K , a * H + h, a * K + k)
    while q != 1:
        a = (m + Vd) // q
        m = a * q - m; q = (d - m * m) // q
        (h, k, H, K) = (H, K , a * H + h, a * K + k)
    return (H, K,   H * H - d * K * K)

def fermat(d):
    """input : d > 1 ,not a perfect square
       output: least solution of x^2 - d y^2 = 1
    """
    (p, q, epsilon) = pell_fermat(d)
    if epsilon == -1:
        (p, q) = (p * p + d * q * q, 2 * p * q)
    return (p, q)

if __name__ == ’__main__’:
    d = 29
    print(pell_fermat(d))
    print(fermat(d))

(70, 13, -1)
(9801, 1820)

So the least solution of $x^2-29y^2=-1$ is $(70,13)$, and the least solution of $x^2-29y^2=1$ is $(9801,1820)$: