# 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 ` `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. "````); ` `        ``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 "````); ` `        ``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 "````); ` `        ``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 "````); ` `    ``else` `        ``printf``(````"Modular square root of %d and %d is %d "````, ` `                ``n, p, x); ` `} `

## Java

 `// 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 "````);  ` `    ``else` `        ``System.out.println(``"Modular square root of "` `+ ` `                            ``n + ``" and "` `+ p + ``" is "` `+ ` `                            ``x + ````" "````);  ` `}  ` `} ` ` `  `// 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. ”);
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 ”);
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 ”);
return -1;

# expressing p – 1, in terms of s * 2^e,
# where s is odd number
ar = convertx2e(p – 1);
s = ar;
e = ar;

# 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 "); else: print("Modular square root of", n, "and", p, "is", x); # This code is contributed by mits [tabby title="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 "````);  ` `    ``else` `        ``Console.WriteLine(``"Modular square root of"` `+  ` `                        ````"{0} and {1} is {2} "````, n, p, x);  ` `}  ` `} ` ` `  `// This code is contributed by mits `

## PHP

 ` 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. "````);  ` `        ``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 "````);  ` `        ``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 "````);  ` `        ``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 "````);  ` `else` `    ``print``(``"Modular square root of "` `.  ` `          ````"\$n and \$p is \$x "````);  ` `     `  `// This code is contributed by mits ` `?> `

Output :

`Modular square root of 2 and 113 is 62`

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