#### Solving the diophantine equation x^{2} - Dy^{2}=N, where D > 0 and not a perfect square

This finds least positive *primitive* solutions for the diophantine equation x^{2} - Dy^{2} = N, D > 0, D not a perfect square in each equivalence class of solutions.

Then all positive solutions in a class are obtained by multiplying each least solution by η^{n}, n ≥ 0, where η is the least solution of Pell's equation x^{2} - Dy^{2} = 1.
E=1 is verbose and prints partial and complete quotients, as well as convergents, up to the end of a period (or double period, in the case of odd period).

The algorithm goes back to Lagrange 1770 and should be better known. Summarizing from a slide-talk by Keith Matthews:

Suppose (x,y) is a solution with gcd(x,y)=1 and x > 0, y > 0. Let x ≡ yP (mod |N|). Then P^{2} ≡ D (mod |N|). Write x = Py + |N|X. Then

|N|X^{2} + 2PXy +(P^{2} - D)y^{2}/|N| = N/|N|

and it can be proved that X/y is a convergent A_{n-1}/B_{n-1} to ω = (-P + √D)/|N| and Q_{n} = (-1)^{n}N/|N|.

Then x = |N|A_{n-1} + PB_{n-1} = G_{n-1}, say.

We have G_{n-1}^{2} - DB_{n-1}^{2} = (-1)^{n}|N|Q_{n} = x^{2} - D y^{2} = N.

Hence Q_{n} = (-1)^{n}N/|N|.

The algorithm: For each congruence solution P, let (P_{n} + √D)/Q_{n} denote the nth complete quotient for ω. We search the continued fraction expansion of ω for the least n such that Q_{n} = (-1)^{n}N/|N| and G_{n-1} ≥ 0. Then (G_{n-1}, B_{n-1}) will be the least positive primitive solution of x^{2} - Dy^{2} = N for the class determined by P.

*Last modified 6th January 2014*

Return to main page