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

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

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

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