Offtonic Pell Solver has been released!

(Part 2)

For my inaugural instructional post, I want to tackle a mathematical problem: solving Pell-type equations (with a computer, not by hand).  What do I mean?  I mean equations that look like this:

x^2 – Dy^2 = N

where D and N are integers and x and y are nonnegative integers.  This makes it a quadratic Diophantine equation, and as it turns out, they’re not so easy to solve.  I came across some of these critters in a Project Euler problem a while ago (no, I won’t tell you which one), and I looked up how to solve them.  Well, maybe my Google-fu is not strong, but I didn’t find anything easy to understand.  There were a few incomplete references, a few references that didn’t explain much, and some papers with quite a bit of theory.  I eventually found this paper (Solving the generalized Pell equation x^2 − Dy^2 = N, John P. Robertson, 2004) that made sense, and I distilled it into an algorithm.  This method is mostly based on Robertson’s.

Part 1 will only cover the trivial cases.  The interesting case lives over at Part 2.

First, if N is 0, the equation reduces to x^2 = Dy^2.  (0,0) is always a solution.  If D is a perfect square, with d^2 = D, (dt,t) is a solution for all nonnegative integers t; if D is not a perfect square (including if it’s negative), (0,0) is the only solution.

Second, if D is 0, then the equation is x^2 = N; if N is a perfect square, with n^2 = N, then (n,t) is a solution for all nonnegative integers t, and if N is not a perfect square, there are no solutions.

So far, this is all very trivial, but I present it for completeness.

Third, if D is negative, the sign changes: x^2 + |D|y^2 = N.  There’s a fairly simple way to find all solutions here, and it runs in O(sqrt(N/|D|)) time: x^2 = N – |D|y^2.  Set y = 0; check if N – |D|y^2 is a perfect square, and if it is, add (x,y) to the list of solutions; increment y, and stop when |D|y^2 > N.  There’s probably a faster way to do this, but brute force works just fine here.  It goes without saying that for negative N there are no solutions!

Fourth, if D is a perfect square d^2, the equation factors into (x – dy)(x + dy) = N.  If N is positive and factors into p*q, with p ≤ q, if we can find x and y such that x – dy = p and x + dy = q has an integer solution, we’ve found a solution to the equation.  Solving, x = (q + p)/2 and y = (q – p)/(2d).  Factoring an arbitrary number N is a difficult problem when N is very large, but it’s at least doable, and most likely we’ll want to focus on not so large numbers anyway.  So we factor N, come up with all its factor pairs, and check whether p and q have the same parity and q – p is divisible by 2d.  If N is odd, p and q have the same parity already, so we don’t need to check that.  Otherwise, p and q need to both be even, which means N must be divisible by 4.  The same parity condition is satisfied when p = 2p’ and q = 2q’, where p’*q’ = N/4, so you can factor N/4 instead and check whether the difference of its factor pairs is divisible by d (because (q – p)/(2d) = (2q’ – 2p’)/(2d) = (q’ – p’)/d).

To summarize, if N is odd, factor it and obtain all of its factor pairs (p,q) where p ≤ q.  If (q – p) is divisible by 2d, ((q + p)/2,(q – p)/(2d)) is a solution.  If N is divisible by 4, let N’ = N/4; factor it and obtain all of its factor pairs (p’,q’) where p’ ≤ q’.  If (q’ – p’) is divisible by d, (q’ + p’,(q’ – p’)/d) is a solution.  If N is even but not divisible by 4, there are no solutions.

If N is negative, things are a bit reversed: (dy – x)(dy + x) = |N|.  The same method still works, almost.  If N is odd, factor |N| and obtain all of its factor pairs (p,q) where p ≤ q.  If (q + p) is divisible by 2d, ((q – p)/2,(q + p)/(2d)) is a solution.  If N is divisible by 4, let N’ = |N|/4 and obtain all of its factor pairs (p’,q’), where p’ ≤ q’.  If (q’ + p’) is divisible by d, (q’ – p’,(q’ + p’)/d) is a solution.  If N is even but not divisible by 4, there are no solutions.

Cases 3 and 4 were less trivial than cases 1 and 2, but still quite trivial.  This, by the way, is a good time to talk about factoring.  Factoring a number doesn’t have to be so hard.  Here’s a relatively fast algorithm to do it.  If you’re not interested in the factoring algorithm, you can skip right on to Part 2!

***

We start by computing a list of prime numbers.  (If we’re factoring many numbers, it would be good to keep the list somewhere so that it doesn’t have to be generated every time.)  If we want to find all the primes up to some number E (E for Eratosthenes), we start with a boolean array e[] of length E + 1 (bonus points if you can make each number take up a single bit), with every bit set to 0.  We set e = 1 and e = 1.  Now, starting from the beginning, we find the next index p such that e[p] = 0 and add that to our list of primes.  Then, we set e[p*p] = 1, e[(p + 1)*p] = 1, and so on until the end of the array, and find the next p with e[p] = 0.  When p*p ≤ E is no longer true, we don’t need to do that second step of setting things to 1; we just read 0’s and add them to our list of primes.  This finds all primes less than E rather quickly.

Now we test each of them against our number N that we’re trying to factor.  There’s a fairly quick way of doing this, too, that yields a prime factorization.  Let n = N.  If n % p is 0, we add 1 to our count of p factors of N, change n to n/p, and continue.  If n % p is not 0, then p becomes the next prime.  Stop when p^2 > n; then n is either prime or 1.  Note that to account for repeated primes, we don’t increment p if it’s a factor of n; instead, we divide n by p.  This way, we end up with a list of primes factors and their exponents.  For example, if the number N is 225225, we end up with ((3,2),(5,2),(7,1),(11,1),(13,1)).  We also get the benefit of n becoming smaller while p becomes bigger as we find more factors.

From here, we’d like to be able to get a list of all factors in order to get the factor pairs (note that is N = n^2, then n appears only once as a factor but the factor pair is (n,n)).  If N = p1^a1 p2^a2 …, then a factor r of N will be p1^b1 p2^b2 …, where b1 ≤ a1, b2 ≤ a2, … .  There’s a simple way to generate these recursively, but recursion is annoying(ly slow), so there’s a clever way to generate these iteratively instead.  Start with an array of 0’s where the length is the number of distinct prime factors; these are the prime exponents of the factors you’re finding.  When you want to step the array, start with a pointer to the first element.  Increment the element at the pointer.  If it’s ≤ the exponent in the original number, return the array.  If it’s bigger than the exponent in the original number, set that element to 0, advance the pointer, and repeat.  If you can’t advance the pointer anymore, return a null array or something because you’re done.  Now, with each array of exponents, you can get a factor.  With the number 225225 above, the array (0,1,0,1,0) would correspond to 3^0 5^1 7^0 11^1 13^0 = 55.  Here are the first few steps:

0 0 0 0 0 = 1

1 0 0 0 0 = 3

2 0 0 0 0 = 9

0 1 0 0 0 = 5

1 1 0 0 0 = 15

2 1 0 0 0 = 45

0 2 0 0 0 = 25

1 2 0 0 0 = 75

2 2 0 0 0 = 225

0 0 1 0 0 = 7

1 0 1 0 0 = 21

…and so on.  Anyway, with each of these we have a factor, so we sort that list of factors.  To find each factor pair, pick an element p from the factor list and let q = N/p.  Stop when q < p, and you’ve found all the factor pairs.

Case 5 is where it gets interesting…

…in the next part.  Feel free to ask any questions in the comments!

(Part 2)