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