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 (2mn,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.

About these ads

Pages: 1 2

6 Responses to “ABC Conjecture”

  1. Paul said

    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)
    
  2. DGel said

    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]
    
  3. [...] 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 [...]

  4. 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 + "]";
        }
    }
    
  5. JP said

    I went with quick and easy, using Racket’s excellent for*/list macro 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.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 600 other followers

%d bloggers like this: