diff options
Diffstat (limited to 'determ.scm')
| -rw-r--r-- | determ.scm | 157 |
1 files changed, 145 insertions, 12 deletions
@@ -1,14 +1,147 @@ -;"determ.scm" Determinant +;;; "determ.scm" Matrix Algebra +;Copyright 2002 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. -(define (determinant m) +(require 'array) + +;;@code{(require 'determinant)} +;;@ftindex determinant + +;;@noindent +;;A Matrix can be either a list of lists (rows) or an array. +;;As with linear-algebra texts, this package uses 1-based coordinates. + +;;; Internal conversion routines +(define (matrix2array matrix prototype) + (let* ((shp (list (list 1 (length matrix)) + (list 1 (length (car matrix))))) + (mat (apply create-array '#() shp))) + (do ((idx 1 (+ 1 idx)) + (rows matrix (cdr rows))) + ((> idx (cadar shp)) rows) + (do ((jdx 1 (+ 1 jdx)) + (row (car rows) (cdr row))) + ((> jdx (cadadr shp))) + (array-set! mat (car row) idx jdx))) + mat)) +(define (matrix2lists matrix) + (let ((shp (array-shape matrix))) + (do ((idx (cadar shp) (+ -1 idx)) + (rows '() + (cons (do ((jdx (cadadr shp) (+ -1 jdx)) + (row '() (cons (array-ref matrix idx jdx) row))) + ((< jdx (caadr shp)) row)) + rows))) + ((< idx (caar shp)) 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 (ones-based) 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 (zero? n) (cdr lst) (cons (car lst) (butnth (+ -1 n) (cdr lst))))) - (define (minor m i j) - (map (lambda (x) (butnth j x)) (butnth i m))) - (define (cofactor m i j) - (* (if (odd? (+ i j)) -1 1) (determinant (minor m i j)))) - (define n (length m)) - (if (eqv? 1 n) (caar m) - (do ((j (+ -1 n) (+ -1 j)) - (ans 0 (+ ans (* (list-ref (car m) j) (cofactor m 0 j))))) - ((negative? j) ans)))) + (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 product of matrices @1 and @2. +(define (matrix:product m1 m2) + (define mat1 (matrix->lists m1)) + (define 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 singlar, 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)))))) |