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.
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)