The Diophantine equation x² − Ny² = 1, where N is a constant, is called Pell's equation, named after John Pell. Pell never actually worked on what is now called Pell's equation, but Euler thought he did, so there you go. Perhaps a better name for the equation might be "Brahmagupta's equation", after the Indian mathematician who studied the equation about a millennium before Pell didn't.
Some of the ancient Greeks may have worked on Pell's equation, but the first significant progress on the equation was made in India by Brahmagupta. In the 12th century, Bhaskara found a general solution to Pell's equation based on Brahmagupta's work. In the 17th century, European mathematicians rediscovered the equation, and in 1766 Lagrange proved that the equation always has a solution in natural numbers when N is not a square number (if N is a square number, solving the equation would require finding two squares that differ by 1; however, the only two such squares are 0 and 1, and 0 isn't a natural number).
The following illustrates how to solve Pell's equation, based on Bhaskara's method. For now, I've omitted proofs that the method works, but I may add this later. The method is called the Chakravala Method (also spelled Cakravâla method), also known as the Cyclic Method (if you're familiar with yoga and Eastern spirituality, you may know that "chakra" is Hindi for "circle" or "wheel"). This is a beautiful algorithm, perhaps one of the most beautiful discovered until relatively recently.
To solve the equation x² − Ny² = 1 (which we'll write as Ny² + 1 = x² from here on), Brahmagupta would first find an auxilliary equation Na² + k = b² where a, b, and k are small enough to be found by inspection and where the absolute value of k is not too big. For example, if you were trying to solve
Now, Brahmagupta would use what is called the Principle of Composition or Brahmagupta's Lemma to compose the auxilliary equation with another auxilliary equation or with itself. The principle of composition states that, if you have a1, b1, and k1 such that
Now, if k = −1, ±2, or ±4, composing the auxilliary equation with itself will allow you to find a solution to the original equation right away. For solving 60y² + 1 = x² we found the auxilliary equation 60 · 1² + 4 = 8². Composing that with itself, we get:
Each term is a multiple of 16, so divide each term by 16:
We now have a solution to the original equation: x = 31, y = 4. This is the smallest solution in natural numbers; to find other solutions, you can apply the principle of composition to the equation immediately above.
Now, this process can be hit-and-miss, and Brahmagupta wasn't able to use it to solve all such equations. However, Bhaskara's method allows us to find a solution all the time. Using Brahmagupta's Lemma, Bhaskara first showed that, if
Choose m such that am + b⁄k is an integer (in which case the other two fractions will be as well) and so that m² − N has the smallest absolute value possible. This will result in new auxilliary equation; if the new value of k is 1, you're done; otherwise, if it's −1, ±2, or ±4, you can use the Brahmagupta's method to find the solution; otherwise, repeat this process with the new equation,
Here's an example of the process for x² − 19y² = 1. Given this equation, one auxilliary equation we could form would be:
Now, we'll create the new auxilliary equation
Now, we need to find m such that 1 · m + 5⁄6 is an integer and |m² − 19| is as small as possible. By inspection, choose m = 1. Our auxilliary equation becomes:
Now, create a new auxilliary equation as above. Find m such that 1 · m + 4⁄−3 is an integer and |m² − 19| is as small as possible. A good choice is m = 5. This gives us the new auxilliary equation:
Since the value of k in this equation is −2, we can compose this equation with itself, as in the previous example. This results in:
Each term can be divided by 4, resulting in:
This gives us a solution: x = 170, y = 39.
The solutions to Pell's equation can sometimes be surprising; for example, the smallest solution in natural numbers for x² − 62y² = 1 is x = 63, y = 8; however, for x² − 61y² = 1 the smallest solutions are x = 1,766,319,049, y = 226,153,980. The smallest solutions in natural numbers for x² − Ny² = 1 for every value of N between 2 and 99 are listed in the table below; you'll notice a few other unusually large numbers there as well.
Sources used (see bibliography page for titles corresponding to numbers): 20.