FLINT3.7.0-dev
/
Algorithm

Square roots modulo a prime

Cost
O(\log n) modular multiplications if n \not\equiv 1 \pmod 8. Otherwise, with n - 1 = q 2^s and q odd, O(\log n + s^2) multiplications for Tonelli–Shanks, or O(\log n) multiplications and a few Jacobi symbols for Müller's method.

Given a prime n and a \in \mathbb{Z}/n\mathbb{Z}, we want r with r^2 = a. Half of the nonzero residues are squares. For a a square, the two roots are r and n - r, and nmod_sqrt returns the one in [0, n/2]. The choice of method depends on n \bmod 8 and on the 2-adic valuation s of n - 1. Standard references are [Coh1996] and [CraPom2005].

Overview#

nmod_sqrt(res, a, mod) proceeds as follows.

  1. If a \le 1, return a. If n is even, return failure.
  2. If n < 128 (NMOD_SQRT_SMALL_CUTOFF), search t = 1, 2, \ldots, \lfloor n/2 \rfloor with t^2 updated incrementally (_nmod_sqrt_small). This is exact for any odd n.
  3. Otherwise choose the arithmetic: half-word lazy Montgomery arithmetic if n < 2^{\mathtt{FLINT\_BITS}/2 - 2}, full-word lazy Montgomery arithmetic if n < 2^{\mathtt{FLINT\_BITS} - 2}, and plain Montgomery arithmetic otherwise. All remaining steps work in Montgomery form, and the three versions are compiled separately so that each is fully inlined.

  4. Compute a candidate root r:

  5. Return r if r^2 = a, and failure otherwise.

No Jacobi symbol of a is computed in advance. For prime n, every method below returns a correct root when a is a square, so the final check is a complete test of whether a is a square. For composite n, the check guarantees that a returned root is genuine, but a failure proves nothing.

The branches in the inner loops are replaced by masks and selections, since their conditions depend on random data and would be mispredicted. For the same reason, exponentiation multiplies by 1 instead of skipping a multiplication when an exponent bit is 0.

The case n \equiv 3 \pmod 4#

Algorithm 1

Input: prime n \equiv 3 \pmod 4 and a.

Output: a candidate root r.

  1. Return r = a^{(n+1)/4}.

If a is a square, Euler's criterion gives a^{(n-1)/2} = 1, so r^2 = a^{(n+1)/2} = a \cdot a^{(n-1)/2} = a.

The case n \equiv 5 \pmod 8#

Here -1 is a square and 2 is not. Atkin's method [Atk1992] uses a square root of -1 built from 2a.

Algorithm 2: Atkin

Input: prime n \equiv 5 \pmod 8 and a \ne 0.

Output: a candidate root r.

  1. t \gets 2a and b \gets t^{(n-5)/8}.
  2. i \gets t b^2.
  3. Return r = ab(i - 1).

Theorem 1. If a is a nonzero square modulo the prime n \equiv 5 \pmod 8, then Algorithm 2 returns a square root of a.

Proof. We have i = t^{(n-1)/4}, so i^2 = (2a)^{(n-1)/2} = -1, because 2 is a nonsquare and a is a square. Hence (i - 1)^2 = i^2 - 2i + 1 = -2i and

r^2 = a^2 b^2 (-2i) = -a i \cdot (2a b^2) = -a i^2 = a.

Tonelli–Shanks#

For n \equiv 1 \pmod 8 write n - 1 = q 2^s with q odd, so s \ge 3. The method of Tonelli [Ton1891] and Shanks [Sha1973] works in the subgroup of order 2^s of (\mathbb{Z}/n\mathbb{Z})^\times, which is cyclic.

Nonresidue. It needs some k that is not a square. FLINT takes the least odd prime k \le 61 with Jacobi symbol (k \mid n) = -1, and otherwise the least odd k \ge 67 with (k \mid n) = -1. Only odd primes need to be tried among the small values, because the least nonresidue is prime and 2 is a square for n \equiv 1 \pmod 8. By quadratic reciprocity (k \mid n) = (n \mid k), and (n \mid k) is read from a bitmask indexed by n \bmod k, computed with two multiplications instead of a division: with m = n k^{-1} \bmod B, the high word h of mk satisfies mk = n + hB, so n \equiv -hB \pmod k, and since B is a square, (n \mid k) = (-h \mid k). The search gives up, which makes nmod_sqrt fail, if n is found to be a perfect square or no nonresidue is found below 2^{20}; neither happens for prime n.

Setup. With u = a^{(q-1)/2} and c = k^{(q-1)/2}, both computed by one exponentiation loop that interleaves the two independent squaring chains, set

