Finding Integers With A Certain Property - Project Euler Problem 221

Question

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?

Answer 1

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.

Answer 2

This article about Chinese remainder, fast implementation, can help you : www.codeproject.com/KB/recipes/CRP.aspx

This is more links for tools and libraries :

Tools:

Maxima http://maxima.sourceforge.net/ Maxima is a system for the manipulation of symbolic and numerical expressions, including differentiation, integration, Taylor series, Laplace transforms, ordinary differential equations, systems of linear equations, polynomials, and sets, lists, vectors, matrices, and tensors. Maxima yields high precision numeric results by using exact fractions, arbitrary precision integers, and variable precision floating point numbers. Maxima can plot functions and data in two and three dimensions.

Mathomatic http://mathomatic.org/math/ Mathomatic is a free, portable, general-purpose CAS (Computer Algebra System) and calculator software that can symbolically solve, simplify, combine, and compare equations, perform complex number and polynomial arithmetic, etc. It does some calculus and is very easy to use.

Scilab www.scilab.org/download/index_download.php Scilab is a numerical computation system similiar to Matlab or Simulink. Scilab includes hundreds of mathematical functions, and programs from various languages (such as C or Fortran) can be added interactively.

mathstudio mathstudio.sourceforge.net An interactive equation editor and step-by-step solver.

Library:

Armadillo C++ Library http://arma.sourceforge.net/ The Armadillo C++ Library aims to provide an efficient base for linear algebra operations (matrix and vector maths) while having a straightforward and easy to use interface.

Blitz++ http://www.oonumerics.org/blitz/ Blitz++ is a C++ class library for scientific computing

BigInteger C# http://msdn.microsoft.com/pt-br/magazine/cc163441.aspx

libapmath http://freshmeat.net/projects/libapmath Welcome to the homepage of the APMath-project. Aim of this project is the implementation of an arbitrary precision C++-library, that is the most convenient in use, this means all operations are implemented as operator-overloadings, naming is mostly the same as that of .

libmat http://freshmeat.net/projects/libmat MAT is a C++ mathematical template class library. Use this library for various matrix operations, finding roots of polynomials, solving equations, etc. The library contains only C++ header files, so no compilation is necessary.

animath http://www.yonsen.bz/animath/animath.html Animath is a Finite Element Method library entirely implemented in C++. It is suited for fluid-structure interaction simulation, and it is mathematically based on higher-order tetrahedral elements.

Answer 3

OK. Here's some more playing with my Chinese Remainder Theorem solution. It turns out that a+b cannot be the product of any prime, p, unless p = 1 (mod 4). This allows faster computation as we only have to check a+b which are multiples of primes such as 2, 5, 13, 17, 29, 37...

So here is a sample of possible a+b values:

[5, 8, 10, 13, 16, 17, 20, 25, 26, 29, 32, 34, 37, 40, 41, 50, 52, 53, 58, 61, 64, 65, 68, 73, 74, 80, 82, 85, 89, 97, 100]

And here is the full program using the Chinese Remainder Theorem:

cachef = {}
def factors(n):
    if n in cachef: cachef[n]

    i = 2
    while i*i <= n:
        if n%i == 0:
            r = set([i])|factors(n//i)
            cachef[n] = r
            return r

        i += 1

    r = set([n])
    cachef[n] = r
    return r

cachet = {}
def table(n):
    if n == 2: return 1
    if n%4 != 1: return

    if n in cachet: return cachet[n]

    a1 = n-1
    for a in range(1, n//2+1):
        if (a*a)%n == a1:
            cachet[n] = a
            return a

cacheg = {}
def extended(a, b):
    if a%b == 0:
        return (0, 1)
    else:
        if (a, b) in cacheg: return cacheg[(a, b)]

        x, y = extended(b, a%b)
        x, y = y, x-y*(a//b)

        cacheg[(a, b)] = (x, y)
        return (x, y)

def go(n):
    f = [a for a in factors(n)]
    m = [table(a) for a in f]

    N = 1
    for a in f: N *= a

    x = 0
    for i in range(len(f)):
        if not m[i]: return 0

        s, t = extended(f[i], N//f[i])
        x += t*m[i]*N//f[i]

    x %= N
    a = x
    while a < n:
        b = n-a
        if (a*b-1)%(a+b) == 0: return a*b*(a*b-1)//(a+b)
        a += N




li = [5, 8, 10, 13, 16, 17, 20, 25, 26, 29, 32, 34, 37, 40, 41, 50, 52, 53, 58, 61, 64, 65, 68, 73, 74, 80, 82, 85, 89, 97, 100]

r = set([6])
find = 6

for a in li:
    g = go(a)
    if g:
        r.add(g)
        #print(g)
    else:
        pass#print(a)


r = list(r)
r.sort()

print(r)
print(len(r), 'with', len(li), 'iterations')

This is better but I hope to improve it further (for example a+b = 2^n seem to never be solutions).

I've also started considering basic substitutions such as:

a = u+1 and b = v+1

ab = 1 (mod a+b)

uv+u+v = 0 (mod u+v+2)

However, I can't see much improvement with that...