[plt-scheme] Matrix Determinants

So, when I get bored doing math assignments, I write short programs to
double-check my answers. It's somewhat of a bad habit actually, as it
often ends up taking long enough that I don't have enough time to
finish the assignment...

Anyways, during my last assignment I came up with the following for
taking matrix determinants and multiplying matrices. Matrices are
lists of rows, rows are lists of numbers, all rows of a matrix are the
same length. I wasn't really worrying about numerical error, since the
matrices I gave it consisted of exact numbers. I'm pretty happy with
how simple matrix-mult is ((apply map list m) is neat!), but taking
the determinant of a matrix proved messy. I add multiples of the first
row to the rest to zero out the rest of the first column, so I only
need to do ~n recursive calls. It gets messy when I need to deal with
the case that the first row has a zero as the first element, swapping
list elements is tricky!

Can anyone come up with a simple matrix determinant implementation?
(other than the obvious O(n!) one of course)

    Henk



(module math mzscheme

  (require (lib "list.ss" "srfi" "1"))

  (provide det matrix-mult)

  (define (find-index pred l)
    (let loop ((index 0) (l l))
      (cond ((null? l) #f)
            ((pred (car l)) index)
            (else (loop (add1 index) (cdr l))))))

  (define (swap l i1 i2)
    (let loop ((l l) (i1 (min i1 i2)) (i2 (max i1 i2)))
      (if (= 0 i1)
        (cons (list-ref l i2)
              (replace-index (cdr l) (sub1 i2) (car l)))
        (cons (car l)
              (loop (cdr l) (sub1 i1) (sub1 i2))))))

  (define (replace-index l i v)
    (if (= i 0)
      (cons v (cdr l))
      (cons (car l) (replace-index (cdr l) (add1 i) v))))

  ; O(n^3) where n is the height of m
  (define (det m)
    (let ((index (find-index (lambda (r) (not (= 0 (car r)))) m)))
      (if (not index)
        0
        (let ((swap-multiplier (if (= index 0) 1 -1))
              (m-swapped (if (= index 0) m (swap m 0 index))))
          (let ((submatrix (gaussian-submatrix m-swapped)))
            (if (null? submatrix)
              (caar m-swapped)
              (* swap-multiplier
                 (caar m-swapped)
                 (det submatrix))))))))

  ; Zeros entries under the top-left entry, by adding multiples of the first
  ; row. Returns the submatrix ommitting the first row and column. The top left
  ; entry must be nonzero. I couldn't really think of a good name, "gaussian
  ; submatrix" doesn't actually mean anything as far as I know.
  ; O(w*h) where w and h are the width and height of m
  (define (gaussian-submatrix m)
    (let ((first-row (car m)))
      (map (lambda (row)
             (let ((mult (/ (car row)
                            (car first-row))))
               (map (lambda (above current)
                      (- current (* mult above)))
                    (cdr first-row)
                    (cdr row))))
           (cdr m))))

  (define (matrix-mult x y)
    (let ((y-cols (apply map list y)))
      (map (lambda (x-row)
             (map (lambda (y-col) (fold + 0 (map * x-row y-col)))
                  y-cols))
           x)))

  ) ; end module


Examples:
> (det '((1 2 3) (4 5 6) (7 8 9)))
0
> (det '((1)))
1
> (det '((10 20) (30 41)))
-190