Diophantine equations and elliptic curves

在上一篇 post 中，我们介绍了在椭圆曲线上的 weil pairing 和 tate pairing。在这一篇中我们考虑椭圆曲线在丢番图方程上的一些应用。

问题零

我们在初中数学中学过，直角三角形满足勾股定理，设两个直角边为 $a$ , $b$ ，斜边为 $c$ 。则公式为：

a^{2} + b^{2} = c^{2}

那么满足三边均为整数的直角三角形有无数个吗？

这便是一个丢番图方程的经典例子。一般而言，丢番图方程（Diophantine equations）是有一个或者几个变量的整系数方程，它们的求解仅仅在整数范围内进行。

我们学到的最简单的方程应该是线性丢番图方程，即形如：
$a x + b y = c$
的方程，其中 $a, b, c$ 为整数， $x, y$ 为未知数。该方程可以使用扩展欧几里得算法来解决。

当然，上面这道题并不需要用到 elliptic curves 的内容来进行分析（所以才叫“问题零”）。

下面说一下问题零的解答，注意到：

(\frac{x^{2} - y^{2}}{x^{2} + y^{2}})^{2} + (\frac{2 x y}{x^{2} + y^{2}})^{2} = 1

这个其实是初中数学老师的某一堂课注意到的。

因此满足三边长度均为整数的直角三角形有无数个，令 $x > y > 0, x, y \in Z$ ，只需要取 $a = x^{2} - y^{2}, b = 2 x y, c = x^{2} + y^{2}$ 即可满足条件。

例如: $x = 2, y = 1 \to (3, 4, 5)$

$x = 3, y = 1 \to (6, 8, 10)$

$x = 3, y = 2 \to (5, 12, 13)$

问题一

我们在小学数学时就学过直角三角形的面积为

S = \frac{1}{2} a b

那么问题是，可以找到面积 $S = 5$ ，边长都为有理数的直角三角形吗？如果存在，那么这样的直角三角形是否有无穷多个？

这道题相当于我们有两个条件的 diophantine equation：

