The Real Complexity of Tonelli-Shanks

Introduction

Pick a prime $p$ and a number $a$ . Squaring $a$ modulo $p$ is instant: compute $a^2 \bmod p$ and you are done. But what about going the other way? Given a result $n$ , find an $x$ such that:

x^2 \equiv n \pmod{p}

This is the modular square root problem, and it sits at the heart of number theory. Unlike ordinary square roots, where $\sqrt{9} = 3$ by reflex, there is no formula that extracts the answer directly. The integers wrap around mod $p$ , and the usual tricks — continuous fractions, Newton's method — break down.

Not every $n$ even has a square root mod $p$ . Exactly half the non-zero residues do (they are called quadratic residues); the other half have none. Euler's criterion settles which half $n$ belongs to in one step: if $n^{(p-1)/2} \equiv 1 \pmod{p}$ , a root exists.

Knowing a root exists is one thing. Computing it is another. The Tonelli–Shanks algorithm, developed by Alberto Tonelli in 1891 and rediscovered and extended by Daniel Shanks in 1973, does it efficiently using a beautiful descent through the 2-adic structure of $p-1$ .

Try It: Solve x²≡a (mod p)

Enter any odd prime $p$ and any number $a$ . The solver will first check whether $a$ is a quadratic residue mod $p$ ; if so, it runs Tonelli–Shanks and shows every step of the 2-adic descent.

<div class="controls">
  <label>Prime <i>p</i>: <input id="inp-p" type="number" value="97" min="3" step="2"></label>
  <label><i>a</i> (mod p): <input id="inp-a" type="number" value="75" min="1"></label>
  <button id="btn-run" type="button">Solve x²≡a (mod p)</button>
