Miller-Rabin-Test

miller-rabin-test.scm

; Prüft, ob eine Zahl eine Primzahl ist laut Miller-Rabin-Test.
;   http://mathworld.wolfram.com/Rabin-MillerStrongPseudoprimeTest.html
; > Given an odd integer n [1], let n=2^rs+1 with s odd [2,3].
; > Then choose a random integer a with 1<=a<=n-1 [4].
; > If a^s=1 (mod n) [5] or a^(2^js)=-1 (mod n) [6] for some 0<=j<=r-1, then n passes the test.
(define (prim? n)
  (cond
    ((< n 2) #f)
    ((>= n 2152302898747) (error "prim?" "Zu große Zahl" n)) ; '(2, 3) sind ausreichend für alle Zahlen < dieser
    ((memq n '(2 3 5 7 11)) #t)
    ((odd? n) ; [1]
      (letrec*
        (
          (calcR
            (lambda (n)
              (if (odd? n)
                0
                (+ 1 (calcR (/fx n 2)))
              )
            )
          )
          (r (calcR (- n 1))) ; [2]
          (s (/ (- n 1) (expt 2 r))) ; [3]

          ; [6] a^(2js) = -1 mod n
          (testPotenzen
            (lambda (a r)
              (cond ((< r 0) #f)
                    ((= (- n 1) (pow a (* (expt 2 r) s) n)) #t)
                    (else (testPotenzen a (- r 1)))
              )
            )
          )

          ; [5] a^s = 1 mod n
          (testEinzeln
            (lambda (a r)
              ;(print "einz " a " " r)
              (if (= 1 (pow a s n))
                #t
                (testPotenzen a r)
              )
            )
          )

          (test
            (lambda (r werte)
              ;(print r " " werte)
              (if (null? werte)
                #t
                (and (testEinzeln (car werte) r) (test r (cdr werte)))
              )
            )
          )
        )
        (test r
          '(2 3 5 7 11) ; [4]
        )
      )
    )
    (else #f)
  )
)