## 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 *x*^{3} + *y*^{3} < *n*, increase *x*, if *x*^{3} + *y*^{3} > *n*, decrease *y*, or if *x*^{3} + *y*^{3} = *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.

[…] 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):

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

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.

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.

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.

Messed up the links. Here they (hopefully) are

David Wilson’s page on Taxicab numbers

Segmented code

Counting triples

Another Python solution:

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

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

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

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

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? […]