Taxicab Numbers

November 9, 2012

Since the two numbers x and y whose cubes are summed must both be between 0 and the cube root of n, one solution is an exhaustive search of all combinations of x and y. A better solution starts with x = 0 and y the cube root of n, when repeatedly makes a three-way decision: if x3 + y3 < n, increase x, if x3 + y3 > n, decrease y, or if x3 + y3 = n, report the success and continue the search for more:

(define (cubes n)
  (let loop ((x 0) (y (root3 n)) (zs (list)))
    (let ((s (+ (* x x x) (* y y y))))
      (cond ((< y x) (reverse zs))
            ((< s n) (loop (+ x 1) y zs))
            ((< n s) (loop x (- y 1) zs))
            (else (loop (+ x 1) (- y 1)
              (cons (list x y) zs)))))))

We compute the cube root using Newton’s method:

(define (root3 n)
  (let loop ((u n))
    (let* ((s u)
           (t (+ (+ s s) (quotient n (* s s))))
           (u (quotient t 3)))
      (if (< u s) (loop u) s))))

Here are all the “taxicab pairs” less than a hundred thousand:

> (let loop ((n 1))
    (when (< 1 (length (cubes n)))
    (display n) (display ": ")
    (display (cubes n)) (newline))
    (when (< n 100000) (loop (+ n 1))))
1729: ((1 12) (9 10))
4104: ((2 16) (9 15))
13832: ((2 24) (18 20))
20683: ((10 27) (19 24))
32832: ((4 32) (18 30))
39312: ((2 34) (15 33))
40033: ((9 34) (16 33))
46683: ((3 36) (27 30))
64232: ((17 39) (26 36))
65728: ((12 40) (31 33))

Here’s another example, which finds taxicab triples; this one ran all night before I interrupted it:

> (let loop ((n 1))
    (when (< 2 (length (cubes n)))
      (display n) (display ": ")
      (display (cubes n)) (newline))
    (loop (+ n 1)))
87539319: ((167 436) (228 423) (255 414))
119824488: ((11 493) (90 492) (346 428))
143604279: ((111 522) (359 460) (408 423))
175959000: ((70 560) (198 552) (315 525))

You can run the program at http://programmingpraxis.codepad.org/I6l4ojtJ, or learn more about taxicab numbers at A004999 and A001235.

About these ads

Pages: 1 2

