Find Square Root under Modulo p | Set 2 (Shanks Tonelli algorithm)

Last Updated : 20 Dec, 2022

Given a number ‘n’ and a prime ‘p’, find square root of n under modulo p if it exists.
Examples:

Input: n = 2, p = 113 Output: 62 62^2 = 3844 and 3844 % 113 = 2 Input: n = 2, p = 7 Output: 3 or 4 3 and 4 both are square roots of 2 under modulo 7 because (3*3) % 7 = 2 and (4*4) % 7 = 2 Input: n = 2, p = 5 Output: Square root doesn't exist

We have discussed Euler’s criterion to check if square root exists or not. We have also discussed a solution that works only when p is in form of 4*i + 3
In this post, Shank Tonelli’s algorithm is discussed that works for all types of inputs.
Algorithm steps to find modular square root using shank Tonelli’s algorithm :
1) Calculate n ^ ((p - 1) / 2) (mod p), it must be 1 or p-1, if it is p-1, then modular square root is not possible.
2) Then after write p-1 as (s * 2^e) for some integer s and e, where s must be an odd number and both s and e should be positive.
3) Then find a number q such that q ^ ((p - 1) / 2) (mod p) = -1
4) Initialize variable x, b, g and r by following values

 x = n ^ ((s + 1) / 2 (first guess of square root) b = n ^ s g = q ^ s r = e (exponent e will decrease after each updation)

5) Now loop until m > 0 and update value of x, which will be our final answer.

 Find least integer m such that b^(2^m) = 1(mod p) and 0 <= m <= r – 1 If m = 0, then we found correct answer and return x as result Else update x, b, g, r as below x = x * g ^ (2 ^ (r – m - 1)) b = b * g ^(2 ^ (r - m)) g = g ^ (2 ^ (r - m)) r = m

so if m becomes 0 or b becomes 1, we terminate and print the result. This loop guarantees to terminate because value of m is decreased each time after updation.
Following is the implementation of above algorithm.

C++

// C++ program to implement Shanks Tonelli algorithm for // finding Modular Square Roots #include <bits/stdc++.h> using namespace std; // utility function to find pow(base, exponent) % modulus int pow(int base, int exponent, int modulus) {  int result = 1;  base = base % modulus;  while (exponent > 0)  {  if (exponent % 2 == 1)  result = (result * base)% modulus;  exponent = exponent >> 1;  base = (base * base) % modulus;  }  return result; } // utility function to find gcd int gcd(int a, int b) {  if (b == 0)  return a;  else  return gcd(b, a % b); } // Returns k such that b^k = 1 (mod p) int order(int p, int b) {  if (gcd(p, b) != 1)  {  printf("p and b are not co-prime.\n");  return -1;  }  // Initializing k with first odd prime number  int k = 3;  while (1)  {  if (pow(b, k, p) == 1)  return k;  k++;  } } // function return p - 1 (= x argument) as x * 2^e, // where x will be odd sending e as reference because // updation is needed in actual e int convertx2e(int x, int& e) {  e = 0;  while (x % 2 == 0)  {  x /= 2;  e++;  }  return x; } // Main function for finding the modular square root int STonelli(int n, int p) {  // a and p should be coprime for finding the modular  // square root  if (gcd(n, p) != 1)  {  printf("a and p are not coprime\n");  return -1;  }  // If below expression return (p - 1) then modular  // square root is not possible  if (pow(n, (p - 1) / 2, p) == (p - 1))  {  printf("no sqrt possible\n");  return -1;  }  // expressing p - 1, in terms of s * 2^e, where s  // is odd number  int s, e;  s = convertx2e(p - 1, e);  // finding smallest q such that q ^ ((p - 1) / 2)  // (mod p) = p - 1  int q;  for (q = 2; ; q++)  {  // q - 1 is in place of (-1 % p)  if (pow(q, (p - 1) / 2, p) == (p - 1))  break;  }  // Initializing variable x, b and g  int x = pow(n, (s + 1) / 2, p);  int b = pow(n, s, p);  int g = pow(q, s, p);  int r = e;  // keep looping until b become 1 or m becomes 0  while (1)  {  int m;  for (m = 0; m < r; m++)  {  if (order(p, b) == -1)  return -1;  // finding m such that b^ (2^m) = 1  if (order(p, b) == pow(2, m))  break;  }  if (m == 0)  return x;  // updating value of x, g and b according to  // algorithm  x = (x * pow(g, pow(2, r - m - 1), p)) % p;  g = pow(g, pow(2, r - m), p);  b = (b * g) % p;  if (b == 1)  return x;  r = m;  } } // driver program to test above function int main() {  int n = 2;  // p should be prime  int p = 113;  int x = STonelli(n, p);  if (x == -1)  printf("Modular square root is not exist\n");  else  printf("Modular square root of %d and %d is %d\n",  n, p, x); }

