Shanks-Tonelli algorithm

The Shanks-Tonelli algorithm is used within modular arithmetic to solve a congruence of the form $x^2 \equiv n \mod p$, where n is a quadratic residue (mod p), and typically $p \equiv 1 \mod 4$.

When $p \equiv 3 \mod 4$, it is much more efficient to use the following identity: $x \equiv n^{\frac{p+1}{4}} \mod p$.

The algorithm

Inputs: p, an odd prime. n, an integer which is a quadratic residue (mod p), meaning that the Legendre symbol] (n/p) = 1.

Outputs: R, an integer satisfying $R^2 \equiv n \mod p$.

1. Factor out powers of 2 from (p-1), defining Q and S as: $p-1 = Q2^S$.
2. Select a W such that the Legendre symbol (W/p) = -1 (that is, W should be a quadratic non-residue modulo p).
3. Let $R = n^{\frac{Q+1}{2}} \mod p$.
4. Let $V = W^Q \mod p$.
5. Loop:
1. Find the lowest i, $0 \leq i \leq n-1$, such that $(R^2n^{-1})^{2^{i}} \equiv 1 \mod p$. This can be done efficiently by starting with $R^2n^{-1}$ and squaring (mod p) until the result is 1.
2. If $i = 0$, return R.
3. Otherwise, let $R' = RV^{2^{S-i-1}} (mod p)$ and repeat the loop with R' as the new R.

Implementations

C

 This code requires GNU Multiple Precision Arithmetic Library, an external library, to run.
#include <gmp.h>

#define mpz_rshift(A,B,l) mpz_tdiv_q_2exp(A, B, l)

int root(mpz_t result, const mpz_t arg, const mpz_t prime)
{
mpz_t y, b, t;
unsigned int r, m;

if (mpz_divisible_p(arg, prime)) {
mpz_set_ui(result, 0);
return 1;
}
if (mpz_legendre(arg, prime) == -1)
return -1;

mpz_init(b);
mpz_init(t);
mpz_init_set_ui(y, 2);
while(mpz_legendre(y, prime) != -1)

mpz_sub_ui(result, prime, 1);
r = mpz_scan1(result, 0);
mpz_rshift(result, result, r);

mpz_powm(y, y, result, prime);

mpz_rshift(result, result, 1);
mpz_powm(b, arg, result, prime);

mpz_mul(result, arg, b);
mpz_mod(result, result, prime);

mpz_mul(b, result, b);
mpz_mod(b, b, prime);

while(mpz_cmp_ui(b, 1)) {

mpz_mul(t, b, b);
mpz_mod(t, t, prime);
for(m = 1; mpz_cmp_ui(t, 1); m++) {
mpz_mul(t, t, t);
mpz_mod(t, t, prime);
}

mpz_set_ui(t, 0);
mpz_setbit(t, r - m - 1);
mpz_powm(t, y, t, prime);

mpz_mul(y, t, t);

r = m;

mpz_mul(result, result, t);
mpz_mod(result, result, prime);

mpz_mul(b, b, y);
mpz_mod(b, b, prime);
}

mpz_clear(y);
mpz_clear(b);
mpz_clear(t);
return 1;
}