title: Finding Integers With A Certain Property - Project Euler Problem 221
body: <p>I've become very addicted to Project Euler recently and am trying to do <a href="http://projecteuler.net/index.php?section=problems&amp;id=221" rel="nofollow noreferrer">this</a> one next! I've started some analysis on it and have reduced the problem down substantially already. Here's my working:</p>
  
  <blockquote>
    <p>A = pqr and</p>
    
    <p>1/A = 1/p + 1/q + 1/r    so    pqr/A =
    pq + pr + qr</p>
    
    <p>And because of the first equation:</p>
    
    <p>pq+pr+qr = 1</p>
    
    <p>Since exactly two of p, q and r have
    to be negative, we can simplify the
    equation down to finding:</p>
    
    <p>abc for which ab = ac+bc+1</p>
    
    <p>Solving for c we get:</p>
    
    <p>ab-1 = (a+b)c</p>
    
    <p>c = (ab-1)/(a+b)</p>
    
    <hr>
    
    <p>This means we need to find a and b for
    which:</p>
    
    <p>ab = 1 (mod a+b)</p>
    
    <p>And then our A value with those a and
    b is:</p>
    
    <p>A = abc = ab(ab-1)/(a+b)</p>
  </blockquote>
  
  <p>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.</p>
  
  <p>For ease I assumed a &lt; b and I have also noticed that gcd(a, b) = 1.</p>
  
  <p>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 <em>smallest</em> solutions. Here's the code anyway:</p>
  
  <pre><code>n = 150000
  seen = set()
  
  a = 3
  while len(seen) &lt; 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
  </code></pre>
  
  <p>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...</p>
  
  <p>So I'm absolutely out of ideas! Anyone got any ideas?</p>