Finding Integers With A Certain Property - Project Euler Problem 221

I've become very addicted to Project Euler recently and am trying to do this one next! I've started some analysis on it and have reduced the problem down substantially already. Here's my working:

A = pqr and

1/A = 1/p + 1/q + 1/r so pqr/A = pq + pr + qr

And because of the first equation:

pq+pr+qr = 1

Since exactly two of p, q and r have to be negative, we can simplify the equation down to finding:

abc for which ab = ac+bc+1

Solving for c we get:

ab-1 = (a+b)c

c = (ab-1)/(a+b)

---

This means we need to find a and b for which:

ab = 1 (mod a+b)

And then our A value with those a and b is:

A = abc = ab(ab-1)/(a+b)

Sorry if that's a lot of math! But now all we have to deal with is one condition and two equations. Now since I need to find the 150,000th smallest integer written as ab(ab-1)/(a+b) with ab = 1 (mod a+b), ideally I want to search (a, b) for which A is as small as possible.

For ease I assumed a < b and I have also noticed that gcd(a, b) = 1.

My first implementation is straight forward and even finds 150,000 solutions fast enough. However, it takes far too long to find the 150,000 smallest solutions. Here's the code anyway:

n = 150000
seen = set()

a = 3
while len(seen) < n:
    for b in range(2, a):
        if (a*b)%(a+b) != 1: continue

        seen.add(a*b*(a*b-1)//(a+b))
        print(len(seen), (a, b), a*b*(a*b-1)//(a+b))

    a += 1

My next thought was to use Stern-Brocot trees but that is just too slow to find solutions. My final algorithm was to use the Chinese remainder theorem to check if different values of a+b yield solutions. That code is complicated and although faster, it isn't fast enough...

So I'm absolutely out of ideas! Anyone got any ideas?

As with many of the Project Euler problems, the trick is to find a technique that reduces the brute force solution into something more straight forward:

A = pqr and 
1/A = 1/p + 1/q + 1/r

So,

pq + qr + rp = 1  or  -r = (pq - 1)/(p + q)

Without loss of generality, 0 < p < -q < -r

There exists k , 1 <= k <= p

-q = p + k
-r = (-p(p + k) – 1) / (p + -p – k)  = (p^2 + 1)/k + p

But r is an integer, so k divides p^2 + 1

pqr = p(p + q)((p^2 + 1)/k + p)

So to compute A we need to iterate over p, and where k can only take values which are divisors of p squared plus 1.

Adding each solution to a set, we can stop when we find the required 150000th Alexandrian integer.