ABC Conjecture
September 18, 2012
Our rad function combines factorization by trial division with a uniq auxiliary function that only adds a factor if it is distinct:
(define (rad n)
(define (uniq x xs)
(if (null? xs) (list x)
(if (= (car xs) x) xs
(cons x xs))))
(let twos ((n n) (fs (list)))
(if (even? n) (twos (/ n 2) (uniq 2 fs))
(if (= n 1) 2
(let odds ((n n) (f 3) (fs fs))
(cond ((< n (* f f))
(apply * (uniq n fs)))
((zero? (modulo n f))
(odds (/ n f) f (uniq f fs)))
(else (odds n (+ f 2) fs))))))))
The q function is simple:
(define (q a b c)
(if (and (= (gcd a b) 1)
(= (gcd b c) 1)
(= (gcd c a) 1))
(/ (log c) (log (rad (* a b c))))
(error 'q "must be coprime")))
To find all the (a,b,c) triples less than a limit, we can either consider all possible pairs of a and b, testing each for coprimality, or use a clever method to generate coprime pairs directly, eliminating the need for a test. We’ll be clever. If you start with the two pairs (2,1) and (3,1), then form two infinite ternary trees starting from each (m,n), each child of the form (2m−n,m), (2m+n,m) and (m+2n,n) forms a new coprime pair; the two trees thus formed are complete and distinct. We return the (a,b,c) triples sorted by descending q:
(define (abc limit)
(define (q a b c)
(/ (log c) (log (rad (* a b c)))))
(define (f a b)
(let* ((c (+ a b)) (qabc (q a b c)))
(if (< limit c) (list)
(append (if (< qabc 1) (list)
(list (list qabc b a c)))
(f (- (* 2 a) b) a)
(f (+ (* 2 a) b) a)
(f (+ a (* 2 b)) b)))))
(sort (lambda (xs ys) (< (car ys) (car xs)))
(append (f 2 1) (f 3 1))))
There are 31 triples q>1 with c less than a thousand:
> (abc 1000)
((1.4265653296335432 3 125 128)
(1.3175706029138485 1 512 513)
(1.311101219926227 1 242 243)
(1.292030029884618 1 80 81)
(1.2727904629543532 13 243 256)
(1.226294385530917 1 8 9)
(1.2251805398372944 1 288 289)
(1.2039689893561185 49 576 625)
(1.198754152359422 169 343 512)
(1.1757189916348774 32 49 81)
(1.1366732909394883 25 704 729)
(1.1126941404922133 1 63 64)
(1.1084359145429683 32 343 375)
(1.1048460623308174 104 625 729)
(1.0921945706283556 1 675 676)
(1.0917548251330267 100 243 343)
(1.0790468171894105 1 624 625)
(1.0458626417646626 1 728 729)
(1.0456203807611666 5 507 512)
(1.0412424573518235 1 48 49)
(1.0370424407259895 81 175 256)
(1.0344309549548174 343 625 968)
(1.0326159011595437 81 544 625)
(1.0326070471078377 7 243 250)
(1.0288287974277142 2 243 245)
(1.027195810121916 4 121 125)
(1.0251241218312794 27 512 539)
(1.018975235452531 5 27 32)
(1.0129028397298183 1 224 225)
(1.0084113092374762 200 529 729)
(1.004797211020379 1 960 961))
You can run the program at http://programmingpraxis.codepad.org/fmbvssfd.
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