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?

Actinouranium

@PythonPower: please see my answer below, which shows how to parameterise A via p and the the divisors of p^2 + 1 – Glenda

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.

Glenda

Unfortunately, pqr doesn't increase monotonically with the increase of p. This means that we should probably generate much more than 150000 numbers (how many?), then sort them and finally take the 150000th from the sorted numbers. Is there any way to generate pqr in a monotonically increasing order? – Isotope

"Unfortunately, pqr doesn't increase monotonically with the increase of p" - could you explain please? – Glenda

Forexample: 1806: x=6, y=7, z=43. But 1428: x=7, y=12, z=17. x1 < x2, but x1y1z1 > x2y2z2. So, if we generate the first 150000 Alexandr. numbers by increasing x, there could exist a bigger x such that xyz is less than the 150000th number we calculated, thus the real 150000th num is something else. – Isotope

The point is that I have decomposed the problem so you only have to iterate over p – Glenda

The point is that one doesn't know how long to iterate over p -- obviously just generating 150000 numbers will not contain the 150000th number based on pqr (not based on p). And the problem requires to produce exactly the 150000th number based on the value of pqr. – Isotope

In N numbers when increasing p, that may/may-not have the 150000th number of pqr -s, regardless of how big N is, and even the 150000th pqr number is among those N p-numbers, how do we now which exactly it is? Even the p-numbers sorted, the 150000th pqr num may still not be one of the N p-numbers – Isotope

clever solution – Electrophilic

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.

Cantilena

Thanks for interesting libraries but my question was really on a specific problem... The algorithm and math rather than the tools. – Actinouranium

Still, the list of links seems to be carefully crafted and may be of interest to some. So I don't see why this should be -1, hence I will upvote. People seem too trigger-happy these days. – Euniceeunuch

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...

Actinouranium

Just out of curiousity: is your program too slow to find a solution at all, or are you just trying to beat the 1 minute rule? (I haven't even looked at this one yet :P) – Swearingen

I haven't solved it yet... I can find enough solutions but they aren't always the smallest 150,000. – Actinouranium

Have a look at your previous question I have posted a much quicker version of table() it runs about 3000 times faster than the one you used in the above code – Thymelaeaceous