An equivalent, but slightly more redundant version of this algorithm was developed by Alberto Tonelli in 1891. The version discussed here was developed independently by Daniel Shanks in 1973, who explained:
(Note: All are taken to mean , unless indicated otherwise).
Inputs: p, an odd prime. n, an integer which is a quadratic residue (mod p), meaning that the Legendre symbol .
Outputs: R, an integer satisfying .
- Factor out powers of 2 from p − 1, defining Q and S as: with Q odd. Note that if , i.e. , then solutions are given directly by .
- Select a z such that the Legendre symbol (that is, z should be a quadratic non-residue modulo p), and set .
- If , return R.
- Otherwise, find the lowest i, , such that ; e.g. via repeated squaring.
- Let , and set and .
Once you have solved the congruence with R the second solution is p − R.
Solving the congruence . It is clear that is odd, and since , 10 is a quadratic residue (by Euler's criterion).
- Step 1: Observe so , .
- Step 2: Take as the quadratic nonresidue (2 is a quadratic nonresidue since (again, Euler's criterion)). Set
- Step 3:
- Step 4: Now we start the loop: so ; i.e.
- Let , so .
- Set . Set , and
- We restart the loop, and since we are done, returning
Indeed, observe that and naturally also . So the algorithm yields two solutions to our congruence.
First write . Now write and , observing that . This latter congruence will be true after every iteration of the algorithm's main loop. If at any point, then and the algorithm terminates with .
If , then consider , a quadratic non-residue of . Let . Then and , which shows that the order of is .
Similarly we have , so the order of divides . Suppose the order of is . Since is a square modulo , is also a square, and hence .
Now we set and with this , and . As before, holds; however with this construction both and have order . This implies that has order with .
If then , and the algorithm stops, returning . Else, we restart the loop with analogous definitions of , , and until we arrive at an that equals 0. Since the sequence of S is strictly decreasing the algorithm terminates.
Speed of the algorithm
The Tonelli–Shanks algorithm requires (on average over all possible input (quadratic residues and quadratic nonresidues))
modular multiplications, where is the number of digits in the binary representation of and is the number of ones in the binary representation of . If the required quadratic nonresidue is to be found by checking if a randomly taken number is a quadratic nonresidue, it requires (on average) computations of the Legendre symbol. The average of two computations of the Legendre symbol are explained as follows: is a quadratic residue with chance , which is smaller than but , so we will on average need to check if a is a quadratic residue two times.
This shows essentially that the Tonelli–Shanks algorithm works very well if the modulus is random, that is, if is not particularly large with respect to the number of digits in the binary representation of . As written above, Cipolla's algorithm works better than Tonelli–Shanks if (and only if) . However, if one instead uses Sutherland's algorithm to perform the discrete logarithm computation in the 2-Sylow subgroup of , one may replace with an expression that is asymptotically bounded by . Explicitly, one computes such that and then satisfies (note that is a multiple of 2 because is a quadratic residue).
The algorithm requires us to find a quadratic nonresidue . There is no known deterministic algorithm that runs in polynomial time for finding such a . However, if the generalized Riemann hypothesis is true, there exists a quadratic nonresidue , making it possible to check every up to that limit and find a suitable within polynomial time. Keep in mind, however, that this is a worst-case scenario; in general, is found in on average 2 trials as stated above.
The Tonelli–Shanks algorithm can (naturally) be used for any process in which square roots modulo a prime are necessary. For example, it can be used for finding points on elliptic curves. It is also useful for the computations in the Rabin cryptosystem.
If many square-roots must be done in the same cyclic group and S is not too large, a table of square-roots of the elements of 2-power order can be prepared in advance and the algorithm simplified and sped up as follows.
- Factor out powers of 2 from p − 1, defining Q and S as: with Q odd.
- Find from the table such that and set
- return R.
- Oded Goldreich, Computational complexity: a conceptual perspective, Cambridge University Press, 2008, p. 588.
- Daniel Shanks. Five Number-theoretic Algorithms. Proceedings of the Second Manitoba Conference on Numerical Mathematics. Pp. 51–70. 1973.
- Gonzalo Tornaria - Square roots modulo p, page 2 http://www.springerlink.com/content/xgxe68edy03la96p/fulltext.pdf
- Sutherland, Andrew V. (2011), "Structure computation and discrete logarithms in finite abelian p-groups", Mathematics of Computation 80: 477–500, doi:10.1090/s0025-5718-10-02356-2
- Bach, Eric (1990), "Explicit bounds for primality testing and related problems", Mathematics of Computation 55 (191): 355–380, doi:10.2307/2008811, JSTOR 2008811
- Adleman, L. M., K. Manders, and G. Miller: 1977, `On taking roots in finite fields'. In: 18th IEEE Symposium on Foundations of Computer Science. pp. 175-177
- Niven, Ivan; Herbert S. Zuckerman, Hugh L. Montgomery (1991). An Introduction to the Theory of Numbers (5th ed.). Wiley. ISBN 0-471-62546-9.
Pages 110–115 describe the algorithm and explain the group theory behind it.
- Daniel Shanks. Five Number Theoretic Algorithms. Proceedings of the Second Manitoba Conference on Numerical Mathematics. Pp. 51–70. 1973.
- Alberto Tonelli, Bemerkung über die Auflösung quadratischer Congruenzen. Nachrichten von der Königlichen Gesellschaft der Wissenschaften und der Georg-Augusts-Universität zu Göttingen. Pp. 344–346. 1891. 
- Gagan Tara Nanda - Mathematics 115: The RESSOL Algorithm 
- Implementation in C# http://shankstonelli.blogspot.com/2010/12/shanks-tonelli-algorithm-in-c.html
- Implementation in Python http://eli.thegreenplace.net/2009/03/07/computing-modular-square-roots-in-python