{\begin{cases} a^{2} + b^{2} = c^{2} \\ a b = 10 \end{cases}

要求其有理数解。如果我们如法炮制用 $x, y$ 来替换 $a, b$ ，结果就变成求：

(x^{2} - y^{2}) x y = 5

的有理数解，方程的次数变成了 $4$ ，看上去似乎比之前的方程更难解。于是我们换一个思路：

我们注意到：

(\frac{a \pm b}{2})^{2} = \frac{a^{2} \pm 2 a b + b^{2}}{4} = (\frac{c}{2})^{2} \pm 5

这个是书上说注意到的，这个应该比上一个显然。

因此我们设 $x = (\frac{c}{2})^{2}$ ，问题可以转换成：

只要找到一组解 $(x - 5, x, x + 5)$ 都是有理数的平方，就可以找到对应满足条件的 $(a, b, c)$ .

如果 $x - 5, x, x + 5$ 都是有理数的平方，那么这三者之积 $x^{3} - 25 x$ 也是有理数的平方。因此我们可以引出这道题的主角，即椭圆曲线：

E : y^{2} = x^{3} - 25 x

来寻找里面 non-trivial 的有理点。

找到第一个解

我们首先可以找到三个有理点 $(0, 0)$ 与 $(\pm 5, 0)$ ，但三角形的面积不能为零，所以这组解是 trivial 的。（并且 $5$ 和 $- 5$ 也不是有理数的平方）

通过简单的搜索，我们可以找到一个非平凡的有理点 $(- 4, 6)$ ，这个点虽然是 non-trivial 的，但这里找不到 $c$ 使得 $(\frac{c}{2})^{2} = - 4$ ，因此不存在这样的实数 $c$ 。

但我们可以从这个点出发做切线，来与 $E$ 相交找到其他的有理点（其实也就是椭圆曲线的倍点操作）。

我们尝试做点 $(- 4, 6)$ 沿 $E$ 的切线。首先我们使用隐函数求导：

2 y y^{'} = 3 x^{2} - 25

y^{'} = \frac{3 x^{2} - 25}{2 y}

代入点 $(- 4, 6)$ 后可得 $y^{'} = \frac{23}{12}$ ，因此我们有切线函数：

y = \frac{23}{12} (x + 4) + 6

将该函数与 $E$ 的方程相联立，有：

(\frac{23}{12} (x + 4) + 6)^{2} = x^{3} - 25 x

展开后可得常数项为 $- \frac{41^{2}}{9}$ 。根据韦达定理， $x_{1} x_{2} x_{3} = - C = \frac{41^{2}}{9}$ （这里 $C$ 指一元三次方程的常数项）。因为该切线切 $E$ 于 $x = - 4$ ，因此该方程有两个重根 $x_{1, 2} = - 4$ ，因此 $x_{3} = \frac{41^{2}}{9} \div (- 4)^{2} = (\frac{41}{12})^{2}$ 。

我们得到了对应的 $x_{3}$ ，就可以算出 $c = \frac{41}{6}$ ，而我们先前有：

(\frac{a \pm b}{2})^{2} = \frac{a^{2} \pm 2 a b + b^{2}}{4} = (\frac{c}{2})^{2} \pm 5

因此：

(\frac{a \pm b}{2})^{2} = (\frac{41}{12})^{2} \pm 5

可以发现 $(\frac{41}{12})^{2} \pm 5$ 也是有理数的平方，进而可以解得 $a = \frac{20}{3}, b = \frac{3}{2}$ 。

综上，我们得到了一个解 $(a, b, c) = (\frac{20}{3}, \frac{3}{2}, \frac{41}{6})$ ，可以验算一下确实满足勾股定理并且面积为 $5$ 。

证明解有无穷多个

下一个问题是，存在无数个这样的三元组 $(a, b, c)$ 吗？答案是肯定的，我们只需重复这样的过程即可（这里手算就不方便了）：

因为 $y = ((x - 5) x (x + 5))^{1 / 2} = \frac{(a - b) c (a + b)}{8} = \frac{(a^{2} - b^{2}) c}{8}$ ，所以可得 $a^{2} - b^{2} = 8 y / c$ . 联合 $a^{2} + b^{2} = c^{2}$ 可得 $a = \sqrt{(8 y / c + c^{2}) / 2}, b = \sqrt{(- 8 y / c + c^{2}) / 2}$ .

因此得到下一个三元组 $(\frac{4920}{1519}, \frac{1519}{492}, \frac{3344161}{747348})$ . 我们可以这样继续通过作切线生成其他符合条件的三元组。这里简单说明一下这种迭代方法始终可以保证 $a, b, c$ 均为有理数：

我们讨论类似这样的一般方程：

E : y^{2} = x^{3} - n^{2} x

中有非平凡有理数点 $(x_{0}, y_{0})$ （即 $x_{0} \neq 0, \pm n$ ），证明沿这个点做切线，于 $E$ 的另外一个交点 $(x_{1}, y_{1})$ 满足 $x_{1} - n, x_{1}, x_{1} + n$ 都是有理数的平方：

首先我们求出切线方程：

y = (\frac{3 x_{0}^{2} - n^{2}}{2 y_{0}}) (x - x_{0}) + y_{0}

根据韦达定理，可得：

x_{a} x_{b} x_{c} = (y_{0} - \frac{(3 x_{0}^{2} - n^{2}) x_{0}}{2 y_{0}})^{2}

= (\frac{2 y_{0}^{2} - 3 x_{0}^{3} + n^{2} x_{0}}{2 y_{0}})^{2}

= (\frac{2 x_{0}^{3} - 2 n^{2} x_{0} - 3 x_{0}^{3} + n^{2} x_{0}}{2 y_{0}})^{2}

= (\frac{x_{0}^{3} + n^{2} x_{0}}{2 y_{0}})^{2}

因为 $x_{a, b}$ 是重根 $x_{0}$ ，因此 $x_{1} = x_{c} = (\frac{x_{0}^{2} + n^{2}}{2 y_{0}})^{2}$ ，于是我们证明了 $x_{1}$ 是有理数的平方，因此 $c$ 是有理数。

注意到：

x_{1} \pm n = (\frac{(x_{0}^{2} + n^{2})^{2} \pm 4 n y_{0}^{2}}{4 y_{0}^{2}}) = (\frac{x_{0}^{2} \pm 2 n x_{0} - n^{2}}{2 y_{0}})^{2}

这个是 deepseek R1 注意到的，我注意力不行（

因此 $x_{1} \pm n$ 也是有理数的平方，这样的话因为前文提到：

(\frac{a \pm b}{2})^{2} = (\frac{c}{2})^{2} \pm n = x \pm n

即 $\frac{a \pm b}{2}$ 也都是有理数，最终我们得到了 $a, b$ 也都是有理数。

我们可以按照这种迭代方法生成 $x_{2}, x_{3}, \dots$ ，使得其对应的 $a, b, c$ 均为有理数。因此如果 $n$ 对应的椭圆曲线 $E$ 存在非平凡有理点，那么边为有理数、面积为 $n$ 的直角三角形就有无数个。Q.E.D.

问题二

这是我在初中时同学给我的一道钓鱼题：

我当时试了一下，发现这个问题能解决的人数占比应该没有 5%。事后同学告诉了我那三个巨长的答案，于是我意识到了一些数论问题的魅力：

看上去十分简单的问题陈述，其答案和解法却非常复杂。

化为 Weierstrass 形式

首先这个方程是齐次的，这意味着只要我们找到一组解 $(a, b, c)$ ，则对于任意 $t \in N_{+}$ ， $(t a, t b, t c)$ 都是这个方程的解。

按照问题一的思路，我们首先应该考虑如何将题目所给的方程：

\frac{a}{b + c} + \frac{b}{a + c} + \frac{c}{a + b} = 4, a, b, c \in N_{+}

转化为椭圆曲线的形式——这个可比问题一的难度大多了。

首先尝试通分：

a^{3} + b^{3} + c^{3} - 3 a^{2} b - 3 a^{2} c - 3 a b^{2} - 3 b^{2} c - 3 a c^{2} - 3 b c^{2} - 5 a b c = 0

然后我们需要用 $a, b, c$ 的线性组合来表示 $x, y$ ，将原式 $F$ 化为长 Weierstrass 形式：

y^{2} + a_{1} x y + a_{3} y = x^{3} + a_{2} x^{2} + a_{4} x + a_{6}

可以通过穷举找到 $F$ 的一个零点 $P = (1, - 1, 0)$ ,于是我们可以过 $P$ 做其切线。先求 $F$ 的偏导数：

{\begin{cases} \frac{\partial F}{\partial a} = 3 a^{2} - 6 a b - 6 a c - 3 b^{2} - 3 c^{2} - 5 b c \\ \frac{\partial F}{\partial b} = 3 b^{2} - 6 a b - 6 b c - 3 a^{2} - 3 c^{2} - 5 a c \\ \frac{\partial F}{\partial c} = 3 c^{2} - 6 a c - 6 b c - 3 a^{2} - 3 b^{2} - 5 a b \end{cases}

然后代入 $(1, - 1, 0)$ 求得了具体的偏导向量 $(6, 6, - 1)$ ，得到了具体的切线 $Z = 6 a + 6 b - c$ .

令 $Z = 0$ ，将 $c = 6 a + 6 b$ 代入 $F$ ，可得 $F = 91 (a + b)^{3}$ ，因为点 $P$ 满足 $a + b = 0$ ，因此零点 $P$ 的重数为 $3$ 。根据参考文献六，我们需要找到一个变换矩阵 $M_{β}$ 满足：

M_{β} = (\begin{matrix} Q_{x} & P_{x} & . \\ Q_{y} & P_{y} & . \\ Q_{z} & P_{z} & . \end{matrix})

是可逆矩阵。一个简单的做法是令点 $Q$ 在 $Z$ 上的另一个点，第三列为 $(1, 0, 0), (0, 1, 0)$ 或 $(0, 0, 1)$ 即可。

例如取 $M_{β}$ 为：

M_{β} = (\begin{matrix} 1 & 1 & 0 \\ 0 & - 1 & 0 \\ 6 & 0 & 1 \end{matrix})

令 $a = u + v, b = - v, c = 6 u + w$ ，我们将 $F (a, b, c)$ 变换为 $F (u, v, w)$ :

91 u^{3} + 69 u^{2} w - u v w - v^{2} w + 15 u w^{2} + w^{3} = 0

v^{2} w + u v w = 91 u^{3} + 69 u^{2} w + 15 u w^{2} + w^{3}

\frac{1}{91} v^{2} w + \frac{1}{91} u v w = u^{3} + \frac{69}{91} u^{2} w + \frac{15}{91} u w^{2} + \frac{1}{91} w^{3}

令 $w^{'} = w / 91$ ，即可化成整系数的长 weierstrass 形式：

v^{2} w^{'} + u v w^{'} = u^{3} + 69 u^{2} w^{'} + 1365 u w^{' 2} + 8281 w^{' 3}

两边同时除以 $w^{' 3}$ ，令 $X = u / w^{'}, Y = v / w^{'}$ ，可得：

Y^{2} + X Y = X^{3} + 69 X^{2} + 1365 X + 8281

这两者便是参考文献六一开始提到的投影形式和仿射形式。

这个时候我们尝试将 $(u, v, w)$ 用 $(a, b, c)$ 来表示：

{\begin{cases} u = a + b \\ v = - b \\ w = - 6 a - 6 b + c \end{cases}

然后就能得到 $(X, Y)$ 的 $(a, b, c)$ 表示了：

{\begin{cases} X = \frac{- 91 (a + b)}{6 a + 6 b - c} \\ Y = \frac{91 b}{6 a + 6 b - c} \end{cases}

同样我们也可以算出其逆变换（这里省略了自由变量 $k$ ）：

{\begin{cases} a = - X - Y \\ b = Y \\ c = - 6 X - 91 \end{cases}

找到一组正整数解

尝试暴力穷举寻找新的方程的整数解，可以通过脚本枚举 $X, Y \in [- 100, 100]$ ，总共找到以下几组：

(-39, 52)
(-39, -13)
(-28, 63)
(-28, -35)
(-15, 11)
(-15, 4)
(-14, 7)
(-13, 13)
(-13, 0)
(0, 91)
(0, -91)

把解带回 $(a, b, c)$ ，排除原式分母为零的情况，有以下几组整数解：

(x, y) = (-39, 52) -> (a, b, c) = (-13, 52, 143) [Valid]
(x, y) = (-39, -13) -> (a, b, c) = (52, -13, 143) [Valid]
(x, y) = (-28, 63) -> (a, b, c) = (-35, 63, 77) [Valid]
(x, y) = (-28, -35) -> (a, b, c) = (63, -35, 77) [Valid]
(x, y) = (-15, 11) -> (a, b, c) = (4, 11, -1) [Valid]
(x, y) = (-15, 4) -> (a, b, c) = (11, 4, -1) [Valid]

其中第五、六行就是参考资料四中使用的特殊解。

但此时对应的 $(a, b, c)$ 并不都是正整数，我们取最后的 $(x, y) = (- 15, 4)$ ，使用问题一中的切线法尝试寻找其他有理数点，尝试反解出 $(a, b, c)$ 来判断是不是正整数。

这里的过程就比较简单了，只需暴力循环无脑判断就行。这里我也懒得造轮子，直接让 GPT 写 sagemath 的椭圆曲线库搞定：

E = EllipticCurve(QQ, [1, 69, 0, 1365, 8281])  # [a1, a2, a3, a4, a6]

def compute_abc(x, y):
    a = -x-y
    b = y
    c = -6*x-91
    return a, b, c

def is_positive_integer(a, b, c):
    return (a > 0 and b > 0 and c > 0)

def find_positive_abc(P, max_iterations):
    iteration = 1
    
    while iteration < max_iterations:
        try:
            Q = (iteration * P)
            x, y = Q.xy()
        except Exception as e:
            print(f"Error in point doubling: {e}. Stopping.")
            break
        
        a, b, c = compute_abc(x, y)

        if is_positive_integer(a, b, c):
            print(f"Found positive integer solution when iteration = {iteration}")
            return Q, a, b, c
        
        iteration += 1
    
    print("No positive integer (a, b, c) found within max iterations.")
    return None


point = E([-15, 4])
result = find_positive_abc(point, max_iterations=10)
if result:
    Q, a, b, c = result
    print(f"Solution: (a, b, c) = ({a}, {b}, {c})")

与参考资料四相同，当迭代次数 iteration=9 时，输出为：

a = \frac{154476802108746166441951315019919837485664325669565431700026634898253202035277999}{12568549348783827476450917658335170150704695563440250577574035123362243258771752}

b = \frac{36875131794129999827197811565225474825492979968971970996283137471637224634055579}{12568549348783827476450917658335170150704695563440250577574035123362243258771752}

c = \frac{9405501114660716625534518421595417184664937929417781}{27028800585432577024141815890005886705385156000813842}

同时乘他们最小公倍数即可得到正整数解：

a = 154476802108746166441951315019919837485664325669565431700026634898253202035277999

b = 36875131794129999827197811565225474825492979968971970996283137471637224634055579

c = 4373612677928697257861252602371390152816537558161613618621437993378423467772036

问题三

有好事者在问题二梗图的启发下，把一个看上去更简单的问题搞成了这种人畜无害的形式：

这就是著名的“三立方和问题”，即求：

x^{3} + y^{3} + z^{3} = k

的整数解 $(x, y, z)$ 。

而图片中的问题，即 $k = 33$ 的情况，直到 2019 年才被 Andrew Booker 解决（问题二是在 2014 年，其实也不是很早）。在随后的一年，Andrew Sutherland 与 Booker 合作，解决了和为 42 的丢番图方程。这样一来，100 及以下的正整数存在解的都找到至少一组解了。

由于这个方程不是齐次方程，所以问题二的倍点方法又行不通了。Andrew 发的论文思路是通过降低暴力枚举的时间复杂度（从 $O (N^{1 + o (1)})$ 降到 $O_{k} (N (\log \log N) (\log \log \log N))$ ），再通过超算算了一个月（实际上是 23 core-years 的计算量）得到了 33 的一个解，数量级在 $N = 10^{16}$ 。所以这里就不太可能把代码写出来，讲解下大致思路即可。

$O (N^{1 + o (1)})$ 的算法

首先考虑问题的弱化版——两立方和问题，一个比较简单的结论是对于质数 $p$ 满足：

x^{3} + y^{3} = p

有整数解的 $p$ 为 $2$ 或形如 $3 n^{2} - 3 n + 1$ 的形式，证明如下：

因为 $x^{3} + y^{3} = (x + y) (x^{2} - x y + y^{2})$ ，因为 $p$ 为质数，故：
$x + y = 1 \lor x^{2} - x y + y^{2} = 1$
先考虑右式，只有 $x = 1, y = 1$ 满足条件。对应的 $p$ 为 $2$ .

再考虑左式，令 $y = 1 - x$ ，右边就成了 $3 x^{2} - 3 x + 1$ 了（当然得必须是质数）。

我们尝试将 $p$ 推广为任意整数 $k$ ，令 $k = r s$ ，于是我们有 $r = x y, s = x^{2} - x y + y^{2}$ 。令 $y = r - x$ ，可得：

s = 3 x^{2} - 3 r x + r^{2}

是个标准的一元二次方程。则 $Δ = (- 3 r)^{2} - 4 * 3 * (r^{2} - s) = 12 s - 3 r^{2}$ ， $Δ$ 为完全平方数是解为整数的必要条件。设 $\sqrt{Δ} = t$ ，则可得 $x = (3 r + t) / 6, y = (3 r - t) / 6$ ，检验 $x, y$ 是否为整数即可。

回到先前的三立方和问题，可以将 $z^{3}$ 移到右边：

x^{3} + y^{3} = k - z^{3}

然后应用之前的两立方和问题判断算法，只需枚举 $z$ 并做质因数分解，即可达到 $O (N^{1 + o (1)})$ 的时间复杂度。

转化为椭圆曲线

假设 $x^{3} + y^{3} + z^{3} = k > 0, | x | > | y | > | z | \geq \sqrt{k}$ 。我们令 $d = | x + y |$ ，则 $z$ 是 $k$ 模 $d$ 意义下的立方根，可以解出：

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

这里 $sgn (x)$ 表示 $x$ 的符号函数，满足 $sgn (0) = 0, sgn (x) = \frac{x}{| x |}$ 。

则

Δ = \frac{4 | k - z^{3} | - d^{3}}{3 d} = ◻

这里 $◻$ 符号表示该式算出来的值必须是完全平方数，因此等价于：

3 d (4 sgn (z) (z^{3} - k) - d^{3}) = ◻ \land z^{3} \equiv k \mod d

对于给定的搜索范围 $N$ ，例如本题的 $N = 10^{16}$ ，我们需要枚举 $d \in [0, (\sqrt[3]{2} - 1) N]$ 中与 $3$ 互质的数。对于每个 $d$ ，我们也需要枚举 $z \in [- N, N]$ ，满足 $z^{3} \equiv k \mod d$ 。只要我们找到一组解 $(d, z)$ ，便可映射到原来的解 $(x, y, z)$ 。

对于指定的整数 $k$ ，可以构造对应的椭圆曲线 $E_{d}$ :

E_{d} : Y^{2} = X^{3} - 2 (6 d)^{3} (d^{3} + 4 sgn (z) k)

使得对应的整点恰好是上式的整数解，具体的映射为

{\begin{cases} X = 12 d | z | \\ Y = (6 d)^{2} | x - y | \end{cases}

于是这个问题又转化为了在 $Q$ 上找椭圆曲线整点的问题。

这种代数变换对于较小的 $d$ 有一定作用。例如在 2000 年以前，尚未找到整数解 $(x, y, z)$ 的 $k$ 有：

33, 42, 114, 165, 390, 579, 627, 633, 732, 795, 906, 921, 975.

Booker 使用了代数软件 Magma 的椭圆曲线求整点功能，验证了在 $d \leq 40$ 时，整数 $(k, d)$ 对中只有 ${(579, 29), (579, 34), (975, 22)}$ 能在 $E_{d}$ 上找到整点。换句话说，就是除了 $k = 579, 975$ ，其余所有上面列出的 $k$ ，都不能在 $d \leq 40$ 时找到整数解。

但是椭圆曲线的参数 $d$ 为 $10^{16}$ 级别，并且这样的椭圆曲线 $E_{d}$ 本身又有 $10^{16}$ 个。所以单纯将问题转化为椭圆曲线的整点问题并没有降低问题本身的复杂度，我们还需进一步的优化。

进一步优化

回到上一节，为了求得满足

3 d (4 sgn (z) (z^{3} - k) - d^{3}) = ◻

的 $(d, z)$ 对，这里 $d \in [0, α N]$ （ $α$ 与前文相同，即 $\sqrt[3]{2} - 1$ ），存在以下两个时间复杂度瓶颈：

为了计算满足 $z^{3} \equiv k \mod d$ 的 $z$ ，我们需要 $d$ 的质因数分解。

可能的 $(d, z)$ 对有 $Ω (N \log N)$ 个。

论文除了使用椭圆曲线排除了一部分整数 $d$ 之外，为了减少需要遍历的 $(d, z)$ 对数量，更主要的方法还是 Montgomery's Batch Inversion Trick 批量计算三次剩余与使用 Legendre symbol 空间换时间的筛法。

首先说一下 batch inversion，给定 $k$ 个数 $a_{1}, a_{2}, \dots a_{k} \mod n$ ，我们需要同时计算 $a_{1}^{- 1}, a_{2}^{- 1}, \dots, a_{k}^{- 1} \mod n$ 。如果分开进行扩展欧几里得，那么时间复杂度就是 $O (k \log n)$ . Montgomery 的做法如下：

计算 $p_{i} = a_{1} a_{2} \dots a_{i} \mod n$ 。从 $p_{1}$ 计算到 $p_{k}$ 只需 $k - 1$ 次模乘法运算。
使用扩展欧几里得计算 $p_{k}^{- 1}$ ，时间复杂度 $O (\log n)$ 。
当 $i = k$ 时，计算 $a_{i}^{- 1} = p_{i - 1} p_{i}^{- 1}$ ，随后计算 $p_{i - 1}^{- 1} = p_{i}^{- 1} a_{i}$ 。
令 $p_{0} = 1$ ，对于 $i = k - 1, k - 2, \dots, 2$ 进行同样操作。

于是计算的时间复杂度从 $O (k \log n)$ 下降到 $O (k + \log n)$ .

回到瓶颈一，这里将 $d$ 进行质因数分解，得到 $d = p_{1}^{e_{1}} p_{2}^{e_{2}} \dots p_{m}^{e_{m}}$ ，然后求解每个三次剩余 $z_{i}^{3} = k \mod p_{i}$ 。

如果 $p = 2 \mod 3$ ，根只有一个，则 $z = k^{(2 p - 1) / 3} \mod p$ ，使用快速幂可以搞定，时间复杂度 $O (\log p)$ 。
如果 $p = 1 \mod 3$ ：
- 令 $p = 3^{w} m + 1, m ∤ 3, b = k^{m} \mod p$ 。
- 迭代计算 $b^{3}, b^{3^{2}}, \dots, b^{3^{v}}$ ，直到 $b^{3^{v}} = 1 \mod p$ ，当且仅当 $v < w$ 时立方根存在。
- 任取 $x$ ，使得 $a = x^{m} \mod p$ 的阶为 $3^{w}$ ，计算 $n = \log_{a} b$ ，则 $3 | n$ ， $a^{n / 3}$ 是 $b$ 的立方根。
- 取 $e \in {1, 2}$ 满足 $3 | (p - e m)$ ，则 $z = a^{n e} k^{(p - e m) / 3}$ . 另外两根为 $a^{3^{w - 1}} z, a^{2 \times 3^{w - 1}} z$ .

然后使用hensel 提升法扩展到 $z_{i}^{3} = k \mod p_{i}^{e_{i}}$ ，使用 CRT 扩展到 $z^{3} \mod d$ .

如果先前用 batch inversion 的时候把逆元存下来，并且使用各个质数的“线性组合”来列举 $d$ ，之后用 CRT 的时候就不用重复计算逆元了，这样计算 $d \leq α B$ 的三次根的总时间复杂度就降到了 $O (N)$ 。

考虑瓶颈二，令

Δ = 3 d (4 sgn (z) (z^{3} - k) - d^{3})

选取 $P = 3 (\log \log N) (\log \log \log N)$ ，计算 $[5, P]$ 之间所有质数乘积 $M$ 。如果 $Δ$ 是完全平方数，对于任意质数 $p ∣ M$ （这里的括号是勒让德记号）：

(\frac{Δ}{p}) = (\frac{3 d}{p}) (\frac{| z |^{3} - sgn (z) k - d^{3} / 4}{p})

如果事先计算 $| z |$ 模 $M$ 的剩余类，然后选择对任意质数 $p ∣ M$ ，以上等式都成立的那些剩余类。可以利用 Hasse 界证明，剩余类的数量不超过 $N (\log \log N) (\log \log \log N)$ ，因此最终的复杂度也降到了 $O_{k} (N (\log \log N) (\log \log \log N))$ .

除了复杂度的降低外，Andrew 在后续工作 Sums of three cubes 提到针对特定 $k$ （例如 $k = 33$ ），使用三次互反律作为筛法，使得剩余类只剩原来的 $\frac{1}{63.6}$ ，从而只需要遍历 $5.5 \times 10^{9}$ 个 $z$ ，大大减少了遍历的工作量。