Java

// Java program to implement Shanks // Tonelli algorithm for finding  // Modular Square Roots  import java.util.*; class GFG {  static int z = 0;   // utility function to find // pow(base, exponent) % modulus  static int pow1(int base1,   int exponent, int modulus)  {   int result = 1;   base1 = base1 % modulus;   while (exponent > 0)   {   if (exponent % 2 == 1)   result = (result * base1) % modulus;   exponent = exponent >> 1;   base1 = (base1 * base1) % modulus;   }   return result;  }  // utility function to find gcd  static int gcd(int a, int b)  {   if (b == 0)   return a;   else  return gcd(b, a % b);  }  // Returns k such that b^k = 1 (mod p)  static int order(int p, int b)  {   if (gcd(p, b) != 1)   {   System.out.println("p and b are" +   "not co-prime.");   return -1;   }   // Initializing k with first  // odd prime number   int k = 3;   while (true)   {   if (pow1(b, k, p) == 1)   return k;   k++;   }  }  // function return p - 1 (= x argument) // as x * 2^e, where x will be odd  // sending e as reference because  // updation is needed in actual e  static int convertx2e(int x)  {   z = 0;  while (x % 2 == 0)   {   x /= 2;   z++;   }   return x;  }  // Main function for finding  // the modular square root  static int STonelli(int n, int p)  {   // a and p should be coprime for   // finding the modular square root   if (gcd(n, p) != 1)   {   System.out.println("a and p are not coprime");   return -1;   }   // If below expression return (p - 1) then modular   // square root is not possible   if (pow1(n, (p - 1) / 2, p) == (p - 1))   {   System.out.println("no sqrt possible");   return -1;   }   // expressing p - 1, in terms of   // s * 2^e, where s is odd number   int s, e;   s = convertx2e(p - 1);  e = z;  // finding smallest q such that q ^ ((p - 1) / 2)   // (mod p) = p - 1   int q;   for (q = 2; ; q++)   {   // q - 1 is in place of (-1 % p)   if (pow1(q, (p - 1) / 2, p) == (p - 1))   break;   }   // Initializing variable x, b and g   int x = pow1(n, (s + 1) / 2, p);   int b = pow1(n, s, p);   int g = pow1(q, s, p);   int r = e;   // keep looping until b   // become 1 or m becomes 0   while (true)   {   int m;   for (m = 0; m < r; m++)   {   if (order(p, b) == -1)   return -1;   // finding m such that b^ (2^m) = 1   if (order(p, b) == Math.pow(2, m))   break;   }   if (m == 0)   return x;   // updating value of x, g and b  // according to algorithm   x = (x * pow1(g, (int)Math.pow(2,   r - m - 1), p)) % p;   g = pow1(g, (int)Math.pow(2, r - m), p);   b = (b * g) % p;   if (b == 1)   return x;   r = m;   }  }  // Driver code  public static void main (String[] args)  {  int n = 2;   // p should be prime   int p = 113;   int x = STonelli(n, p);   if (x == -1)  System.out.println("Modular square" +   "root is not exist\n");   else  System.out.println("Modular square root of " +  n + " and " + p + " is " +  x + "\n");  }  } // This code is contributed by mits

Python3

# Python3 program to implement Shanks Tonelli # algorithm for finding Modular Square Roots  # utility function to find pow(base,  # exponent) % modulus  def pow1(base, exponent, modulus): result = 1; base = base % modulus; while (exponent > 0): if (exponent % 2 == 1): result = (result * base) % modulus; exponent = int(exponent) >> 1; base = (base * base) % modulus; return result; # utility function to find gcd  def gcd(a, b): if (b == 0): return a; else: return gcd(b, a % b); # Returns k such that b^k = 1 (mod p)  def order(p, b): if (gcd(p, b) != 1): print("p and b are not co-prime.\n"); return -1; # Initializing k with first  # odd prime number  k = 3; while (True): if (pow1(b, k, p) == 1): return k; k += 1; # function return p - 1 (= x argument) as  # x * 2^e, where x will be odd sending e  # as reference because updation is needed # in actual e  def convertx2e(x): z = 0; while (x % 2 == 0): x = x / 2; z += 1; return [x, z]; # Main function for finding the # modular square root  def STonelli(n, p): # a and p should be coprime for # finding the modular square root  if (gcd(n, p) != 1): print("a and p are not coprime\n"); return -1; # If below expression return (p - 1) then # modular square root is not possible  if (pow1(n, (p - 1) / 2, p) == (p - 1)): print("no sqrt possible\n"); return -1; # expressing p - 1, in terms of s * 2^e,  # where s is odd number  ar = convertx2e(p - 1); s = ar[0]; e = ar[1]; # finding smallest q such that  # q ^ ((p - 1) / 2) (mod p) = p - 1  q = 2; while (True): # q - 1 is in place of (-1 % p)  if (pow1(q, (p - 1) / 2, p) == (p - 1)): break; q += 1; # Initializing variable x, b and g  x = pow1(n, (s + 1) / 2, p); b = pow1(n, s, p); g = pow1(q, s, p); r = e; # keep looping until b become  # 1 or m becomes 0  while (True): m = 0; while (m < r): if (order(p, b) == -1): return -1; # finding m such that b^ (2^m) = 1  if (order(p, b) == pow(2, m)): break; m += 1; if (m == 0): return x; # updating value of x, g and b  # according to algorithm  x = (x * pow1(g, pow(2, r - m - 1), p)) % p; g = pow1(g, pow(2, r - m), p); b = (b * g) % p; if (b == 1): return x; r = m; # Driver Code n = 2; # p should be prime  p = 113; x = STonelli(n, p); if (x == -1): print("Modular square root is not exist\n"); else: print("Modular square root of", n, "and", p, "is", x); # This code is contributed by mits

