Diophantine equations and elliptic curves

In the previous post, we introduced the Weil pairing and Tate pairing on elliptic curves. In this article, we explore some applications of elliptic curves in Diophantine equations.

Problem Zero

In middle school mathematics, we learned that right triangles satisfy the Pythagorean theorem. Let the two legs be $a$ and $b$ , and the hypotenuse be $c$ . The formula is:

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

So, are there infinitely many right triangles with all three sides being integers?

This is a classic example of a Diophantine equation. Generally, Diophantine equations are equations with one or more variables and integer coefficients, and their solutions are sought only within the integer domain.

The simplest equations we learn are linear Diophantine equations, which take the form:
$a x + b y = c$
where $a, b, c$ are integers, and $x, y$ are unknowns. This equation can be solved using the extended Euclidean algorithm.

Of course, the above problem does not require the use of elliptic curves for analysis (which is why it's called "Problem Zero").

Below is the solution to Problem Zero. Note that:

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

This was actually noticed in one of my middle school math teacher's classes.

Therefore, there are infinitely many right triangles with all three sides being integers. Let $x > y > 0, x, y \in Z$ , and simply take $a = x^{2} - y^{2}$ , $b = 2 x y$ , $c = x^{2} + y^{2}$ to satisfy the condition.

For example: $x = 2, y = 1 \to (3, 4, 5)$

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

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

Question 1

In elementary school mathematics, we learned that the area of a right triangle is:

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

So the question is: can we find a right triangle with area $S = 5$ and all side lengths rational? If such a triangle exists, are there infinitely many such right triangles?

This problem is equivalent to solving a Diophantine equation with two conditions:

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

We are to find its rational solutions. If we try to substitute $a$ and $b$ with $x$ and $y$ as before, the equation becomes:

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

The degree of the equation becomes $4$ , which seems more difficult to solve than the previous one. So let's try a different approach:

We observe that:

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

This is noted in the textbook, and it should be more obvious than the previous one.

Therefore, let $x = {(\frac{c}{2})}^{2}$ . The problem can be transformed into:

As long as we find a solution where $(x - 5, x, x + 5)$ are all rational squares, we can find the corresponding $(a, b, c)$ that satisfy the conditions.

If $x - 5$ , $x$ , and $x + 5$ are all rational squares, then their product $x^{3} - 25 x$ is also a rational square. This leads us to the main character of this problem—the elliptic curve:

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

We aim to find non-trivial rational points on this curve.

Finding the First Solution

We can first identify three rational points: $(0, 0)$ and $(\pm 5, 0)$ . However, the area of the triangle cannot be zero, so this set of solutions is trivial. (Moreover, $5$ and $- 5$ are not squares of rational numbers.)

Through a simple search, we can find a non-trivial rational point $(- 4, 6)$ . Although this point is non-trivial, there is no real number $c$ such that $(\frac{c}{2})^{2} = - 4$ , so no such real $c$ exists.

But we can use this point to draw a tangent to $E$ and find other rational points (essentially performing a point doubling operation on the elliptic curve).

We attempt to draw the tangent to $E$ at the point $(- 4, 6)$ . First, we use implicit differentiation:

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

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

Substituting the point $(- 4, 6)$ gives $y^{'} = \frac{23}{12}$ , so the tangent function is:

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

Substituting this function into the equation of $E$ , we get:

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

Expanding this yields a constant term of $- \frac{41^{2}}{9}$ . By Vieta's formulas, $x_{1} x_{2} x_{3} = - C = \frac{41^{2}}{9}$ (where $C$ is the constant term of the cubic equation). Since the tangent touches $E$ at $x = - 4$ , the equation has a double root $x_{1, 2} = - 4$ , so $x_{3} = \frac{41^{2}}{9} \div (- 4)^{2} = (\frac{41}{12})^{2}$ .

With the corresponding $x_{3}$ , we can compute $c = \frac{41}{6}$ . Previously, we had:

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

Therefore:

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

We observe that $(\frac{41}{12})^{2} \pm 5$ are also squares of rational numbers, allowing us to solve for $a = \frac{20}{3}, b = \frac{3}{2}$ .

In summary, we obtain a solution $(a, b, c) = (\frac{20}{3}, \frac{3}{2}, \frac{41}{6})$ , which can be verified to satisfy the Pythagorean theorem and yield an area of $5$ .

Proving There Are Infinitely Many Solutions

The next question is: are there infinitely many such triples $(a, b, c)$ ? The answer is yes, and we can simply repeat this process (though it becomes inconvenient to do by hand):

Since $y = ((x - 5) x (x + 5))^{1 / 2} = \frac{(a - b) c (a + b)}{8} = \frac{(a^{2} - b^{2}) c}{8}$ , we have $a^{2} - b^{2} = 8 y / c$ . Combined with $a^{2} + b^{2} = c^{2}$ , we get $a = \sqrt{(8 y / c + c^{2}) / 2}, b = \sqrt{(- 8 y / c + c^{2}) / 2}$ .

Thus, we obtain the next triple $(\frac{4920}{1519}, \frac{1519}{492}, \frac{3344161}{747348})$ . We can continue generating other valid triples by drawing tangents in this way. Here is a brief explanation of why this iterative method always ensures that $a, b, c$ are rational numbers:

We consider a general equation of the form:

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

with a non-trivial rational point $(x_{0}, y_{0})$ (i.e., $x_{0} \neq 0, \pm n$ ). We prove that by drawing the tangent at this point, the other intersection point $(x_{1}, y_{1})$ with $E$ satisfies that $x_{1} - n, x_{1}, x_{1} + n$ are all rational squares:

First, we find the tangent line:

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

By Vieta's formulas, we have:

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}

Since $x_{a}$ and $x_{b}$ are both $x_{0}$ (a double root), we have $x_{1} = x_{c} = {(\frac{x_{0}^{2} + n^{2}}{2 y_{0}})}^{2}$ , proving that $x_{1}$ is a rational square, and thus $c$ is rational.

Note that:

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}

