views:

2371

answers:

4

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?

+3  A: 

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.

lsalamon
Thanks for interesting libraries but my question was really on a specific problem... The algorithm and math rather than the tools.
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.
mstrobl
A: 

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

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)
Dana
I haven't solved it yet... I can find enough solutions but they aren't always the smallest 150,000.
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
Rex Logan
A: 

I think you need a criterion that if you have checked all "small" pairs of numbers you have found all Alexandrian numbers up to a certain bound. Then you need to check enough pairs so that you have 150000 Alexandrian numbers below the bound.

A concrete bound is the following: Suppose you have checked all pairs a,b with ab <= n. Then for a given unchecked pair a,b with ab>n we get A = ab(ab-1)/(a+b) > n(n-1)/(n+1).

starblue
+2  A: 

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.

Mitch Wheat
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?
Dimitre Novatchev
"Unfortunately, pqr doesn't increase monotonically with the increase of p" - could you explain please?
Mitch Wheat
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.
Dimitre Novatchev
The point is that I have decomposed the problem so you only have to iterate over p
Mitch Wheat
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 p*q*r (not based on p). And the problem requires to produce *exactly* the 150000th number based on the value of p*q*r.
Dimitre Novatchev
In N numbers when increasing p, that may/may-not have the 150000th number of p*q*r -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
Dimitre Novatchev