Quadratic Residues and Tonelli-Shanks, in depth: square roots modulo p, Cipolla, Hensel lifting, CRT and point decompression

Residues and the Euler test

Work in the integers modulo an odd prime p. Squaring is two-to-one on the nonzero elements, because x and p - x have the same square and a polynomial of degree two has at most two roots in a field. So exactly (p - 1)/2 nonzero values are squares, called quadratic residues, and the other (p - 1)/2 are non-residues. For p = 113 the residues start 1, 2, 4, 7, 8, 9, 11, 13, 14, 15, 16, 18, while 3 and 5 are non-residues.

The test is the Euler criterion: for a not divisible by p, a^((p-1)/2) mod p is 1 when a is a residue and p - 1 (that is, -1) when it is not. The reason is that the multiplicative group modulo p is cyclic of order p - 1. Write a = g^e for a generator g; a is a square exactly when e is even, and a^((p-1)/2) = g^(e(p-1)/2) is 1 exactly when e is even. This is the Legendre symbol, treated fully in Legendre and Jacobi symbols in depth; here it is a gate costing one modular exponentiation. The gate matters: every algorithm below loops forever or returns garbage on a non-residue.

The easy primes

If p = 3 (mod 4), the answer is a single exponentiation: x = a^((p+1)/4) mod p. Check it: x^2 = a^((p+1)/2) = a * a^((p-1)/2) = a * 1 = a, using the Euler criterion in the last step. The NIST P-256 field prime is 3 mod 4, so its point decompression is one exponentiation.

There is a similar closed form when p = 5 (mod 8), usually credited to Atkin: compute v = (2a)^((p-5)/8), i = 2a*v^2, and x = a*v*(i - 1). The field prime of Curve25519, 2^255 - 19, is 5 mod 8.

Primes with p = 1 (mod 8) have no such formula, and they are common: NTT primes such as 998244353 = 119 * 2^23 + 1 are built with a large power of two in p - 1, and the NIST P-224 prime 2^224 - 2^96 + 1 has p - 1 = 2^96 times an odd number. Those are the cases Tonelli-Shanks and Cipolla exist for.

Tonelli-Shanks from first principles

Write p - 1 = Q * 2^S with Q odd. The whole difficulty lives in the subgroup of order 2^S: Tonelli-Shanks starts with a guess that is right up to an error in that subgroup and cancels the error one bit at a time. Set R = a^((Q+1)/2). Then R^2 = a * a^Q, so R would be the root if t = a^Q were 1. In general t is not 1, but it is a 2^S-th root of unity, so it lives in the 2-power subgroup. Find a non-residue z (try 2, 3, 4 and so on) and set c = z^Q. Because z is a non-residue, c has order exactly 2^S: it generates the whole 2-power subgroup. The algorithm keeps the invariant R^2 = a * t and repeatedly multiplies R by a power of c that strictly reduces the order of t, until t = 1 and R is the answer.

