Tonelli-Shanks Algorithm

November 23, 2012

Here is our version of the function:

```(define (tonelli a p)   (define (expp b e) (expm b e p))   (define (modp n) (modulo n p))   (let* ((p-1 (- p 1)) (p-1/2 (/ p-1 2)))     (if (even? p) (error 'tonelli "must be odd"))     (if (not (prime? p)) (error 'tonelli "must be prime"))     (if (< p 3) (error 'tonelli "must be greater than two"))     (if (< 1 (gcd a p)) (error 'tonelli "must be coprime"))     (if (= (expp a p-1/2) -1)       (error 'tonelli "must be quadratic residue"))     (let loop ((s p-1) (e 0))       (if (even? s) (loop (/ s 2) (+ e 1))         (let loop ((n 2))           (if (not (= (expp n p-1/2) p-1)) (loop (+ n 1))             (let loop1 ((x (expp a (/ (+ s 1) 2)))                         (b (expp a s)) (g (expp n s)) (r e))               (let loop2 ((m 0) (mm 1))                 (if (not (= (expp b mm) 1)) (loop2 (+ m 1) (* mm 2))                   (if (zero? m) (list x (- p x)) ; found result                     (loop1 (modp (* x (expp g (expt 2 (- r m 1)))))                            (modp (* b (expp g (expt 2 (- r m)))))                            (expp g (expt 2 (- r m))) m)))))))))))```

We defined functions `expp` and `modp`, and save the values of p−1 and (p−1)/2, to save redundancies later in the function. After the error returns, the first `loop` computes s and e and the second loop computes n. Then `loop1` is the main processing loop, with variables x, b, g and r, `loop2` calculates m, and then the program exits when m = 0 or loops with new values of x, b, g and r. It’s all much simpler than the description on the previous page. Here’s an example:

```> (tonelli 2 113) (62 51)```

You can check that by squaring: 62 · 62 = 3844 = 34 · 113 + 2, and 51 · 51 = 23 · 113 + 2.

We use `expm` from the Standard Prelude and `prime?` from a previous exercise. You can run the program at http://programmingpraxis.codepad.org/CAPseLXK.

Pages: 1 2

5 Responses to “Tonelli-Shanks Algorithm”

1. […] Pages: 1 2 […]

2. Paul said

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 = m
```
3. sushma said

when i try to execute, i get the error as No module named ma.primee

4. fossjon said

This is great, was looking for an implementation of this instead of just the theory on Wikipedia, thank you!

5. vishal said

awesome thanks alot…!! (y)