Tonelli-Shanks Algorithm, in depth: the cost model in S, one shared exponentiation, the constant-time RFC 9380 loop and table-driven roots for high-2-adicity fields

S is the number that matters

Write p - 1 = Q * 2^S with Q odd. S is the 2-adicity of the field, and it decides almost everything about the cost. If S = 1 (p = 3 mod 4) the root is one exponentiation, x^((p+1)/4). If S = 2 (p = 5 mod 8) a fixed formula needs one exponentiation and a correction by sqrt(-1). Above that you need a loop that works inside the subgroup of order 2^S, and its length grows with S. The fields people actually use span the whole range:

The large values are not accidents. Proof systems choose fields with big S because a large power-of-two subgroup gives the roots of unity that a radix-2 number-theoretic transform needs. The fields where fast polynomial arithmetic is easiest are therefore the fields where square roots cost most.

Share the exponentiation

Textbook statements compute two exponentiations: R = x^((Q+1)/2) and t = x^Q. Both are derived from one. Let w = x^((Q-1)/2); then R = w * x and t = w * w * x, which costs two extra multiplications instead of a second exponentiation. The loop then keeps the invariant R^2 = x * t and multiplies t by powers of g = z^Q, where z is any non-residue, until t = 1.

def sqrt_ts(x, p, Q, S, g):
    """x a nonzero square mod p; p - 1 = Q * 2^S; g = z^Q for a non-residue z (precomputed)."""
    w = pow(x, (Q - 1) // 2, p)
    R, t = w * x % p, w * w % p * x % p
    M, c = S, g
    while t != 1:
        i, t2 = 0, t
        while t2 != 1:                     # least i with t^(2^i) = 1
            t2 = t2 * t2 % p
            i += 1
        b = c
        for _ in range(M - i - 1):
            b = b * b % p
        M, c = i, b * b % p
        t, R = t * c % p, R * b % p
    return R

On secp256k1, where the loop never runs, sharing halves the work: 1,005 multiplications against 505.

The cost model, measured

The loop has two parts. Finding i takes up to M squarings, and raising c to 2^(M-i-1) takes the rest of the round. The number of rounds and their lengths depend on the order of t, which depends on x, so the cost is data-dependent. Mean multiplications over 200 random squares:

Two facts stand out. For S = 32 on a 255-bit field the loop costs roughly as much as the exponentiation. For S = 96 the loop is most of the work, and the constant-time version, which always runs the worst case, costs almost twice the classic mean.

The constant-time loop from RFC 9380

If x is secret, the classic loop leaks it: the number of rounds and squarings depends on the order of x^Q, and timing reveals it. RFC 9380, the hash-to-curve standard, says constant-time implementations are required when encoding inputs are secret, and its Appendix I.4 gives a constant-time Tonelli-Shanks credited to Bowe, Grigg and Ogilvie-Wigley and optimised by Michael Scott. It always runs S - 1 rounds with a fixed number of squarings each, and replaces the branch with CMOV(a, b, e), which returns a when e is false and b when e is true. Transcribed:

def sqrt_ts_ct(x, p, S, Q, c5):
    """RFC 9380 Appendix I.4. c5 = z^Q for a non-square z. Shows the operation
    sequence only: Python integers and the if-expression are NOT constant-time."""
    c3 = (Q - 1) // 2
    z = pow(x, c3, p)
    t = z * z % p * x % p
    z = z * x % p
    b, c = t, c5
    for i in range(S, 1, -1):              # i = S, S-1, ..., 2
        for _ in range(i - 2):
            b = b * b % p
        e = (b == 1)
        zt = z * c % p
        z = z if e else zt                 # z = CMOV(zt, z, e)
        c = c * c % p
        tt = t * c % p
        t = t if e else tt                 # t = CMOV(tt, t, e)
        b = t
    return z

The squaring count is fixed at the sum of (i - 2) for i from 2 to S, which is (S - 1)(S - 2)/2, plus three multiplications per round. For P-224 that is 4,465 + 285 = 4,750, and adding the exponentiation gives the 5,007 measured on every single input. In a real library, CMOV is a masked select over limbs and the comparison with one has no early exit.

Worked example: the constant-time loop modulo 97

Take p = 97, so p - 1 = 96 = 3 * 2^5: Q = 3, S = 5. The smallest non-residue is z = 5, so c5 = 5^3 mod 97 = 28. Find the root of x = 2, which is a square because 2^48 mod 97 = 1. Setup gives z = x^1 * x = 4 and t = 2 * 2 * 2 = 8. The value 8 has order 16, the largest a square can give t when S = 5, so every round corrects something.

The result is 83, and 83^2 = 6889 = 71 * 97 + 2. The other root is 97 - 83 = 14. Notice that the loop shape would be identical for an input whose t was already 1: the same squarings and multiplications run, and only the CMOV selections differ.

Tables for a fixed field

For a fixed field with large S, and inputs that are public, there is a faster route. Since t = x^Q lies in the subgroup generated by g, t = g^e for an even e, and R * g^(-e/2) is the root. So the loop is really a discrete logarithm in a group of order 2^S, and that can be solved w bits at a time with precomputed tables, the same Pohlig-Hellman idea described in the primitive roots and discrete logarithms article. Recover the low w bits of e by squaring t up into the subgroup of order 2^w and looking the result up; strip those bits with a table entry; repeat.

class WindowSqrt:
    def __init__(self, p, w):
        self.p, self.w = p, w
        self.Q, self.S = split(p)                       # p - 1 = Q * 2^S
        g = pow(nonresidue(p), self.Q, p)
        self.ks, self.look, self.inv, self.half = range(0, self.S, w), {}, {}, {}
        for k in self.ks:
            ww = min(w, self.S - k)
            base = pow(g, 2 ** (self.S - ww), p)
            for d in range(2 ** ww):
                self.look[(ww, pow(base, d, p))] = d
                self.inv[(k, d)] = pow(g, -(d << k) % (p - 1), p)
                self.half[(k, d)] = pow(g, -((d << k) // 2) % (p - 1), p)

    def sqrt(self, x):
        p, S, w = self.p, self.S, self.w
        x0 = pow(x, (self.Q - 1) // 2, p)
        R, t = x0 * x % p, x0 * x0 % p * x % p
        for k in self.ks:
            ww = min(w, S - k)
            y = t
            for _ in range(S - k - ww):
                y = y * y % p
            d = self.look[(ww, y)]                      # next ww bits of e
            if d:
                t = t * self.inv[(k, d)] % p
                R = R * self.half[(k, d)] % p
        return R

The lowest window always yields an even digit, because e is even for a square, so the halving table is exact. Costs for P-224: w = 4 needs 784 table entries and 1,406 multiplications on average; w = 8 needs 6,400 entries, about 179 KB at 28 bytes each, and 809 multiplications, 3.5 times fewer than the classic version. Goldilocks with w = 8 drops from 406 to 121. As written, the lookups and the branch on d depend on the input, so this is for public data, such as decompressing public keys or verifying, and not for secret inputs; a constant-time version would scan whole tables with masked selects and lose much of the gain.

Failure modes

• Non-residue input. The classic loop is only correct for squares. Given a non-residue, t has order 2^S, i reaches M and the exponent M - i - 1 goes negative. The constant-time loop quietly returns a wrong value. Always check r * r = x at the end; RFC 9380's sqrt_ratio returns an is-square flag with its result.
• Zero. With x = 0 the classic inner loop squares 0 forever. Handle 0 first.
• A residue mistaken for z. If z is a square, g does not generate the 2-power subgroup and the loop fails for about half of all inputs. Derive z with an Euler test, never by assuming 2 or 3 works; for P-224 the first non-residue is 11.
• Constants from another field. Libraries hardcode S, Q, (Q - 1)/2 and z^Q per field. Copying them between the BLS12-381 base and scalar fields compiles and then fails.
• Unnormalised output. Which of r and p - r a variant returns is arbitrary. Protocols that need a canonical root define one, such as the sgn0 rule in RFC 9380; apply it after the square root, never inside it.
• Representation slips. In Montgomery form the test b = 1 must compare against the Montgomery representation of one. A test against the integer 1 fails every round.

Testing and trade-offs

Run every variant exhaustively over all squares for all primes below 2,000, which catches off-by-one errors in loop bounds; the listings above pass that test. Then run thousands of random squares in each production field and require all variants to agree up to sign. Feed known non-residues and require rejection. Finally, recompute the constants from p in the test suite, by factoring out powers of 2 and searching for z, and compare them with the hardcoded ones. For the constant-time variant, add a timing harness that compares inputs whose t has order 1 with inputs whose t has order 2^(S-1), and read the side-channel article for what such tests can and cannot prove.

What to do next

1. Compute S for every field your code uses, with the factoring loop above, and record it next to the field definition.
2. Make sure your exponentiation is shared: one power, then R = w * x and t = w * w * x. Count multiplications with a wrapper class before and after.
3. If any input to the root can be secret, use the constant-time loop with a real CMOV and confirm with a timing test that the loop length does not vary.
4. For a public-input hot path in a field with S of 28 or more, prototype the windowed table and measure multiplications against memory for w = 4 and w = 8.
5. Add the exhaustive small-prime test, the non-residue rejection test and the recomputed-constants test to CI. Then read Legendre and Jacobi symbols for a faster residuosity test and modular exponentiation for the addition chains that shrink the setup cost.