(defun expt-mod (base exponent modulus)
  (declare (optimize (speed 3) (safety 1))
           (type integer base exponent modulus))
  (cond ((= modulus 1) 0)
        ((= base 2) (mod (ash 1 exponent) modulus))
        (t (let ((b (mod base modulus))
                 (e exponent)
                 (acc 1))
             (declare (type integer b e acc))
             (loop while (plusp e) do
               (when (logbitp 0 e)
                 (setf acc (mod (* acc b) modulus)))
               (setf b (mod (* b b) modulus)
                     e (ash e -1)))
             acc))))



(defun miller-rabin-base (n a d s)
  (declare (optimize (speed 3) (safety 1))
	   (type integer n a d s))
  (let ((x (expt-mod a d n)))
    (cond ((or (= x 1) (= x (1- n))) t)
          (t (loop repeat s
		   do (setf x (expt-mod x 2 n))	      
                   when (= x (1- n)) return t
		     finally (return nil))))))

(defun miller-rabin (n &optional (k 10))
  (declare (optimize (speed 3) (safety 1))
           (type integer n))
  (cond ((<= n 1) nil)
        ((<= n 3) t)
        ((evenp n) nil)
        (t (let ((d (1- n))
                 (s 0))
             (loop while (evenp d) do
               (setf d (ash d -1)
                     s (1+ s)))
             (loop repeat k
                   for a = (+ 2 (random (- n 2)))
                   unless (miller-rabin-base n a d s)
                     return nil
                   finally (return t))))))


(defun primep (n &optional (k 10))
  (declare (optimize (speed 3) (safety 1))
	   (type integer n k))
  (when (miller-rabin n k) n))




(defun rhoff (n)
  (declare (optimize (speed 3) (safety 1))
           (type (integer 2 *) n))
  (cond ((evenp n) 2)
        ((primep n) n)
        (t
         (let ((c (1+ (random (1- n))))
               (x 2)
               (y 2)
               (d 1))
           (declare (type integer c x y d))
           (flet ((f (z)
                    (declare (type integer z))
                    (mod (+ (* z z) c) n)))
             (loop while (= d 1) do
               (setf x (f x)
                     y (f (f y))
                     d (gcd (abs (- x y)) n))
		   finally (return (and (< 1 d n) d))))))))

(defun rho (n)
  (declare (optimize (speed 3) (safety 1))
           (type integer n))
  (when (> n 1)
    (let ((pending (make-array 1
                               :element-type 'integer
                               :adjustable t
                               :fill-pointer 1
                               :initial-contents (list n)))
          (results (make-array 0
                               :element-type 'integer
                               :adjustable t
                               :fill-pointer 0)))
      (declare (type vector pending results))
      (loop while (plusp (fill-pointer pending)) do
        (let ((m (vector-pop pending)))
          (declare (type integer m))
          (cond
            ((= m 1)
             nil)
            ((primep m)
             (vector-push-extend m results))
            ((evenp m)
             (vector-push-extend 2 results)
             (vector-push-extend (ash m -1) pending))
            (t
             (loop for d = (rhoff m)
                   until d
                   finally
		      (let ((q (/ m d)))
			(declare (type integer q))
			(vector-push-extend d pending)
			(vector-push-extend q pending)))))))
      (sort results #'<))))


(defun factor (n)
  (declare (optimize (speed 3) (safety 1))
	   (type integer n))
  (rho n))



(defun repunit-value (n &key (digit 1) (base 2))
  (declare (type (integer 2) base)
	   (type (integer 1) digit)
	   (optimize (speed 3) (safety 2)))
  (when (> base digit)
  (loop for i from 0 to (1- n)
	for j = digit then (+ j (* digit (expt base i)))
	finally (return j))))


(defun mersenne-number (p)
  (declare (type (integer 0 *) p)
           (optimize (speed 3) (safety 1)))
  (1- (ash 1 p)))

(defun lucas-lehmer-residue (p)
  (declare (type (integer 3 *) p)
           (optimize (speed 3) (safety 1)))
  (let ((m (mersenne-number p)))
    (loop repeat (- p 1)
          for s of-type integer = 4 then (mod (- (* s s) 2) m)
          finally (return s))))

(defun lucas-lehmer-primep (p)
  (declare (type (integer 3 *) p)
           (optimize (speed 3) (safety 1)))
  (and (primep p)
       (zerop (lucas-lehmer-residue p))))





(defgeneric digits (number base))

(defmethod digits ((number integer) base)
  (declare (type (integer 0 *) number)
           (type (integer 2 *) base)
           (optimize (speed 3) (safety 1) (debug 0)))
  (cond
    ((= base 2)
     (let* ((len (max 1 (integer-length number)))
            (result (make-array len :element-type 'bit)))
       (loop for i from 0 below len
             do (setf (sbit result (- len 1 i))
                      (if (logbitp i number) 1 0)))
       result))
    ((< number base)
     (vector number))
    (t
     (let ((result (make-array 0
                               :element-type `(integer 0 ,(1- base))
                               :adjustable t
                               :fill-pointer 0)))
       (loop while (plusp number)
             do (multiple-value-bind (q r) (floor number base)
                  (vector-push-extend r result)
                  (setf number q)))
       (nreverse result)))))

(defmethod digits ((number vector) base)
  (declare (type (integer 2 *) base)
	   (optimize (speed 3) (safety 1) (debug 0)))
  (loop with result = 0
	for digit across number do
	  (setf result (+ (* result base) digit))
	finally (return result)))

(defmethod digits ((seq sequence) base)
  (digits (coerce seq 'vector) base))








(defun of-n-bits (n)
  (declare (optimize (speed 3) (safety 1) (debug 0))
	   (type integer n))
  (when (>= n 2)
    (loop for i from (1- n) downto 0
	  for sum = (expt 2 i) then (+ sum (* (random 2) (expt 2 i)))
	  finally (return sum))))

(defun prime-of-n-bits (n)
  (declare (optimize (speed 3) (safety 1) (debug 0))
	   (type integer n))
  (when (>= n 2)
    (loop for i = (of-n-bits n)
	  when (primep i) return i)))