guile-project-euler-solutions/lib/rabin-miller-test.scm

(library (lib rabin-miller-test)
  (export
   rabin-miller-test)

  (import
   (except (rnrs base) let-values)
   (only (guile)
         lambda* λ
         remainder)
   (srfi srfi-1))

  ;; It is sufficient to check numbers with the following prime numbers as witnesses up to:
  ;; n < 3.317.044.064.679.887.385.961.981
  ;; This number is larger than 64bit integer range.
  ;; https://en.wikipedia.org/wiki/Miller–Rabin_primality_test
  (define SIGNIFICANT-PRIMES (list 2 3 5 7 11 13 17 19 23 29 31 37 41))

  (define square
    (λ (x)
      (* x x)))

  (define halve
    (λ (x)
      (/ x 2)))

  (define divides?
    (λ (a b)
      (= (remainder b a) 0)))

  (define even?
    (λ (n)
      (divides? 2 n)))

  (define (next n)
    (cond
     ((< n 2) 2)
     ((even? n) (+ n 1))
     (else (+ n 2))))

  (define one-rabin-miller-check
    (λ (base exp m)
      "This is the actual test for non-trivial square roots of
1 congruent regarding m."
      (cond
       [(= exp 0) 1]
       [(even? exp)
        ;; If the exponent is an even number, we can halve
        ;; it, and the squaring happens outside!  Why this
        ;; is possible? Check the rabin-miller-test.pdf /
        ;; .odt.
        (remainder (square (one-rabin-miller-check base
                                                   (halve exp)
                                                   m))
                   m)]
       [else
        ;; in the else branch we take a factor outside instead:
        (remainder (* base (one-rabin-miller-check base
                                                   (- exp 1)
                                                   m))
                   m)])))

  (define rabin-miller-test-single-significant-prime
    (λ (number-to-check a)
      "This procedure is only for calling
one-rabin-miller-check with the correct parameters."
      (cond
       [(> 1 number-to-check) #f]  ; safeguard against 0 or lower parameter
       [else
        (= (one-rabin-miller-check a
                                   (- number-to-check 1)
                                   number-to-check)
           1)])))

  (define rabin-miller-test
    (λ (potential-prime significant-primes)
      (let iter ([number-to-check potential-prime]
                 [remaining-significant-primes significant-primes])
        (cond
         [(null? remaining-significant-primes) #t]
         [else
          (let ([current-significant-prime
                 (first remaining-significant-primes)])
            (cond
             [(>= current-significant-prime number-to-check) #t]
             [(rabin-miller-test-single-significant-prime
               number-to-check
               current-significant-prime)
              (iter number-to-check
                    (drop remaining-significant-primes 1))]
             [else #f]))])))))


  ;; (define find-primes-limited
  ;;   (λ (min max)
  ;;     (cond
  ;;      [(< min max)
  ;;       (cond
  ;;        [(rabin-miller-test min)
  ;;         ;; test successful, number is prime according to
  ;;         ;; test guaranteed to be true up to:
  ;;         ;; n < 3.317.044.064.679.887.385.961.981
  ;;         (display "prime:")(displayln min)
  ;;         (find-primes-limited (next min) max)]
  ;;        [else
  ;;         (find-primes-limited (next min) max)])]
  ;;      [else
  ;;       (displayln "finished")])))


  ;; (time (find-primes-limited 0 1000000)))