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