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.
Programming Praxis 
[…] today’s Programming Praxis exercise, our goal is to prove that 1729 is the smallest number that can be […]
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[…] Pages: 1 2 […]
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)]
[…] had another programming puzzle by Programming Praxis. This time, we are looking for a very special sort of number, a Taxicab […]
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…
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; }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.
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 FalseMessed up the links. Here they (hopefully) are
David Wilson’s page on Taxicab numbers
Segmented code
Counting triples
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))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))]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))]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): ....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!
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;
}
}
}
}
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; } }Here’s a solution in Go
http://play.golang.org/p/oXbLYjL-J2
[…] Determine if a positive number can be expressed as a sum of two cubes? […]
def taxicab_numbers(n) (1..Math.cbrt(n)) .to_a .combination(2) .select { |a, b| ((a ** 3) + (b ** 3)) == n } end puts taxicab_numbers(ARGV[0].to_i).inspect