Tonelli-Shanks Algorithm

Purpose: Compute a modular square root — find $x$ with $x^2\equiv a(\mathrm{mod}p)$ for an odd prime $p$ — in $O({\log }^{2}p)$ expected time.

Algorithm

1. Check solvability with the Euler criterion: $a^{(p-1)/2}\equiv 1(\mathrm{mod}p)$, otherwise no root exists.
2. Easy case $p\equiv 3(\mathrm{mod}4)$: the answer is $x=a^{(p+1)/4}$. Done in one modpow.
3. General case. Write $p-1=q\cdot 2^s$ with $q$ odd.
4. Find any quadratic non-residue $z$ (try $z=2,3,\dots$; half of all residues work, so this takes $O(1)$ tries on average). Set $c=z^q$.
5. Initialise $x=a^{(q+1)/2}$, $t=a^q$, $m=s$.
6. While $t\ne 1$:
   • find the least $i\in [1,m)$ with $t^{2^i}=1$;
   • set $b=c^{2^{m-i-1}}$, then $x\leftarrow xb$, $t\leftarrow t{b}^2$, $c\leftarrow b^2$, $m\leftarrow i$.
7. Return $x$ (and $p-x$, the other root).

Code

long long tonelli(long long a, long long p) {
    a %= p; if (a < 0) a += p;
    if (a == 0) return 0;
    if (powmod(a, (p - 1) / 2, p) != 1) return -1;      // not a residue
    if (p % 4 == 3) return powmod(a, (p + 1) / 4, p);
 
    long long q = p - 1; int s = 0;
    while (q % 2 == 0) { q /= 2; s++; }
 
    long long z = 2;
    while (powmod(z, (p - 1) / 2, p) != p - 1) z++;      // find a non-residue
 
    long long c = powmod(z, q, p);
    long long x = powmod(a, (q + 1) / 2, p);
    long long t = powmod(a, q, p);
    int m = s;
 
    while (t != 1) {
        int i = 0;
        long long tt = t;
        while (tt != 1) { tt = tt * tt % p; i++; }
        long long b = powmod(c, 1LL << (m - i - 1), p);
        x = x * b % p;
        c = b * b % p;
        t = t * c % p;
        m = i;
    }
    return x;
}

Paradigm

Group-theoretic descent. Each iteration halves the order of $t$ inside the 2-Sylow subgroup of ${\mathbb{Z}}_p^{\ast }$, so the loop terminates in at most $s$ rounds.

Complexity

• Time: $O({\log }^{2}p)$ expected — $O(\log p)$ for the non-residue search (constant tries × one modpow) plus $O(s\log p)\le O({\log }^{2}p)$ for the descent
• Space: $O(1)$

Why It Works

Throughout the loop the invariant $x^2\equiv at(\mathrm{mod}p)$ holds, and $t$ lies in the subgroup of order $2^m$. The value $i$ is the exact order of $t$ (as a power of two), and $b$ is chosen so that $b^2$ has the same order as $t$ but is its inverse there — multiplying kills the top bit. Since $m$ strictly decreases, we reach $t=1$, at which point the invariant reads $x^2\equiv a$. ∎

Variants / Use Cases

• Cipolla’s algorithm — an alternative with $O(\log p)$ worst case, better when $p-1$ is divisible by a huge power of two
• Quadratic residues — the theory this implements
• Discrete root — solving $x^k\equiv a$ generalises this via primitive roots
• Prime powers and composites — lift with Hensel lifting, then recombine with CRT
• Counting points on curves, Cornacchia — need a square root as a subroutine