Project Euler 127: Investigating the number of abc-hits below a given limit.

Problem 127 of Project Euler is a really awesome number theoretic problem using concepts of GCD and radicals. The problem reads

The radical of n, rad(n), is the product of distinct prime factors of n. For example, 504 = 2³ x 3² x 7, so rad(504) = 2 x 3 x 7 = 42.

We shall define the triplet of positive integers (a, b, c) to be an abc-hit if:

1. GCD(a, b) = GCD(a, c) = GCD(b, c) = 1
2. a < b
3. a + b = c
4. rad(abc) < c

For example, (5, 27, 32) is an abc-hit, because:

1. GCD(5, 27) = GCD(5, 32) = GCD(27, 32) = 1
2. 5 < 27
3. 5 + 27 = 32
4. rad(4320) = 30 < 32

It turns out that abc-hits are quite rare and there are only thirty-one abc-hits for c < 1000, with ∑c = 12523.

Find ∑c for c < 120000.

I will come up with two solutions for this problem. Where the second is actually a refinement of the first. But lets jump right to some a few key insights into the problem

Key insights

The first insight which is not really the key insight but one that eases my latter writing

We have defined the following relationship a + b = c. Then If and only if GCD(a,b) = 1 the it follows that GCD(a,c) = GCD(b, c) = 1

And we can actually prove this

Assume that GCD(a,b) = d. Then we have a positive integer d that divides both a and b such that we can write a = da’ and b = db’ and thus we can write c = a + b = da’ + db’ = d*(a’+b’). And this we will have that GCD(a,c) and GCD(b,c) will have at least a common factor of d.

If GCD(a,b) = 1 then we have no factors d which divides both numbers and GCD(a,c) =1 since otherwise we would have a factor d that divides a and c ,  c / d = (a+b) / d = a/d + b/d.  However since GCD(a,b) = 1 this factor d cannot exist. The same goes for GCD(b,c)

Hope you can follow me this long. Now to the key insight of the this problem

If GCD(a,b) = 1 then we have that rad(a) * rad(b) * rad(c) = rad(abc).

The radical of a number is the product of unique prime factors. Since a,b and c doesn’t have any common prime factors then each prime factor in each of the components of abc is unique and thus none of that is removed when we take rad(abc).

With this insight we can rewrite property 4 to rad(a) * rad(b) * rad(c) < c if GCD(a,b) = 1 and thus limit the size of radicals we need significantly.

First solution

In problem 124 we made a sieve to calculate all the radicals. We will use that again here in a bit simpler version since we don’t need to sort them afterwards.

Besides that we can rather easily fulfill property 2 and 3 by making two nested loops, and then we just need to check the two last properties. The code for all of this looks like

long result = 0;
long limit = 120000;

//Sieve all radicals
long[] radicals = new long[limit+1];
for (long i = 0; i < radicals.Length; i++) {
    radicals&#91;i&#93; = 1;
}

for (long i = 2; i < limit; i++) {
    if (radicals&#91;i&#93; == 1) {
        radicals&#91;i&#93; = i;

        for (long j = i + i; j < limit; j += i) {
            radicals&#91;j&#93; *= i;
        }
    }
}

//Check for the properties
for (long a = 1; a < limit; a++) {
    for (long b = a + 1; b &lt; limit-a; b++) {

        if (radicals&#91;a&#93; * radicals&#91;b&#93; * radicals&#91;a + b&#93; >= a+b) continue;
        if (GCD(a, b) != 1) continue;
        result += a + b;
    }
}

Pretty simple I think. The main problem is that it takes more than 30 seconds to execute the code. The problem wields the correct result

sum of c = 18407904

A bit more clever solution

We are checking a whole lot of useless combinations, since we are discarding all solutions where rad(a)rad(b)rad(c) > c.

If we have a sorted list of radicals we could list of radicals (by coincidence that is exactly what we have form problem 124) we could say something clever about that.

So in this version I will turn everything upside down. We will loop over c. Since we have the sorted radicals we will loop over them and use them as a.

So now that we have a and c all we need is to find b. However, the radical function is rad(1) = 1, but since a < b we know that rad(b) = 2 as a bare minimum. That means if rad(a) * rad(c) >c/2 we don’t really need to continue the loop over a, since the resulting b will never yield a valid abc hit.

The code looks like

Radical[] sortedRadicals = new Radical[limit];
Array.Copy(radicals, 1, sortedRadicals, 0, limit);
Array.Sort(sortedRadicals);

for (long c = 3; c <= limit; c++) { long radc = radicals[c].rad; long chalf = c / 2; foreach (var a in sortedRadicals){ if (a.rad*radicals[c].rad > chalf) break;
if (a.number >= chalf) continue;

long b = c-a.number;
if (a.rad * radicals[b].rad * radc < c && GCD(a.rad, radicals[b].rad) == 1) result += c; } } [/csharp] Where the radicals is a list of radicals generated as in 124. This is only a part of the code, you can see the whole source code here. I can execute this code in 228ms, which is much more feasible than before.