This was noticed by DeepSeek R1—my attention wasn't sharp enough (

Thus, $x_{1} \pm n$ are also rational squares. As mentioned earlier:

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

That is, $\frac{a \pm b}{2}$ are also rational numbers. Ultimately, we conclude that $a$ and $b$ are rational.

We can generate $x_{2}, x_{3}, \dots$ in this iterative manner, ensuring that the corresponding $a, b, c$ are all rational. Therefore, if the elliptic curve $E$ corresponding to $n$ has a non-trivial rational point, then there are infinitely many right triangles with rational sides and area $n$ . Q.E.D.

Question Two

This is a fishing question a classmate gave me in middle school:

I tried it back then and realized that probably less than 5% of people could solve it. Later, my classmate told me those three extremely long answers, and I began to appreciate the charm of some number theory problems:

A seemingly simple problem statement can have incredibly complex solutions and answers.

Transforming into Weierstrass Form

First, this equation is homogeneous, meaning that once we find a solution $(a, b, c)$ , then for any $t \in N_{+}$ , $(t a, t b, t c)$ is also a solution.

Following the idea from Problem 1, we should first consider how to transform the given equation:

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

into the form of an elliptic curve—this is significantly more challenging than Problem 1.

First, attempt to combine the fractions:

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

Then, we need to express $x$ and $y$ using linear combinations of $a, b, c$ to transform the original expression $F$ into the long Weierstrass form:

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

A zero point $P = (1, - 1, 0)$ of $F$ can be found through enumeration. We then draw the tangent line to $F$ at $P$ . First, compute the partial derivatives of $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}

Substitute $(1, - 1, 0)$ to obtain the specific gradient vector $(6, 6, - 1)$ , which gives the tangent line $Z = 6 a + 6 b - c$ .

Set $Z = 0$ , substitute $c = 6 a + 6 b$ into $F$ , and obtain $F = 91 (a + b)^{3}$ . Since point $P$ satisfies $a + b = 0$ , the multiplicity of the zero point $P$ is $3$ . According to Reference 6, we need to find a transformation matrix $M_{β}$ such that:

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

is an invertible matrix. A simple approach is to take point $Q$ as another point on $Z$ , and the third column as $(1, 0, 0)$ , $(0, 1, 0)$ , or $(0, 0, 1)$ .

For example, take $M_{β}$ as:

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

Let $a = u + v$ , $b = - v$ , $c = 6 u + w$ . Transform $F (a, b, c)$ into $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}

Let $w^{'} = w / 91$ to obtain an integer-coefficient long Weierstrass form:

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

Divide both sides by $w^{' 3}$ , and let $X = u / w^{'}$ , $Y = v / w^{'}$ to obtain:

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

These are the projective and affine forms mentioned at the beginning of Reference 6.

Now, express $(u, v, w)$ in terms of $(a, b, c)$ :

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

Then, we can express $(X, Y)$ in terms of $(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}

Similarly, we can compute the inverse transformation (omitting the free variable $k$ ):

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

Finding a Set of Positive Integer Solutions

Attempting a brute-force enumeration to find integer solutions for the new equation, a script can iterate over $X, Y \in [- 100, 100]$ , yielding the following sets:

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

Substituting these solutions back into $(a, b, c)$ and excluding cases where the denominator of the original expression is zero, the following integer solutions are obtained:

(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]

The fifth and sixth rows correspond to the special solutions used in Reference 4.

However, the corresponding $(a, b, c)$ values are not all positive integers. Taking the last solution $(x, y) = (- 15, 4)$ , we attempt to use the tangent method from Problem 1 to find other rational points and then solve back for $(a, b, c)$ to determine if they are positive integers.

This process is relatively straightforward, involving a brute-force loop for judgment. To avoid reinventing the wheel, I directly used GPT to write a script leveraging SageMath's elliptic curve library:

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})")

As in Reference 4, when the iteration count iteration=9, the output is:

a = \frac{154476802108746166441951315019919837485664325669565431700026634898253202035277999}{12568549348783827476450917658335170150704695563440250577574035123362243258771752}

b = \frac{36875131794129999827197811565225474825492979968971970996283137471637224634055579}{12568549348783827476450917658335170150704695563440250577574035123362243258771752}

c = \frac{9405501114660716625534518421595417184664937929417781}{27028800585432577024141815890005886705385156000813842}

Multiplying by their least common multiple yields the positive integer solutions:

a = 154476802108746166441951315019919837485664325669565431700026634898253202035277999

b = 36875131794129999827197811565225474825492979968971970996283137471637224634055579

c = 4373612677928697257861252602371390152816537558161613618621437993378423467772036

Question 3

Inspired by the meme in Question 2, a curious individual turned an apparently simpler problem into this seemingly harmless form:

This is the famous "sums of three cubes problem," which seeks integer solutions $(x, y, z)$ to the equation:

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

The problem depicted in the image corresponds to the case where $k = 33$ , which was not solved until 2019 by Andrew Booker (Question 2 was solved in 2014, which is also relatively recent). In the following year, Andrew Sutherland collaborated with Booker to solve the Diophantine equation for $k = 42$ . With this, all positive integers up to 100 that admit solutions have been found with at least one set of solutions.

Since this equation is not homogeneous, the doubling method from Question 2 is not applicable. The approach described in Andrew's paper involves reducing the time complexity of brute-force enumeration (from $O (N^{1 + o (1)})$ down to $O_{k} (N (\log \log N) (\log \log \log N))$ ), followed by a month of supercomputing (equivalent to 23 core-years of computation) to obtain a solution for 33, with a magnitude on the order of $N = 10^{16}$ . Therefore, it is impractical to write code for this problem; instead, a general explanation of the idea will suffice.

Algorithm with $O (N^{1 + o (1)})$ Complexity

First, consider a simplified version of the problem—the two-cube-sum problem. A relatively straightforward conclusion is that for a prime $p$ satisfying:

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

the integer solutions exist only when $p = 2$ or when $p$ takes the form $3 n^{2} - 3 n + 1$ . The proof is as follows:

Since $x^{3} + y^{3} = (x + y) (x^{2} - x y + y^{2})$ and $p$ is prime, we have:
$x + y = 1 \lor x^{2} - x y + y^{2} = 1.$
First, consider the right-hand condition: only $x = 1, y = 1$ satisfies it, corresponding to $p = 2$ .

Next, consider the left-hand condition: let $y = 1 - x$ , then the expression becomes $3 x^{2} - 3 x + 1$ (which must also be prime).

We now attempt to generalize $p$ to any integer $k$ . Let $k = r s$ , so we have $r = x + y$ and $s = x^{2} - x y + y^{2}$ . Let $y = r - x$ , which gives:

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

This is a standard quadratic equation in $x$ . The discriminant $Δ = (- 3 r)^{2} - 4 \cdot 3 \cdot (r^{2} - s) = 12 s - 3 r^{2}$ must be a perfect square for the solution to be an integer. Let $\sqrt{Δ} = t$ , then we obtain $x = (3 r + t) / 6$ and $y = (3 r - t) / 6$ . It suffices to verify whether $x$ and $y$ are integers.

Returning to the original three-cube-sum problem, we can move $z^{3}$ to the right-hand side:

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

We can then apply the algorithm for the two-cube-sum problem described above. By enumerating $z$ and performing prime factorization, we achieve a time complexity of $O (N^{1 + o (1)})$ .

Transformation to an Elliptic Curve

Assume $x^{3} + y^{3} + z^{3} = k > 0$ , with $| x | > | y | > | z | \geq \sqrt{k}$ . Let $d = | x + y |$ . Then $z$ is a cube root of $k$ modulo $d$ , and we can solve for:

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

Here, $sgn (x)$ denotes the sign function of $x$ , satisfying $sgn (0) = 0$ and $sgn (x) = \frac{x}{| x |}$ .

Then,

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

Here, the symbol $◻$ indicates that the computed value must be a perfect square. Therefore, this is equivalent to:

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

For a given search range $N$ , such as $N = 10^{16}$ in this problem, we need to enumerate $d \in [0, (\sqrt[3]{2} - 1) N]$ such that $d$ is coprime with $3$ . For each $d$ , we also need to enumerate $z \in [- N, N]$ satisfying $z^{3} \equiv k \mod d$ . Once we find a solution pair $(d, z)$ , it can be mapped back to the original solution $(x, y, z)$ .

For a specified integer $k$ , we can construct the corresponding elliptic curve $E_{d}$ :

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

such that the integer points on it correspond exactly to the integer solutions of the above equation. The specific mapping is:

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

Thus, the problem is transformed into finding integer points on an elliptic curve over $Q$ .

This algebraic transformation is effective for relatively small $d$ . For example, before the year 2000, the values of $k$ for which no integer solution $(x, y, z)$ had been found included:

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

Booker used the integer point finding functionality for elliptic curves in the algebraic software Magma and verified that for $d \leq 40$ , among all integer pairs $(k, d)$ , only ${(579, 29), (579, 34), (975, 22)}$ yield integer points on $E_{d}$ . In other words, except for $k = 579$ and $k = 975$ , none of the other listed $k$ values have integer solutions for $d \leq 40$ .

However, the parameter $d$ for the elliptic curve is on the order of $10^{16}$ , and there are $10^{16}$ such elliptic curves $E_{d}$ . Therefore, simply transforming the problem into finding integer points on an elliptic curve does not reduce the inherent complexity of the problem. Further optimization is still required.

Further Optimization

Returning to the previous section, in order to find pairs $(d, z)$ satisfying:

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

where $d \in [0, α N]$ ( $α$ is the same as before, i.e., $\sqrt[3]{2} - 1$ ), there are two time complexity bottlenecks:

To compute $z$ such that $z^{3} \equiv k \mod d$ , we need the prime factorization of $d$ .

The number of possible $(d, z)$ pairs is $Ω (N \log N)$ .

Apart from using elliptic curves to exclude a portion of integers $d$ , the paper primarily employs Montgomery's Batch Inversion Trick for batch computing cubic residues and a space-time trade-off sieve using the Legendre symbol to reduce the number of $(d, z)$ pairs that need to be traversed.