r = ua = a^{(q+1)/2}, \qquad t = ur = a^q, \qquad g = c^2 k = k^q.

Then r^2 = a t, and g has order exactly 2^s since g^{2^{s-1}} = k^{(n-1)/2} = -1. If a is a square, then t^{2^{s-1}} = a^{(n-1)/2} = 1. It remains to find x with t g^x = 1; then x is even when a is a square, and r g^{x/2} is a root.

Algorithm 3: Tonelli–Shanks descent

Input: r, t, g as above, with r^2 = at and g of order 2^s.

Output: a root r of a, or failure.

  1. m \gets s.

  2. While t \ne 1:

    1. Find the least i with 0 < i < m and t^{2^i} = 1, by repeated squaring. If there is none, return failure.
    2. b \gets g^{2^{m - i - 1}}, g \gets b^2, t \gets t g, r \gets r b, m \gets i.
  3. Return r.

Theorem 2. If a is a nonzero square modulo the prime n, then Algorithm 3 terminates and returns a square root of a. If a is not a square, it returns failure.

Proof. The loop maintains r^2 = at, g of order exactly 2^m, and t^{2^{m-1}} = 1. This holds initially. In step 2.1, t \ne 1 and t^{2^{m-1}} = 1, so i exists and t has order exactly 2^i. Then b has order 2^{i+1}, so b^2 has order 2^i. In a cyclic group, the elements of order exactly 2^i all satisfy x^{2^{i-1}} = -1, so (t b^2)^{2^{i-1}} = 1. The updates give (rb)^2 = at b^2 = a t_{\text{new}}, and the invariant holds with m = i. Since m decreases strictly, the loop terminates, and on exit t = 1 and r^2 = a. If a is not a square, then t^{2^{s-1}} = -1, so no i < s exists in the first round.

Each round costs up to m dependent squarings, so the total is O(s^2) multiplications in one long dependency chain. This is cheap for small s, which is the common case: asymptotically, half of the primes n \equiv 1 \pmod 8 have s = 3. FLINT uses it for s < 12.

Bitwise variant#

For larger s the latency of the dependent squarings dominates. The second form finds the bits x_0, x_1, \ldots of x from the bottom, as in Pohlig–Hellman (see Discrete logarithms modulo a prime (Pohlig–Hellman)), using products that do not depend on each other.

Algorithm 4: Tonelli–Shanks by bits

Input: r, t, g as above.

Output: a root r of a, or failure.

  1. For 0 \le i < s: D_i \gets t^{2^i} and G_i \gets g^{2^i}.
  2. If D_{s-1} \ne 1, return failure.
  3. For j = 1, \ldots, s - 1: if D_{s-1-j} \ne 1, set D_i \gets D_i G_{i+j} for 0 \le i < s - 1 - j, and r \gets r G_{j-1}.
  4. Return r.

Theorem 3. If a is a nonzero square modulo the prime n, then Algorithm 4 returns a square root of a. If a is not a square, it returns failure.

Proof. Let x \in [0, 2^s) with t g^x = 1, which exists because t lies in the cyclic group generated by g. Let X_j be the part of x below bit j. We claim that before iteration j of step 3, D_i = (t g^{X_j})^{2^i} for i \le s - 1 - j, and r = a^{(q+1)/2} g^{X_j/2}.

Since t g^{X_j} = g^{-(x - X_j)} and 2^j divides x - X_j, we get D_{s-1-j} = g^{-2^{s-1} (x - X_j)/2^j} = (-1)^{x_j}, so the test in step 3 reads off bit x_j. Setting the bit multiplies t g^{X_j} by g^{2^j}, hence each D_i by G_{i+j}, and r by g^{2^{j-1}} = G_{j-1}. Only D_i with i < s - 1 - j are used later. Step 2 is the case j = 0: D_{s-1} = a^{(n-1)/2}, which is 1 if and only if a is a square, and then x_0 = 0. At the end X = x and r^2 = a^{q+1} g^{x} = a \cdot t g^x = a.

Step 1 is two independent chains of s - 1 squarings, and step 3 performs at most s(s-1)/2 multiplications. The work is still O(s^2), but the multiplications of each iteration are independent, so their cost is bounded by multiplication throughput instead of latency.

Müller's method#

The cost of Tonelli–Shanks grows with s. The Cipolla–Lehmer method works in \mathbb{F}_{n^2} instead, with a cost that does not depend on s. Müller [Mue2004] gave a version for n \equiv 1 \pmod 4 that needs a single Lucas sequence. For P \in \mathbb{Z}/n\mathbb{Z}, let V_0 = 2, V_1 = P and V_{k+1} = P V_k - V_{k-1}, so that V_k = \alpha^k + \alpha^{-k} where \alpha, \alpha^{-1} are the roots of X^2 - PX + 1.

