ABC Conjecture
September 18, 2012
The ABC Conjecture has recently been in the news on math blogs because of the claim that it has been proved by Shinichi Mochizuki. Though the proof is being taken seriously, due to Mochizuki’s reputation, it is five hundred pages long, and confirmation will take several months. The conjecture states that given two positive integers a and b and their sum c with no common factors, the product of the distinct prime factors of abc is rarely much smaller than c.
The radical of a number n is the product of the distinct prime factors of n; for instance, rad(18) = 6 because 18 = 2 × 3 × 3 and, eliminating the duplicate occurrence of 3, 2 × 3 = 6. The quality of an (a,b,c) triple is given by q(a,b,c) = log(c) / log(rad(abc)). The precise statement of the ABC conjecture is that for every ε > 0 there are only finitely many triples (a,b,c) with a and b coprime positive integers and a + b = c such that rad(abc)1+ε < c, or equivalently, such that q(a,b,c) > 1.
Your task is to write the functions rad and q and find all the triples with c < 1000 for which q(a,b,c) > 1. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
A python version. I left out code for gcd and prime factors
def rad(n): return reduce(operator.mul, set(ma.prime.factors3(n))) def q(a, b, c): return math.log(c) / math.log(rad(a * b * c)) def find_abc(limit=1000): pfac = [0] + [set(ma.prime.factors3(n)) for n in range(1, limit)] num = 0 for a in range(1, limit - 2): for b in range(a + 1, limit - a): c = a + b if c < limit and not pfac[a] & pfac[b]: qabc = q(a, b, c) if qabc > 1: num += 1 print "{0:3d} {1:5.3f} {2:3d} {3:3d} {4:3d}".format(num, qabc, a, b, c) find_abc(1000)Haskell version of the suggested solution. Uses primes and factorization algorithm taken from the haskellwiki
import Data.Function (fix) import Data.List (nub, sortBy, sort) primesTMEf = 2 : g (fix g) where g xs = 3 : gaps 5 (joinT [[x*x, x*x+2*x..] | x <- xs]) joinT ((x:xs):t) = x : union xs (joinT (pairs t)) pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t gaps k s@(x:xs) | k<x = k:gaps (k+2) s -- equivalent to | True = gaps (k+2) xs -- [k,k+2..]`minus`s, k<=head s union (x:xs) (y:ys) = case (compare x y) of LT -> x : union xs (y:ys) EQ -> x : union xs ys GT -> y : union (x:xs) ys union xs [] = xs union [] ys = ys primeFactors n = factor primesTMEf n where factor ps@(p:pt) n | p*p > n = [n] | rem n p == 0 = p : factor ps (quot n p) | otherwise = factor pt n rad = product . nub . primeFactors q a b c = log (fromIntegral c) / log (fromIntegral $ rad (a * b * c)) triples :: Int -> [(Int, Int, Int)] triples limit = triples' [(2,1),(3,1)] where triples' [] = [] triples' ((a,b):xs) | c <= limit = (a,b,c) : triples' ([(2 * a - b, a), (2 * a + b, a), (a + 2 * b, b)] ++ xs) | otherwise = triples' xs where c = a + b main = putStr $ unlines . map show . sort $ [(qval, a, b, c) | (a, b, c) <- triples 1000, let qval = q a b c, qval > 1][…] make a long story shorter, there was a challenge on Programming Praxis that intrigued me, which was to write code that given a upper bound on c would generate a list of […]
A Java solution:
import java.util.*; public class ABCConjecture { public static void main(String[] args) { List<ABCTriple> triples = abc(1000); Collections.sort(triples, new Comparator<ABCTriple>() { public int compare(ABCTriple t1, ABCTriple t2) { if (t1.q < t2.q) return 1; return -1; } }); for (ABCTriple t : triples) { System.out.println(t); } } public static List<ABCTriple> abc(final int n) { List<ABCTriple> result = new LinkedList<ABCTriple>(); for (int a = 1; a <= n - 1; a++) { for (int b = a + 1; b <= n; b++) { int c = a + b; if (c < n && coPrime(a, b, c)) { double q = q(a, b, c); if (q > 1.0) { result.add(new ABCTriple(a, b, c, q)); } } } } return result; } public static double q(int a, int b, int c) { return Math.log(c) / Math.log(rad(a * b * c)); } public static int rad(int n) { return product(unique(primeFactors(n))); } private static boolean coPrime(Integer... numbers) { for (int i = 0; i < numbers.length - 1; i++) { for (int j = i + 1; j < numbers.length; j++) { if (gcd(numbers[i], numbers[j]) != 1) return false; } } return true; } private static int gcd(int a, int b) { int max = Math.max(a, b); int min = Math.min(a, b); while (min != 0) { int r = max % min; max = min; min = r; } return max; } private static int product(Collection<Integer> numbers) { int result = 1; for (int i : numbers) { result *= i; } return result; } private static Set<Integer> unique(Collection<Integer> numbers) { return new HashSet<Integer>(numbers); } private static List<Integer> primeFactors(int n) { List<Integer> primeFactors = new LinkedList<Integer>(); for (int p : PRIMES) { if (p * p > n) break; while (n % p == 0) { primeFactors.add(p); n /= p; } } if (n > 1) primeFactors.add(n); return primeFactors; } private static List<Integer> PRIMES = sieve(1000); private static List<Integer> sieve(int n) { List<Integer> primes = new LinkedList<Integer>(); BitSet sieve = new BitSet(n + 1); for (int i = 2; i * i <= n; i = sieve.nextClearBit(i + 1)) { primes.add(i); for (int j = i + i; j <= n; j += i) { sieve.set(j); } } return primes; } } class ABCTriple { final int a; final int b; final int c; final double q; ABCTriple(int a, int b, int c, double q) { this.a = a; this.b = b; this.c = c; this.q = q; } public String toString() { return q + " [" + a + ", " + b + ", " + c + "]"; } }I went with quick and easy, using Racket’s excellent
for*/listmacro to generate the final step. Interestingly, although it just brute forces the problem rather than the “smarter” solution posted here, it actually runs about 20% faster on the n=1000 case than the code posted here (on the same machine and the same environment). Not entirely sure why, although I am curious if anyone has any thoughts.In any case, you can see my version on my blog.
This algorithm is still to slow to solve http://projecteuler.net/problem=127