def tonelli_shanks(a, p):
    """Return x with x*x % p == a % p, or None if a is a non-residue. p prime."""
    a %= p
    if a == 0:
        return 0
    if p == 2:
        return a
    if pow(a, (p - 1) // 2, p) != 1:        # Euler criterion: no root exists
        return None
    if p % 4 == 3:                          # one-exponentiation shortcut
        return pow(a, (p + 1) // 4, p)
    Q, S = p - 1, 0
    while Q % 2 == 0:                       # p - 1 = Q * 2^S, Q odd
        Q //= 2
        S += 1
    z = 2
    while pow(z, (p - 1) // 2, p) != p - 1: # any non-residue
        z += 1
    M, c, t, R = S, pow(z, Q, p), pow(a, Q, p), pow(a, (Q + 1) // 2, p)
    while t != 1:
        i, t2 = 0, t                        # least i with t^(2^i) == 1
        while t2 != 1:
            t2 = t2 * t2 % p
            i += 1
        b = pow(c, 1 << (M - i - 1), p)     # b = c^(2^(M-i-1))
        M, c, t, R = i, b * b % p, t * b * b % p, R * b % p
    return R

Why it terminates: at each step t has order exactly 2^i with i below M, and c has order 2^M. Then b = c^(2^(M-i-1)) has order 2^(i+1), so b^2 has order 2^i, and multiplying t by b^2 cancels its top bit of order: the new t has order at most 2^(i-1). Meanwhile R is multiplied by b, so R^2 = a * t is preserved. Since i drops every round, there are at most S rounds, each doing at most S squarings, giving O(S^2) multiplications on top of the exponentiations; negligible for small S, dominant when S is 96 as for P-224. Finding z takes two Euler tests on average; for a fixed field, compute c once and cache it.

Worked trace: the square root of 13 modulo 113

Take p = 113 and a = 13. First the gate: 13^56 mod 113 = 1, so 13 is a residue. Next, 113 is 1 mod 4, so the shortcut does not apply. Factor p - 1 = 112 = 7 * 2^4, so Q = 7 and S = 4. The first non-residue is 3, so z = 3.

t is now 1, so R = 76. Check: 76^2 = 5776 = 51 * 113 + 13. The other root is 113 - 76 = 37. The first round removed the order-8 part of the error, leaving t = -1, and the second removed that. Cipolla on the same input also returns 76, a useful cross-check when porting either routine.

Cipolla algorithm

Cipolla takes a different route. Find r such that d = r^2 - a is a non-residue, then work in the two-dimensional field of numbers u + v*w with w^2 = d. In that field, (r + w)^((p+1)/2) is a square root of a, and it lands back in the ordinary integers modulo p. The proof uses the Frobenius map: (r + w)^p = r - w, so (r + w)^(p+1) = (r + w)(r - w) = r^2 - d = a.

def cipolla(a, p):
    """Square root of a mod prime p via arithmetic in F_p[w]/(w^2 - d)."""
    a %= p
    if a == 0 or p == 2:
        return a
    if pow(a, (p - 1) // 2, p) != 1:
        return None
    r = 0
    while pow((r * r - a) % p, (p - 1) // 2, p) != p - 1:
        r += 1                               # need r*r - a to be a non-residue
    d = (r * r - a) % p
    def mul(x, y):                           # (x0 + x1 w)(y0 + y1 w), w^2 = d
        return ((x[0] * y[0] + x[1] * y[1] * d) % p,
                (x[0] * y[1] + x[1] * y[0]) % p)
    res, base, e = (1, 0), (r, 1), (p + 1) // 2
    while e:
        if e & 1:
            res = mul(res, base)
        base = mul(base, base)
        e >>= 1
    return res[0]                            # the w-coordinate is provably 0

The cost is one exponentiation in the extension field, a few base-field multiplications per bit, regardless of S. So Cipolla wins when S is large, and Tonelli-Shanks wins for small S. Cipolla is a safe default when p arrives at run time and could have any shape.

Prime powers and composite moduli

For a prime power p^k with p odd and a coprime to p, a root modulo p lifts uniquely to a root modulo p^k. The update is a Newton step on f(x) = x^2 - a: x becomes x - f(x)/f'(x), and f'(x) = 2x is invertible at every level, via the extended Euclidean algorithm.

def sqrt_mod_prime_power(a, p, k):
    """Root of a modulo p^k for odd p with gcd(a, p) == 1, via Newton/Hensel."""
    r = tonelli_shanks(a, p)
    if r is None:
        return None
    m = p
    for _ in range(k - 1):                   # linear lifting, one power at a time
        m *= p
        r = (r - (r * r - a) * pow(2 * r, -1, m)) % m
    return r

Worked example: the square root of 2 modulo 7 is 4 (16 = 2 + 14). Lifting once gives 39 modulo 49, since 39^2 = 1521 = 31 * 49 + 2. Lifting again gives 235 modulo 343.

For a general modulus n = p1^k1 * ... * pm^km, find a root modulo each prime power and combine them with the Chinese remainder theorem, as in the Chinese remainder theorem article. Each odd prime-power factor contributes two roots, so there are 2^m roots modulo n. That multiplicity is why square roots modulo n = pq are as hard as factoring: two roots that are not negatives of each other reveal a factor through gcd(x - y, n), the engine of the Rabin cryptosystem proof and the quadratic sieve. Powers of two and a divisible by p need separate handling; refuse them loudly rather than return a wrong answer.

Point decompression in elliptic-curve code

A compressed elliptic-curve point stores x and one bit of y. The decoder computes rhs = x^3 + a*x + b modulo the field prime, takes a square root, and picks the root whose parity matches the stored bit. Three things go wrong:

• Skipping the residue check. If the right-hand side is a non-residue, the x coordinate is not on the curve. A decoder that runs the p = 3 mod 4 formula anyway returns a value whose square is -rhs, producing an invalid point; invalid-point inputs are the basis of real key-recovery attacks. Always square the result and compare.
• Timing leaks. The Tonelli-Shanks loop runs a data-dependent number of rounds and squarings. When the input is secret, use a constant-time variant with a fixed number of iterations and conditional moves, or a field whose shape allows a fixed exponentiation.
• The zero root. When rhs is 0 there is one root; reject an encoding that asks for an odd zero.

Failure modes and testing

The most valuable habit is an oracle test: for every prime below a few thousand, compare against a brute-force table of squares for every a, including 0 and the non-residues. For the cyclic-group structure behind all of this, see primitive roots and discrete logarithms.

Trade-offs

What to do next

1. Implement the Euler test and a brute-force oracle, then make every later routine pass the oracle for all primes below 1,000.
2. Type in the Tonelli-Shanks listing and reproduce the p = 113, a = 13 trace by printing M, c, t and R each round.
3. Implement Cipolla and time both on 998244353 and on the P-224 prime; observe the S dependence.
4. Add Hensel lifting and CRT, then enumerate all four square roots of 4 modulo 15 and check them by hand.
5. Write a point decompressor for a curve you use and add an explicit on-curve check.
6. Read the Legendre and Jacobi article to replace the Euler test with the faster Jacobi algorithm where it matters.