Shanks-Tonelli algorithm

From CodeCodex

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

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

1 The algorithm
2 Implementations
- 2.1 C
3 External links

[edit] 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 <math>R^2 \equiv n \mod p</math>.

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

[edit] Implementations

[edit] C

Image:Info.png

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_add_ui(y, y, 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;
 }