最近 $x^3 + y^3 + z^3 = 33$ 的整数解被找到了：
$$8866128975287528^3 + (-8778405442862239)^3 + (-2736111468807040)^3 = 33,$$
着实让三立方和问题火了一把。到目前为止还没人能从理论上解决这个问题，在三立方和问题上的成果其实都只是各种计算机暴力算法罢了，只不过人们不断地用更多的计算资源，以及尽量不那么「暴力」的暴力算法去尝试增加搜索上界，以找到更多的整数解。
在这次算出整数解33的Andrew Booker之前，已经有许多人在这个问题上作出过尝试了。Elsenhans and Jahnel (2009)的搜索上界 $\max(|x|,|y|,|z|)$ 是 $10^{14}$ ，Huisman (2016)将上界提升到了 $10^{15}$（并找到了 $x^3 + y^3 + z^3 = 74$ 的整数解），这次Brooker则是把上界进一步提高到了 $10^{16}$。Elsenhans and Jahnel与Huisman用的是同一种方法（该方法最早由Noam Elkies提出）。本来直接放上论文链接就结束了，但他们的文章中对具体方法的实现细节着墨很少，我也是费了一番功夫才感觉大致理解了他们的方法，所以这里就简单讲下我个人的理解。
首先，不妨假定 $|x|\leq |y|\leq |z|$ 。问题等价于 $x^3 + y^3 = z^3 + n$，此处 $n$ 要远小于 $|x|,|y|,|z|$（否则就直接穷举了）。两边同时除以 $z^3$，
$$\Big(\frac{x}{z}\Big)^3 + \Big(\frac{y}{z}\Big)^3 = 1 + \frac{n}{z^3}.$$
此处 $n$ 也可以是负数，我们可以不失一般性地假定 $0<x<y<z$。由于 $n$ 要比 $z$ 小的多，因而 $n/z^3$ 会是个很小的数。令 $X=x/z$、$Y=y/z$ ，于是问题就变为了在曲线 $Y = (1 - X^3 )^{ 1/3 }$ 附近找有理数点。因为 $x<y$，只需搜索 $0<X<1/\sqrt[3]{2}$ 之间的值。
对于曲线上的一点 $(X_0,Y_0)$，可以在其四周定义一个很小的平行四边形，其中两条边平行于 $Y$ 轴，另两条边平行于该点处切线。这个平行四边形可以由
$$|X-X_0|\leq\epsilon,\quad |Y-a(X_0)X-b(X_0)|\leq\epsilon'$$
表示。$\epsilon$ 和 $\epsilon'$ 的值是根据搜索的 $n$ 和 $z$ 的范围来确定的。接下来要做的就是在 $0<X<1/\sqrt[3]{2}$ 范围内的一个个小平行四边形中找到所有的有理数点 $(x/z,y/z)$，然后计算对应的 $x^3 + y^3 - z^3$，看看是不是33或-33就行了。
令 $z$ 的搜索上界为 $N$，即 $0<z\leq N$。由上面的两个式子可以得到
$$|x-X_0z|\leq\epsilon z \leq \epsilon N,$$
$$|y-a(X_0)x-b(X_0)z|\leq\epsilon'z\leq\epsilon'N.$$
加上 $0<z\leq N$，一共有三个约束条件。于是，我们其实是要解一个线性方程组
$$\begin{pmatrix} 1 & 0 & -X_0 \\ -a(X_0) & 1 & -b(X_0) \\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} x \\ y \\ z \end{pmatrix}=\begin{pmatrix} \pm \epsilon N \\ \pm \epsilon'N \\ N \end{pmatrix},$$
上面这个方程组共有四个解，加上原点正好组成一个棱锥，然后在棱锥中遍历所有的整数点就好了。但问题在于这个棱锥的高度会远大于另两个锥度，也就是说棱锥的顶点会非常尖。这样直接遍历的效率相当低。为了提高效率，我们可以把上面的方程组改写为
$$\begin{pmatrix} \frac{1}{\epsilon N} & 0 & -\frac{X_0}{\epsilon N} \\ -\frac{a(X_0)}{\epsilon'N} & \frac{1}{\epsilon'N} & -\frac{b(X_0)}{\epsilon'N} \\ 0 & 0 & \frac{1}{N} \end{pmatrix}\begin{pmatrix} x \\ y \\ z \end{pmatrix}=\begin{pmatrix} \pm 1 \\ \pm 1 \\ 1 \end{pmatrix}.$$
令这个矩阵的三个列向量分别为 $\mathbf{b}_1,\,\mathbf{b}_2,\,\mathbf{b}_3$。这三个向量可以看作是一个格（lattice）$\mathcal{L}=\{x\mathbf{b}_1+y\mathbf{b}_2+z\mathbf{b}_3\,|\,x,y,z\in\mathbb{Z}\}$ 的格基。现在问题就转化为了找到 $\mathcal{L}$ 上原点附近的一些格点对应的坐标 $(x,y,z)$ 。这是一个已知的问题。先通过 LLL（Lenstra–Lenstra–Lovász）格基规约算法计算出 $\mathcal{L}$ 上近似正交的一组格基 $\mathbf{b}_1^{*},\, \mathbf{b}_2^{*},\, \mathbf{b}_3^{*}$，再使用Fincke-Pohst算法就能找到所有满足条件的 $(x,y,z)$ 了。
Huisman (2016)提到他将上界提高到 $10^{15}$ 总共用了大约10万个CPU小时。我自己的科研中有些比较大的case 会跑上百万CPU小时，所以按同样的方法把上界提高到 $10^{16}$ 甚至 $10^{17}$ 应该还是很可行的。不过Booker这次用的方法和之前的不太一样。上面说的方法适合于在一定范围的 $n$ （比如 $n\leq 1000$）中找到所有的整数解。但如果只是想针对具体的 $n$ （比如 33）来算的话，Booker的方法更有效率。他计算了 $n$ 等于33和42的情况，在 $\min(|x|,|y|,|z|)\leq 10^{16}$ 的范围内找到了 $n=33$ 的一个解。这样，42就成了100以内仅剩的一个还无解的数了。他的计算用了十几万CPU小时，比Huisman (2016) 稍多一些，但在一个量级上。