18 Responses to “Taxicab Numbers”

  1. [...] today’s Programming Praxis exercise, our goal is to prove that 1729 is the smallest number that can be [...]

  2. My Haskell solution (see http://bonsaicode.wordpress.com/2012/11/09/programming-praxis-taxicab-numbers/ for a version with comments):

    import Text.Printf
    
    cubesums :: Integer -> [(Integer, Integer)]
    cubesums n = f 0 (round $ fromIntegral n ** (1/3)) where
        f x y = if y < x then [] else case compare (x^3 + y^3) n of
                    EQ -> (x,y) : f (x + 1) (y - 1)
                    LT -> f (x + 1) y
                    GT -> f x (y - 1)
    
    taxicab :: [(Integer, (Integer, Integer), (Integer, Integer))]
    taxicab = [(n,p,q) | n <- [1..99999], [p,q] <- [cubesums n]]
    
    main :: IO ()
    main = mapM_ (\(n,p,q) -> printf "%d : %s %s\n" n (show p) (show q)) taxicab
    
  3. Paul said

    My version in Python. This gives the resuls below in a few seconds. It does not scale very well. N cannot be as large as 10000 on my machine due to MemoryError.

    from collections import defaultdict
    
    N =700
    ntriples = 3 # for original Ramanujan problem set to 2
    
    results = defaultdict(list)
    for x in range(1, N):
        x3 = x ** 3
        for y in range(x+1, N):
            results[x3 + y ** 3].append((x, y))
    
    for key in sorted(results.keys()):
        if len(results[key]) == ntriples:
            print key, results[key]
    

    87539319 [(167, 436), (228, 423), (255, 414)]
    119824488 [(11, 493), (90, 492), (346, 428)]
    143604279 [(111, 522), (359, 460), (408, 423)]
    175959000 [(70, 560), (198, 552), (315, 525)]
    327763000 [(300, 670), (339, 661), (510, 580)]

  4. [...] had another programming puzzle by Programming Praxis. This time, we are looking for a very special sort of number, a Taxicab [...]

  5. jpverkamp said

    Here’s my crack at it in Racket: Taxicab numbers

    Had some floating point woes with (expt n 1/3), but once I worked those out it works great.

    Now it’d be interested to optimize a bit and see how long it takes to calculate Ta(6). Or better yet, Ta(7). At least according to the Wikipedia page, it hasn’t been done…

  6. danaj said

    As usual, a few examples in Perl. I have to say I enjoy seing the Haskell solutions, as they’re remarkably concise. My code is much less elegant, but it’s fun to work on. I hadn’t planned on spending time on this, but it kept nagging at me today, since it looked similar to work I did for Goldbach and additive variable length codes. Would it work just to plug in cubes instead of primes into a modified find_pairs function? Indeed it does.

    #!/usr/bin/env perl
    use strict;
    use warnings;
    
    my $lim   = $ARGV[0] || 100000;    # Search this far
    my $paira = $ARGV[1] || 2;         # display all numbers with this many pairs
    my @basis = map { $_*$_*$_ } (1 .. int($lim ** (1.0/3.0) + 1));
    
    foreach my $p (1 .. $lim) {
      my @pairs = _find_pairs(\@basis, $p, sub { sprintf("%d^3 + %d^3", $_[0], $_[1]) } );
      next unless scalar @pairs == $paira;
      print "$p  =  ", join(" = ", @pairs), "\n";
    }
    
    
    sub _find_pairs {
      my($p, $val, $formatsub) = @_;
      return unless $#$p >= 0;
    
      # Determine how far to look in the basis.
      my $maxbasis = 0;
      if ($#$p > 16 && $val > $p->[-1]) {   # binary search
        my $i = 0;
        my $j = scalar $#$p;
        while ($i < $j) {
          my $mid = int( ($i + $j) / 2 );
          if ($p->[$mid] < $val)  { $i = $mid+1; }
          else                    { $j = $mid;   }
        }
        $maxbasis = $i-1;
      } else {                              # linear search
        $maxbasis+=4 while exists $p->[$maxbasis+5] && $val > $p->[$maxbasis+4];
        $maxbasis++  while exists $p->[$maxbasis+1] && $val > $p->[$maxbasis];
      }
    
      my @pairs;
      my $i = 0;
      my $j = $maxbasis;
      my $pi = $p->[$i];
      my $pj = $p->[$j];
      while ($i <= $j) {
        my $sum = $pi + $pj;
        if    ($sum < $val) {  $pi = $p->[++$i];  }
        elsif ($sum > $val) {  $pj = $p->[--$j];  }
        else {
          push @pairs, $formatsub->($i+1, $j+1);
          $pi = $p->[++$i];
        }
      }
      @pairs;
    }
    

    Given a basis and a target value, _find_pairs returns a list of all the pairs that sum to the value. What I normally use this for is to find the smallest encoding of a pair for a given arbitrary number. In this example it works reasonably well for Ta(2), but running it to 500k or more shows it’s obviously slowing down, and it would be a long wait for Ta(3). Timing on my system for Ta(2) to 500k is basically the same as the two-liner Pari code on the OEIS page — about 11s.

    OK, what can we do to make this work? Similar to how Paul did it, we can flip the problem and calculate all the values instead of doing all the work for just one. My first solution was a simple double loop, using a hash for the summations to control memory with large inputs. It dropped the time for Ta(2)/500k from 11s to 0.02s. Ta(3) found in 0.2s. Ta(3) to 1000000000 in only a little over a second, so let’s try for Ta(4)! Oh. Way too much memory use. I went through a few iterations trimming memory use and speeding it up some, but it was still hitting ~32GB before getting to Ta(4). Obvious answer: segment. This is rather a brute force solution, but it manages to get Ta(4) in 19.5 minutes (searching to 7000000000000) while staying at a constant ~3.3GB of memory. Dropping the segment size can make it use much less memory at a little speed expense. Still too slow for practically looking for larger values, but far, far faster than what I started out with.

    #!/usr/bin/env perl
    use strict;
    use warnings;
    use Data::Dumper;
    
    my $lim   = $ARGV[0] || 100000;    # Search this far
    my $paira = $ARGV[1] || 2;         # display all numbers with this many pairs
    my @basis = map { $_*$_*$_ } (1 .. int($lim ** (1.0/3.0) + 1));
    
    {
      # 100B => 3.3GB, 10B => 800MB, 1B => 270MB
      my $segsize = 10_000_000_000;
      my $low;
      my $high = 0;
      do {
        $low = $high+1;
        $high = $low + $segsize - 1;
        $high = $lim if $high > $lim;
        foreach my $p (_find_pairs_segment(\@basis, $paira, $low, $high, sub { sprintf("%d^3 + %d^3", $_[0], $_[1]) })) {
          my $n = shift @$p;
          print "$n = ", join(" = ", @$p), "\n";
        }
        #printf "  $high  %3d%%\n", int(100*$high/$lim);
      } while ($high < $lim);
    }
    
    sub _find_pairs_segment {
      my($p, $len, $start, $end, $formatsub) = @_;
      my $plen = $#$p;
    
      my %allpairs;
      foreach my $i (0 .. $plen) {
        my $pi = $p->[$i];
        next if ($pi+$p->[$plen]) < $start;
        last if (2*$pi) > $end;
        foreach my $j ($i .. $plen) {
          my $sum = $pi + $p->[$j];
          next if $sum < $start;
          last if $sum > $end;
          push @{ $allpairs{$sum} }, $i, $j;
        }
        # If we wanted to save more memory, we could filter and delete every entry
        # where $n < 2 * $p->[$i+1].  This can cut memory use in half, but is slow.
      }
    
      my @retlist;
      foreach my $list (grep { scalar @$_ == $len*2 } values %allpairs) {
        my $n = $p->[$list->[0]] + $p->[$list->[1]];
        my @pairlist;
        while (@$list) {
          push @pairlist, $formatsub->(1 + shift @$list, 1 + shift @$list);
        }
        push @retlist, [$n, @pairlist];
      }
      @retlist = sort { $a->[0] <=> $b->[0] } @retlist;
      return @retlist;
    }
    
  7. Ankur Dhama said

    A bit of logic programming in core.logic:

    (defn cubo [a r] (fresh [a1] (*fd a a a1) (*fd a a1 r)))

    (defn taxi-cab [n]
    (run* [q] (fresh [a ac b bc] (infd a b (interval 1 (int (java.lang.Math/cbrt n)) )) (infd ac bc (interval 0 n)) (cubo a ac) (cubo b bc) (+fd ac bc n) (== q [a b]))))

    NOTE: FYI, could be made more efficient.

  8. Paul said

    I have had a lot of fun with this exercise. Like danaj I created a segmented version of my first entry. It can be found here. Also I created a fast method to count triples, based on David Wilson’s page on Taxicab numbers and used it to generate numbers with 4 triples using his combination and magnification methods. Maybe a topic for a next exercise?
    Here a method to generate the triples more in line with this exercise. The nice thing about this version is that it is possible to step up with x with a step of 2, as long you make sure that the sum of the triples of x and y are odd (even) if n is odd (even). Therefore this version is about 2 times faster compared to a version with a step of 1. The full version is here.

    def invpow(x, n):
        """ x , n > 0 integers
            return y, such that y ** n <= x and (y+1) ** n > x
        """
        return int((x + 0.5) ** (1.0/3))
    
    def cubes(n):
        """generator for all sums of 3rd powers to get n
           start with pair x, y, where y is the highest number, such that y**3<=n
           If n is odd, (x,y) can only be an odd-even pair. If n is even, (x,y) can
           only be odd-odd or even-even. While x^3+y^3 is smaller than n, x can
           advance with 2 steps i.s.o one step, as a one step increase cannot be a
           solution. If y is decreased x must increase to make sure that the sum can
           be a solution (is odd or even, like n).
           This version is about a factor 2 faster than cubes0 (with step of 1)
        """ 
        y = invpow(n, 3)
        y3 = y ** 3
        x = (n - y) & 1 # make x odd if n - y^3 is odd
        while x <= y:
            test = x ** 3 + y3
            if test < n:
                x += 2
            else:
                if test == n:
                    yield x, y
                x += 1
                y -= 1
                y3 = y ** 3
                
    def check_cubes(n, ncubes):
        """check if number of cubes is ncubes for number n"""
        cub = list(itertools.islice(cubes(n), ncubes+1))
        return cub if len(cub) == ncubes else False
    
  9. Another Python solution:

    rom math import pow
    
    def cubic_root(x):
        return int(round(pow(x, 1.0 / 3)))
    
    def sum_of_cubes(x):
        result = []
        i = 1
        j = cubic_root(x)
        while i < j:
            s = i ** 3 + j ** 3
            if s == x:
                result.append((i, j))
                i += 1
                j -= 1
            elif s < x:
                i += 1
            else:
                j -= 1
        return result
    
    print(sum_of_cubes(1729))
    
  10. Mike Zraly said

    Another python solution. This one precomputes cubes and cube roots of interest, but only calls math.pow() once.

    import math
    import sys
    
    def cube_pairs(n):
        """Yield pairs of integers whose cubes sum to n."""
    
        cubes = []
        cbrts = {}
        for i in xrange(int(math.floor(math.pow(n, 1.0/3.0)))):
    	root = i + 1
    	cube = root ** 3
    	cubes.append(cube)
    	cbrts[cube] = root
    
        for cube in cubes:
    	if cube + cube > n:
    	    continue
    	diff = n - cube
    	if diff in cbrts:
    	    yield(cbrts[cube], cbrts[diff])
    
    if __name__ == "__main__":
        for arg in sys.argv[1:]:
    	print
    	print arg, ":", [pair for pair in cube_pairs(int(arg))]
    
  11. Mike Zraly said

    Whoops! Let’s remove tabs first for proper formatting:

    import math
    import sys
    
    def cube_pairs(n):
        """Yield pairs of integers whose cubes sum to n."""
    
        cubes = []
        cbrts = {}
        for i in xrange(int(math.floor(math.pow(n, 1.0/3.0)))):
            root = i + 1
            cube = root ** 3
            cubes.append(cube)
            cbrts[cube] = root
    
        for cube in cubes:
            if cube + cube > n:
                continue
            diff = n - cube
            if diff in cbrts:
                yield(cbrts[cube], cbrts[diff])
    
    if __name__ == "__main__":
        for arg in sys.argv[1:]:
            print
            print arg, ":", [pair for pair in cube_pairs(int(arg))]
    
  12. Mike Zraly said

    Oops, lines 9 and 10 in my previous comment can be collapsed into one line by using a two-argument version of xrange.

    for root in xrange(1, xrange(int(math.floor(math.pow(n, 1.0/3.0)))) + 1):
        ....
    
  13. Here’s a PHP script that checks all the numbers up to 2000 thus confirming the postulation.

    http://pastebin.com/kZZj1jnD

    I also tried it up to 1 million.

    1729 has 2 pairs found
    4104 has 2 pairs found
    13832 has 2 pairs found
    20683 has 2 pairs found
    32832 has 2 pairs found
    39312 has 2 pairs found
    40033 has 2 pairs found
    46683 has 2 pairs found
    64232 has 2 pairs found
    65728 has 2 pairs found

    This is my first attempt on here so be gentle!

  14. Maithily said

    Here is a Java Solution:

    public static void isTaxicabNumber(int inputNumber, int order)
    {
    String s;

    int upperLimit = inputNumber-1;

    for ( int i = 1; i inputNumber)
    {
    break;
    }
    for (int j = i+1; i inputNumber)
    {
    break;
    }
    if (inputNumber == Math.pow(i,3) + Math.pow(j,3))
    {
    System.out.println(inputNumber +" is a Taxicab Number. Pair :(" +i+","+j+")");
    break;
    }

    }
    }
    }

  15. Colt Jones said

    Here is a PHP script to test a given number. Took me a few iterations to come up with this. I had a few less than optimal solutions that were O(n^2) type solutions. This script can take up a good amount of memory for large numbers. I.E. Ta(5)

    Pastebin Version: taxicabNumberTest

    #!/usr/bin/php
    <?php
    
    if ($argc != 2) {
       echo "Usage cube.php <number>".PHP_EOL;
       exit;
    }
    $num = $argv[1];
    
    echo "Finding cubes less than ".$num.PHP_EOL;
    
    $factor = 1;
    $arr = array();
    do {
       $cube = pow($factor,3);
       if ($cube <= $num) {
          $arr[$factor] = $cube;
       }
       $factor++;
    } while ($cube <= $num);
    
    //Find the "top" items
    $bottom = array();
    $top = array();
    $half = floor($num/2);
    foreach ($arr as $factor => $product) {
       if ($product <= $half) {
          $bottom[$product] = $factor;
       } else {
          $top[$factor] = $product;
       }
    }
    
    echo "Found ".count($top)." cubes > lower half of list".PHP_EOL;
    echo "Found ".count($bottom)." cubes <= lower half of list".PHP_EOL;
    
    foreach ($top as $factor => $product) {
       $remainder = $num-$product;
       if (isset($bottom[$remainder])) {
          echo "Found ".$factor."^3 + ".$bottom[$remainder]."^3 = ".$num.PHP_EOL;
       }
    }
    

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 634 other followers

%d bloggers like this: