Реализация схемы Миллера-Рабина непредсказуемым выходом - PullRequest
3 голосов
/ 10 февраля 2010

Я новичок в Схеме. Я попробовал и реализовал вероятностный вариант алгоритма Рабина-Миллера с использованием схемы PLT. Я знаю, что это вероятностно и все такое, но большую часть времени я получаю неправильные результаты. Я реализовал то же самое с помощью C, и это сработало хорошо (никогда не подводил попытки). Я получаю ожидаемый результат при отладке, но при запуске он почти всегда возвращает неверный результат. Я использовал алгоритм из Википедии .

(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))))

Кроме того, правильно ли я программирую на схеме? (Я не уверен насчет части цикла, где я использую call/cc. Я нашел ее на каком-то сайте и с тех пор использую ее.)

Заранее спасибо.

1 Ответ

6 голосов
/ 10 февраля 2010

в общем, вы программируете слишком «императивно»; более элегантный экспорт будет

(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))

, что позволяет избежать использования набора! и просто рекурсивно реализует правила

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

Точно так же вещь rab_mil запрограммирована очень не по схеме. Вот альтернативная реализация. Обратите внимание, что нет «разрыва» циклов и нет вызова / cc; вместо этого прерывание реализовано в виде хвостового рекурсивного вызова, который действительно соответствует 'goto' на схеме:

(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)))
...