; matrix determinant and inverse
(define (make-matrix rows columns . value)
(do ((m (make-vector rows)) (i 0 (+ i 1)))
((= i rows) m)
(if (null? value)
(vector-set! m i (make-vector columns))
(vector-set! m i (make-vector columns (car value))))))
(define (matrix-rows x) (vector-length x))
(define (matrix-cols x) (vector-length (vector-ref x 0)))
(define (matrix-ref m i j) (vector-ref (vector-ref m i) j))
(define (matrix-set! m i j x) (vector-set! (vector-ref m i) j x))
(define-syntax for
(syntax-rules ()
((for (var first past step) body ...)
(let ((ge? (if (< first past) >= <=)))
(do ((var first (+ var step)))
((ge? var past))
body ...)))
((for (var first past) body ...)
(let* ((f first) (p past) (s (if (< first past) 1 -1)))
(for (var f p s) body ...)))
((for (var past) body ...)
(let* ((p past)) (for (var 0 p) body ...)))))
(define (matrix-add a b)
(let ((ar (matrix-rows a)) (ac (matrix-cols a))
(br (matrix-rows b)) (bc (matrix-cols b)))
(if (or (not (= ar br)) (not (= ac bc)))
(error 'matrix-add "incompatible matrices")
(let ((c (make-matrix ar ac)))
(for (i ar)
(for (j ac)
(matrix-set! c i j
(+ (matrix-ref a i j)
(matrix-ref b i j)))))
c))))
(define (matrix-scalar-multiply n a)
(let* ((ar (matrix-rows a))
(ac (matrix-cols a))
(c (make-matrix ar ac)))
(for (i ar)
(for (j ac)
(matrix-set! c i j
(* n (matrix-ref a i j)))))
c))
(define (matrix-multiply a b)
(let ((ar (matrix-rows a)) (ac (matrix-cols a))
(br (matrix-rows b)) (bc (matrix-cols b)))
(if (not (= ac br))
(error 'matrix-multiply "incompatible matrices")
(let ((c (make-matrix ar bc 0)))
(for (i ar)
(for (j bc)
(for (k ac)
(matrix-set! c i j
(+ (matrix-ref c i j)
(* (matrix-ref a i k)
(matrix-ref b k j)))))))
c))))
(define (matrix-transpose a)
(let* ((ar (matrix-rows a))
(ac (matrix-cols a))
(c (make-matrix ac ar)))
(for (i ar)
(for (j ac)
(matrix-set! c j i
(matrix-ref a i j))))
c))
(define (sub-matrix a i j)
(let ((r (matrix-rows a)) (c (matrix-cols a)))
(let ((m (make-matrix (- r 1) (- c 1))) (new-i -1))
(for (old-i c)
(when (not (= old-i i))
(set! new-i (+ new-i 1))
(let ((new-j -1))
(for (old-j r)
(when (not (= old-j j))
(set! new-j (+ new-j 1))
(matrix-set! m new-i new-j
(matrix-ref a old-i old-j)))))))
m)))
(define (matrix-determinant a) ; assume a is square
(let ((n (matrix-rows a)))
(if (= n 2)
(- (* (matrix-ref a 0 0) (matrix-ref a 1 1))
(* (matrix-ref a 1 0) (matrix-ref a 0 1)))
(let loop ((j 0) (k 1) (d 0))
(if (= j n) d
(loop (+ j 1) (* k -1)
(+ d (* k (matrix-ref a 0 j)
(matrix-determinant
(sub-matrix a 0 j))))))))))
(define (matrix-cofactors a) ; assume a is square
(let* ((n (matrix-rows a)) (cof (make-matrix n n)))
(if (= n 2)
(for (i n)
(for (j n)
(matrix-set! cof i j
(* (expt -1 (+ i j))
(matrix-ref a (- 1 i) (- 1 j))))))
(for (i n)
(for (j n)
(matrix-set! cof i j
(* (expt -1 (+ i j))
(matrix-determinant (sub-matrix a i j)))))))
cof))
(define (matrix-adjugate a)
(matrix-transpose (matrix-cofactors a)))
(define (matrix-inverse a)
(matrix-scalar-multiply
(/ (matrix-determinant a))
(matrix-adjugate a)))
(define a '#( #(6 24 1) #(13 16 10) #(20 17 15)))
(display (matrix-multiply a (matrix-inverse a))) (newline)