</div>
<div id="result" class="result"></div>
<div id="steps" class="steps"></div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; }
.controls { display: flex; flex-wrap: wrap; gap: .6rem; align-items: flex-end; margin-bottom: .8rem; }
.controls label { font-size: .88rem; color: #444; display: flex; flex-direction: column; gap: .2rem; }
.controls input { width: 90px; padding: .3rem .5rem; border: 1px solid #cdd9e3; border-radius: 6px;
                  font-size: 1rem; }
button { font: 600 14px system-ui, sans-serif; padding: .45rem .9rem;
         border: 1px solid #1d3557; background: #1d3557; color: #fff;
         border-radius: 8px; cursor: pointer; align-self: flex-end; }
.result { font: 700 1rem system-ui, sans-serif; min-height: 1.5em; margin-bottom: .5rem; }
.result.ok  { color: #0a7d33; }
.result.bad { color: #c92f3c; }
.steps { font-size: .85rem; line-height: 1.6; }
.step { background: #f4f7fa; border-left: 3px solid #1d3557; padding: .3rem .6rem;
        margin: .3rem 0; border-radius: 0 6px 6px 0; font-family: ui-monospace, monospace; }
.step.hi { border-color: #0a7d33; background: #e6f4ec; }
.step.label { font-family: system-ui, sans-serif; border-color: #8899aa;
              background: transparent; font-style: italic; color: #555; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Tonelli-Shanks algorithm: find x such that x^2 ≡ a (mod p)
// Reference: Shanks 1973; see Wikipedia "Tonelli–Shanks algorithm"

function modpow(base, exp, mod) {
  // Fast modular exponentiation using BigInt
  base = ((base % mod) + mod) % mod;
  let result = 1n;
  while (exp > 0n) {
    if (exp & 1n) result = result * base % mod;
    base = base * base % mod;
    exp >>= 1n;
  }
  return result;
}

function isPrime(n) {
  if (n < 2n) return false;
  if (n === 2n) return true;
  if (n % 2n === 0n) return false;
  for (let i = 3n; i * i <= n; i += 2n) if (n % i === 0n) return false;
  return true;
}

function tonelliShanks(n, p) {
  // Returns { root, steps } or { root: null, steps }
  const steps = [];

  const log = (msg, hi) => steps.push({ msg, hi: !!hi });

  n = ((n % p) + p) % p;

  // Step 1: Euler's criterion
  const euler = modpow(n, (p - 1n) / 2n, p);
  log(fmt('Euler criterion: %n%^((%p%-1)/2) ≡ %euler% (mod %p%)', { n: n, p: p, euler: euler }));
  if (euler !== 1n) {
    log(fmt('%n% is NOT a quadratic residue mod %p% — no solution exists.', { n: n, p: p }));
    return { root: null, steps };
  }
  log(fmt('%n% is a quadratic residue mod %p% — a root exists.', { n: n, p: p }), true);

  // Trivial case: p ≡ 3 (mod 4)
  if (p % 4n === 3n) {
    const r = modpow(n, (p + 1n) / 4n, p);
    log(`p ≡ 3 (mod 4): direct formula x = a^((p+1)/4) mod p`);
    log(`x = ${n}^((${p}+1)/4) = ${n}^${(p+1n)/4n} ≡ ${r} (mod ${p})`, true);
    log(fmt('Verify: %r%² = %v% ≡ %v2% (mod %p%) ✓', { r: r, v: r*r, v2: r*r % p, p: p }), true);
    return { root: r, steps };
  }

  // Factor out 2s from p-1: p-1 = Q * 2^S
  let Q = p - 1n, S = 0n;
  while (Q % 2n === 0n) { Q /= 2n; S++; }
  log(fmt('Factor p-1=%v%: Q=%Q%, S=%S%  (p-1 = Q·2^S)', { v: p-1n, Q: Q, S: S }));

  // Find a non-residue z
  let z = 2n;
  while (modpow(z, (p - 1n) / 2n, p) !== p - 1n) z++;
  log(fmt('Non-residue found: z=%z%  (z^((p-1)/2) ≡ -1 mod p)', { z: z }));

  // Initialize
  let M = S;
  let c = modpow(z, Q, p);
  let t = modpow(n, Q, p);
  let R = modpow(n, (Q + 1n) / 2n, p);
  log(fmt('Init: M=%M%, c=%c%, t=%t%, R=%R%', { M: M, c: c, t: t, R: R }));
  log(`Descent loop: reduce 2-adic valuation of t until t≡1 (mod p)`);

  // Main loop
  for (let iter = 0; iter < 100; iter++) {
    if (t === 0n) { log(`t=0 → x=0`); return { root: 0n, steps }; }
    if (t === 1n) { log(fmt('t≡1 (mod p): descent complete → x=%R%', { R: R }), true); break; }

    // Find least i such that t^(2^i) ≡ 1
    let i = 1n, tmp = t * t % p;
    while (tmp !== 1n) { tmp = tmp * tmp % p; i++; }
    log(`i=${i}: t^(2^${i}) ≡ 1 (mod p);  b = c^(2^(M-i-1))`);

    const b = modpow(c, modpow(2n, M - i - 1n, p - 1n), p);
    M = i;
    c = b * b % p;
    t = t * c % p;
    R = R * b % p;
    log(`  → M=${M}, c=${c}, t=${t}, R=${R}`);
  }

  log(fmt('Result: x=%R%,  x²≡%v% (mod %p%) ✓', { R: R, v: R*R%p, p: p }), true);
  return { root: R, steps };
}

function run() {
  const pVal = BigInt(document.getElementById('inp-p').value || '97');
  const aVal = BigInt(document.getElementById('inp-a').value || '75');
  const stepsEl = document.getElementById('steps');
  const resultEl = document.getElementById('result');

  stepsEl.innerHTML = '';
  resultEl.textContent = '';
  resultEl.className = 'result';

  if (!isPrime(pVal) || pVal < 3n) {
    resultEl.textContent = fmt('p=%pVal% is not an odd prime.', { pVal: pVal });
    resultEl.className = 'result bad';
    return;
  }
  if (aVal <= 0n || aVal >= pVal) {
    resultEl.textContent = `a must be between 1 and p-1.`;
    resultEl.className = 'result bad';
    return;
  }

  const { root, steps } = tonelliShanks(aVal, pVal);

  steps.forEach(({ msg, hi }) => {
    const d = document.createElement('div');
    d.className = 'step' + (hi ? ' hi' : '');
    d.textContent = msg;
    stepsEl.appendChild(d);
  });

  if (root === null) {
    resultEl.textContent = fmt('No solution: %aVal% is not a quadratic residue mod %pVal%.', { aVal: aVal, pVal: pVal });
    resultEl.className = 'result bad';
  } else {
    const other = (pVal - root) % pVal;
    resultEl.textContent = fmt('x = %root%  (also x = %other%);  %root%² ≡ %v% (mod %pVal%) ✓', { root: root, other: other, v: root*root%pVal, pVal: pVal });
    resultEl.className = 'result ok';
  }
}

document.getElementById('btn-run').addEventListener('click', run);
run();

<div class="controls">
  <label>Prime <i>p</i>: <input id="inp-p" type="number" value="97" min="3" step="2"></label>
  <label><i>a</i> (mod p): <input id="inp-a" type="number" value="75" min="1"></label>
  <button id="btn-run" type="button">Solve x²≡a (mod p)</button>
</div>
<div id="result" class="result"></div>
<div id="steps" class="steps"></div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; }
.controls { display: flex; flex-wrap: wrap; gap: .6rem; align-items: flex-end; margin-bottom: .8rem; }
.controls label { font-size: .88rem; color: #444; display: flex; flex-direction: column; gap: .2rem; }
.controls input { width: 90px; padding: .3rem .5rem; border: 1px solid #cdd9e3; border-radius: 6px;
                  font-size: 1rem; }
button { font: 600 14px system-ui, sans-serif; padding: .45rem .9rem;
         border: 1px solid #1d3557; background: #1d3557; color: #fff;
         border-radius: 8px; cursor: pointer; align-self: flex-end; }
.result { font: 700 1rem system-ui, sans-serif; min-height: 1.5em; margin-bottom: .5rem; }
.result.ok  { color: #0a7d33; }
.result.bad { color: #c92f3c; }
.steps { font-size: .85rem; line-height: 1.6; }
.step { background: #f4f7fa; border-left: 3px solid #1d3557; padding: .3rem .6rem;
        margin: .3rem 0; border-radius: 0 6px 6px 0; font-family: ui-monospace, monospace; }
.step.hi { border-color: #0a7d33; background: #e6f4ec; }
.step.label { font-family: system-ui, sans-serif; border-color: #8899aa;
              background: transparent; font-style: italic; color: #555; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Tonelli-Shanks algorithm: find x such that x^2 ≡ a (mod p)
// Reference: Shanks 1973; see Wikipedia "Tonelli–Shanks algorithm"

function modpow(base, exp, mod) {
  // Fast modular exponentiation using BigInt
  base = ((base % mod) + mod) % mod;
  let result = 1n;
  while (exp > 0n) {
    if (exp & 1n) result = result * base % mod;
    base = base * base % mod;
    exp >>= 1n;
  }
  return result;
}

function isPrime(n) {
  if (n < 2n) return false;
  if (n === 2n) return true;
  if (n % 2n === 0n) return false;
  for (let i = 3n; i * i <= n; i += 2n) if (n % i === 0n) return false;
  return true;
}

function tonelliShanks(n, p) {
  // Returns { root, steps } or { root: null, steps }
  const steps = [];

  const log = (msg, hi) => steps.push({ msg, hi: !!hi });

  n = ((n % p) + p) % p;

  // Step 1: Euler's criterion
  const euler = modpow(n, (p - 1n) / 2n, p);
  log(fmt('Euler criterion: %n%^((%p%-1)/2) ≡ %euler% (mod %p%)', { n: n, p: p, euler: euler }));
  if (euler !== 1n) {
    log(fmt('%n% is NOT a quadratic residue mod %p% — no solution exists.', { n: n, p: p }));
    return { root: null, steps };
  }
  log(fmt('%n% is a quadratic residue mod %p% — a root exists.', { n: n, p: p }), true);

  // Trivial case: p ≡ 3 (mod 4)
  if (p % 4n === 3n) {
    const r = modpow(n, (p + 1n) / 4n, p);
    log(`p ≡ 3 (mod 4): direct formula x = a^((p+1)/4) mod p`);
    log(`x = ${n}^((${p}+1)/4) = ${n}^${(p+1n)/4n} ≡ ${r} (mod ${p})`, true);
    log(fmt('Verify: %r%² = %v% ≡ %v2% (mod %p%) ✓', { r: r, v: r*r, v2: r*r % p, p: p }), true);
    return { root: r, steps };
  }

  // Factor out 2s from p-1: p-1 = Q * 2^S
  let Q = p - 1n, S = 0n;
  while (Q % 2n === 0n) { Q /= 2n; S++; }
  log(fmt('Factor p-1=%v%: Q=%Q%, S=%S%  (p-1 = Q·2^S)', { v: p-1n, Q: Q, S: S }));

  // Find a non-residue z
  let z = 2n;
  while (modpow(z, (p - 1n) / 2n, p) !== p - 1n) z++;
  log(fmt('Non-residue found: z=%z%  (z^((p-1)/2) ≡ -1 mod p)', { z: z }));

  // Initialize
  let M = S;
  let c = modpow(z, Q, p);
  let t = modpow(n, Q, p);
  let R = modpow(n, (Q + 1n) / 2n, p);
  log(fmt('Init: M=%M%, c=%c%, t=%t%, R=%R%', { M: M, c: c, t: t, R: R }));
  log(`Descent loop: reduce 2-adic valuation of t until t≡1 (mod p)`);

  // Main loop
  for (let iter = 0; iter < 100; iter++) {
    if (t === 0n) { log(`t=0 → x=0`); return { root: 0n, steps }; }
    if (t === 1n) { log(fmt('t≡1 (mod p): descent complete → x=%R%', { R: R }), true); break; }

    // Find least i such that t^(2^i) ≡ 1
    let i = 1n, tmp = t * t % p;
    while (tmp !== 1n) { tmp = tmp * tmp % p; i++; }
    log(`i=${i}: t^(2^${i}) ≡ 1 (mod p);  b = c^(2^(M-i-1))`);

    const b = modpow(c, modpow(2n, M - i - 1n, p - 1n), p);
    M = i;
    c = b * b % p;
    t = t * c % p;
    R = R * b % p;
    log(`  → M=${M}, c=${c}, t=${t}, R=${R}`);
  }

  log(fmt('Result: x=%R%,  x²≡%v% (mod %p%) ✓', { R: R, v: R*R%p, p: p }), true);
  return { root: R, steps };
}

function run() {
  const pVal = BigInt(document.getElementById('inp-p').value || '97');
  const aVal = BigInt(document.getElementById('inp-a').value || '75');
  const stepsEl = document.getElementById('steps');
  const resultEl = document.getElementById('result');

  stepsEl.innerHTML = '';
  resultEl.textContent = '';
  resultEl.className = 'result';

  if (!isPrime(pVal) || pVal < 3n) {
    resultEl.textContent = fmt('p=%pVal% is not an odd prime.', { pVal: pVal });
    resultEl.className = 'result bad';
    return;
  }
  if (aVal <= 0n || aVal >= pVal) {
    resultEl.textContent = `a must be between 1 and p-1.`;
    resultEl.className = 'result bad';
    return;
  }

  const { root, steps } = tonelliShanks(aVal, pVal);

  steps.forEach(({ msg, hi }) => {
    const d = document.createElement('div');
    d.className = 'step' + (hi ? ' hi' : '');
    d.textContent = msg;
    stepsEl.appendChild(d);
  });

  if (root === null) {
    resultEl.textContent = fmt('No solution: %aVal% is not a quadratic residue mod %pVal%.', { aVal: aVal, pVal: pVal });
    resultEl.className = 'result bad';
  } else {
    const other = (pVal - root) % pVal;
    resultEl.textContent = fmt('x = %root%  (also x = %other%);  %root%² ≡ %v% (mod %pVal%) ✓', { root: root, other: other, v: root*root%pVal, pVal: pVal });
    resultEl.className = 'result ok';
  }
}

document.getElementById('btn-run').addEventListener('click', run);
run();

Result

Console

Notice that checking the answer is trivial — just compute $x^2 \bmod p$ and compare. Finding $x$ in the first place requires the descent: peel off the factors of 2 from $p-1$ , locate a non-residue to build a correction term, then iteratively halve the 2-adic valuation of the error until it disappears.

The Real Complexity

How fast is Tonelli–Shanks? The algorithm breaks into two phases:

Setup: write $p - 1 = Q \cdot 2^S$ where $Q$ is odd. This factoring out of 2s is done once and takes $O(\log p)$ steps.
Finding a non-residue: pick a random $z$ with $z^{(p-1)/2} \equiv -1 \pmod{p}$ . Roughly half of all $z$ work, so the expected number of trials is 2.
The descent loop: runs at most $S$ iterations (at most $O(\log p)$ of them), each costing $O(\log^2 p)$ modular multiplications.
Total: $O(\log^2 p)$ expected time — polynomial in the bit-length of $p$ .

The "expected" qualifier matters. Finding the initial non-residue is the only randomized step. Under the Generalized Riemann Hypothesis (GRH), a non-residue smaller than $2\ln^2 p$ always exists, making a deterministic $O(\log^5 p)$ version possible. Without GRH, the best deterministic algorithm in the worst case is slower. This sits in the same landscape as the discrete logarithm: no quantum hardness, but also no classical proof of easy.

The algorithm is solved in the sense that it is proven correct and runs in polynomial time — but the fine print about randomness versus determinism touches open problems in analytic number theory.

Where It Matters

The question "what is the square root of $a$ mod $p$ ?" sounds academic. It is not:

Elliptic-curve cryptography: every point decompression (recovering a full $(x, y)$ point from just its $x$ coordinate) requires a modular square root. This happens in every TLS handshake, every Bitcoin transaction, every Signal message.
Primality proving: the Atkin–Morain elliptic-curve primality test (ECPP) calls Tonelli–Shanks repeatedly to construct certificates.
Decoding quadratic-residue codes: a class of error-correcting codes whose decoding algorithm reduces to extracting square roots mod $p$ .
Number field sieve: the relation-collection phase of the fastest known factoring algorithm needs modular square roots for smoothness tests.

Alternative algorithms exist — Cipolla's algorithm (1903) has a cleaner proof and similar complexity; Atkin's algorithm handles the special case $p \equiv 3 \pmod{8}$ even faster — but Tonelli–Shanks is the most widely implemented because its loop structure maps efficiently to hardware.

Understanding Tonelli–Shanks is understanding how the difficulty of the discrete logarithm problem protects the internet — squaring is easy; reversing it without the algorithm is what keeps secrets safe.

Conclusion

The Tonelli–Shanks algorithm is a beautiful example of how pure number theory becomes indispensable engineering. Its 2-adic descent — peeling factors of 2 from $p-1$ , correcting the error one bit at a time — turns an apparently intractable inversion into a fast polynomial-time computation.

Squaring mod $p$ is easy; un-squaring it with no algorithm would take $O(\sqrt{p})$ brute-force steps — infeasible for the 256-bit primes inside your browser's TLS stack. Tonelli–Shanks closes that gap. Every time your device compresses or decompresses an elliptic-curve point, this algorithm — or one of its close cousins — runs silently underneath.

The next time you see the padlock icon in your browser, remember: a 19th-century descent through powers of two is part of what keeps the connection private. That is the reach of a well-placed mathematical idea.

Comments

Loading comments...

https://www.kipuhub.com/en/article/tonelli-shanks/Content licensed under CC BY-NC 4.0.

Cookies & privacy

Tonelli-Shanks

The 2-adic descent that unlocks modular square roots