Theorem 4 (Müller). Let n \equiv 1 \pmod 4 be prime, let a be a nonzero square, and let t \ne 0 with (a t^2 - 4 \mid n) = -1. With P = a t^2 - 2,

\left(\frac{V_{(n-1)/4}(P)}{t}\right)^2 = a.

Proof. Let u = t \sqrt a \in \mathbb{Z}/n\mathbb{Z}, and let w be a square root of u^2 - 4 = a t^2 - 4 in \mathbb{F}_{n^2}; it is not in \mathbb{Z}/n\mathbb{Z}. Put \beta = (u + w)/2 and \beta' = (u - w)/2. Then \beta \beta' = 1, \beta + \beta' = u and \beta^2 + \beta'^2 = u^2 - 2 = P, so \alpha = \beta^2 is a root of X^2 - PX + 1. The Frobenius map fixes u and sends w to -w, so \beta^n = \beta' = \beta^{-1} and \beta^{n+1} = 1. Hence \varepsilon = \beta^{(n+1)/2} satisfies \varepsilon^2 = 1 and \varepsilon = \pm 1. Now

V_{(n-1)/4} = \beta^{(n-1)/2} + \beta^{-(n-1)/2} = \varepsilon \beta^{-1} + \varepsilon \beta = \varepsilon u = \pm t \sqrt a.
Algorithm 5: Müller

Input: prime n \equiv 1 \pmod 4 and a.

Output: a candidate root r, or no answer.

  1. For j = 0, 1, \ldots, 32, let t = 2^j and test whether (a t^2 - 4 \mid n) = -1; stop at the first success. If none succeeds, give no answer.
  2. P \gets a t^2 - 2.
  3. Compute V = V_{(n-1)/4}(P) with the ladder (V_m, V_{m+1}) \mapsto (V_{2m}, V_{2m+1}) or (V_{2m+1}, V_{2m+2}) over the bits of (n-1)/4, using V_{2m} = V_m^2 - 2 and V_{2m+1} = V_m V_{m+1} - P.
  4. Return r = V / t.

Choosing t = 2^j makes the division by t in step 4 a sequence of halvings, each a shift and a conditional addition of (n+1)/2. The Jacobi symbols in step 1 are computed for a - 4 (j = 0) and for c_j = 4^{j-1} a - 1 (j \ge 1), because a 4^j - 4 = 4 c_j and 4 is a square; they satisfy c_1 = a - 1 and c_{j+1} = 4 c_j + 3. Each candidate fails with probability about 1/2, so giving no answer after 33 candidates is very unlikely; nmod_sqrt then falls back to Tonelli–Shanks.

Each step of the ladder costs one squaring and one independent multiplication, so the latency is about that of one exponentiation. The fixed cost of the Jacobi symbols makes the method pay off only for large s. FLINT uses it when 6s \ge b + 120, where b is the bit length of n. Since s < b, this happens only for n \ge 2^{25}, and for 64-bit n it requires s \ge 31.

Use in FLINT#

int nmod_sqrt(ulong * res, ulong a, nmod_t mod)
int _nmod_sqrt_small(ulong * res, ulong a, ulong n)
ulong n_sqrtmod(ulong a, ulong p)

Referenced from#

References#

  • [Coh1996] Cohen, Henri : A course in computational algebraic number theory, Springer, 1996
  • [CraPom2005] R. Crandall and C. Pomerance, Prime Numbers: A Computational Perspective, second edition, Springer (2005).
  • [Atk1992] A. O. L. Atkin, "Probabilistic primality testing", summary by F. Morain, in Algorithms Seminar 1991–1992, INRIA Research Report 1779 (1992), 159–163
  • [Ton1891] A. Tonelli, "Bemerkung über die Auflösung quadratischer Congruenzen", Nachrichten von der Königlichen Gesellschaft der Wissenschaften und der Georg-Augusts-Universität zu Göttingen (1891), 344–346
  • [Sha1973] D. Shanks, "Five number-theoretic algorithms", Proceedings of the Second Manitoba Conference on Numerical Mathematics, Congressus Numerantium 7 (1973), 51–70
  • [Mue2004] S. Müller, "On the computation of square roots in finite fields", Designs, Codes and Cryptography 31 (2004), 301–312, https://doi.org/10.1023/B:DESI.0000015890.44831.e2
Cite this page

The FLINT team, “Square roots modulo a prime,” FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev, 2026. https://flintlib.org/doc/algorithms/number-theory/sqrt-mod-p.html (generated 2026-10-20).