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.

About these ads

Pages: 1 2

4 Responses to “Tonelli-Shanks Algorithm”

  1. 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    
    
  2. sushma said

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

  3. fossjon said

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 576 other followers

%d bloggers like this: