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

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++

filter_none

edit
close

play_arrow

link
brightness_4
code

// 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);
}

chevron_right


Java

filter_none

edit
close

play_arrow

link
brightness_4
code

// Java program to implement Shanks
// Tonelli algorithm for finding 
// Modular Square Roots 
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

chevron_right


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 [tabby title="C#"]

filter_none

edit
close

play_arrow

link
brightness_4
code

// 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

chevron_right


PHP

filter_none

edit
close

play_arrow

link
brightness_4
code

<?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
?>

chevron_right



Output :

Modular square root of 2 and 113 is 62

For more detail about above algorithm please visit :
http://cs.indstate.edu/~sgali1/Shanks_Tonelli.pdf

For detail of example (2, 113) see :
http://www.math.vt.edu/people/brown/class_homepages/shanks_tonelli.pdf

This article is contributed by Utkarsh Trivedi. Please write comments if you find anything incorrect, or you want to share more information about the topic discussed above



My Personal Notes arrow_drop_up

Improved By : Mithun Kumar