[plt-scheme] Implementing Rabin-Miller
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!