;;; "determ.scm" Matrix Algebra
;Copyright 2002, 2004 Aubrey Jaffer
;
;Permission to copy this software, to modify it, to redistribute it,
;to distribute modified versions, and to use it for any purpose is
;granted, subject to the following restrictions and understandings.
;
;1. Any copy made of this software must include this copyright notice
;in full.
;
;2. I have made no warranty or representation that the operation of
;this software will be error-free, and I am under no obligation to
;provide any services, by way of maintenance, update, or otherwise.
;
;3. In conjunction with products arising from the use of this
;material, there shall be no use of my name in any advertising,
;promotional, or sales literature without prior written consent in
;each case.
(require 'array)
(require 'multiarg-apply)
;;@code{(require 'determinant)}
;;@ftindex determinant
;;@noindent
;;A Matrix can be either a list of lists (rows) or an array.
;;Unlike linear-algebra texts, this package uses 0-based coordinates.
;;; Internal conversion routines
(define (matrix2array matrix prototype)
(let* ((dim1 (length matrix))
(dim2 (length (car matrix)))
(mat (make-array '#() dim1 dim2)))
(do ((idx 0 (+ 1 idx))
(rows matrix (cdr rows)))
((>= idx dim1) rows)
(do ((jdx 0 (+ 1 jdx))
(row (car rows) (cdr row)))
((>= jdx dim2))
(array-set! mat (car row) idx jdx)))
mat))
(define (matrix2lists matrix)
(let ((dims (array-dimensions matrix)))
(do ((idx (+ -1 (car dims)) (+ -1 idx))
(rows '()
(cons (do ((jdx (+ -1 (cadr dims)) (+ -1 jdx))
(row '() (cons (array-ref matrix idx jdx) row)))
((< jdx 0) row))
rows)))
((< idx 0) rows))))
(define (coerce-like-arg matrix arg)
(cond ((array? arg) (matrix2array matrix arg))
(else matrix)))
;;@body
;;Returns the list-of-lists form of @1.
(define (matrix->lists matrix)
(cond ((array? matrix)
(if (not (eqv? 2 (array-rank matrix)))
(slib:error 'not 'matrix matrix))
(matrix2lists matrix))
((and (pair? matrix) (list? (car matrix))) matrix)
((vector? matrix) (list (vector->list matrix)))
(else (slib:error 'not 'matrix matrix))))
;;@body
;;Returns the array form of @1.
(define (matrix->array matrix)
(cond ((array? matrix)
(if (not (eqv? 2 (array-rank matrix)))
(slib:error 'not 'matrix matrix))
matrix)
((and (pair? matrix) (list? (car matrix)))
(matrix2array matrix '#()))
((vector? matrix) matrix)
(else (slib:error 'not 'matrix matrix))))
(define (matrix:cofactor matrix i j)
(define mat (matrix->lists matrix))
(define (butnth n lst)
(if (<= n 1) (cdr lst) (cons (car lst) (butnth (+ -1 n) (cdr lst)))))
(define (minor matrix i j)
(map (lambda (x) (butnth j x)) (butnth i mat)))
(coerce-like-arg
(* (if (odd? (+ i j)) -1 1) (determinant (minor mat i j)))
matrix))
;;@body
;;@1 must be a square matrix.
;;@0 returns the determinant of @1.
;;
;;@example
;;(require 'determinant)
;;(determinant '((1 2) (3 4))) @result{} -2
;;(determinant '((1 2 3) (4 5 6) (7 8 9))) @result{} 0
;;@end example
(define (determinant matrix)
(define mat (matrix->lists matrix))
(let ((n (length mat)))
(if (eqv? 1 n) (caar mat)
(do ((j n (+ -1 j))
(ans 0 (+ ans (* (list-ref (car mat) (+ -1 j))
(matrix:cofactor mat 1 j)))))
((<= j 0) ans)))))
;;@body
;;Returns a copy of @1 flipped over the diagonal containing the 1,1
;;element.
(define (transpose matrix)
(if (number? matrix)
matrix
(let ((mat (matrix->lists matrix)))
(coerce-like-arg (apply map list mat)
matrix))))
;;@body
;;Returns the element-wise sum of matricies @1 and @2.
(define (matrix:sum m1 m2)
(define mat1 (matrix->lists m1))
(define mat2 (matrix->lists m2))
(coerce-like-arg (map (lambda (row1 row2) (map + row1 row2)) mat1 mat2)
m1))
;;@body
;;Returns the element-wise difference of matricies @1 and @2.
(define (matrix:difference m1 m2)
(define mat1 (matrix->lists m1))
(define mat2 (matrix->lists m2))
(coerce-like-arg (map (lambda (row1 row2) (map - row1 row2)) mat1 mat2)
m1))
(define (matrix:scale m1 scl)
(coerce-like-arg (map (lambda (row1) (map (lambda (x) (* scl x)) row1))
(matrix->lists m1))
m1))
;;@args m1 m2
;;Returns the product of matrices @1 and @2.
;;@args m1 z
;;Returns matrix @var{m1} times scalar @var{z}.
;;@args z m1
;;Returns matrix @var{m1} times scalar @var{z}.
(define (matrix:product m1 m2)
(cond ((number? m1) (matrix:scale m2 m1))
((number? m2) (matrix:scale m1 m2))
(else
(let ((mat1 (matrix->lists m1))
(mat2 (matrix->lists m2)))
(define (dot-product v1 v2) (apply + (map * v1 v2)))
(coerce-like-arg
(map (lambda (arow)
(apply map
(lambda bcol (dot-product bcol arow))
mat2))
mat1)
m1)))))
;;@body
;;@1 must be a square matrix.
;;If @1 is singular, then @0 returns #f; otherwise @0 returns the
;;@code{matrix:product} inverse of @1.
(define (matrix:inverse matrix)
(let* ((mat (matrix->lists matrix))
(det (determinant mat))
(rank (length mat)))
(and (not (zero? det))
(do ((i rank (+ -1 i))
(inv '() (cons
(do ((j rank (+ -1 j))
(row '()
(cons (/ (matrix:cofactor mat j i) det) row)))
((<= j 0) row))
inv)))
((<= i 0)
(coerce-like-arg inv matrix))))))