Main reference is an answer to the following question:
According to a comment by the same author, the method in this answer can be employed as well for solving the following problem:
-
What is the probability that a natural number is a sum of three cubes?
The reason is that $1/n = (1/n^{1/3})^3$, so an analogous method should work for three cubes.
We define an initial segment of the naturals with length $N$ and count all sums of three cubes
in that segment. The idea is that if $N\to\infty$, then we have the desired probability. (Oh yeah,
it assumed that zero is part of the naturals as well) The accompanying program is in Pascal:
program kubussen;
procedure test(N : integer; show : boolean);
var
i,j,k,t,L : integer;
memo : array of boolean;
begin
t := 0;
i := 0;
SetLength(memo,N);
while true do
begin
if sqr(i)*i > N-1 then Break;
j := i;
while true do
begin
if sqr(i)*i+sqr(j)*j > N-1 then Break;
k := j;
while true do
begin
L := sqr(i)*i+sqr(j)*j+sqr(k)*k;
if L > N-1 then Break;
if show then
begin
Writeln(L:3,' = ',i,'^3 + ',j,'^3 + ',k,'^3');
if memo[L] then Writeln('duplicate');
end;
memo[L] := true;
k := k + 1;
t := t + 1;
end;
j := j + 1;
end;
i := i + 1;
end;
Writeln(t,'/',N);
end;
begin
test(10,FALSE);
test(100,TRUE);
test(1000,FALSE);
test(100000,FALSE);
test(1000000,FALSE);
test(10000000,FALSE);
test(100000000,FALSE);
test(1000000000,FALSE);
end.
Sample output, details (show=TRUE) for $N=100$ :
6/10
0 = 0^3 + 0^3 + 0^3
1 = 0^3 + 0^3 + 1^3
8 = 0^3 + 0^3 + 2^3
27 = 0^3 + 0^3 + 3^3
64 = 0^3 + 0^3 + 4^3
2 = 0^3 + 1^3 + 1^3
9 = 0^3 + 1^3 + 2^3
28 = 0^3 + 1^3 + 3^3
65 = 0^3 + 1^3 + 4^3
16 = 0^3 + 2^3 + 2^3
35 = 0^3 + 2^3 + 3^3
72 = 0^3 + 2^3 + 4^3
54 = 0^3 + 3^3 + 3^3
91 = 0^3 + 3^3 + 4^3
3 = 1^3 + 1^3 + 1^3
10 = 1^3 + 1^3 + 2^3
29 = 1^3 + 1^3 + 3^3
66 = 1^3 + 1^3 + 4^3
17 = 1^3 + 2^3 + 2^3
36 = 1^3 + 2^3 + 3^3
73 = 1^3 + 2^3 + 4^3
55 = 1^3 + 3^3 + 3^3
92 = 1^3 + 3^3 + 4^3
24 = 2^3 + 2^3 + 2^3
43 = 2^3 + 2^3 + 3^3
80 = 2^3 + 2^3 + 4^3
62 = 2^3 + 3^3 + 3^3
99 = 2^3 + 3^3 + 4^3
81 = 3^3 + 3^3 + 3^3
29/100
178/1000
13142/100000
124398/1000000
1213524/10000000
11991346/100000000
119249735/1000000000
Trivial duplicates are:
$$
i^3+j^3+k^3 = i^3+k^3+j^3 = k^3+i^3+j^3 = j^3+i^3+k^3 = j^3+k^3+i^3 = k^3+j^3+i^3
$$
They are eliminated in the program by starting loops accordingly : $i\ge 0,j\ge i,k\ge j$ .
If these duplicates are not eliminated, then the outcome must be divided by six.
Non-trivial duplicates are retained, though. For example:
$$
855 = 1^3 + 5^3 + 9^3 = 0^3 + 7^3 + 8^3 \\
729 = 1^3 + 6^3 + 8^3 = 0^3 + 0^3 + 9^3 \\
251 = 2^3 + 3^3 + 6^3 = 1^3 + 5^3 + 5^3 \\
216 = 3^3 + 4^3 + 5^3 = 0^3 + 0^3 + 6^3 \\
344 = 4^3 + 4^3 + 6^3 = 0^3 + 1^3 + 7^3
$$
With help of the above, it is conjectured that, in the limit for $N\to\infty$, the desired probability
is equal to the following integral, over the domain $x\ge 0,y\ge 0,z\ge 0, x^3+y^3+z^3 \le 1$ :
$$
\frac{1}{6} \iiint dx\,dy\,dz
$$
But here I'm stuck. My knowledge about multiple integrals has become rusty and I have no idea
how to calculate other than numerically (and the latter is rather prohibitive as well).
It is suggested by the sample output that the outcome must be somewhat close to $0.119249735$ ,
but how good has it converged? So the question is: does there exist an exact solution for the integral?
The integral is a little bit simpler than I first thought:
$$
\frac{1}{6} \iiint dx\,dy\,dz =
\int_{x=0}^{x=1} \left[ \int_{y=0}^{y=(1-x^3)^{1/3}} (1-x^3-y^3)^{1/3}\, dy \right] dx
$$
This is what MAPLE
says about the (indefinite) inner integral:
> int((1-x^3-y^3)^(1/3),y);
3
3 (1/3) y
(1 - x ) y hypergeom([-1/3, 1/3], [4/3], - -------)
3
-1 + x
The "hypergeom" is abacedabra for me. And if not: we still have the outer integral, which
seemingly cannot be evaluated exactly. So let's evaluate the integral numerically, with help of MAPLE again:
> evalf(Int(Int((1-x^3-y^3)^(1/3),y=0..(1-x^3)^(1/3)),x=0..1))/6;
0.1186788238
Which I find satisfactory close to the last iterand ($0.119249735$) of the program.