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