First, let's discuss batch inversion. Given $k$ numbers $a_{1}, a_{2}, \dots a_{k} \mod n$ , we need to compute $a_{1}^{- 1}, a_{2}^{- 1}, \dots, a_{k}^{- 1} \mod n$ simultaneously. If done separately using the extended Euclidean algorithm, the time complexity would be $O (k \log n)$ . Montgomery's approach is as follows:

Compute $p_{i} = a_{1} a_{2} \dots a_{i} \mod n$ . Calculating from $p_{1}$ to $p_{k}$ requires only $k - 1$ modular multiplications.
Use the extended Euclidean algorithm to compute $p_{k}^{- 1}$ , with time complexity $O (\log n)$ .
When $i = k$ , compute $a_{i}^{- 1} = p_{i - 1} p_{i}^{- 1}$ , then compute $p_{i - 1}^{- 1} = p_{i}^{- 1} a_{i}$ .
Let $p_{0} = 1$ , and perform the same operation for $i = k - 1, k - 2, \dots, 2$ .

Thus, the time complexity is reduced from $O (k \log n)$ to $O (k + \log n)$ .

Returning to bottleneck one, factorize $d$ to obtain $d = p_{1}^{e_{1}} p_{2}^{e_{2}} \dots p_{m}^{e_{m}}$ , then solve each cubic residue $z_{i}^{3} = k \mod p_{i}$ .

If $p = 2 \mod 3$ , there is only one root, so $z = k^{(2 p - 1) / 3} \mod p$ , which can be computed using fast exponentiation with time complexity $O (\log p)$ .
If $p = 1 \mod 3$ :
- Let $p = 3^{w} m + 1, m ∤ 3, b = k^{m} \mod p$ .
- Iteratively compute $b^{3}, b^{3^{2}}, \dots, b^{3^{v}}$ until $b^{3^{v}} = 1 \mod p$ . A cubic root exists if and only if $v < w$ .
- Arbitrarily choose $x$ such that $a = x^{m} \mod p$ has order $3^{w}$ , then compute $n = \log_{a} b$ . If $3 | n$ , then $a^{n / 3}$ is a cubic root of $b$ .
- Choose $e \in {1, 2}$ such that $3 | (p - e m)$ , then $z = a^{n e} k^{(p - e m) / 3}$ . The other two roots are $a^{3^{w - 1}} z$ and $a^{2 \times 3^{w - 1}} z$ .

Then, use Hensel's lemma to extend to $z_{i}^{3} = k \mod p_{i}^{e_{i}}$ , and use the Chinese Remainder Theorem (CRT) to extend to $z^{3} \mod d$ .

If the modular inverses are stored during batch inversion and $d$ is enumerated using "linear combinations" of the primes, the inverses do not need to be recalculated during CRT. This reduces the total time complexity for computing cubic roots for $d \leq α B$ to $O (N)$ .

Considering bottleneck two, let:

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

Select $P = 3 (\log \log N) (\log \log \log N)$ , and compute $M$ as the product of all primes in $[5, P]$ . If $Δ$ is a perfect square, then for any prime $p ∣ M$ (here the parentheses denote the Legendre symbol):

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

If the residue classes of $| z | \mod M$ are precomputed, and those residue classes for which the above equality holds for all primes $p ∣ M$ are selected, the Hasse bound can be used to prove that the number of residue classes does not exceed $N (\log \log N) (\log \log \log N)$ . Thus, the final complexity is reduced to $O_{k} (N (\log \log N) (\log \log \log N))$ .

In addition to the reduction in complexity, Andrew mentioned in subsequent work Sums of three cubes that for specific $k$ (e.g., $k = 33$ ), using the cubic reciprocity law as a sieve reduces the number of residue classes to $\frac{1}{63.6}$ of the original. This means only $5.5 \times 10^{9}$ values of $z$ need to be traversed, significantly reducing the workload.

In terms of implementation, the paper used Intel intrinsics and leveraged supercomputing parallelism (with $10^{6}$ subtasks) to enumerate $d$ . After a month of computation (equivalent to $23$ core-years), the paper found the first solution for $k = 33$ :