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

Posted by programmingpraxis
Filed in Exercises
5 Comments »

5 Responses to “Tonelli-Shanks Algorithm”

  2. Paul said
    November 25, 2012 at 1:10 PM

    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
    October 23, 2013 at 8:46 AM

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

  4. fossjon said
    April 4, 2014 at 7:31 PM

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

  5. vishal said
    June 15, 2014 at 3:29 PM

    awesome thanks alot…!! (y)

Leave a comment