// C# program to implement Shanks // Tonelli algorithm for finding  // Modular Square Roots  using System; class GFG {   static int z=0;   // utility function to find  // pow(base, exponent) % modulus  static int pow1(int base1,   int exponent, int modulus)  {   int result = 1;   base1 = base1 % modulus;   while (exponent > 0)   {   if (exponent % 2 == 1)   result = (result * base1) % modulus;   exponent = exponent >> 1;   base1 = (base1 * base1) % modulus;   }   return result;  }  // utility function to find gcd  static int gcd(int a, int b)  {   if (b == 0)   return a;   else  return gcd(b, a % b);  }  // Returns k such that b^k = 1 (mod p)  static int order(int p, int b)  {   if (gcd(p, b) != 1)   {   Console.WriteLine("p and b are" +  "not co-prime.");   return -1;   }   // Initializing k with   // first odd prime number   int k = 3;   while (true)   {   if (pow1(b, k, p) == 1)   return k;   k++;   }  }  // function return p - 1 (= x argument)  // as x * 2^e, where x will be odd sending // e as reference because updation is  // needed in actual e  static int convertx2e(int x)  {   z = 0;  while (x % 2 == 0)   {   x /= 2;   z++;   }   return x;  }  // Main function for finding // the modular square root  static int STonelli(int n, int p)  {   // a and p should be coprime for   // finding the modular square root   if (gcd(n, p) != 1)   {   Console.WriteLine("a and p are not coprime");   return -1;   }   // If below expression return (p - 1) then   // modular square root is not possible   if (pow1(n, (p - 1) / 2, p) == (p - 1))   {   Console.WriteLine("no sqrt possible");   return -1;   }   // expressing p - 1, in terms of s * 2^e,  // where s is odd number   int s, e;   s = convertx2e(p - 1);  e=z;  // finding smallest q such that q ^ ((p - 1) / 2)   // (mod p) = p - 1   int q;   for (q = 2; ; q++)   {   // q - 1 is in place of (-1 % p)   if (pow1(q, (p - 1) / 2, p) == (p - 1))   break;   }   // Initializing variable x, b and g   int x = pow1(n, (s + 1) / 2, p);   int b = pow1(n, s, p);   int g = pow1(q, s, p);   int r = e;   // keep looping until b become   // 1 or m becomes 0   while (true)   {   int m;   for (m = 0; m < r; m++)   {   if (order(p, b) == -1)   return -1;   // finding m such that b^ (2^m) = 1   if (order(p, b) == Math.Pow(2, m))   break;   }   if (m == 0)   return x;   // updating value of x, g and b   // according to algorithm   x = (x * pow1(g, (int)Math.Pow(2, r - m - 1), p)) % p;   g = pow1(g, (int)Math.Pow(2, r - m), p);   b = (b * g) % p;   if (b == 1)   return x;   r = m;   }  }  // Driver code  static void Main()  {   int n = 2;   // p should be prime   int p = 113;   int x = STonelli(n, p);   if (x == -1)  Console.WriteLine("Modular square root" +  "is not exist\n");   else  Console.WriteLine("Modular square root of" +   "{0} and {1} is {2}\n", n, p, x);  }  } // This code is contributed by mits

PHP

<?php // PHP program to implement Shanks Tonelli  // algorithm for finding Modular Square Roots  // utility function to find pow(base, // exponent) % modulus  function pow1($base, $exponent, $modulus) { $result = 1; $base = $base % $modulus; while ($exponent > 0) { if ($exponent % 2 == 1) $result = ($result * $base) % $modulus; $exponent = $exponent >> 1; $base = ($base * $base) % $modulus; } return $result; } // utility function to find gcd  function gcd($a, $b) { if ($b == 0) return $a; else return gcd($b, $a % $b); } // Returns k such that b^k = 1 (mod p)  function order($p, $b) { if (gcd($p, $b) != 1) { print("p and b are not co-prime.\n"); return -1; } // Initializing k with first // odd prime number  $k = 3; while (1) { if (pow1($b, $k, $p) == 1) return $k; $k++; } } // function return p - 1 (= x argument) as // x * 2^e, where x will be odd sending e  // as reference because updation is needed // in actual e  function convertx2e($x, &$e) { $e = 0; while ($x % 2 == 0) { $x = (int)($x / 2); $e++; } return $x; } // Main function for finding the  // modular square root  function STonelli($n, $p) { // a and p should be coprime for  // finding the modular square root  if (gcd($n, $p) != 1) { print("a and p are not coprime\n"); return -1; } // If below expression return (p - 1) then  // modular square root is not possible  if (pow1($n, ($p - 1) / 2, $p) == ($p - 1)) { printf("no sqrt possible\n"); return -1; } // expressing p - 1, in terms of s * 2^e,  // where s is odd number  $e = 0; $s = convertx2e($p - 1, $e); // finding smallest q such that  // q ^ ((p - 1) / 2) (mod p) = p - 1  $q = 2; for (; ; $q++) { // q - 1 is in place of (-1 % p)  if (pow1($q, ($p - 1) / 2, $p) == ($p - 1)) break; } // Initializing variable x, b and g  $x = pow1($n, ($s + 1) / 2, $p); $b = pow1($n, $s, $p); $g = pow1($q, $s, $p); $r = $e; // keep looping until b become // 1 or m becomes 0  while (1) { $m = 0; for (; $m < $r; $m++) { if (order($p, $b) == -1) return -1; // finding m such that b^ (2^m) = 1  if (order($p, $b) == pow(2, $m)) break; } if ($m == 0) return $x; // updating value of x, g and b  // according to algorithm  $x = ($x * pow1($g, pow(2, $r - $m - 1), $p)) % $p; $g = pow1($g, pow(2, $r - $m), $p); $b = ($b * $g) % $p; if ($b == 1) return $x; $r = $m; } } // Driver Code $n = 2; // p should be prime  $p = 113; $x = STonelli($n, $p); if ($x == -1) print("Modular square root is not exist\n"); else print("Modular square root of " . "$n and $p is $x\n"); // This code is contributed by mits ?>

