; In: even nat x
; Out: (values q m), m odd such that x = m*(2^q)
(define (decompose x)
  (let loop ((q 1))
    (let ((m (/ x (expt 2 q))))
      (if (and (integer? m) (odd? m))
          (values q (inexact->exact m))
          (loop (+ q 1))))))

; In: odd nat n, n > 1
; Out: #f if x is composite; #t if inconclusive.
(define (rabin-miller-round n)
  ; b is random in [1, n-1]
  (let ((b (+ 1 (random (- n 1)))))
    (let-values (((q m) (decompose (- n 1)))) ; m odd, n-1=m*(2^q)
      ; n may be prime if b^m == 1 (mod n)
      (if (= (modulo (expt b m) n) 1)
          #t
          ; n may be prime if there exists some i in [0, q-1] such that
          ; n | b^(m*2^i)+1
          (let loop ((i 0))
            (cond ((> i (- q 1)) #f)
                  ((= 0 (modulo (+ (expt b (* m (expt 2 i))) 1) n)) #t)
                  (else (loop (+ i 1)))))))))

; In: nat n, n>0, nat rounds, rounds>0
; Out: #f if in is composite, #t if n is probably prime
;      (n is "probably prime" to certainty at least: 1 - (1/4)^rounds)
(define (rabin-miller-prime? n rounds)
  (cond ((= n 1) #t)
        ((even? n) #f)
        (else (let loop ((rounds rounds))
                (cond ((= rounds 0) #t)
                      ((rabin-miller-round n) (loop (- rounds 1)))
                      (else #f))))))

; In: nat n, n>0
; Out: #f if in is composite, #t if n is almost definitely prime
(define (probably-prime? n)
  (rabin-miller-prime? n 100))