Tonelli-Shanks Algorithm
November 23, 2012
One of the operations of modular arithmetic, and an important step in many algorithms of number theory, is finding modular square roots. In the expression x2 ≡ a (mod p), x is a modular square root of a; for instance, 62 is a square root of 2, mod 113, because 62 * 62 = 3844 = 36 * 113 + 2. Like real numbers, modular square roots come in pairs, so -62 ≡ 51 (mod 113) is also a square root of 2.
The usual method for finding modular square roots involves an algorithm developed by Alberto Tonelli in the 1890s; there are several variants, including one due to Daniel Shanks in 1973, one of which we saw in a previous exercise. The particular variant that we will look at today is described at http://www.math.vt.edu/people/brown/doc/sqrts.pdf.
We will assume that the modulus p = s · 2e + 1 is an odd prime, that the number a for which we seek a square root is coprime to p, and that the Legendre symbol (a/p) is 1, which indicates that the square root exists. Find an n such that n(p−1)/2 ≡ −1 (mod p); it won’t take long if you start from n = 2 and increase n by 1 until you find one. Then set x ≡ a(s+1)/2 (mod p), b ≡ as (mod p), g ≡ ns (mod p), and r = e. Next find the smallest non-negative integer m such that b2m ≡ 1 (mod p), which is done by repeated squarings and reductions. If m is zero, you’re finished; report the value of x and halt. Otherwise, set x to x · g2r−m−1 (mod p), set b to b · g2r−m (mod p), set g to g2r−m (mod p), and set r = m, then loop back and compute a new m, continuing until m is zero.
As an example we find the square root of 2 mod 113. First compute p = 113 = 7 · 24 + 1 and n = 3, then x = 16, b = 15, g = 40, and r = 4. The first computation of m is 2, so set x = 62, b = 1, g = 98, and r = 2. Then the second computation of m is 0, so x = 62 is a square root of 2 (mod 113), and so is 113 – 62 = 51.
Your task is to write a function that computes modular square roots according to the algorithm described above. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
[…] Pages: 1 2 […]
A python version.
from fractions import gcd import itertools import ma.primee as PR def ord_r(r, n): """the multiplicative order of n modulo r""" if gcd(r, n) != 1: raise ValueError("n and r are not coprime") k = 3 while 1: if pow(n, k, r) == 1: return k k += 1 def ifrexp(x): e = 0 while x % 2 == 0: x //= 2 e += 1 return x, e def tonelli(a, p): if not PR.is_prime(p): raise ValueError("p must be prime") if gcd(a, p) != 1: raise ValueError("a and p are not coprime") if pow(a, (p - 1) // 2, p) == (-1 % p): raise ValueError("no sqrt possible") s, e = ifrexp(p - 1) for n in itertools.count(2): if pow(n, (p - 1) // 2, p) == (-1 % p): break else: raise ValueError("Not found") # just to be sure x = pow(a, (s + 1) // 2, p) b = pow(a, s, p) g = pow(n, s, p) r = e while 1: for m in xrange(r): if ord_r(p, b) == 2 ** m: break if m == 0: return x, -x % p x = (x * pow(g, 2 ** (r - m - 1), p)) % p g = pow(g, 2 ** (r - m), p) b = (b * g) % p if b == 1: return x, -x % p r = mwhen i try to execute, i get the error as No module named ma.primee
This is great, was looking for an implementation of this instead of just the theory on Wikipedia, thank you!
awesome thanks alot…!! (y)