Taxicab numbers

From Rosetta Code
Revision as of 13:12, 27 March 2014 by rosettacode>Bearophile (+ intermediate D version)
Taxicab numbers is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

A taxicab number (the definition that is being used here) is a positive integer that can be expressed as the sum of two positive cubes in more than one way.

The first taxicab number is :

Task requirements
  • Compute and display the lowest 25 taxicab numbers (in numeric order, and in a human-readable format).
    • For each of the taxicab numbers, show the number as well as it's constituent cubes.
  • Show the 2,000th taxicab number and a half dozen more (extra credit)
See also

C

Using a priority queue to emit sum of two cubs in order. It's reasonably fast and doesn't use excessive amount of memory (the heap is only at 245 length upon the 2006th taxi). <lang c>#include <stdio.h>

  1. include <stdlib.h>

typedef unsigned long long xint; typedef unsigned uint; typedef struct { uint x, y; // x > y always xint value; } sum_t;

xint *cube; uint n_cubes;

sum_t *pq; uint pq_len, pq_cap;

void add_cube(void) { uint x = n_cubes++; cube = realloc(cube, sizeof(xint) * (n_cubes + 1)); cube[n_cubes] = (xint) n_cubes*n_cubes*n_cubes; if (x < 2) return; // x = 0 or 1 is useless

if (++pq_len >= pq_cap) { if (!(pq_cap *= 2)) pq_cap = 2; pq = realloc(pq, sizeof(*pq) * pq_cap); }

sum_t tmp = (sum_t) { x, 1, cube[x] + 1 }; // upheap uint i, j; for (i = pq_len; i >= 1 && pq[j = i>>1].value > tmp.value; i = j) pq[i] = pq[j];

pq[i] = tmp; }

void next_sum(void) { redo: while (!pq_len || pq[1].value >= cube[n_cubes]) add_cube();

sum_t tmp = pq[0] = pq[1]; // pq[0] always stores last seen value if (++tmp.y >= tmp.x) { // done with this x; throw it away tmp = pq[pq_len--]; if (!pq_len) goto redo; // refill empty heap } else tmp.value += cube[tmp.y] - cube[tmp.y-1];

uint i, j; // downheap for (i = 1; (j = i<<1) <= pq_len; pq[i] = pq[j], i = j) { if (j < pq_len && pq[j+1].value < pq[j].value) ++j; if (pq[j].value >= tmp.value) break; } pq[i] = tmp; }

uint next_taxi(sum_t *hist) { do next_sum(); while (pq[0].value != pq[1].value);

uint len = 1; hist[0] = pq[0]; do { hist[len++] = pq[1]; next_sum(); } while (pq[0].value == pq[1].value);

return len; }

int main(void) { uint i, l; sum_t x[10]; for (i = 1; i <= 2006; i++) { l = next_taxi(x); if (25 < i && i < 2000) continue; printf("%4u:%10llu", i, x[0].value); while (l--) printf(" = %4u^3 + %4u^3", x[l].x, x[l].y); putchar('\n'); } return 0; }</lang>

Output:
   1:      1729 =   12^3 +    1^3 =   10^3 +    9^3
   2:      4104 =   15^3 +    9^3 =   16^3 +    2^3
   3:     13832 =   20^3 +   18^3 =   24^3 +    2^3
   4:     20683 =   27^3 +   10^3 =   24^3 +   19^3
   5:     32832 =   30^3 +   18^3 =   32^3 +    4^3
   6:     39312 =   33^3 +   15^3 =   34^3 +    2^3
   7:     40033 =   33^3 +   16^3 =   34^3 +    9^3
   8:     46683 =   30^3 +   27^3 =   36^3 +    3^3
   9:     64232 =   36^3 +   26^3 =   39^3 +   17^3
  10:     65728 =   33^3 +   31^3 =   40^3 +   12^3
  11:    110656 =   40^3 +   36^3 =   48^3 +    4^3
  12:    110808 =   45^3 +   27^3 =   48^3 +    6^3
  13:    134379 =   43^3 +   38^3 =   51^3 +   12^3
  14:    149389 =   50^3 +   29^3 =   53^3 +    8^3
  15:    165464 =   48^3 +   38^3 =   54^3 +   20^3
  16:    171288 =   54^3 +   24^3 =   55^3 +   17^3
  17:    195841 =   57^3 +   22^3 =   58^3 +    9^3
  18:    216027 =   59^3 +   22^3 =   60^3 +    3^3
  19:    216125 =   50^3 +   45^3 =   60^3 +    5^3
  20:    262656 =   60^3 +   36^3 =   64^3 +    8^3
  21:    314496 =   66^3 +   30^3 =   68^3 +    4^3
  22:    320264 =   66^3 +   32^3 =   68^3 +   18^3
  23:    327763 =   58^3 +   51^3 =   67^3 +   30^3
  24:    373464 =   60^3 +   54^3 =   72^3 +    6^3
  25:    402597 =   61^3 +   56^3 =   69^3 +   42^3
2000:1671816384 = 1168^3 +  428^3 =  944^3 +  940^3
2001:1672470592 = 1124^3 +  632^3 = 1187^3 +   29^3
2002:1673170856 = 1034^3 +  828^3 = 1164^3 +  458^3
2003:1675045225 = 1153^3 +  522^3 = 1081^3 +  744^3
2004:1675958167 = 1096^3 +  711^3 = 1159^3 +  492^3
2005:1676926719 = 1188^3 +   63^3 = 1095^3 +  714^3
2006:1677646971 =  990^3 +  891^3 = 1188^3 +   99^3

D

High Level Version

Translation of: Python

<lang d>void main() {

   import std.stdio, std.range, std.algorithm,std.typecons,std.string;
   const cubes = iota(1, 1201).map!(x => tuple(x, ulong(x)^^3)).array;
   bool[Tuple!(const int, const int)][ulong] sum2cubes;
   foreach (const p1; cubes)
       foreach (const p2; cubes[p1[0] .. $])
           sum2cubes[p1[1] + p2[1]][tuple(p2[0], p1[0])] = true;
   const taxied = zip(sum2cubes.byKey, sum2cubes.byValue)
                  .filter!(kv => kv[1].length > 1).array
                  .schwartzSort!(kv => kv[0]).release;
   foreach (const r; [[0, 25], [2000-1, 2000+6]]) {
       foreach (immutable i, const t; taxied[r[0] .. r[1]])
           writefln("%4d: %10d = %-(%s = %)", i + r[0] + 1, t[0],
                    t[1].byKey.map!q{"%4d^3 + %4d^3".format(a[])});
       writeln;
   }

}</lang>

Output:
   1:       1729 =   12^3 +    1^3 =   10^3 +    9^3
   2:       4104 =   16^3 +    2^3 =   15^3 +    9^3
   3:      13832 =   24^3 +    2^3 =   20^3 +   18^3
   4:      20683 =   27^3 +   10^3 =   24^3 +   19^3
   5:      32832 =   30^3 +   18^3 =   32^3 +    4^3
   6:      39312 =   34^3 +    2^3 =   33^3 +   15^3
   7:      40033 =   34^3 +    9^3 =   33^3 +   16^3
   8:      46683 =   30^3 +   27^3 =   36^3 +    3^3
   9:      64232 =   36^3 +   26^3 =   39^3 +   17^3
  10:      65728 =   40^3 +   12^3 =   33^3 +   31^3
  11:     110656 =   48^3 +    4^3 =   40^3 +   36^3
  12:     110808 =   48^3 +    6^3 =   45^3 +   27^3
  13:     134379 =   43^3 +   38^3 =   51^3 +   12^3
  14:     149389 =   53^3 +    8^3 =   50^3 +   29^3
  15:     165464 =   54^3 +   20^3 =   48^3 +   38^3
  16:     171288 =   54^3 +   24^3 =   55^3 +   17^3
  17:     195841 =   58^3 +    9^3 =   57^3 +   22^3
  18:     216027 =   59^3 +   22^3 =   60^3 +    3^3
  19:     216125 =   50^3 +   45^3 =   60^3 +    5^3
  20:     262656 =   60^3 +   36^3 =   64^3 +    8^3
  21:     314496 =   68^3 +    4^3 =   66^3 +   30^3
  22:     320264 =   68^3 +   18^3 =   66^3 +   32^3
  23:     327763 =   67^3 +   30^3 =   58^3 +   51^3
  24:     373464 =   72^3 +    6^3 =   60^3 +   54^3
  25:     402597 =   69^3 +   42^3 =   61^3 +   56^3

2000: 1671816384 =  944^3 +  940^3 = 1168^3 +  428^3
2001: 1672470592 = 1187^3 +   29^3 = 1124^3 +  632^3
2002: 1673170856 = 1164^3 +  458^3 = 1034^3 +  828^3
2003: 1675045225 = 1153^3 +  522^3 = 1081^3 +  744^3
2004: 1675958167 = 1159^3 +  492^3 = 1096^3 +  711^3
2005: 1676926719 = 1188^3 +   63^3 = 1095^3 +  714^3
2006: 1677646971 = 1188^3 +   99^3 =  990^3 +  891^3

Heap-Based Version

Translation of: Java

<lang d>import std.stdio, std.string, std.container;

struct CubeSum {

   ulong x, y, value;
   this(in ulong x_, in ulong y_) pure nothrow {
       this.x = x_;
       this.y = y_;
       this.value = x_ ^^ 3 + y_ ^^ 3;
   }

}

final class Taxi {

   BinaryHeap!(Array!CubeSum, "a.value > b.value") pq;
   CubeSum last;
   ulong n = 0;
   this() {
       last = nextSum;
   }
   CubeSum nextSum() {
       while (pq.empty || pq.front.value >= n ^^ 3)
           pq.insert(CubeSum(++n, 1));
       auto s = pq.front;
       pq.removeFront;
       if (s.x > s.y + 1)
           pq.insert(CubeSum(s.x, s.y + 1));
       return s;
   }
   CubeSum[] nextTaxi() {
       CubeSum s;
       typeof(return) train;
       while ((s = nextSum).value != last.value)
           last = s;
       train ~= last;
       do {
           train ~= s;
       } while ((s = nextSum).value == last.value);
       last = s;
       return train;
   }

}

void main() {

   auto taxi = new Taxi;
   foreach (immutable i; 1 .. 2007) {
       const t = taxi.nextTaxi;
       if (i > 25 && i < 2000)
           continue;
       writef("%4d: %10d", i, t[0].value);
       foreach (const s; t)
           writef(" = %4d^3 + %4d^3", s.x, s.y);
       writeln;
   }

}</lang>

Output:
   1:       1729 =   10^3 +    9^3 =   12^3 +    1^3
   2:       4104 =   15^3 +    9^3 =   16^3 +    2^3
   3:      13832 =   20^3 +   18^3 =   24^3 +    2^3
   4:      20683 =   24^3 +   19^3 =   27^3 +   10^3
   5:      32832 =   30^3 +   18^3 =   32^3 +    4^3
   6:      39312 =   33^3 +   15^3 =   34^3 +    2^3
   7:      40033 =   33^3 +   16^3 =   34^3 +    9^3
   8:      46683 =   30^3 +   27^3 =   36^3 +    3^3
   9:      64232 =   39^3 +   17^3 =   36^3 +   26^3
  10:      65728 =   40^3 +   12^3 =   33^3 +   31^3
  11:     110656 =   40^3 +   36^3 =   48^3 +    4^3
  12:     110808 =   45^3 +   27^3 =   48^3 +    6^3
  13:     134379 =   51^3 +   12^3 =   43^3 +   38^3
  14:     149389 =   50^3 +   29^3 =   53^3 +    8^3
  15:     165464 =   48^3 +   38^3 =   54^3 +   20^3
  16:     171288 =   54^3 +   24^3 =   55^3 +   17^3
  17:     195841 =   57^3 +   22^3 =   58^3 +    9^3
  18:     216027 =   59^3 +   22^3 =   60^3 +    3^3
  19:     216125 =   50^3 +   45^3 =   60^3 +    5^3
  20:     262656 =   60^3 +   36^3 =   64^3 +    8^3
  21:     314496 =   66^3 +   30^3 =   68^3 +    4^3
  22:     320264 =   68^3 +   18^3 =   66^3 +   32^3
  23:     327763 =   67^3 +   30^3 =   58^3 +   51^3
  24:     373464 =   60^3 +   54^3 =   72^3 +    6^3
  25:     402597 =   69^3 +   42^3 =   61^3 +   56^3
2000: 1671816384 = 1168^3 +  428^3 =  944^3 +  940^3
2001: 1672470592 = 1124^3 +  632^3 = 1187^3 +   29^3
2002: 1673170856 = 1164^3 +  458^3 = 1034^3 +  828^3
2003: 1675045225 = 1153^3 +  522^3 = 1081^3 +  744^3
2004: 1675958167 = 1159^3 +  492^3 = 1096^3 +  711^3
2005: 1676926719 = 1095^3 +  714^3 = 1188^3 +   63^3
2006: 1677646971 =  990^3 +  891^3 = 1188^3 +   99^3

Run-time: about 0.67 seconds with ldc2 compiler.

Low Level Heap-Based Version

Translation of: C

<lang d>struct Taxicabs {

   alias CubesSumT = uint; // Or ulong.
   static struct Sum {
       CubesSumT value;
       uint x, y;
   }
   // The cubes can be pre-computed if CubesSumT is a BigInt.
   private uint nCubes;
   private Sum[] pq;
   private uint pq_len;
   private void addCube() pure nothrow {
       nCubes = nCubes ? nCubes + 1 : 2;
       if (nCubes < 2)
           return; // 0 or 1 is useless.
       pq_len++;
       if (pq_len >= pq.length)
           pq.length = (pq.length == 0) ? 2 : (pq.length * 2);
       immutable tmp = Sum(CubesSumT(nCubes - 2) ^^ 3 + 1,
                           nCubes - 2, 1);
       // Upheap.
       uint i = pq_len;
       for (; i >= 1 && pq[i >> 1].value > tmp.value; i >>= 1)
           pq[i] = pq[i >> 1];
       pq[i] = tmp;
   }


   private void nextSum() pure nothrow {
       while (!pq_len || pq[1].value >= (nCubes - 1) ^^ 3)
           addCube();
       Sum tmp = pq[0] = pq[1]; //pq[0] always stores last seen value.
       tmp.y++;
       if (tmp.y >= tmp.x) { // Done with this x; throw it away.
           tmp = pq[pq_len];
           pq_len--;
           if (!pq_len)
               return nextSum(); // Refill empty heap.
       } else
           tmp.value += tmp.y ^^ 3 - (tmp.y - 1) ^^ 3;
       // Downheap.
       uint i = 1;
       while (true) {
           uint j = i << 1;
           if (j > pq_len)
               break;
           if (j < pq_len && pq[j + 1].value < pq[j].value)
               j++;
           if (pq[j].value >= tmp.value)
               break;
           pq[i] = pq[j];
           i = j;
       }
       pq[i] = tmp;
   }


   Sum[] nextTaxi(size_t N)(ref Sum[N] hist) pure nothrow {
       do {
           nextSum();
       } while (pq[0].value != pq[1].value);
       uint len = 1;
       hist[0] = pq[0];
       do {
           hist[len] = pq[1];
           len++;
           nextSum();
       } while (pq[0].value == pq[1].value);
       return hist[0 .. len];
   }

}


void main() nothrow {

   import core.stdc.stdio;
   Taxicabs t;
   Taxicabs.Sum[3] x;
   foreach (immutable uint i; 1 .. 2007) {
       const triples = t.nextTaxi(x);
       if (i > 25 && i < 2000)
           continue;
       printf("%4u: %10lu", i, triples[0].value);
       foreach_reverse (const s; triples)
           printf(" = %4u^3 + %4u^3", s.x, s.y);
       putchar('\n');
   }

}</lang>

Output:
   1:       1729 =   12^3 +    1^3 =   10^3 +    9^3
   2:       4104 =   15^3 +    9^3 =   16^3 +    2^3
   3:      13832 =   20^3 +   18^3 =   24^3 +    2^3
   4:      20683 =   27^3 +   10^3 =   24^3 +   19^3
   5:      32832 =   30^3 +   18^3 =   32^3 +    4^3
   6:      39312 =   33^3 +   15^3 =   34^3 +    2^3
   7:      40033 =   33^3 +   16^3 =   34^3 +    9^3
   8:      46683 =   30^3 +   27^3 =   36^3 +    3^3
   9:      64232 =   36^3 +   26^3 =   39^3 +   17^3
  10:      65728 =   33^3 +   31^3 =   40^3 +   12^3
  11:     110656 =   40^3 +   36^3 =   48^3 +    4^3
  12:     110808 =   45^3 +   27^3 =   48^3 +    6^3
  13:     134379 =   43^3 +   38^3 =   51^3 +   12^3
  14:     149389 =   50^3 +   29^3 =   53^3 +    8^3
  15:     165464 =   48^3 +   38^3 =   54^3 +   20^3
  16:     171288 =   54^3 +   24^3 =   55^3 +   17^3
  17:     195841 =   57^3 +   22^3 =   58^3 +    9^3
  18:     216027 =   59^3 +   22^3 =   60^3 +    3^3
  19:     216125 =   50^3 +   45^3 =   60^3 +    5^3
  20:     262656 =   60^3 +   36^3 =   64^3 +    8^3
  21:     314496 =   66^3 +   30^3 =   68^3 +    4^3
  22:     320264 =   66^3 +   32^3 =   68^3 +   18^3
  23:     327763 =   58^3 +   51^3 =   67^3 +   30^3
  24:     373464 =   60^3 +   54^3 =   72^3 +    6^3
  25:     402597 =   61^3 +   56^3 =   69^3 +   42^3
2000: 1671816384 = 1168^3 +  428^3 =  944^3 +  940^3
2001: 1672470592 = 1124^3 +  632^3 = 1187^3 +   29^3
2002: 1673170856 = 1034^3 +  828^3 = 1164^3 +  458^3
2003: 1675045225 = 1153^3 +  522^3 = 1081^3 +  744^3
2004: 1675958167 = 1096^3 +  711^3 = 1159^3 +  492^3
2005: 1676926719 = 1188^3 +   63^3 = 1095^3 +  714^3
2006: 1677646971 =  990^3 +  891^3 = 1188^3 +   99^3

Run-time: about 0.08 seconds with ldc2 compiler.

J

<lang J> 25 {. ;({."#. <@(0&#`({.@{.(;,)<@}."1)@.(1<#))/. ])/:~~.,/(+,/:~@,)"0/~3^~1+i.100 ┌──────┬────────────┬─────────────┐ │1729 │1 1728 │729 1000 │ ├──────┼────────────┼─────────────┤ │4104 │8 4096 │729 3375 │ ├──────┼────────────┼─────────────┤ │13832 │8 13824 │5832 8000 │ ├──────┼────────────┼─────────────┤ │20683 │1000 19683 │6859 13824 │ ├──────┼────────────┼─────────────┤ │32832 │64 32768 │5832 27000 │ ├──────┼────────────┼─────────────┤ │39312 │8 39304 │3375 35937 │ ├──────┼────────────┼─────────────┤ │40033 │729 39304 │4096 35937 │ ├──────┼────────────┼─────────────┤ │46683 │27 46656 │19683 27000 │ ├──────┼────────────┼─────────────┤ │64232 │4913 59319 │17576 46656 │ ├──────┼────────────┼─────────────┤ │65728 │1728 64000 │29791 35937 │ ├──────┼────────────┼─────────────┤ │110656│64 110592 │46656 64000 │ ├──────┼────────────┼─────────────┤ │110808│216 110592 │19683 91125 │ ├──────┼────────────┼─────────────┤ │134379│1728 132651 │54872 79507 │ ├──────┼────────────┼─────────────┤ │149389│512 148877 │24389 125000 │ ├──────┼────────────┼─────────────┤ │165464│8000 157464 │54872 110592 │ ├──────┼────────────┼─────────────┤ │171288│4913 166375 │13824 157464 │ ├──────┼────────────┼─────────────┤ │195841│729 195112 │10648 185193 │ ├──────┼────────────┼─────────────┤ │216027│27 216000 │10648 205379 │ ├──────┼────────────┼─────────────┤ │216125│125 216000 │91125 125000 │ ├──────┼────────────┼─────────────┤ │262656│512 262144 │46656 216000 │ ├──────┼────────────┼─────────────┤ │314496│64 314432 │27000 287496 │ ├──────┼────────────┼─────────────┤ │320264│5832 314432 │32768 287496 │ ├──────┼────────────┼─────────────┤ │327763│27000 300763│132651 195112│ ├──────┼────────────┼─────────────┤ │373464│216 373248 │157464 216000│ ├──────┼────────────┼─────────────┤ │402597│74088 328509│175616 226981│ └──────┴────────────┴─────────────┘</lang>

Explanation:

First, generate 100 cubes.

Then, form a 3 column table of unique rows: sum, small cube, large cube

Then, gather rows where the first entry is the same. Keep the ones with at least two such entries.

Note that the cube root of the 25th entry is slightly smaller than 74, so testing against the first 100 cubes is more than sufficient.

Java

<lang java>import java.util.PriorityQueue; import java.util.ArrayList;

class CubeSum implements Comparable<CubeSum> { public long x, y, value;

public CubeSum(long x, long y) { this.x = x; this.y = y; this.value = x*x*x + y*y*y; }

public String toString() { return String.format("%4d^3 + %4d^3", x, y); }

public int compareTo(CubeSum that) { return value < that.value ? -1 : value > that.value ? 1 : 0; } }

public class Taxi { PriorityQueue<CubeSum> pq; CubeSum last; long n;

public Taxi() { pq = new PriorityQueue<CubeSum>(); n = 0; last = NextSum(); }

CubeSum NextSum() { while (pq.size() == 0 || pq.peek().value >= n*n*n) pq.add(new CubeSum(++n, 1));

CubeSum s = pq.remove(); if (s.x > s.y + 1) pq.add(new CubeSum(s.x, s.y+1));

return s; }

ArrayList<CubeSum> NextTaxi() { CubeSum s; ArrayList<CubeSum> train = new ArrayList<CubeSum>();

while ((s = NextSum()).value != last.value) last = s;

train.add(last);

do { train.add(s); } while ((s = NextSum()).value == last.value); last = s;

return train; }

public static final void main(String[] args) { Taxi taxi = new Taxi();

for (int i = 1; i <= 2006; i++) { ArrayList<CubeSum> t = taxi.NextTaxi(); if (i > 25 && i < 2000) continue;

System.out.printf("%4d: %10d", i, t.get(0).value); for (CubeSum s: t) System.out.print(" = " + s); System.out.println(); } } }</lang>

Output:
   1:       1729 =   10^3 +    9^3 =   12^3 +    1^3
   2:       4104 =   15^3 +    9^3 =   16^3 +    2^3
   3:      13832 =   20^3 +   18^3 =   24^3 +    2^3
   4:      20683 =   24^3 +   19^3 =   27^3 +   10^3
   5:      32832 =   30^3 +   18^3 =   32^3 +    4^3
   6:      39312 =   33^3 +   15^3 =   34^3 +    2^3
   7:      40033 =   34^3 +    9^3 =   33^3 +   16^3
   8:      46683 =   30^3 +   27^3 =   36^3 +    3^3
   9:      64232 =   36^3 +   26^3 =   39^3 +   17^3
  10:      65728 =   33^3 +   31^3 =   40^3 +   12^3
  11:     110656 =   40^3 +   36^3 =   48^3 +    4^3
  12:     110808 =   45^3 +   27^3 =   48^3 +    6^3
  13:     134379 =   43^3 +   38^3 =   51^3 +   12^3
  14:     149389 =   50^3 +   29^3 =   53^3 +    8^3
  15:     165464 =   48^3 +   38^3 =   54^3 +   20^3
  16:     171288 =   54^3 +   24^3 =   55^3 +   17^3
  17:     195841 =   57^3 +   22^3 =   58^3 +    9^3
  18:     216027 =   59^3 +   22^3 =   60^3 +    3^3
  19:     216125 =   50^3 +   45^3 =   60^3 +    5^3
  20:     262656 =   60^3 +   36^3 =   64^3 +    8^3
  21:     314496 =   66^3 +   30^3 =   68^3 +    4^3
  22:     320264 =   66^3 +   32^3 =   68^3 +   18^3
  23:     327763 =   58^3 +   51^3 =   67^3 +   30^3
  24:     373464 =   60^3 +   54^3 =   72^3 +    6^3
  25:     402597 =   61^3 +   56^3 =   69^3 +   42^3
2000: 1671816384 = 1168^3 +  428^3 =  944^3 +  940^3
2001: 1672470592 = 1124^3 +  632^3 = 1187^3 +   29^3
2002: 1673170856 = 1164^3 +  458^3 = 1034^3 +  828^3
2003: 1675045225 = 1153^3 +  522^3 = 1081^3 +  744^3
2004: 1675958167 = 1159^3 +  492^3 = 1096^3 +  711^3
2005: 1676926719 = 1095^3 +  714^3 = 1188^3 +   63^3
2006: 1677646971 =  990^3 +  891^3 = 1188^3 +   99^3

PARI/GP

<lang parigp>taxicab(n)=my(t); for(k=sqrtnint((n-1)\2,3)+1, sqrtnint(n,3), if(ispower(n-k^3, 3), if(t, return(1), t=1))); 0; cubes(n)=my(t); for(k=sqrtnint((n-1)\2,3)+1, sqrtnint(n,3), if(ispower(n-k^3, 3, &t), print(n" = \t"k"^3\t+ "t"^3"))) select(taxicab, [1..402597]) apply(cubes, %);</lang>

Output:
%1 = [1729, 4104, 13832, 20683, 32832, 39312, 40033, 46683, 64232, 65728, 110656, 110808, 134379, 149389, 165464, 171288, 195841, 216027, 216125, 262656, 314496, 320264, 327763, 373464, 402597]
1729 =          10^3    + 9^3
1729 =          12^3    + 1^3
4104 =          15^3    + 9^3
4104 =          16^3    + 2^3
13832 =         20^3    + 18^3
13832 =         24^3    + 2^3
20683 =         24^3    + 19^3
20683 =         27^3    + 10^3
32832 =         30^3    + 18^3
32832 =         32^3    + 4^3
39312 =         33^3    + 15^3
39312 =         34^3    + 2^3
40033 =         33^3    + 16^3
40033 =         34^3    + 9^3
46683 =         30^3    + 27^3
46683 =         36^3    + 3^3
64232 =         36^3    + 26^3
64232 =         39^3    + 17^3
65728 =         33^3    + 31^3
65728 =         40^3    + 12^3
110656 =        40^3    + 36^3
110656 =        48^3    + 4^3
110808 =        45^3    + 27^3
110808 =        48^3    + 6^3
134379 =        43^3    + 38^3
134379 =        51^3    + 12^3
149389 =        50^3    + 29^3
149389 =        53^3    + 8^3
165464 =        48^3    + 38^3
165464 =        54^3    + 20^3
171288 =        54^3    + 24^3
171288 =        55^3    + 17^3
195841 =        57^3    + 22^3
195841 =        58^3    + 9^3
216027 =        59^3    + 22^3
216027 =        60^3    + 3^3
216125 =        50^3    + 45^3
216125 =        60^3    + 5^3
262656 =        60^3    + 36^3
262656 =        64^3    + 8^3
314496 =        66^3    + 30^3
314496 =        68^3    + 4^3
320264 =        66^3    + 32^3
320264 =        68^3    + 18^3
327763 =        58^3    + 51^3
327763 =        67^3    + 30^3
373464 =        60^3    + 54^3
373464 =        72^3    + 6^3
402597 =        61^3    + 56^3
402597 =        69^3    + 42^3

Perl 6

This uses a pretty simple search algorithm that doesn't necessarily return the Taxicab numbers in order. Assuming we want all the Taxicab numbers within some range S to N, we'll search until we find N values. When we find the Nth value, we continue to search up to the cube root of the largest Taxicab number found up to that point. That ensures we will find all of them inside the desired range without needing to search arbitrarily or use magic numbers. Defaults to returning the Taxicab numbers from 1 to 25. Pass in a different start and end value if you want some other range. <lang perl6>sub MAIN ($start = 1, $end = 25) {

   my %taxi;
   my $taxis = 0;
   my $terminate = 0;
   for 1 .. * -> $c1 {
       display( %taxi, $start, $end ) and exit if 0 < $terminate < $c1;
       my $c = $c1 ** 3;
       for 1 .. $c1 -> $c2 {
           my $this = $c2 ** 3 + $c;
           %taxi{$this}.push: [$c2, $c1];
           $taxis++ if %taxi{$this}.elems == 2;
   	    $terminate = %taxi.grep( { $_.value.elems > 1 } ).sort( +*.key )[*-1].key**(1/3)
               if $taxis == $end and !$terminate;
       }
   }

}

sub display (%this_stuff, $start, $end) {

   my $i = $start; 
   printf "%4d %10d  =>\t%s\n", $i++, $_.key, 
       ($_.value.map({ sprintf "%4d³ + %-s", $_[0], "$_[1]³" })).join: ",\t"
       for %this_stuff.grep( { $_.value.elems > 1 } ).sort( +*.key )[$start-1..$end-1];
   1;

}</lang>

Output:

With no passed parameters (default)

   1       1729  =>	   9³ + 10³,	   1³ + 12³
   2       4104  =>	   9³ + 15³,	   2³ + 16³
   3      13832  =>	  18³ + 20³,	   2³ + 24³
   4      20683  =>	  19³ + 24³,	  10³ + 27³
   5      32832  =>	  18³ + 30³,	   4³ + 32³
   6      39312  =>	  15³ + 33³,	   2³ + 34³
   7      40033  =>	  16³ + 33³,	   9³ + 34³
   8      46683  =>	  27³ + 30³,	   3³ + 36³
   9      64232  =>	  26³ + 36³,	  17³ + 39³
  10      65728  =>	  31³ + 33³,	  12³ + 40³
  11     110656  =>	  36³ + 40³,	   4³ + 48³
  12     110808  =>	  27³ + 45³,	   6³ + 48³
  13     134379  =>	  38³ + 43³,	  12³ + 51³
  14     149389  =>	  29³ + 50³,	   8³ + 53³
  15     165464  =>	  38³ + 48³,	  20³ + 54³
  16     171288  =>	  24³ + 54³,	  17³ + 55³
  17     195841  =>	  22³ + 57³,	   9³ + 58³
  18     216027  =>	  22³ + 59³,	   3³ + 60³
  19     216125  =>	  45³ + 50³,	   5³ + 60³
  20     262656  =>	  36³ + 60³,	   8³ + 64³
  21     314496  =>	  30³ + 66³,	   4³ + 68³
  22     320264  =>	  32³ + 66³,	  18³ + 68³
  23     327763  =>	  51³ + 58³,	  30³ + 67³
  24     373464  =>	  54³ + 60³,	   6³ + 72³
  25     402597  =>	  56³ + 61³,	  42³ + 69³

With passed parameters 2000 2006:

2000 1671816384  =>	 940³ + 944³,	 428³ + 1168³
2001 1672470592  =>	 632³ + 1124³,	  29³ + 1187³
2002 1673170856  =>	 828³ + 1034³,	 458³ + 1164³
2003 1675045225  =>	 744³ + 1081³,	 522³ + 1153³
2004 1675958167  =>	 711³ + 1096³,	 492³ + 1159³
2005 1676926719  =>	 714³ + 1095³,	  63³ + 1188³
2006 1677646971  =>	 891³ + 990³,	  99³ + 1188³

Python

(Magic number 1201 found by trial and error) <lang python>from collections import defaultdict from itertools import product from pprint import pprint as pp

cube2n = {x**3:x for x in range(1, 1201)} sum2cubes = defaultdict(set) for c1, c2 in product(cube2n, cube2n): if c1 >= c2: sum2cubes[c1 + c2].add((cube2n[c1], cube2n[c2]))

taxied = sorted((k, v) for k,v in sum2cubes.items() if len(v) >= 2)

  1. pp(len(taxied)) # 2068

for t in enumerate(taxied[:25], 1):

   pp(t)

print('...') for t in enumerate(taxied[2000-1:2000+6], 2000):

   pp(t)</lang>
Output:
(1, (1729, {(12, 1), (10, 9)}))
(2, (4104, {(16, 2), (15, 9)}))
(3, (13832, {(20, 18), (24, 2)}))
(4, (20683, {(27, 10), (24, 19)}))
(5, (32832, {(30, 18), (32, 4)}))
(6, (39312, {(33, 15), (34, 2)}))
(7, (40033, {(33, 16), (34, 9)}))
(8, (46683, {(30, 27), (36, 3)}))
(9, (64232, {(36, 26), (39, 17)}))
(10, (65728, {(33, 31), (40, 12)}))
(11, (110656, {(48, 4), (40, 36)}))
(12, (110808, {(48, 6), (45, 27)}))
(13, (134379, {(51, 12), (43, 38)}))
(14, (149389, {(50, 29), (53, 8)}))
(15, (165464, {(54, 20), (48, 38)}))
(16, (171288, {(54, 24), (55, 17)}))
(17, (195841, {(57, 22), (58, 9)}))
(18, (216027, {(60, 3), (59, 22)}))
(19, (216125, {(60, 5), (50, 45)}))
(20, (262656, {(64, 8), (60, 36)}))
(21, (314496, {(66, 30), (68, 4)}))
(22, (320264, {(66, 32), (68, 18)}))
(23, (327763, {(58, 51), (67, 30)}))
(24, (373464, {(72, 6), (60, 54)}))
(25, (402597, {(69, 42), (61, 56)}))
...
(2000, (1671816384, {(1168, 428), (944, 940)}))
(2001, (1672470592, {(1187, 29), (1124, 632)}))
(2002, (1673170856, {(1164, 458), (1034, 828)}))
(2003, (1675045225, {(1153, 522), (1081, 744)}))
(2004, (1675958167, {(1159, 492), (1096, 711)}))
(2005, (1676926719, {(1188, 63), (1095, 714)}))
(2006, (1677646971, {(990, 891), (1188, 99)}))

Although, for this task it's simply faster to look up the cubes in the sum when we need to print them, because we can now store and sort only the sums: <lang python>cubes, crev = [x**3 for x in range(1,1200)], {}

  1. for cube root lookup

for x,x3 in enumerate(cubes): crev[x3] = x + 1

sums = sorted([x+y for x in cubes for y in cubes if y < x])

idx = 0 for i in range(1, len(sums)-1):

   if sums[i-1] != sums[i] and sums[i] == sums[i+1]:
       idx += 1
       if idx > 25 and idx < 2000 or idx > 2006: continue
       n,p = sums[i],[]
       for x in cubes:
           if n-x < x: break
           if n-x in crev:
               p.append((crev[x], crev[n-x]))
       print "%4d: %10d"%(idx,n),
       for x in p: print " = %4d^3 + %4d^3"%x,
       print</lang>
Output:

Output trimmed to reduce clutter.

   1:       1729  =    1^3 +   12^3  =    9^3 +   10^3
   2:       4104  =    2^3 +   16^3  =    9^3 +   15^3
   3:      13832  =    2^3 +   24^3  =   18^3 +   20^3
   4:      20683  =   10^3 +   27^3  =   19^3 +   24^3
   5:      32832  =    4^3 +   32^3  =   18^3 +   30^3
...
2004: 1675958167  =  492^3 + 1159^3  =  711^3 + 1096^3
2005: 1676926719  =   63^3 + 1188^3  =  714^3 + 1095^3
2006: 1677646971  =   99^3 + 1188^3  =  891^3 +  990^3

Using heapq module

A priority queue that holds cube sums. When consecutive sums come out with the same value, they are taxis. <lang python>from heapq import heappush, heappop

def cubesum():

   h,n = [],1
   while True:
       while not h or h[0][0] > n**3: # could also pre-calculate cubes
           heappush(h, (n**3 + 1, n, 1))
           n += 1
       (s, x, y) = heappop(h)
       yield((s, x, y))
       y += 1
       if y < x:    # should be y <= x?
           heappush(h, (x**3 + y**3, x, y))

def taxis():

   out = [(0,0,0)]
   for s in cubesum():
       if s[0] == out[-1][0]:
           out.append(s)
       else:
           if len(out) > 1: yield(out)
           out = [s]

n = 0 for x in taxis():

   n += 1
   if n >= 2006: break
   if n <= 25 or n >= 2000:
       print(n, x)</lang>
Output:
(1, [(1729, 10, 9), (1729, 12, 1)])
(2, [(4104, 15, 9), (4104, 16, 2)])
(3, [(13832, 20, 18), (13832, 24, 2)])
(4, [(20683, 24, 19), (20683, 27, 10)])
(5, [(32832, 30, 18), (32832, 32, 4)])
(6, [(39312, 33, 15), (39312, 34, 2)])
(7, [(40033, 33, 16), (40033, 34, 9)])
(8, [(46683, 30, 27), (46683, 36, 3)])
(9, [(64232, 36, 26), (64232, 39, 17)])
(10, [(65728, 33, 31), (65728, 40, 12)])
(11, [(110656, 40, 36), (110656, 48, 4)])
(12, [(110808, 45, 27), (110808, 48, 6)])
(13, [(134379, 43, 38), (134379, 51, 12)])
(14, [(149389, 50, 29), (149389, 53, 8)])
(15, [(165464, 48, 38), (165464, 54, 20)])
(16, [(171288, 54, 24), (171288, 55, 17)])
(17, [(195841, 57, 22), (195841, 58, 9)])
(18, [(216027, 59, 22), (216027, 60, 3)])
(19, [(216125, 50, 45), (216125, 60, 5)])
(20, [(262656, 60, 36), (262656, 64, 8)])
(21, [(314496, 66, 30), (314496, 68, 4)])
(22, [(320264, 66, 32), (320264, 68, 18)])
(23, [(327763, 58, 51), (327763, 67, 30)])
(24, [(373464, 60, 54), (373464, 72, 6)])
(25, [(402597, 61, 56), (402597, 69, 42)])
(2000, [(1671816384, 944, 940), (1671816384, 1168, 428)])
(2001, [(1672470592, 1124, 632), (1672470592, 1187, 29)])
(2002, [(1673170856, 1034, 828), (1673170856, 1164, 458)])
(2003, [(1675045225, 1081, 744), (1675045225, 1153, 522)])
(2004, [(1675958167, 1096, 711), (1675958167, 1159, 492)])
(2005, [(1676926719, 1095, 714), (1676926719, 1188, 63)])

REXX

<lang rexx>/*REXX program displays the first (lowest) taxicab numbers. */ parse arg L1 H1 L2 H2 . /*obtain the optional numbers. */ if L1== | L1==',' then L1=1 /*L1 " " " " " */ if H1== | H1==',' then H1=25 /*H1 " " " " " */ if L2== | L2==',' then L2=2000 /*L2 " " " " " */ if H2== | H2==',' then H2=2007 /*H2 " " " " " */ mx=max(L1, H1, L2, H2) /*find how many taxicab #s needed*/ mx=mx+mx%10 /*cushion, compensate for triples*/ w=length(mx) /*width is used formatting output*/ numeric digits 20 /*prepare to use larger numbers. */

  1. =0; @.=0; $.=; b='■'; p='**3' /*initialize some REXX variables.*/
                                      /* [↓]  generate extra taxicab #s*/
do j=1   until  #>=mx                 /*taxicab nums won't be in order.*/
jjj=j**3                              /*might as well calculate a cube.*/
         do k=1  for j-1; s=jjj+k**3  /*define a whole bunch of cubes. */
         if @.s==0  then do           /*if cube not defined, then do it*/
                         @.s = "────►"right(j,6,b)p"■■■+"right(k,6,b)p
                         iterate      /* ··· and then go and do another*/
                         end          /* [↑] define one cube at a time.*/
         comma=pos(',',@.s)\==0       /*has it has been defined before?*/
         @.s=@.s","right(j,9,b)p"■■■+"right(k,6,b)p   /*build the text.*/
         $.s=right(s,15,b)'■■■'@.s    /*define the rest of taxicab #s. */
         if comma  then iterate       /*S  is a triple (or better).    */
         #=#+1                        /*bump the taxicab number count. */
         #.#=s                        /*define a   #.   taxicab number.*/
         end   /*k*/                  /* [↑]  build cubes one-at-a-time*/
end   /*j*/                           /* [↑]  complete with overage #s.*/

h=mx /*H= ½─way point for pivot sort. */

     do  while  h>1;  h=h%2;   do i=1  for mx-h; j=i; k=h+i    /*sort. */
     do  while  #.k<#.j;       _=#.j;  #.j=#.k;  #.k=_         /*sort. */
         if h>=j  then leave;  j=j-h;  k=k-h;  end;  end       /*sort. */
     end   /*while h>1*/              /* [↑]  sort the taxicab # array.*/

call show L1,H1 /*show 1st range of taxicab #s. */ call show L2,H2 /*show 2nd range of taxicab #s. */ exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────SHOW subroutine─────────────────────*/ show: parse arg low,high; if low==0 then return /*show bunch taxicab#s,*/

    do t=low to high;  _=#.t          /*get single taxicab # at a time.*/
    say right(t,w)  translate($._,,b) /*display taxicab # (with blanks)*/
    end   /*t*/                       /* [↑] ■ are translated to blanks*/

return</lang> output   using the default input:

 
   1            1729            10**3   +     9**3,       12**3   +     1**3
   2            4104            15**3   +     9**3,       16**3   +     2**3
   3           13832            20**3   +    18**3,       24**3   +     2**3
   4           20683            24**3   +    19**3,       27**3   +    10**3
   5           32832            30**3   +    18**3,       32**3   +     4**3
   6           39312            33**3   +    15**3,       34**3   +     2**3
   7           40033            33**3   +    16**3,       34**3   +     9**3
   8           46683            30**3   +    27**3,       36**3   +     3**3
   9           64232            36**3   +    26**3,       39**3   +    17**3
  10           65728            33**3   +    31**3,       40**3   +    12**3
  11          110656            40**3   +    36**3,       48**3   +     4**3
  12          110808            45**3   +    27**3,       48**3   +     6**3
  13          134379            43**3   +    38**3,       51**3   +    12**3
  14          149389            50**3   +    29**3,       53**3   +     8**3
  15          165464            48**3   +    38**3,       54**3   +    20**3
  16          171288            54**3   +    24**3,       55**3   +    17**3
  17          195841            57**3   +    22**3,       58**3   +     9**3
  18          216027            59**3   +    22**3,       60**3   +     3**3
  19          216125            50**3   +    45**3,       60**3   +     5**3
  20          262656            60**3   +    36**3,       64**3   +     8**3
  21          314496            66**3   +    30**3,       68**3   +     4**3
  22          320264            66**3   +    32**3,       68**3   +    18**3
  23          327763            58**3   +    51**3,       67**3   +    30**3
  24          373464            60**3   +    54**3,       72**3   +     6**3
  25          402597            61**3   +    56**3,       69**3   +    42**3

2000      1671816384           944**3   +   940**3,     1168**3   +   428**3
2001      1672470592          1124**3   +   632**3,     1187**3   +    29**3
2002      1673170856          1034**3   +   828**3,     1164**3   +   458**3
2003      1675045225          1081**3   +   744**3,     1153**3   +   522**3
2004      1675958167          1096**3   +   711**3,     1159**3   +   492**3
2005      1676926719          1095**3   +   714**3,     1188**3   +    63**3
2006      1677646971           990**3   +   891**3,     1188**3   +    99**3
2007      1680918365          1132**3   +   613**3,     1189**3   +    16**3