summaryrefslogtreecommitdiffstats
path: root/determ.scm
diff options
context:
space:
mode:
Diffstat (limited to 'determ.scm')
-rw-r--r--determ.scm157
1 files changed, 145 insertions, 12 deletions
diff --git a/determ.scm b/determ.scm
index 4b53e5f..1078750 100644
--- a/determ.scm
+++ b/determ.scm
@@ -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))))))