A latent bug in optimized Fp2 square-root algorithms

Jul 20, 2026 16 min read

TL;DR

Three widely-deployed $F_{p^2}$ square-root algorithms — the Adj–Rodríguez “complex method” (ePrint 2012/685, §V-A / Algorithm 8), Scott’s Tricks of the Trade §6.3 (ePrint 2020/1497), and the SQIsign-optimized variant of Aardal et al. (ePrint 2024/1563, Algorithm 3) — silently return $(0, 0)$ on inputs of shape $(a, 0)$ where $a \in F_p$ is a quadratic non-residue. The input is, in fact, always a square in $F_{p^2}$ for any nonzero $a \in F_p$ (Euler’s criterion: $a^{(q^2 - 1)/2} = (a^{q-1})^{(q+1)/2} = 1$), and the correct root is $(0, \sqrt{a/\beta})$ where $\beta = u^2$ is the QNR defining $F_{p^2}$. The algorithms just don’t find it.

The omission is at the paper level, not in any particular library: the proofs in §V-A of Adj–Rodríguez carry an unstated precondition $b \neq 0$. Scott’s exposition inherits it. Aardal et al.’s further optimization fuses away the explicit $\alpha = -1$ branch from a different, correct algorithm in the same Adj–Rodríguez paper (Algorithm 9), losing the protection.

Any library implementing any of the three algorithms verbatim inherits the bug.

1. Square roots in $F_p$

Let $p$ be a large odd prime and $F_p = \mathbb{Z}/p\mathbb{Z}$. For $a \in F_p^*$, computing $\sqrt{a}$ reduces to a single exponentiation in $F_p$ plus an $O(1)$ correction, with the exact form depending on $p \bmod 4$ (or $\bmod 8$, $\bmod 16$):

1.1. The easy case: $p \equiv 3 \pmod 4$

$$ y = a^{(p+1)/4} $$

Check:

$$ y^2 = a^{(p+1)/2} = a \cdot a^{(p-1)/2} = a \cdot \left(\tfrac{a}{p}\right). $$

If $a$ is a QR, $\left(\frac{a}{p}\right) = 1$ and $y^2 = a$. If $a$ is a QNR, $y^2 = -a$ (so the function “computes a sqrt of $-a$ instead”; the caller must check by squaring).

Cost: one $F_p$ exponentiation of bitlength roughly $\log_2 p$, i.e. about $1.5 \log_2 p$ $F_p$ multiplications via square-and-multiply, or fewer via addition chains.

1.2. The general case: Tonelli–Shanks

For $p \equiv 1 \pmod 4$, write $p - 1 = 2^s \cdot q$ with $q$ odd. Tonelli–Shanks finds $\sqrt{a}$ via a precomputed non-residue and $O(s^2)$ $F_p$ multiplications on top of one $F_p$ exponentiation. It’s still dominated by the exp, just with a $(\log_2 p)$-bit exponent and a small bookkeeping tail.

Variants: Atkin’s algorithm (for $p \equiv 5 \pmod 8$), Müller’s algorithm, and Cipolla–Lehmer (variable time). All amount to “one $F_p$ exp + bookkeeping”.

1.3. Cost summary

For any of these algorithms, the dominant cost is one $F_p$ exponentiation with an exponent the size of $p$. We’ll write this as a single unit $M_q$ (one exponentiation in the field of order $q$). Everything else is in the noise. Note the algorithm also tells us whether $a$ is a QR — either by inspecting the post-exp Legendre symbol or by squaring the result and comparing.

2. Square roots in $F_{p^2}$ — the naive approach

The same algorithms transfer directly: if $q = p^2$ is the order of $F_{p^2}$, then for $a \in F_{p^2}^*$, $\sqrt{a}$ is a single $F_{p^2}$ exponentiation away (with the same $\bmod 4 / \bmod 8 / \bmod 16$ casework on $q$, which depends on $p \bmod \cdot$).

Naive cost: one $F_{p^2}$ exponentiation $\approx \log_2(q^2) \cdot 3 = 6 \log_2 p$ $F_p$ multiplications (square-and-multiply with Karatsuba $F_{p^2}$ mul costing 3 $F_p$ muls). That’s a $\sim 3\times$ slowdown vs. $F_p$ Sqrt, because the exponent doubled in bitlength ($\log_2 q^2 = 2 \log_2 p$) and each step cost $\sim 3\times$ more (each $F_{p^2}$ mul $\approx 3$ $F_p$ muls).

So $F_{p^2}$ Sqrt by the naive method is feasible but costly — and pairing libraries do it a lot (G2 point decompression, hash-to-curve, twist arithmetic).

3. The descent: Adj–Rodríguez’s complex method (ePrint 2012/685, §V-A)

Adj and Rodríguez-Henríquez (2012) observed that $\sqrt{a}$ in $F_{p^2}$ can be reduced to two square roots in $F_p$ plus a handful of $F_p$ operations, by exploiting the structure of $F_{p^2} = F_p[u]/(u^2 - \beta)$ for a fixed QNR $\beta \in F_p$. The construction is the “complex method”, named because it generalizes the classical

$$ \sqrt{a + bi} = \sqrt{\frac{a + \sqrt{a^2 + b^2}}{2}} ;+; i \cdot \frac{b}{2 \sqrt{(a + \sqrt{a^2 + b^2})/2}} $$

formula from $\mathbb{C} = \mathbb{R}[i]/(i^2 + 1)$.

3.1. The setup

Let $x = a + b u \in F_{p^2}$ with $a, b \in F_p$ and $u^2 = \beta$ (the defining QNR). Suppose $y = c + d u$ is a square root of $x$, i.e. $y^2 = x$. Expanding:

$$ y^2 = (c + du)^2 = c^2 + 2cd,u + d^2 u^2 = (c^2 + \beta d^2) + (2cd),u. $$

Matching coefficients with $x$:

$$ c^2 + \beta d^2 = a, \qquad 2cd = b. \tag{$\star$} $$

3.2. The norm trick

Take the field norm of both sides of $y^2 = x$:

$$ N(y^2) = N(y)^2 = N(x), $$

where $N(z) = z \cdot z^p = z \cdot \bar z = (\Re z)^2 - \beta (\Im z)^2 \in F_p$ is the $F_p$-valued norm map. So if we compute

$$ s := \sqrt{N(x)} \quad \text{(a single $F_p$ Sqrt — one $F_p$ exp)} $$

then $s = \pm N(y) = \pm(c^2 - \beta d^2)$ (note the sign — minus $\beta$, not plus $\beta$; this is the norm of $y$, not of $y^2$).

3.3. Recovering $c$ and $d$

Now we have two equations involving $c$, $d$:

$$ \begin{aligned} c^2 + \beta d^2 &= a \quad \text{(from $\Re$ of $y^2 = x$)} \ c^2 - \beta d^2 &= \pm s \quad \text{(from $N(y) = \pm s$)} \end{aligned} $$

Adding: $2c^2 = a \pm s$, so $c^2 = (a \pm s)/2$. Subtracting: $2\beta d^2 = a \mp s$, so $d^2 = (a \mp s)/(2\beta)$.

Either of $c$ or $d$ can be recovered with one more $F_p$ Sqrt:

$$ c = \sqrt{(a + s)/2} \quad \text{(one $F_p$ Sqrt — the second $F_p$ exp)}, \qquad d = \frac{b}{2c} \quad \text{(one $F_p$ inv + one $F_p$ mul, from $2cd = b$).} $$

The sign of $s$ is resolved by checking which of $(a + s)/2$ or $(a - s)/2$ is a QR in $F_p$ — one of them is guaranteed to be, by a parity argument: their product is

$$ \frac{(a + s)(a - s)}{4} = \frac{a^2 - s^2}{4} = \frac{a^2 - N(x)}{4} = \frac{a^2 - (a^2 - \beta b^2)}{4} = \frac{\beta b^2}{4}, $$

and $\beta$ is a QNR so $\beta b^2 / 4$ is (QNR)$\cdot$(square) = QNR. Hence one of the two factors must be QR and the other QNR.

3.4. The algorithm (Adj–Rodríguez §V-A “Algorithm 8”)

