Solving the diophantine equation x2 - Dy2=N, where D > 0 and not a perfect square

This finds least positive primitive solutions for the diophantine equation x2 - Dy2 = 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 x2 - Dy2 = 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 P2 ≡ D (mod |N|). Write x = Py + |N|X. Then
|N|X2 + 2PXy +(P2 - D)y2/|N| = N/|N|
and it can be proved that X/y is a convergent An-1/Bn-1 to ω = (-P + √D)/|N| and Qn = (-1)nN/|N|.
Then x = |N|An-1 + PBn-1 = Gn-1, say.
We have Gn-12 - DBn-12 = (-1)n|N|Qn = x2 - D y2 = N.
Hence Qn = (-1)nN/|N|.
The algorithm: For each congruence solution P, let (Pn + √D)/Qn denote the nth complete quotient for ω. We search the continued fraction expansion of ω for the least n such that Qn = (-1)nN/|N| and Gn-1 ≥ 0. Then (Gn-1, Bn-1) will be the least positive primitive solution of x2 - Dy2 = N for the class determined by P.


Enter D (< 1020):
Enter N (non-zero):
Enter E (0 or 1):

Last modified 6th January 2014
Return to main page