The Tonelli-Shanks Algorithm, Explained

04 Jun, 2025

Modular arithmetic — perhaps surprisingly — has tons of practical applications, ranging from hashing strings to forming cryptosystems. And — maybe not so surprisingly — these involve solving equations which would normally be trivial in real numbers, but aren't nearly as simple in the discrete case. Such is the case with $x^2=a$, aka finding the square root of $a$. That is exactly where the Tonelli-Shanks algorithm comes in. It solves $x^2=a$ in — pardon my math — any cyclic group $G$ of even order, in $\Theta ({\log }^2\left|G\right|)$ group operations.

Sadly, the way it's often presented is convoluted both in terms of the pseudocode provided and the actual reasoning behind its correctness. As such, this article is my attempt at making sense of Tonelli-Shanks, in words myself from two years ago would understand.

The basics

I'll try not to overload the reader with too much mathematical notation and formalisms. That said, we do first need to introduce some basics. Feel free to skip over what you already know.

Modular arithmetic

Throughout this article, we'll be working with integers modulo a prime $p$. In simple terms, this means that after every operation we divide the result by $p$ and leave just the remainder, effectively limiting ourselves to the numbers $\{0,\dots ,p-1\}$. For example, $2·2=4\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}3)$. That last part just means we did this taking $p=3$. From now on I'll omit the special symbols and write $2·2=1$ where $p$ is known from the context. Furthermore, in our example $-4=3+3-4=2$, that is negative numbers get reduced as well.

An important property of multiplication modulo $p$ is that for any non-zero $a$ and $b$, it holds that $a·b\ne 0$. This is somewhat crucial, so we will assume the number $a$ in $x^2=a$ is never zero. We're not losing out on much anyway since the only solution to $x^2=0$ is zero itself.

Generators and square roots

It turns out there always exists a number $g$ such that if we keep multiplying it by itself (aka take $g^1,g^2,g^3,\dots$), we will eventually get all the numbers $\{1,\dots ,p-1\}$. This $g$ is called a generator.

It follows from the above that every number from $\{1,\dots ,p-1\}$ can be written as $g^i$ for some $i$ in $\{1,\dots ,p-1\}$. Additionally, the number of elements in $\{1,\dots ,p-1\}$ is $p-1$ which is even, unless $p=2$ (again, a trivial case: $\sqrt{0}=0$ and $\sqrt{1}=1$).

And these powers of $g$ cycle around. That is, $g^{p-1}$ becomes $1$ again, $g^p$ becomes $g$, etc. So the exponents sort of act like numbers modulo $p-1$. An important realization stems from this fact: only $a=g^k$ for even $k$ have square roots. This is because if $k$ is odd, then were any integer $i$ to satisfy $x^2=(g^i{)}^2=g^{2i}=g^k=a$, it would mean that $2i\equiv k(\mathrm{m}\mathrm{o}\mathrm{d}p-1)$, which can't be true since $p-1$ and $2i$ are even, and $k$ is not. We call these numbers quadratic residues, in contrast with quadratic nonresidues which do not have a square root.

Oh, and one more thing. You know how in real numbers $\sqrt{4}=2$? Well if you think about it, $-2$ is a perfectly valid solution too. We sort of arbitrarily choose the positive one. Among (non-zero) integers modulo $p$, every quadratic residue also has exactly two square roots, each being the negative of the other. Our algorithm will find one of them, and it's not really ever clear which one. If you badly need consistency, simply choosing the smaller one is an option.

Fast exponentiation

One last thing we'll need is an efficient way to calculate $a^n$ modulo $p$. This is accomplished using a simple recursive algorithm often referred to as fast exponentiation or square-and-multiply:

• If $n=1$, return $a$

• If $n$ is odd, return $a^{n-1}·a$

• Otherwise, set $r=a^{\frac{n}{2}}$ and return $r^2$

Clearly, this yields the correct value of $a^n$. And because every (at most) two iterations the value of $n$ halves, the runtime of this algorithm is $𝒪(\log n)$.

The Tonelli-Shanks algorithm

With the preliminaries out of the way, let's get down to business. The Tonelli-Shanks algorithm is as follows:

• Let $p$ be an odd prime. We will solve $x^2=a$ for a quadratic residue $a$.

• First, we somehow find a quadratic nonresidue $z$.

• We will keep two numbers $s$ and $t$, such that $a^s{z}^t=1$ at all times.

• Also, $t$ will always stay strictly "more even" than $s$, that is if some $2^k$ divides $s$ then $2^{k+1}$ divides $t$.

• To start out, $s=\frac{p-1}{2}$ and $t=p-1$.

• As long as $s$ is even:

  • If $a^{\frac{s}{2}}{z}^{\frac{t}{2}}=1$, set $s\leftarrow \frac{s}{2}$ and $t\leftarrow \frac{t}{2}$.

  • Otherwise, set $s\leftarrow \frac{s}{2}$ and $t\leftarrow \frac{t}{2}+\frac{p-1}{2}$.

• The solution to $x^2=a$ is $a^{\frac{s+1}{2}}{z}^{\frac{t}{2}}$.

So why on earth does this work? Let's go through the algorithm step by step.

Finding the nonresidue

For Tonelli-Shanks to work, we need a quadratic nonresidue $z$. Thankfully, finding one is easy. We can make use of a theorem called Euler's criterion:

If $z^{\frac{p-1}{2}}=1$ then $z$ is a quadratic residue.
Otherwise, $z^{\frac{p-1}{2}}=-1$ and $z$ is a nonresidue.

I won't show it here, but the proof is not that elaborate and you can easily find it online.

So we can check if any given $z$ is good in $𝒪(\log p)$ using the fast exponentiation algorithm described above. And as it turns out, we can just keep choosing a random $z$ until one passes our test. This is because exactly half of the numbers in $\{1,\dots ,p-1\}$ are nonresidues, so the expected number of times we'll have to repeat this is $\Theta (\log p)$. And again, I'll leave the proof as an exercise to the reader.

In total, this part of the algorithm takes $\Theta ({\log }^{2}p)$ time.

Initial state

We start with $s=\frac{p-1}{2}$ and $t=p-1$, which as we can see satisfies our invariants:

And obviously, $\frac{p-1}{2}$ is divisible by one less $2$ than $p-1$.

The loop

First of all, if $a^{\frac{s}{2}}{z}^{\frac{t}{2}}=1$ then we clearly can use these new values for $s$ and $t$, keeping both our invariants.

The second case is more interesting. We can't just do the same thing, because $a^{\frac{s}{2}}{z}^{\frac{t}{2}}\ne 1$ which would break one of our invariants. But, $a^s{z}^t=1$ is currently true. And if we realize that the only square roots of $1$ are $1$ and $-1$ (as there can only be two), then we know $a^{\frac{s}{2}}{z}^{\frac{t}{2}}=-1$, since dividing the exponent by two is exactly what taking a square root is. So if we multiply this number by another $-1$, we will get $1$ again. And this is what adding $\frac{p-1}{2}$ to $t/2$ does, as $z^{\frac{p-1}{2}}=-1$ per Euler's criterion:

How about the other invariant? Can we be sure that adding $\frac{p-1}{2}$ won't cause the new $t$ to be "less even" than $s$? Yes, because both $\frac{t}{2}$ and $\frac{p-1}{2}$ are multiples of $2·\frac{s}{2}$, so their sum is still divisible by at least one more $2$ than $\frac{s}{2}$.

The result

The final answer looks just as mysterious. But you can verify it's correct by simply taking its square:

Runtime complexity

This part is relatively simple. Each iteration involves two fast exponentiations, and because after each iteration $s$ becomes $\frac{s}{2}$, there can't be more than $𝒪(\log p)$ iterations before all the twos get divided out. In total, this nets us a $𝒪({\log }^{2}p)$ algorithm.

There is also the initial search for $z$, which is the only randomized part of Tonelli-Shanks. But unless we get hugely unlucky this will take a comparable amount of time, and for any given $p$ we only need to find $z$ once.

Finishing thoughts

That's it! As you can (hopefully) see, there is no magic involved. I wouldn't go as far as to call this algorithm intuitive, but you can surely make sense of it, even without any advanced math knowledge. And the implementation isn't hard either. Let's end with coding one up in Python:

import random

# These work for prime p > 2.
# Python's 3-argument pow performs fast exponentiation mod p.

def tonelli_shanks(a, z, p):
    s = (p - 1) // 2
    t = p - 1
    while s % 2 == 0:
        apow = pow(a, s // 2, p)
        zpow = pow(z, t // 2, p)
        azpow = apow * zpow % p
        if azpow == 1:
            s = s // 2
            t = t // 2
        else:
            s = s // 2
            t = t // 2 + (p - 1) // 2
    apow = pow(a, (s + 1) // 2, p)
    zpow = pow(z, t // 2, p)
    return apow * zpow % p

def is_nonresidue(x, p):
    return pow(x, (p - 1) // 2, p) == p - 1

def find_nonresidue(p):
    while True:
        z = random.randint(2, p - 1)
        if is_nonresidue(z, p):
            return z

Special thanks to tfkls for helpful feedback, as well as to dr. Lech Duraj for the original explanation given during his lecture.

#algebra #algorithms #discrete root #number theory