Lisp code for the Frobenius problem
The following code was used in support of a class
assignment for CS523, March 2009.
It goes along with the paper
Calculating Frobenius numbers
with Boolean Toeplitz matrix multiplication
;;
;; Code to accompany the paper
;; "Calculating Frobenius Numbers with Boolean Toeplitz
;; Matrix Multiplication"
;;
;; Chris Bogart 3/16/2009 CS523
;;
;; First, the final algorithm (described as "Toeplitz shortcut" in
;; the paper). Other methods, for comparison, are further
;; down in the file.
;;
;; Structure "toemat" holds a toeplitz matrix
(defstruct toemat size contents)
;; tmat-ref checks a value in a Toeplitz matrix,
;; given an index between -n and n,
;; indexed as follows:
;;
;; a_0 a_{-1} a_{-2} ... a_{-n}
;; a_1 a_0 a_{-1} ... a_{-n+1}
;; a_2 a_1 a_0 ...
;; a_3
;; ...
;; a_n
;;
(defun tmat-ref (tm i) (aref (toemat-contents tm)
(+ i (toemat-size tm) )))
;; tmat-set sets a value in a Toeplitz matrix.
(defun tmat-set (tm i v)
(setf (aref (toemat-contents tm)
(+ i (toemat-size tm) )) v))
;; Create a toeplitz matrix
(defun empty-tmat (n)
(make-toemat :size n
:contents (make-array (+ (* 2 n) 1)
:initial-element 0)))
;; Make a minimal Frobenius matrix
;; given a list of coin denominations
(defun make-tmat (coins)
(let ((tm (empty-tmat (apply #'max coins))))
(dolist (c coins) (tmat-set tm (- c 1) 1))
(tmat-set tm -1 1)
tm))
;; Print a toeplitz matrix
;; in full square matrix form
;;
(defun print-tmat (tm)
(let ((n (toemat-size tm)))
(loop for j from 0 to (- n 1) do
(loop for i from 0 to (- n 1) do
(format t "~d " (tmat-ref tm (- j i))))
(format t "~%"))))
;; Multiply two Toeplitz matrices. Danger: this
;; assumes the result will be another Toeplitz matrix.
;; This is not true in general; only use this function
;; when you know the result will be another Boolean Toeplitz
;; matrix. I suspect it is a valid assumption when
;; everything above the superdiagonal is zero.
;;
(defun tmatmult (a b)
(let* ((n (toemat-size a))
(c (empty-tmat n)))
(loop for r from 0 to n do
(loop for i from 0 to (- n 1) until (> (tmat-ref c (- r)) 0) do
(if (and (> (tmat-ref a (- i)) 0) (> (tmat-ref b (- i r)) 0))
(tmat-set c (- r) 1))))
(loop for r from 1 to n do
(loop for i from 0 to (- n 1) until (> (tmat-ref c r) 0) do
(if (and (> (tmat-ref a (- r i)) 0) (> (tmat-ref b i) 0))
(tmat-set c r 1))))
c))
;; Check: are all entries in this toeplitz matrix
;; nonzero?
(defun tpos-matrix (tmat)
(let ((n (toemat-size tmat))
(found 't))
(loop for i from (- n) to n do
(if (= 0 (tmat-ref tmat i))
(setf found '())))
found))
;; Find the index of primitivity of a Toeplitz matrix.
;; Again, this assumes that powers of a toeplitz matrix
;; are toeplitz matrices, which is not true in general,
;; but works for this application
;;
(defun tfind-ix-prim (mx)
(let* ((n (toemat-size mx))
(k '1)
(k1 '1)
(ans '())
(done '())
(ps (list mx)))
(loop while (and (not done)
(< k (+ (expt (- n 1) 2) 1))) do
(setf ps (cons (tsquare (car ps)) ps))
(setf k (* 2 k))
(if (tpos-matrix (car ps)) (setf done 't)))
(if (not done) -1
;; else
(progn
(setf ps (cdr ps))
(setf ans (pop ps))
(setf k (/ k 2))
(setf k1 (/ k 2))
(loop while ps do
(let ((ans2 (tmatmult (car ps) ans)))
(if (not (tpos-matrix ans2))
(progn
(setf ans ans2)
(setf k (+ k k1)) ))
(pop ps)
(setf k1 (/ k1 2))))
(+ 1 k)))))
;; Square a toeplitz matrix, with the same assumptions.
;;
(defun tsquare (a) (tmatmult a a))
;; Take an exponent of a toeplitz matrix; same
;; assumptions
(defun tpower (a k)
(if (= k 1)
a
(if (evenp k)
(tsquare (tpower a (/ k 2)))
(tmatmult a (tpower a (- k 1))))))
;; Find the frobenius number from a list of coin denominations.
;;
(defun tfrob (cs)
(- (tfind-ix-prim (make-tmat cs)) (apply #'max cs)))
;;
;; =========================================
;;
;; Next are similar functions which compute the same
;; results using fully elaborated matrices. This has
;; neither the dangerous assumptions nor the comuptational
;; speed of the Toeplitz stuff above.
;; Find the frobenius number from a list of coin denoms
(defun frob (cs)
(- (find-ix-prim (make-min-frob cs)) (apply #'max cs)))
;; Turn a list of coin denoms into a minimal Frobenius
;; matrix.
(defun make-min-frob (cs)
(let* ((n (apply #'max cs))
(lcs (cons 0 cs))
(r (make-array (list n n) :initial-element 0)))
(loop for j from 0 to (- n 1) do
(loop for i from 0 to (- n 1) do
(if (member (- j i -1) lcs)
(setf (aref r i j) 1))))
r))
;; Print a matrix (represented as a 2-d array)
;;
(defun p (mx)
(let ((n (array-dimension mx 0)))
(loop for j from 0 to (- n 1) do
(loop for i from 0 to (- n 1) do
(format t " ~d" (aref mx i j)))
(format t "~%"))))
;; Query the size of a matrix
(defun msize (mx) (array-dimension mx 0))
;; Are all entries in this matrix nonzero?
(defun pos-matrix (mx)
(let ((n (msize mx)) (found 't))
(loop for i from 0 to (- n 1) do
(loop for j from 0 to (- n 1) do
(if (= 0 (aref mx i j)) (setf found '()))))
found))
;; Semi-Ring Matrix Multiplication
;; In each dot product step, STOP as soon as you have
;; a non-zero product. All we care about is if they
;; are non-zero, not how large they are.
;; Note that all we care about is non-zeroness;
;; multiplying 54 * 17 just gives us 1 here
;;
(defun matmult (a b)
(let* ((n (array-dimension a 0))
(c (make-array (list n n) :initial-element 0)))
(loop for i from 0 to (- n 1) do
(loop for k from 0 to (- n 1) do
(loop for j from 0 to (- n 1) until (> 0 (aref c i k)) do
(if (> (* (aref a i j) (aref b j k)) 0) (setf (aref c i k) 1)))))
c))
;; Real matrix multiplication; this is not currently
;; used, but you can swap it in for verification
(defun matmult-traditional (a b)
(let* ((n (array-dimension a 0))
(c (make-array (list n n) :initial-element 0)))
(loop for i from 0 to (- n 1) do
(loop for j from 0 to (- n 1) do
(loop for k from 0 to (- n 1) do
(incf (aref c i k)
(* (aref a i j) (aref b j k))))))
c))
;; Square a matrix
(defun square (a) (matmult a a))
;; Take a power of a matrix. if it's an
;; even power, divide and conquer. If it's an
;; odd power, subtract and conquer.
(defun mpower (mx k)
(if (= k 1)
mx
(if (evenp k)
(square (mpower mx (/ k 2)))
(matmult mx (mpower mx (- k 1))))))
;; Find index of primitivity of a boolean
;; matrix. First, square it up to (n-1)^2+1,
;; or until a positive matrix is found.
;; Then, multiply back in old saved squares until
;; the exact exponent is found where the
;; product becomes all positive.
;;
(defun find-ix-prim (mx)
(let* ((n (array-dimension mx 0))
(k '1)
(k1 '1)
(ans '())
(done '())
(ps (list mx)))
(loop while (and (not done)
(< k (+ (expt (- n 1) 2) 1))) do
(setf ps (cons (square (car ps)) ps))
(setf k (* 2 k))
(if (pos-matrix (car ps)) (setf done 't)))
(if (not done) -1
;; else
(progn
(setf ps (cdr ps))
(setf ans (pop ps))
(setf k (/ k 2))
(setf k1 (/ k 2))
(loop while ps do
(let ((ans2 (matmult (car ps) ans)))
(if (not (pos-matrix ans2))
(progn
(setf ans ans2)
(setf k (+ k k1)) ))
(pop ps)
(setf k1 (/ k1 2))))
(+ 1 k)))))
;;;
;;; ==============================
;;;
;;; The following are utility functions for
;;; experimenting with frobenius, primitivity, and
;;; toeplitzity; they're not needed by the stuff above.
;;;
;; return all lists of k integers less than or equal to max
;;
(defun allcombos (max k)
(if (= k 0)
'(())
;;else
(if (= 1 max)
(if (= 1 k) '((1)) '())
;; else
(append
(mapcar (lambda (l) (cons max l))
(allcombos (- max 1) (- k 1)))
(allcombos (- max 1) k)))))
;; All combinations of k integers less than or equal
;; to max that do not include 1.
(defun allnon1s (max k) (mapcar (lambda (x) (cons max x))
(remove-if (lambda (x) (member 1 x))
(allcombos (- max 1) (- k 1)))))
;;
;; Run Frobenius calculation on all possible problems
;; where n=k and a_n=max and there is no penny
;;
(defun rRange (max k) (mapcar #'frob (allnon1s max k)))
;; Is a full matrix a toeplitz matrix?
;;
(defun toeplitzp (mx)
(let* ((n (array-dimension mx 0))
(toe 't))
(loop for i from 1 to (- n 1) do
(loop for j from 1 to (- n 1) do
(if (not (equal (> (aref mx i j) 0) (> (aref mx (- i 1) (- j 1)) 0))) (setf toe '()))))
toe))
;; Detect all cycles in a matrix
;;
(defun cycles (mx)
(let* ((n (array-dimension mx 0))
(cyc '())
(pow mx))
(loop for k from 1 to n do
(if (> (aref pow 1 1) 0) (push k cyc))
(setf pow (matmult pow mx)))
cyc))
;;
;; Produce a "coin matrix" as described in the paper.
;;
(defun coins (cs)
(let* ((n (apply #'max cs))
(r (make-array (list n n) :initial-element 0)))
(loop for i from 0 to (- n 2) do
(setf (aref r (+ i 1) i) 1))
(dolist (c cs) (setf (aref r 0 (- c 1)) 1))
r))
;; Create a full nxn Toeplitz matrix in order to test
;; out theories about powers of them. cs should be a list
;; of nonzero indices, with 0 in the corner, negs along the top,
;; and pos running down. n tells the overall size.
(defun full-toeplitz (cs n)
(let* ((r (make-array (list n n) :initial-element 0)))
(loop for j from 0 to (- n 1) do
(loop for i from 0 to (- n 1) do
(if (member (- j i) cs)
(setf (aref r i j) 1))))
r))
;; Check and see if a toeplitz matrix defined by cs, of size
;; n, can be raised from powers from 1 to k and still remain
;; Toeplitz
(defun powerable-p (r k)
(let* ((up r)
(success 't))
(loop for i from 1 to k do
(if (not (toeplitzp up)) (progn (setf success '()))
) ;; (format t "OK ~d~%" i))
(setf up (matmult r up)))
success))