How do I see if $a^3+b^3=c^3+d^3$ has any solutions where $ 1 \le a,b,c,d \in \mathbb{Z} \le 1000$ and $a \ne b \ne c \ne d$? How do I see if  $a^3+b^3=c^3+d^3$ has any solutions where $ 1 \le a,b,c,d \in \mathbb{Z} \le 1000$ and $a \ne b \ne c \ne d$ ? 
I know I can write a program to brute force this and find out, but is there a way I can determine this through algebra?
I thought I could use the Difference of Powers formula: 
$(a-c)(a^2+ac+c^2)=a^3-c^3=d^3-b^3=(d-b)(d^2+db+b^2)$. 
At this point I am stuck. 
 A: You are asking about Taxicab Numbers, of which the most famous is
$$1729 = 12^3 + 1^3 = 10^3 + 9^3.$$
There is no (known) algebraic method to solve this equation. All of our current understand of these types of numbers come from computer calculation. There is also an OEIS entry dedicated to them.
If you're curious, we also have found numbers which can be written with $3$+ pairs:
\begin{align*}
    87539319 &=167^3 + 436^3 \\
    &= 228^3 + 423^3 \\
    &= 255^3 + 414^3 \\
    6963472309248 &= 2421^3 + 19083^3 \\
    &= 5436^3 + 18948^3 \\
    &= 10200^3 + 18072^3 \\
    &= 13322^3 + 16630^3. \\
\end{align*}
A: If you are looking to generate the list as opposed to using OEIS, a good place to start is to render $n^3\equiv n\bmod 6$. So $a^3+b^3=c^3+d^3\implies d\equiv a+b-c\bmod 6$ reducing your trials by a factor of $6$.
Slightly more sophisticated logic can be used to eliminate more of the chaff, using modulo $5$ and modulo $8$.
For modulo $5$, we develop cases corresponding to whether or not $a,b,c$ are each $\equiv0$. We use the face that for nonzero residues $r$, $r^3\equiv1/r$. Then:

*

*If $a\equiv b\equiv0$, then $d\equiv -c$.


*If $a\equiv c\equiv 0$ and $b\not\equiv0$, then $d\equiv b$.


*In all other cases we find that
$d\equiv(abc)(bc+ac-ab)^{-1}$
where the "reciprocal" function is defined by the mapping $(0,1,2,3,4)\to(0,1,3,2,4)$. In other words, the zero divisor is patched in by defining the quotient to be $0$, but only if no more than one of $a,b,c$ separately has zero residue.
Now look at modulo $8$. Here all cubes are $\equiv0$ if they are even or $\equiv\text{ the cube root }$ if they are odd. We then cibstruct the following cases:

*

*If $a,b,c$ are all even then $d$ is even.


*If $a,b$ are even and $c$ is odd, then $d\equiv -c\bmod 8$.


*If $a,c$ are even and $b$ is odd, then $d\equiv b\bmod 8$.


*If $a$ is even and $b,c$ are odd, we must have $b\equiv c\bmod 8$ and $d$ even, or else there is no solution.


*If $c$ is even and $a,b$ are odd, we must have $b\equiv -a\bmod 8$ and $d$ even, or else there is no solution.


*If $a,b,c$ are all odd, then $d\equiv a+b-c\bmod 8$.
When we implement the conditions above for $\bmod 6,5,8$, we find that each ordered triple $(a,b,c)\bmod 30$ gives zero, one or four residues $\bmod 120$ for $d$, with the average being one residue. Thereby the trials would be reduced by a factor of $120$.
A: One programming point here: it's probably better to find all numbers that are sums of two cubes and then check the list for duplicates, rather than iterating over all quadruples a, b, c, d.
For example, in Python, code to generate solutions in integers less than 20:
from collections import defaultdict
n = 20;
sums_of_cubes = [[x,y, x**3 + y**3] for x in range(1, n) for y in range(1, n) if x < y];

pairs_by_sum = defaultdict(list)

for z in sums_of_cubes:
  pairs_by_sum[z[2]].append([z[0], z[1]])

[[z, pairs_by_sum[z]] for z in pairs_by_sum if len(pairs_by_sum[z]) > 1]

This returns
[[1729, [[1, 12], [9, 10]]], [4104, [[2, 16], [9, 15]]]]

which gives the two smallest solutions.
