[plt-scheme] Implementing Rabin-Miller

From: Paulo Jorge O. C. Matos (pocm at mega.ist.utl.pt)
Date: Tue Oct 15 16:09:36 EDT 2002

Hi all,

I've been sometimes confronted with some messy situations in scheme. Two 
are of special importance. Sometimes while implementing algorithms 
(usually specified in an imperative pseudo-language) there are steps 
like this: If <condition> go to step 3 else return <value>. There is no 
goto blabla or return-from in Scheme. I'm curious how you would 
implement the following simple algorithm in Scheme. I've done it as 
follows (extremely ugly, I agree):
(define isprime?
   ;Rabin-Miller algorithm
   (lambda (p)
     (let* ((b (get-b p))
            (m (/ (- p 1) (expt 2 b)))
            (a (random (if (> p 2147483647) 2147483647 p)))
            (z (modulo (expt a m) p)))
       (printf "z = ~a, a = ~a, m = ~a, b = ~a~n" z a m b)
       (if (or (= z 1)
               (= z (- p 1)))
           (begin
             (printf "Step 1: z = ~a so p may be prime.~n" z)
             #t) ;p MAY BE prime!
           (let loop ((j 0))
             (printf "Step 4, setting j = 0 and z = z^2 mod p.~n")
             (if (and (= z 1)
                      (> j 0))
                 #f ;p is SURELY not prime
                 (begin
                   (set! j (+ j 1))
                   (if (and (< j b)
                            (not (= z (- p 1))))
                       (begin
                         (set! z (modulo (* z z) p))
                         (loop j))
                       (if (= z (- p 1))
                           #t
                           (if (and (= j b)
                                    (not (= z (- p 1))))
                               #f))))))))))

The rabin-miller is as follows:

Choose a random number, p, to test. Calculate b, where b is the number 
of times 2 divides p -  1 (i.e., 2b is the largest power of 2 that 
divides p - 1). Then calculate m, such that p = 1 + 2b *m.

(1)  Choose a random number, a, such that a is less than p.
(2)  Set j = 0 and set z = am mod p.
(3)  If z = 1, or if z = p - 1, then p passes the test and may be prime. 
(4)  If j > 0 and z = 1, then p is not prime.
(5)  Set j = j + 1. If j < b and z ? p - 1, set z = z2 mod p and go back 
to step (4). If z = p - 1, then p passes the test and may be prime.
(6)  If j = b and z ? p - 1, then p is not prime.


Best regards,
-- 
Paulo J. Matos : pocm(_at_)mega.ist.utl.pt
Instituto Superior Tecnico - Lisbon
Software & Computer Engineering - A.I.
  - > http://mega.ist.utl.pt/~pocm
  ---	
	Yes, God had a deadline...
		So, He wrote it all in Lisp!



Posted on the users mailing list.