The site's tenth Number Theory entry, and a companion to a harder problem already
here. Baby-Step Giant-Step recovers an exponent —
given g, h, and p, find x with gx ≡
h (mod p) — and no known algorithm does that in time polynomial in the digits of p.
Tonelli–Shanks recovers something else: given a prime p and a number n that's
a quadratic residue mod p (a square, mod p), find r
with r² ≡ n (mod p) — a modular square root. It sounds like the same shape of
problem, "recover something Modular
Exponentiation would otherwise compute forward," but this one has
no hardness gap at all: Tonelli–Shanks solves it in polynomial time, always, for any odd prime modulus.
The two entries make a useful pair precisely because they look alike and aren't: one recovery problem is
provably easy, the other is the basis of Diffie–Hellman.
Enter a prime p and a quadratic residue n, then step through. p = 19,
n = 6 hits the fast path: 19 ≡ 3 (mod 4), so the answer falls out of one exponentiation, no loop
needed. p = 41, n = 5 needs the general loop, and takes it through two full iterations —
toggle "raw exponent (broken)" below on this preset to see the first pitfall: an
exponentiation bug this page's own draft implementation shipped with, until the verification step below
caught it. p = 41, n = 6 has no solution — 6 fails Euler's criterion, and the demo reports
that cleanly before any loop runs.
loop trail — each entry is one search-and-update round:
Before anything else, check whether n is even a quadratic residue mod p at
all: Euler's criterion says n(p−1)/2 mod p comes out to
1 if n is a residue and p − 1 (that is, −1) if it
isn't — one exponentiation, using the same modPow routine every other entry in this category
shares. If it's not a residue, stop: no square root exists, and there's no point running the rest of the
algorithm to find that out (see the second pitfall for what happens if this check
gets skipped anyway).
If n passes, the modulus splits the work into two shapes. When p ≡ 3 (mod 4),
there's a direct formula: r = n(p+1)/4 mod p. Squaring both sides shows why it
works — r² = n(p+1)/2 = n · n(p−1)/2 = n · 1 = n, using Euler's
criterion's own result for the last step, since n was already confirmed a residue. One
exponentiation, done.
Most primes aren't 3 mod 4, though, and that formula doesn't extend to them: when p ≡ 1 (mod 4),
(p + 1) / 4 isn't even an integer, so there's real work to do. Write p − 1 = Q · 2S
with Q odd (every p − 1 is even, so this factoring always succeeds), find any
quadratic non-residue z (trying z = 2, 3, 4, … in order — roughly half of
all nonzero residues are non-residues, so this takes only a couple of tries in expectation, checked against
Euler's criterion each time), and initialize:
M = S
c = z^Q mod p
t = n^Q mod p
R = n^((Q+1)/2) mod pThen loop: if t == 1, R is the answer. Otherwise, find the least i
with 0 < i < M such that t2^i ≡ 1 (mod p) — that i is
guaranteed to exist as long as everything upstream stayed correct, because t's order always
divides 2M by construction. Set b = c2^(M−i−1) mod p, then
update M = i, c = b² mod p, t = t · b² mod p, R = R · b mod p,
and repeat. Each round strictly shrinks M, so the loop terminates in at most S
rounds — and S is bounded by the number of times 2 divides p − 1, at most
log₂(p − 1).
function modPow(base, exp, mod) {
base = ((base % mod) + mod) % mod;
let result = 1n % mod;
while (exp > 0n) {
if (exp & 1n) result = (result * base) % mod;
base = (base * base) % mod;
exp >>= 1n;
}
return result;
}
function tonelliShanks(n, p) { // p odd prime, n reduced mod p; returns null if n has no square root
n = ((n % p) + p) % p;
if (n === 0n) return 0n;
const euler = modPow(n, (p - 1n) / 2n, p);
if (euler === p - 1n) return null; // Euler's criterion: n is not a quadratic residue
if (p % 4n === 3n) { // fast path
return modPow(n, (p + 1n) / 4n, p);
}
let Q = p - 1n, S = 0n; // general path: factor p-1 = Q * 2^S, Q odd
while (Q % 2n === 0n) { Q /= 2n; S++; }
let z = 2n;
while (modPow(z, (p - 1n) / 2n, p) !== p - 1n) z++; // first quadratic non-residue
let M = S;
let c = modPow(z, Q, p);
let t = modPow(n, Q, p);
let R = modPow(n, (Q + 1n) / 2n, p);
while (t !== 1n) {
let i = 0n, temp = t;
while (temp !== 1n) { temp = (temp * temp) % p; i++; } // least i with t^(2^i) = 1
const b = modPow(c, 1n << (M - i - 1n), p); // c^(2^(M-i-1)) -- the exponent is itself a power of two
M = i;
c = (b * b) % p;
t = (t * c) % p;
R = (R * b) % p;
}
return R;
}The interactive demo above uses an equivalent generator function that yields one event per Euler-criterion
check, per fast-path computation, and per general-loop search-and-update round, purely to drive the
step-through — the arithmetic is identical, including the 1n << (M - i - 1n) bit-shift
that computes 2M−i−1 as the exponent (see the first pitfall for why that shift
matters and can't be simplified away). Verified two ways before writing any of the prose above: 8,000
randomized trials (random odd prime p up to 20 bits, random x, n = x² mod p,
confirming tonelliShanks(n, p) squares back to n) plus an exhaustive sweep of
every (p, n) pair for every prime p < 2000 — 277,048 pairs, covering both
residues and non-residues, 0 failures either way. The three worked presets above were each checked directly:
p = 19, n = 6 → r = 5 (5² mod 19 = 6), p = 41, n = 5 → r = 28
(28² mod 41 = 5), and p = 41, n = 6 confirmed to have no square root by brute-force
squaring every value from 0 to 40.
The exponent in b = c2^(M−i−1) is itself a power of two, not the bare
integer M − i − 1 — and this page's own first draft got it wrong. Writing
modPow(c, M - i - 1n, p) instead of modPow(c, 1n << (M - i - 1n), p) compiles,
runs, and produces plausible-looking BigInts at every step, with no error until the loop hits an input it
can't recover from. Traced on p = 41, n = 5 (toggle the checkbox above on this preset to watch
it live): the correct first round computes i = 2, exponent 2(3−2−1) =
20 = 1, so b = c1 = 38. The broken version computes the same
i = 2 but takes the exponent literally — M − i − 1 = 0 — giving b = c0
= 1. A step that should shrink t toward 1 instead multiplies everything by 1 and changes
nothing real. The second round then searches for a valid i in the now-shrunk range
0 < i < 2 and never finds one — t stays stuck at a value whose order doesn't
fit that range, because the previous round's no-op step never actually reduced it. The two exponents happen
to agree whenever M − i − 1 is exactly 0 or 1 (since 2⁰ = 1 lands close enough to
look right in isolation), which is exactly what let this ship past a first read-through — it took tracing
the second round, not the first, to see the values stop moving.
Skipping Euler's criterion doesn't make the general loop guess a wrong answer for a
non-residue — it makes the inner search exhaust its range and find nothing, every time. Removing
the upfront check and running the general-case loop directly on p = 41, n = 6 (a genuine
non-residue) doesn't silently return a bogus root: the inner "least i" search reaches
i = M without t ever hitting 1, the same out-of-range failure the first pitfall
produces. Checked across all 22 non-residues mod 41: every single one exhausts the search this way, never
once returning a plausible-but-wrong value instead. The reason is structural, not coincidental: t = nQ
always has order dividing 2S by Fermat's little theorem regardless of whether
n is a residue, but only when n genuinely is one does that order divide
2S strictly — for a non-residue, t's order is exactly
2S, the one value the open range 0 < i < M is built to exclude.
The Euler check isn't just an early-exit optimization, in other words — without it, a genuine non-residue
crashes the general path outright rather than reporting "no solution" cleanly the way the fast path does.
The fast-path formula only works when p ≡ 3 (mod 4) — applying it unconditionally
is a common shortcut, and it fails silently rather than with an error. Many minimal Tonelli–Shanks
write-ups only implement the p ≡ 3 (mod 4) case and stop there, sometimes without stating the
restriction prominently. Feeding p = 41 (41 ≡ 1 mod 4) into the fast-path formula anyway:
(p + 1) / 4 = 42 / 4, which isn't an integer — but BigInt division truncates instead of
erroring, silently giving exponent 10 instead of a real fourth root of anything. Checked
directly: n = 5, exponent 10, gives r = 40, and
40² mod 41 = 1, not 5 — a confidently wrong answer with no exception, no NaN,
nothing that looks like failure. The general loop above isn't an optional fallback for exotic moduli; it's
required for the majority of primes, since half of all odd primes are 1 mod 4.
O(log² p) expected. The fast path is one modPow call, O(log p)
modular multiplications. The general path runs at most S ≤ log₂(p − 1) loop rounds, each doing
an inner search of at most M ≤ S squarings plus a couple of modPow calls
(O(log p) each) — O(S · log p) = O(log² p) total, plus the O(log p)
expected cost of finding a non-residue z (a small constant number of Euler's-criterion checks in
practice). That's polynomial in the number of digits of p, the same shape of bound
Miller–Rabin and
Modular Exponentiation get, and a completely different regime from
Baby-Step Giant-Step's O(√p) — exponential
in the digits of p, because no polynomial algorithm for that problem is known. Both entries
answer a "recover X" question about a prime modulus; one happens to be easy and the other is the reason
Diffie–Hellman is considered hard at all. Space is O(1) beyond the BigInts themselves — no table
to build, unlike baby-step giant-step's O(√p) memory cost.
This site's guide, Choosing a Number Theory Algorithm, compares this entry against the other thirteen Number Theory entries side by side.