(defun primep (n)
  "Determines if integer N is prime."
  (assert (integerp n))
  (if (< n 0) (setq n (* n -1)))
  (cond ((<= n 3) t)
	((eq (mod n 2) 0) nil)
	(t (labels ((recursive-prime-p (n x n-sqrt)
				       (cond ((eq (mod n x) 0) nil)
					     ((>= x n-sqrt) t)
					     (t (recursive-prime-p n
								   (+ x 2)
								   n-sqrt)))))
	     (recursive-prime-p n 3 (floor (sqrt n)))))))

(defun factor (n)
  "Returns a list of the primes and their exponents that compose the
integer N.  For example, (factor 12) results in ((2 2) (3 1))."
  (cond ((eq n 0)
	 nil)
;	((<= (abs n) 3)
;	 (list (list n 1)))
	((eq (mod n 2) 0)
	 (do ((exponent 2 (1+ exponent))
	      (x 4 (* x 2)))
	     ((not (eq (mod n x) 0))
	      (append (list (list 2 (1- exponent)))
		      (factor (/ n (/ x 2)))))))
	(t (do (factors
		(x 3 (+ x 2))
		(max-x (floor (sqrt (abs n)))))
	       ((or (eq n 1)
		    (> x max-x))
		(if (not (eq n 1))
		    (append factors (list (list n 1)))
		  factors))
	     (do ((exponent 1 (1+ exponent))
		  (y x (* y x)))
		 ((not (eq (mod n y) 0))
		  (when (> exponent 1)
		    (setq factors
			  (append factors
				  (list (list x (1- exponent))))
			  n (/ n (/ y x))))))))))

(defun unfactor (factors)
  "Takes a list of ((FACTOR EXPONENT) ...) as returned by the factor
function.  Returns the number those factors compose."
  (reduce '* (mapcar (lambda (l) (apply 'expt l)) factors)))

(defun coprime (&rest numbers)
  "Returns t if all arguments are coprime (relative to all other
arguments) or nil otherwise."
  (eq (apply 'gcd numbers) 1))

(defun euler-phi (n)
  "Returns Euler's phi function for the number N."
  (reduce '* (mapcar (lambda (l)
		       (let ((p (first l))
			     (e (first (rest l))))
			 (* (expt p (1- e)) (1- p)))) (factor n))))

(defun z-mod-n (n)
  "Returns a list of the elements of the set of integer mod N."
  (let (set)
    (dotimes (i n)
      (setq set (append set (list i))))
    set))

(defun multiplicative-group (n)
  "Returns the subset of elements from Z mod N that are invertable.
This represents the multplicative group of Z mod N."
  (remove-if (lambda (a) (not (eq (gcd a n) 1))) (z-mod-n n)))

; Secretly, the argument to possible-orders can be a list, which
; is the result of (factor).
(defun possible-orders (n)
  "Return all possible orders for numbers in the set Z sub N."
  (let (f)
    (if (numberp n)
	(setq f (factor (euler-phi n)))
      (setq f n))
    (if (null f)
	nil
      (let ((products
	     (do* ((factor (first f))
		   (p (first factor))
		   (e (first (rest factor)))
		   (i 0 (1+ i))
		   (n 1 (* n p))
		   set)
		 ((> i e) set)
	       (setq set (append set (list n)))))
	    (further-products (possible-orders (rest f)))
	    results)
	(if (null further-products)
	    products
	  (progn
	    (mapcar (lambda (l)
		      (setq results (append results l)))
		    (mapcar (lambda (factor-a)
			      (mapcar (lambda (factor-b)
					(* factor-a factor-b))
				      products))
			    further-products))
	    results))))))

(defun order-of-integer (a n)
  "Find the order of integer A in the multiplicative group of Z mod N."
  (let ((orders (possible-orders n)))
    (do* ((order (car orders) (car orders))
	  (orders (cdr orders) (cdr orders)))
	((or (null order) (eq (mod (expt a order) n) 1)) order))))