Algorithm 8 (Adj-Rodríguez "Complex Method"):
Input:  x = a + b·u ∈ Fp²*
Output: y ∈ Fp² s.t. y² = x, or ⊥

  1: if b = 0:
  2:    return (√a, 0)                         ← the bug lives here
  3: end if
  4: α ← a² − β·b²                              (= N(x))
  5: γ ← χ(α)                                   (Legendre symbol of α in Fp)
  6: if γ = −1:
  7:    return ⊥                                (x is not a square in Fp²)
  8: end if
  9: s ← √α                                     (1st Fp Sqrt)
 10: δ ← (a + s) / 2
 11: if δ is QNR in Fp: δ ← (a − s) / 2
 12: c ← √δ                                     (2nd Fp Sqrt)
 13: d ← b / (2c)                               (1 Fp inv, 1 Fp mul)
 14: return (c, d)

Cost: $2 M_q + O(1)$ $F_p$ muls/inverses $\approx 2$ $F_p$ Sqrts. That’s a $\sim 3\times$ speedup vs. naive $F_{p^2}$ Sqrt (which costs $\sim 6$ $F_p$ exps’ worth).

4. The bug: lines 1–2

The purely-real shortcut “if $b = 0$: return $(\sqrt{a}, 0)$” is wrong when $a$ is a non-residue in $F_p$.

4.1. The math

Why $(a, 0)$ is always a square in $F_{p^2}$ for any nonzero $a \in F_p$:

$$ \begin{aligned} (a, 0)^{(q^2 - 1)/2} &= a^{(q^2 - 1)/2} \ &= a^{(q - 1)(q + 1)/2} \ &= \bigl(a^{q - 1}\bigr)^{(q + 1)/2} \ &= 1^{(q + 1)/2} \ &= 1 \qquad \text{(Fermat: $a^{q - 1} = 1$ in $F_p^*$).} \end{aligned} $$

So $(a, 0)$ must have a square root in $F_{p^2}$ — and we even know where it lives. From the equations $(\star)$ of §3.1 with $b = 0$:

$$ c^2 + \beta d^2 = a, \qquad 2cd = 0 ;\Longrightarrow; c = 0 \text{ or } d = 0. $$

  • $d = 0$ branch: $c^2 = a$. Works iff $a$ is a QR in $F_p$. Yields $(\sqrt{a}, 0)$.
  • $c = 0$ branch: $\beta d^2 = a$, so $d^2 = a/\beta$. Works iff $a/\beta$ is a QR in $F_p$. Yields $(0, \sqrt{a/\beta})$.

The two branches are complementary by Euler:

$$ \left(\tfrac{a/\beta}{p}\right) = \left(\tfrac{a}{p}\right) \cdot \left(\tfrac{1/\beta}{p}\right) = \left(\tfrac{a}{p}\right) \cdot \left(\tfrac{\beta}{p}\right) = \left(\tfrac{a}{p}\right) \cdot (-1) = -\left(\tfrac{a}{p}\right) $$

(since $\beta$ is QNR by construction). So exactly one of ${a, a/\beta}$ is a QR. If $a$ is QR, the $d = 0$ branch gives the real root. If $a$ is QNR, the $c = 0$ branch gives the purely imaginary root $(0, \sqrt{a/\beta})$. The output is always nonzero.

4.2. Why Algorithm 8 misses it

Line 2 unconditionally returns $(\sqrt{a}, 0)$. When $a$ is QNR, $\sqrt{a}$ doesn’t exist in $F_p$. What the $F_p$ Sqrt routine actually returns depends on the implementation:

  • Some return $\bot$ / None, in which case the call site propagates $\bot$, claiming $x$ is not a square — wrong, since it is.
  • Most pairing libraries’ $F_p$ Sqrt returns the exponentiation $a^{(p+1)/4}$ (for $p \equiv 3 \pmod 4$) without checking, yielding a value $r$ such that $r^2 = -a$ rather than $r^2 = a$. The caller sees $(r, 0)$ and may treat it as a sqrt, which it isn’t.
  • Constant-time variants using ccopy / fp_select cascades based on is_zero silently produce $(0, 0)$, also wrong.

The “proof” given on page 15 of Adj–Rodríguez — that one of $(a + s)/2$, $(a - s)/2$ is a QR because their product is the QNR $\beta b^2 / 4$ — fails when $b = 0$: the product collapses to $0$, which is neither QR nor QNR. The proof has an unstated precondition $b \neq 0$ that the algorithm’s “if $b = 0$” shortcut acknowledges in form but mishandles in substance.

4.3. The right way to handle $b = 0$

The correct branch dispatches on $\left(\frac{a}{p}\right)$:

if b = 0:
    if Legendre(a) = +1: return (√_Fp(a), 0)            ← purely-real root
    else:                return (0, √_Fp(a/β))          ← purely-imaginary root

Three $F_p$ operations more than the buggy version (one Legendre symbol, one $F_p$ inv, one $F_p$ Sqrt that runs on a known-QR input). Adding this branch is the fix that gnark-crypto kb8, SQIsign, relic, and (in this writeup’s companion PR) constantine apply.

Adj–Rodríguez had this branch in their own paper — in Algorithm 9, the novel $q \equiv 3 \pmod 4$ algorithm (see §7) — but their Algorithm 8 in §V-A doesn’t, and that’s the one downstream libraries copy.

5. Scott’s inverse-sqrt optimization (ePrint 2020/1497, §6.3)

Michael Scott’s 2020 Tricks of the Trade expands on Adj–Rodríguez’s complex method by fusing the final $F_p$ inversion $d = b/(2c)$ into the second $F_p$ Sqrt via the inverse-square-root trick. Concretely, replace lines 12-13 of Algorithm 8 with:

12': r ← invsqrt(δ)            (= 1/√δ; computed via a single Fp exp)
13': c ← r · δ                  (= √δ)
14': d ← b · r / 2              (= b/(2·√δ) = b/(2c))

The trick: $\mathrm{invsqrt}(\delta) = \delta^{(p-3)/4}$ for $p \equiv 3 \pmod 4$. Computing both $c$ and $d$ from $r$ saves the standalone inversion. Combined with progenitor tricks (Scott’s other contribution), the total cost drops from $\sim 2 M_q + 1$ inv to $\sim 1.5 M_q$.

Scott inherits the bug

Scott’s exposition reads:

If $a + i b$ is a quadratic residue, then one of $(a + \sqrt{a^2 - \beta b^2})/2$ or $(a - \sqrt{a^2 - \beta b^2})/2$ will be a quadratic residue in $F_p$. And indeed this must be the case since their product is equal to $\beta b^2/4$. This is a QNR times a perfect square, and hence a QNR. If a QNR is the product of two other field elements then one of them must be a QR, and the other a QNR. Therefore $\ldots$

Same reasoning, same gap: the “their product is $\beta b^2 / 4$, a QNR” argument requires $b \neq 0$. With $b = 0$ the product is $0$, neither QR nor QNR, and the QNR-times-square argument collapses. Scott §6.3 doesn’t discuss the case at all.

The resulting algorithm has the same shape as Algorithm 8 — same “if $b = 0$” shortcut, same lack of $\left(\frac{a}{p}\right)$-conditional handling. The Scott optimization speeds up the body of Algorithm 8 without revisiting the precondition Algorithm 8 was tacitly resting on.

Result: every library that ports Scott §6.3 has the bug. This includes arkworks-algebra (which patched it locally with an explicit if self.c1.is_zero() { ... } branch citing the Adj–Rodríguez paper’s Algorithm 8), relic’s $p \not\equiv 3 \pmod 4$ branch (which has an explicit purely-real handler), and historically gnark-crypto’s main internal/fptower/e2.go for all six pairing curves (which still has the bug at the time of writing).

6. Aardal et al.’s SQIsign optimization (ePrint 2024/1563, Algorithm 3)

Aardal, Adj, Alblooshi, Aranha, Canales-Martínez, Chávez-Saab, Filho, Reijnders and Rodríguez-Henríquez (2024) further specialized the Scott-optimized algorithm for the SQIsign post-quantum signature scheme, where $p \equiv 3 \pmod 4$. Their Algorithm 3 inlines the $q \equiv 3 \pmod 4$ casework, replaces the final inversion with a partial exponentiation, and reduces the entire computation to two $F_p$ exponentiations with no standalone inversions:

Algorithm 3 (Aardal et al., q ≡ 3 mod 4):
Input:  a = (a₀, a₁) ∈ Fp²*
Output: x s.t. x² = a, deterministic and constant-time

  1: δ ← (a₀² + a₁²)^((p+1)/4)                      (one Fp exp)
  2: x₀ ← a₀ + δ
  3: t₀ ← 2 · x₀
  4: x₁ ← t₀^((p−3)/4)                              (second Fp exp; inverse-sqrt of t₀)
  5: x₀ ← x₀ · x₁                                   (= √x₀ · scale)
  6: x₁ ← a₁ · x₁
  7: t₁ ← (2 · x₀)²
  8: if t₁ = t₀: return (x₀, x₁)
  9: else:        return (x₁, −x₀)

Steps 8-9 are a clever sign-correction trick: only one of the two pairs $(x_0, x_1)$ or $(x_1, -x_0)$ is the true sqrt, and the test $t_1 = t_0$ picks the right one without an extra exponentiation. With $\beta = -1$ (the SQIsign / BN254 / BLS12-381 setting), the algorithm runs entirely in $F_p$ with two $F_p$ exponentiations and a handful of $F_p$ muls. Beautifully tight.

6.1. Where the bug returns

Algorithm 3 has no “if $b = 0$” shortcut at all — Aardal et al. fused the case-split into the generic path. The intent is to make the algorithm branchless and constant-time. The effect is to remove the only check that might have stopped the bug at the front door.

Trace Algorithm 3 with $a = (a_0, 0)$, $a_0$ non-QR in $F_p$, $p \equiv 3 \pmod 4$:

Step 1. Compute $\delta$:

$$ \delta = (a_0^2 + 0)^{(p+1)/4} = a_0^{(p+1)/2} = a_0 \cdot a_0^{(p-1)/2} = a_0 \cdot \left(\tfrac{a_0}{p}\right) = a_0 \cdot (-1) = -a_0. $$

Step 2. Compute $x_0 = a_0 + \delta = a_0 + (-a_0) = 0$. This is where the bug fires.

Steps 3–7. Cascade through zero:

$$ t_0 = 0, \quad x_1 = 0^{(p-3)/4} = 0, \quad x_0 \leftarrow x_0 \cdot x_1 = 0, \quad x_1 \leftarrow a_1 \cdot x_1 = 0, \quad t_1 = (2 \cdot 0)^2 = 0. $$

Step 8. $t_1 = t_0$ (both zero), so return $(x_0, x_1) = (0, 0)$.

The output is $(0, 0)$, but the correct answer is $(0, \sqrt{a_0/\beta}) = (0, \sqrt{-a_0})$ (for $\beta = -1$), which exists in $F_p$ because $-1$ is QNR (so $-a_0 = (-1) \cdot \text{non-QR} = \text{QR}$) and is in the SQIsign codomain.

6.2. Why fusing broke things

Adj–Rodríguez 2012/685 actually has a bug-free algorithm for $q \equiv 3 \pmod 4$ — their Algorithm 9, presented in §V-B. It explicitly tests $\alpha = \left(\frac{x}{F_{p^2}}\right)$ (computed cheaply as $\alpha = a^{(q-1)/2}$) and handles the $\alpha = -1$ case via $x \leftarrow i \cdot x_0$ where $i = \sqrt{-1} \in F_{p^2}$ (which exists because $q \equiv 3 \pmod 4 \Rightarrow -1$ is QNR in $F_p \Rightarrow i \notin F_p$ but $i \in F_{p^2}$).

The $\alpha = -1$ branch is exactly the case where the input has its $F_{p^2}$-Legendre symbol equal to $-1$ only because the imaginary root is in $F_{p^2} \setminus F_p$ — i.e., the purely-real-non-QR-input case. Algorithm 9 catches it; Algorithm 3 fuses it away.

Aardal et al. mention that they’re building on Adj–Rodríguez §V (Algorithms 8, 9, 10), but the “fusion” pattern they pursue ends up moving from the structure of Algorithm 9 (which has the $\alpha = -1$ branch) back toward the structure of Algorithm 8 (which doesn’t). The optimization’s appeal is the absence of branches — but the absent branch is the one keeping correctness.

Result: every library that ports Aardal Algorithm 3 verbatim has the bug. The fix is to add a guard that detects the collapse (e.g., “if $a_1 = 0$ then $\delta \leftarrow a_0$” overrides the bad $\delta = -a_0$) or to switch back to Algorithm 9’s $\alpha = -1$ test. SQIsign’s own reference C implementation does the guard (fp_select(&x0, &x0, &(a->re), fp_is_zero(&(a->im))) at src/gf/ref/lvlx/fp2.c:162). Relic does the same (fp_copy_sec(t[0], a[0], fp_is_zero(a[1]))). gnark-crypto and constantine’s BLS12-377 path inherited the bug raw, and were patched only after independent rediscovery.

7. The fix: Adj–Rodríguez’s own Algorithm 9

The bug-free recipe has been in the same 2012 paper the whole time. Algorithm 9 (Adj–Rodríguez §V-B):

Algorithm 9 (Adj-Rodríguez, q ≡ 3 mod 4):
Input:  a ∈ Fp²*, i ∈ Fp² with i² = −1
Output: x s.t. x² = a, or ⊥

  1: a₁ ← a^((q − 3)/4)              (one Fp² exp)
  2: α  ← a₁ · (a₁ · a)               (= a^((q−1)/2) ∈ Fp²)
  3: a₀ ← α · α^q                     (= α^(q+1); Legendre check)
  4: if a₀ = −1: return ⊥             (a is not a square in Fp²)
  5: x₀ ← a₁ · a                      (= a^((q+1)/4))
  6: if α = −1:                       ★ KEY BRANCH
  7:    return i · x₀                  (purely-imaginary root)
  8: else:
  9:    b ← (1 + α)^((q − 1)/2)        (Tonelli-style correction)
 10:    return b · x₀
 11: end if

The $\alpha = -1$ branch is the structural fix the optimization layers fuse away. Adj–Rodríguez Algorithm 9 is the only one of the algorithms discussed here that handles the purely-real non-QR input correctly without an explicit “if $b = 0$” guard — because the guard is implicit in the $\alpha = -1$ test.

The Rust implementation zkcrypto/bls12_381’s Fp2::sqrt ports Algorithm 9 directly (citing it by name in the source comment: // Algorithm 9, https://eprint.iacr.org/2012/685.pdf, with constant time modifications.). It has been bug-free since the crate’s first release in $\sim 2020$. blst uses a different but equivalent rotation-alignment trick that also catches the case structurally. arkworks-algebra ports Algorithm 8 but adds the if c1.is_zero() branch with explicit Legendre dispatch (the §4.3 fix above).

For libraries that ship Algorithm 8 / Scott §6.3 / Aardal Algorithm 3 in their pristine published form, the bug is silent until a caller happens to ask for $\sqrt{(\text{non-QR}, 0)}$. That hardly ever happens in mainstream pairing protocols (curve constants block it), so the bug has lain dormant. But it’s directly reachable from any code that constructs primitive cube roots of unity in $F_{p^2}$ (via $(\sqrt{-3} - 1)/2$), Frobenius endomorphism constants, or post-quantum isogeny intermediates — that’s how it got rediscovered three times in the last three years.

8. Provenance summary

Adj-Rodríguez 2012/685 Algorithm 8  ──► buggy (purely-real shortcut)
                                       │
                                       ▼
                       Scott 2020/1497 §6.3  ──► buggy (inherits the §V-A proof gap)
                                       │
                                       ▼ (fused via Aardal optimization)
              Aardal 2024/1563 Algorithm 3  ──► buggy (re-introduces by absorbing α-detection)


Adj-Rodríguez 2012/685 Algorithm 9  ──► CORRECT (q ≡ 3 mod 4, has the α = −1 branch)
                                       │
                                       └─► zkcrypto/bls12_381 Fp2::sqrt  ──► clean (since ~2020)


Adj-Rodríguez 2012/685 Algorithm 10 ──► CORRECT (q ≡ 1 mod 4)

9. Impact

It turns out that libraries that implement the faster algorithms are the ones impacted:

10. References