Import Upstream version 3a1upstream/3a1

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))))))