Miller-Rabin scheme implementation unpredictable output

Question

Greetings everyone, I am new to scheme. I have tried and implemented probabilistic variant of rabin-miller algorithm using plt scheme. I know it is probabilistic and all, but I am getting the wrong results most of the time. I have implemented the same thing using C, and it worked well(never failed a try). I get the expected output while debugging :/ but when I run, it almost always returns with an incorrect result. I used the algorithm from http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Algorithm_and_running_time

(define expmod( lambda(b e m)
                 ;(define result 1)
                 (define r 1)
                 (let loop()
                   (if (bitwise-and e 1)
                       (set! r (remainder (* r b) m)))
                   (set! e (arithmetic-shift e -1))
                   (set! b (remainder (* b b) m))
                   (if (> e 0)
                       (loop)))r))



(define rab_mil( lambda(n k)
                  (call/cc (lambda(breakout)
                  (define s 0)
                  (define d 0)
                  (define a 0)
                  (define n1 (- n 1))
                  (define x 0)          
                  (let loop((count 0))
                    (if (=(remainder n1 2) 0)
                        (begin
                          (set! count (+ count 1))
                          (set! s count)
                          (set! n1 (/ n1 2))
                          (loop count))
                        (set! d n1)))
                  (let loop((count k))
                    (set! a (random (- n 3)))
                    (set! a (+ a 2))
                    (set! x (expmod a d n))
                    (set! count (- count 1))

                    (if (or (= x 1) (= x (- n 1)))
                        (begin
                          (if (> count 0)(loop count))))
                    (let innerloop((r 0))
                      (set! r (+ r 1))
                      (if (< r (- s 1)) (innerloop  r))
                      (set! x (expmod x 2 n))
                      (if (= x 1)
                          (begin
                          (breakout #f)))
                      (if (= x (- n 1)) 
                          (if (> count 0)(loop count)))
                      )
                    (if (= x (- s 1)) 
                        (breakout #f))(if (> count 0) (loop count)))#t))))

Any scheme programmers out there who could help? Also, Am I programming the right way in scheme?(I am not sure about the breaking out of loop part where I use call/cc. I found it on some site and been using it ever since) Thanks in advance.

Answer 1

Hi,

in general you are programming in a too "imperative" fashion; a more elegant expmod would be

(define (expmod b e m)
  (define (emod b e)
    (case ((= e 1) (remainder b m))
          ((= (remainder e 2) 1)
           (remainder (* b (emod b (- e 1))) m)
          (else (emod (remainder (* b b) m) (/ e 2)))))))
  (emod b e))

which avoids the use of set! and just implements recursively the rules

b^1 == b (mod m)     
b^k == b b^(k-1) (mod m) [k odd]
b^(2k) == (b^2)^k (mod m)

Similarly the rab_mil thing is programmed in a very non-scheme fashion. Here's an alternative implementation. Note that there is no 'breaking' of the loops and no call/cc; instead the breaking out is implemented as a tail-recursive call which really corresponds to 'goto' in Scheme:

(define (rab_mil n k)
  ;; calculate the number 2 appears as factor of 'n'
  (define (twos-powers n)
     (if (= (remainder n 2) 0)
         (+ 1 (twos-powers (/ n 2)))
         0))
  ;; factor n to 2^s * d where d is odd:
  (let* ((s (twos-powers n 0))
         (d (/ n (expt 2 s))))
    ;; outer loop
    (define (loop k)
      (define (next) (loop (- k 1)))
      (if (= k 0) 'probably-prime
          (let* ((a (+ 2 (random (- n 2))))
                 (x (expmod a d n)))
            (if (or (= x 1) (= x (- n 1)))
                (next)
                (inner x next))))))
    ;; inner loop
    (define (inner x next)
      (define (i r x)
        (if (= r s) (next)
            (let ((x (expmod x 2 n)))
              (case ((= x 1) 'composite)
                    ((= x (- n 1)) (next))
                    (else (i (+ 1 r) x))))
      (i 1 x))
    ;; run the algorithm
    (loop k)))