JavaScript

<script> // JavaScript program to implement Shanks // Tonelli algorithm for finding  // Modular Square Roots   let z = 0;   // utility function to find // pow(base, exponent) % modulus  function pow1(base1,   exponent, modulus)  {   let result = 1;   base1 = base1 % modulus;   while (exponent > 0)   {   if (exponent % 2 == 1)   result = (result * base1) % modulus;   exponent = exponent >> 1;   base1 = (base1 * base1) % modulus;   }   return result;  }    // utility function to find gcd  function gcd(a, b)  {   if (b == 0)   return a;   else  return gcd(b, a % b);  }    // Returns k such that b^k = 1 (mod p)  function order(p, b)  {   if (gcd(p, b) != 1)   {   document.write("p and b are" +   "not co-prime." + "<br/>");   return -1;   }     // Initializing k with first  // odd prime number   let k = 3;   while (true)   {   if (pow1(b, k, p) == 1)   return k;   k++;   }  }    // function return p - 1 (= x argument) // as x * 2^e, where x will be odd  // sending e as reference because  // updation is needed in actual e  function convertx2e(x)  {   z = 0;  while (x % 2 == 0)   {   x /= 2;   z++;   }   return x;  }    // Main function for finding  // the modular square root  function STonelli(n, p)  {   // a and p should be coprime for   // finding the modular square root   if (gcd(n, p) != 1)   {   System.out.prletln("a and p are not coprime");   return -1;   }     // If below expression return (p - 1) then modular   // square root is not possible   if (pow1(n, (p - 1) / 2, p) == (p - 1))   {   document.write("no sqrt possible" + "<br/>");   return -1;   }     // expressing p - 1, in terms of   // s * 2^e, where s is odd number   let s, e;   s = convertx2e(p - 1);  e = z;    // finding smallest q such that q ^ ((p - 1) / 2)   // (mod p) = p - 1   let q;   for (q = 2; ; q++)   {   // q - 1 is in place of (-1 % p)   if (pow1(q, (p - 1) / 2, p) == (p - 1))   break;   }     // Initializing variable x, b and g   let x = pow1(n, (s + 1) / 2, p);   let b = pow1(n, s, p);   let g = pow1(q, s, p);     let r = e;     // keep looping until b   // become 1 or m becomes 0   while (true)   {   let m;   for (m = 0; m < r; m++)   {   if (order(p, b) == -1)   return -1;     // finding m such that b^ (2^m) = 1   if (order(p, b) == Math.pow(2, m))   break;   }   if (m == 0)   return x;     // updating value of x, g and b  // according to algorithm   x = (x * pow1(g, Math.pow(2,   r - m - 1), p)) % p;   g = pow1(g, Math.pow(2, r - m), p);   b = (b * g) % p;     if (b == 1)   return x;   r = m;   }  } // Driver Code  let n = 2;     // p should be prime   let p = 113;     let x = STonelli(n, p);     if (x == -1)  document.write("Modular square" +   "root is not exist\n");   else  document.write("Modular square root of " +  n + " and " + p + " is " +  x + "\n");  </script>