解决33问题──将33写成3个整数的立方和

这篇文章内容翻译自论文 Cracking the problem with 33，论文研究了方程 $x^{3} + y^{3} + z^{3} = k$ 在一些小的 $k$ 值的解，并首次将33写成了3个整数的立方和。完成中文可以查看项目 qiwihui/cracking-the-problem-with-33。截止到目前，100以内的自然数就剩下42还没有找到关于立方和的整数解了！

Answer to the Ultimate Question of Life, the Universe, and Everything. – 42

以下是论文正文翻译：

解决33问题

作者：ANDREW R. BOOKER

摘要受到Tim Browning和Brady Haran的 Numberphile 视频“未解决的33问题“的启发，我们研究了方程 $x^{3} + y^{3} + z^{3} = k$ 在一些小的 $k$ 值的解。我们找到了 $k = 33$ 的第一个已知解。

1. 简介

令 $k$ 为正整数，其中 $k \equiv \pm 4 (mod 9)$ 。然后Heath-Brown [HB92] 推测有无限多的三元组 $(x ， y ， z) \in Z^{3}$ 满足

$k = x^{3} + y^{3} + z^{3} . (1)$

早在1954年就开始对（1）进行各种数值研究 [MW55]；请参阅 [BPTYJ07]，了解截至2000年的这些研究的历史。自那时起进行的计算由于Elkies [Elk00] 而被算法所主导。我们所知道的最新内容是Huisman [Hui16] 的论文，该论文确定了（1）的所有解，其中 $k \leq 1000$ 且 $max {∣ x ∣, ∣ y ∣, ∣ z ∣} \leq 1 0^{1} 5$ 。特别是，Huisman报告说除了13个 $k \leq 1000$ 的值以外的所有解决方案都是已知的：

$33, 42, 114, 165, 390, 579, 627, 633, 732, 795, 906, 921, 975. (2)$

Elkies的算法通过使用格基减少（lattice basis reduction）在Fermat曲线 $X^{3} + Y^{3} = 1$ 附近寻找有理点来工作；它非常适合同时找到许多 $k$ 值的解。在本文中，我们描述了一种在k值确定时更有效的不同方法。它的优点是可以找到所有具有最小坐标界限的解，而不是Elkies算法中的最大坐标。这总是产生搜索范围的非平凡的扩张（nontrivial expansion），因为除了可以单独考虑的有限多个例外之外，还有

$max {∣ x ∣, ∣ y ∣, ∣ z ∣} > 32 min {∣ x ∣, ∣ y ∣, ∣ z ∣}$

此外，根据经验，通常情况是其中一个变量比其他变量小得多，因此我们希望实际上增益更大。

我们的策略类似于一些早期的方法（特别参见 [HBLtR93]，[Bre95]，[KTS97] 和 [BPTYJ07]），并且基于观察： $k - z^{3} = x^{3} + y^{3}$ 的任何解都具有 $x + y$ 作为一个因子。相对于早期研究，我们的主要贡献是注意到，通过一些时间空间权衡，运行时间在高度边界内非常接近线性，并且在现代64位计算机上实现时非常实用。

更详细地说，假设 $(x ， y ， z)$ 是（1）的解，并且不失一般性，假设 $∣ x ∣ \geq ∣ y ∣ \geq ∣ z ∣$ 。然后我们有

$k - z^{3} = x^{3} + y^{3} = (x + y) (x^{2} - x y + y^{2})$

如果 $k - z^{3} = 0$ 则 $y = - x$ ，并且 $x$ 的每个值都产生一个解。否则，设 $d = ∣ x + y ∣ = ∣ x ∣ + y sgn x$ ，我们看到 $d$ 可以除 $∣ k - z^{3} ∣$ 并且

$\frac{k - z ^{3}}{d} = x^{2} - x y + y^{2} = x (2 x - (x + y)) + y^{2} = ∣ x ∣ (2∣ x ∣ - d) + (d - ∣ x ∣)^{2} = 3 x^{2} - 3 d ∣ x ∣ + d^{2}$

得到

${x, y} = {\frac{1}{2} sgn (k - z^{3}) (d \pm \frac{4∣ k - z ^{3} ∣ - d ^{3}}{3 d})}$

因此，给定 $z$ 的候选值，通过遍历 $∣ k - z^{3} ∣$ 的所有除数，有一个有效的程序来查找 $x$ 和 $y$ 的所有相应值。这个基本算法在假设整数分解的时间复杂度的标准启发式（standard heuristics）下，已经能在时间 $O (B^{1 + ε})$ 内找到满足 $min {∣ x ∣, ∣ y ∣, ∣ z ∣} \geq B$ 的所有解。在下一节中，我们将解释如何避免因子分解并更有效地实现相同目的。

感谢感谢Roger Heath-Brown提供了有用的意见和建议。

