(defun random-between (lower upper)
  "Return a random number between LOWER and UPPER, inclusive."
  (+ lower (random (1+ (- upper lower)))))

(defun |in 2^r*d+1 form| (integer)
  "Return the values of r, d where integer = 2^r*d+1."
  (assert (oddp integer))
  ;; Get rid of the + 1 at the end.
  (decf integer)
  (loop for r from 0
        for d = integer then (/ d 2)
        when (oddp d)
          return (values r d)))

(defun expt-mod (b e m)
  "Compute b^e mod m"
  (let ((r 1))
    (setf b (mod b m))
    (loop while (plusp e)
          if (oddp e)
            do (setf r (mod (* r b) m))
          do (setf e (floor e 2)
                   b (mod (expt b 2) m)))
    r))

(defun prime-p (integer &key (rounds 32))
  "Test if an integer is prime, using the Miller-Rabin algorithm.
This has a 4^(-rounds) probability of incorrectly deciding that a composite (non-prime) integer is prime."
  (check-type integer (integer 2 *))
  (when (evenp integer)
    (return-from prime-p (= integer 2)))
  (multiple-value-bind (r d)
      (|in 2^r*d+1 form| integer)
    (loop repeat rounds
          for a = (random-between 2 (- integer 2))
          for x = (expt-mod a d integer)
          unless (or (= x 1) (= x (1- integer)))
            do (block inner-loop
                 (dotimes (i (1- r))
                   (setf x (expt-mod x 2 integer))
                   (when (= x (1- integer))
                     (return-from inner-loop)))
                 (return-from prime-p nil))
          finally (return t))))