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 = 23 x 32 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:
- GCD(a, b) = GCD(a, c) = GCD(b, c) = 1
- a < b
- a + b = c
- rad(abc) < c
For example, (5, 27, 32) is an abc-hit, because:
- GCD(5, 27) = GCD(5, 32) = GCD(27, 32) = 1
- 5 < 27
- 5 + 27 = 32
- 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[i] = 1; } for (long i = 2; i < limit; i++) { if (radicals[i] == 1) { radicals[i] = i; for (long j = i + i; j < limit; j += i) { radicals[j] *= i; } } } //Check for the properties for (long a = 1; a < limit; a++) { for (long b = a + 1; b < limit-a; b++) { if (radicals[a] * radicals[b] * radicals[a + b] >= 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.
The race for the abc-hits is always very active
De beste abc-drietallen van vandaag
http://www.tammo80.nl/java/abc/
___
Table 1
Re: Project Euler 127
Number of abc-hits below 1000
http://img823.imageshack.us/img823/8215/pe127tableone.jpg
___
Sources:
1) Project Euler 127
http://www.mathblog.dk/project-euler-127-abc-hits/
2) Good abc Triples
http://www.phfactor.net/abc/index.php
3) The abc conjecture (also known as Oesterlé–Masser conjecture) is
a conjecture in number theory, first proposed by Joseph Oesterlé (1988)
and David Masser (1985) as an integer analogue of the
Mason–Stothers theorem for polynomials. […]
http://en.wikipedia.org/wiki/Abc_conjecture
4) THE ABC CONJECTURE HOME PAGE
http://www.math.unicaen.fr/~nitaj/abc.html
You can speed it up even more by pre-calculating C/rad(C)
Then, instead of testing rad(A)*rad(C) < C/2, you can test rad(A) < C/rad(C)
Mine executed in 44 ms using this test.