2. 方法

为了便于表示，我们假设 $k \equiv \pm 3 (mod 9)$ ；请注意，这适用于（2）中的所有 $k$ 。由于上述基本算法对于寻找小解是合理的，因此我们将假设 $∣ z ∣ > k$ 。此外，如果我们将（1）专门用于 $y = z$ 的解，那么我们得到Thue方程 $x^{3} + 2 y^{3} = k$ ，这是有效可解的。使用 PARI/GP [The18] 中的Thue求解器，我们验证了（2）中的 $k$ 不存在这样的解。因此，我们可以进一步假设 $y \neq = z$ 。

由于 $∣ z ∣ > k \geq 3 k$ ，我们有

$sgn z = - sgn (k - z^{3}) = - sgn (x^{3} + y^{3}) = - sgn x .$

同样，因为 $x^{3} + z^{3} = k - y^{3}$ 和 $∣ y ∣ \geq ∣ z ∣$ ，我们有 $sgn y = - sgn x = sgn z$ 。将（1）的两边乘以 $- sgn z$ ，我们得到

$∣ x ∣^{3} - ∣ y ∣^{3} - ∣ z ∣^{3} = - k sgn z (4)$

令 $α = 32 - 1$ ，并且 $d = ∣ x + y ∣ = ∣ x ∣ - ∣ y ∣$ 。如果 $d \geq α ∣ z ∣$ 则

$- k sgn z = ∣ x ∣^{3} - ∣ y ∣^{3} - ∣ z ∣^{3} \geq (∣ y ∣ + α ∣ z ∣)^{3} - ∣ y ∣^{3} - ∣ z ∣^{3} = 3 α (α + 2) (∣ y ∣ - ∣ z ∣) z^{2} + 3 α (∣ y ∣ - ∣ z ∣)^{2} ∣ z ∣ \geq 3 α (α + 2) ∣ y - z ∣ z^{2}$

由于 $3 α (α + 2) > 1$ ，这与我们的假设不相容，即 $y \neq = z$ 和 $∣ z ∣ > k$ 。因此我们必然有 $0 < d < α ∣ z ∣$ 。

接下来，减少（4）模3并回想我们的假设 $k \equiv \pm 3 (mod 9)$ ，我们有

$d = ∣ x ∣ - ∣ y ∣ \equiv ∣ z ∣ (mod 3) .$

设 $ϵ \in {\pm 1}$ 使得 $k \equiv 3 ϵ (mod 9)$ 。然后，由于每个立方数都与 $0$ 或 $\pm 1 (m o d 9)$ 相等，我们必然有 $x \equiv y \equiv z \equiv ϵ (mod 3)$ ，因此 $sgn z = ϵ (\frac{∣ z ∣}{3}) = ϵ (\frac{d}{3})$ 。基于（3），当且仅当 $d ∣ z^{3} - k$ 以及 $3 d (4∣ z^{3} - k ∣ - d^{3}) = 3 d (4 ϵ (\frac{d}{3}) (z^{3} - k) - d^{3})$ 是平方数时，我们得到（1）的解。

总之，找到（1）的所有解并且满足 $∣ x ∣ \geq ∣ y ∣ \geq ∣ z ∣ > k$ ， $y \neq = z$ 和 $∣ z ∣ \leq B$ ，对于每个与3互质的 $d \in Z \cap (0, α B)$ ，解决以下系统就足够了：

$\frac{d}{3 2 - 1} < ∣ z ∣ \leq B, sgn z = ϵ (\frac{d}{3}), z^{3} \equiv k (mod d) 3 d (4 ϵ (\frac{d}{3}) (z^{3} - k) - d^{3}) = □ (5)$

我们解决这个问题的方法很简单：我们通过它们的主要因子分解递归地计算 $d$ 的值，并应用中国剩余定理来将 $z^{3} \equiv k (mod d)$ 的解减少到素数模幂的情况下，其中标准算法可以适用。设 $r_{d} (k) = # {z (mod d) : z^{3} \equiv k (mod d)}$ 表示 $k$ 模 $d$ 的立方根数。通过标准分析估计，由于 $k$ 不是立方数，我们有

$d \leq α B \sum r_{d} (k) ≪_{k} B$

启发式地，计算对所有素数 $p \leq α B$ 的 $z^{3} \equiv k (mod p)$ 的解可以用 $[0, α B]$ 上的整数在 $O (B)$ 算术运算来完成；见例如 [[NZM91]，§2.9，练习8]中描述的算法。假设这一点，可以看出，使用Montgomery的批量反转技巧[[Mon87]，§10.3.1]，计算对所有正整数 $p \leq α B$ 的 $z^{3} \equiv k (mod p)$ 的根的剩余工作可以再次用 $O (B)$ 算术运算完成。

因此，我们可以在线性时间内计算满足（5）的第一行的所有 $z$ ，作为算术进展（arithmetic progressions）的并集。为了检测最后一行的解，有一个快速的方法来确定 $Δ := 3 d (4 ϵ (\frac{d}{3}) (z^{3} - k) - d^{3})$ 是一个平方数至关重要。我们首先注意到对于固定 $d$ ，这种情况减少到在椭圆曲线上找到积分点；特别是，令 $X = 12 d ∣ z ∣$ 和 $Y = (6 d^{2} ∣ x - y ∣$ ，从（3）中我们看到（X，Y）位于Mordell曲线上

$Y^{2} = X^{3} - 2 (6 d)^{3} (d^{3} + 4 ϵ (\frac{d}{3}) k) . (6)$

因此，对于固定 $d$ ，存在至多有限多个解，并且它们可以被有效地约束。对于 $d$ 的一些小值，找到（6）上的所有积分点并检查是否产生任何满足（1）的解是切实可行的。例如，使用Magma[[BCFS18]，§128.2.8]中的积分点函数（functionality），我们验证了如（2）中的 $k$ 和 $d \leq 40$ 情况下没有解，除了 $(k, d) \in {(579, 29), (579, 34), (975, 22)}$ 。

接下来我们自然注意到一些同余和可分性约束：

引理设 $z$ 为（5）的解，设 $p$ 为素数，设 $s = or d_{p} d$ ， $t = or d_{p} (z^{3} - k)$ 。则

(i) $z \equiv \frac{4}{3} k (2 - d^{2}) + 9 (k + d) (mod 18)$ ； (ii) 如果 $p \equiv 2 (mod 3)$ 则 $t \leq 3 s$ ； (iii) 如果 $t \leq 3 s$ 则 $s \equiv t (mod 2)$ ； (iv) 如果 $or d_{p} k \in {1, 2}$ 则 $s \in {0, or d_{p} k}$ 。

证明令 $Δ = 3 d (4 ϵ (\frac{d}{3}) (z^{3} - k) - d^{3})$ ，令 $δ = (\frac{d}{3})$ ，我们有 $∣ z ∣ \equiv d \equiv δ (mod 3)$ ，观察到 $(δ + 3 n)^{3} \equiv δ + 9 n (mod 27)$ ，模27，我们有

$\frac{Δ}{3 d} = 4 ϵδ (z^{3} - k) - d^{3} = 4∣ z ∣^{3} - d^{3} - 4 ϵδ k \equiv 4 [δ + 3 (∣ z ∣ - δ)] - [δ + 3 (d - δ)] - 4 ϵδ k = 3 (4∣ z ∣ - d) - δ [18 + 4 (ϵ k - 3)] \equiv 3 (4∣ z ∣ - d) - d [18 + 4 (ϵ k - 3)] = 12∣ z ∣ - 9 d - 4 ϵ d k \equiv 3∣ z ∣ - 4 ϵ d k$

这消失了模9，所以为了使 $Δ$ 成为平方数，它也必须消除mod 27。于是

$z = ϵδ ∣ z ∣ \equiv \frac{4 δ d k}{3} \equiv \frac{4 ( 2 - d ^{2} ) k}{3} (mod 9)$

减少（1）模2我们得到 $z \equiv k + d (mod 2)$ ，这得到（i）。

接下来设 $u = p^{- s} d$ 和 $v = p^{- t} ϵδ (z^{3} - k)$ ，这样就有

$Δ = 3 (4 p^{s + t} uv - p^{4 s} u^{4})$

如果 $3 s < t$ 则 $p^{- 4 s} Δ \equiv - 3 u^{4} (mod 4 p)$ ，但是当 $p \equiv 2 (mod 3)$ 时这是不可能的，因为 $- 3$ 不是 $4 p$ 的平方模。因此，在这种情况下我们必须 $t < 3 s$ 。

接下来假设 $t < 3 s$ 。我们考虑以下情况，涵盖所有可能性：

若 $p = 3$ 则 $s = t = 0$ ，那么 $s \equiv t (mod 2)$ 。
若 $p \neq = 3$ 且 $3 s > t + 2 ord_{p} 2$ ，则 $ord_{p} Δ = s + t + 2 ord_{p} 2$ ，那么 $s \equiv t (mod 2)$ 。
若 $3 s \in {t, t + 2}$ 则 $s \equiv t (mod 2)$ 。
如果 $p = 2$ 且 $3 s = t + 1$ 则 $2^{- 4 s} Δ = 3 (2 uv - u^{4}) \equiv 3 (mod 4)$ ，这是不可能的。

因此，在任何情况我们得出结论 $s \equiv t (mod 2)$ 。

最后，假设 $p ∣ k$ 和 $p \neq ∣3 k$ 。如果 $s = 0$ 则无需证明的，所以假设不然。由于 $d ∣ z^{3} - k$ ，我们必须有 $d ∣ k$ ，因为

$0 < s \leq t = ord_{p} (z^{3} - k) = ord_{p} k < 3 s$

通过部分（iii）得出 $s \equiv ord_{p} k (mod 2)$ ，因此 $s = ord_{p} k$ 。

因此，一旦 $z (mod d)$ 的残差类（residue class）固定，则其残差模 $l c m (d, 18)$ 是确定的。还要注意，条件（ii）和（iii）对于测试 $p = 2$ 是有效的。

然而，即使有这些优化，也有 $≪ B lo g B$ 对 $d, z$ 满足（5）的第一行和引理的结论（i）和（iv）。因此，为了实现比 $O (B lo g B)$ 更好的运行时间，需要从一开始就消除一些 $z$ 值。我们通过标准的时间空间交换来实现这一目标。确切地说，设置 $P = 3 (lo g lo g B) (lo g lo g lo g B)$ ，并且让 $M = \prod_{5 \leq p \leq P} p$ 是区间 $[5, P]$ 之间的素数的乘积。根据素数定理，我们得到 $lo g M = (1 + o (1)) P$ 。如果 $Δ$ 是平方数，那么对于任意素数 $p ∣ M$ 我们有

$(\frac{Δ}{p}) = (\frac{3 d}{p}) (\frac{∣ z ∣ ^{3} - c}{p}) \in {0, 1} (7)$

其中 $c \equiv ϵ (\frac{d}{3}) k + \frac{d ^{3}}{4}$ 。当 $lcm (d, 18) \leq α B / M$ 时，我们首先为每个残差类 $∣ z ∣ (mod M)$ 计算该函数，并且仅选择对于每个 $p ∣ M$ 满足（7）的那些残基。由Hasse约束，允许的残差的数量最多为

$\frac{M}{2 ^{ω (M / (M, d))}} p ∣ \frac{M}{( M , d )} \prod (1 + O (\frac{1}{p})) = \frac{M}{2 ^{ω (M / (M, d))}} e^{O (P / l o g P)}$

因此，要考虑的 $z$ 值的总数最多为

$\sum_{lcm (d, 18) \leq \frac{α B}{M}} r_{d} (k) [M + \frac{e ^{O (P / l o g P)}}{2 ^{ω (M / (M, d))}} \frac{α B}{d}] + \sum_{d \leq α B, l c m (d, 18) \leq \frac{α B}{M}} \frac{r _{d} ( k ) α B}{d} ≪_{k} B lo g M + \frac{e ^{O (P / l o g P)}}{2 ^{ω (M)}} \sum_{g ∣ M} \frac{2 ^{ω (g)} r _{g} ( k )}{g} \sum_{d^{'} \leq \frac{α B}{9 g M}} \frac{r _{d^{'}} ( k ) α B}{d ^{'}} ≪_{k} B lo g M + B lo g B \frac{e ^{O (P / l o g P)}}{2 ^{ω (M)}} \prod_{p ∣ M} (1 + \frac{2 r _{p} ( k )}{p}) ≪ BP + \frac{B l o g B}{2 ^{(1 + o (1)) P / l o g P}} ≪ B (lo g lo g B) (lo g lo g lo g B)$

对于没有以这种方式消除的 $z$ ，我们遵循类似的策略，其中一些其他辅助模 $M^{'}$ 由较大的素数组成，以加速平方测试。我们预先计算模为 $M^{'}$ 的立方数表和Legendre符号模 $p ∣ M^{'}$ ，因此将测试（7）简化为了表查找。只有当所有这些测试都通过时，我们才能在多精度算术中计算 $Δ$ 并应用一般的平方检验，这种情况对于一小部分候选值来说都是如此。事实上，我们期望Legendre测试的数量平均有限，所以总的来说，找到所有解决方案的 $∣ z ∣ \leq B$ 应该要求不超过 $O_{k} (B (lo g lo g B) (lo g lo g lo g B))$ 次表查找和对 $[0, B]$ 中整数的算术运算。

因此，当 $B$ 符合机器字大小时，我们预计运行时间几乎是线性的，这就是我们在实践中观察到的 $B < 2^{64}$ 。

3. 实现

我们在C中实现了上述算法，其中有一些内联汇编程序来源于由Ben Buhrow [Buh19] 编写的Montgomery算法 [Mon85]，以及Kim Walisch的用于枚举素数的 primesieve 库 [Wal19]。

该算法自然地在具有超过 $α B$ 的素因子和具有 $α B$ -平滑的素数的 $d$ 的值之间分配。前一组 $d$ 消耗超过运行时间的三分之二，但更容易并行化。我们在布里斯托大学高级计算研究中心的大规模并行集群Bluecrystal Phase 3上运行了这一部分。对于平滑的 $d$ ，我们使用了一个单独的32核和64核节点的小集群。