home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The World of Computer Software
/
World_Of_Computer_Software-02-387-Vol-3of3.iso
/
j
/
jacal1a0.zip
/
jacal
/
vect.scm
< prev
Wrap
Text File
|
1992-12-23
|
10KB
|
383 lines
;;; JACAL: Symbolic Mathematics System. -*-scheme-*-
;;; Copyright 1989, 1990, 1991, 1992 Aubrey Jaffer.
;;; See the file "COPYING" for terms applying to this program.
(define (bunch_norm x)
(if (null? (cdr x))
(car x)
x))
(define (infinite-list-of . elts)
(let ((lst (copy-list elts)))
(nconc lst lst)))
(define (poly_elim poleqns vars)
(do ((vs vars (cdr vs)) (polys poleqns) (poly #f))
((null? vs) polys)
(do ((var (car vs))
(pl polys (if (null? pl)
(math-error "not enough equations" poleqns vars)
(cdr pl)))
(npl '() (cons (car pl) npl)))
((poly_find-var? (car pl) var)
(set! poly (promote var (car pl)))
(do ((pls (cdr pl) (cdr pls)))
((null? pls) (set! polys npl))
(if (bunch? (car pls)) (math-error "elim bunch?" (car pls)))
(set! npl (cons (poly_resultant poly (car pls) var)
npl))))
(if (bunch? (car pl)) (math-error "elim bunch?" (car pl))))))
;;; This tries to solve the equations no matter what is involved.
;;; It will eliminate variables in vectors of equations.
(define (eliminate eqns vars)
(bunch_norm
(if (some bunch? eqns)
(apply map
(lambda arglist (eliminate arglist vars))
(map (lambda (x)
(if (bunch? x) x (infinite-list-of x)))
eqns))
(poly_elim eqns vars))))
(define (copymatrix m)
(if (matrix? m)
(map copy-list m)
(math-error "Not a matrix: " m)))
(define (row? a)
(and (bunch? a)
(notevery bunch? a)))
(define (matrix? a)
(and (bunch? a)
(bunch? (car a))
(let ((len (length (car a))))
(and (every (lambda (r) (and (bunch? r) (= (length r) len)))
(cdr a))
len))))
(define (matrix . lst)
(if (matrix? lst) lst (matrix lst)))
(define (butnth n lst)
(cond ((null? lst) (math-error "Coordinate out of range: " n " " lst))
((zero? n) (cdr lst))
(else (cons (car lst)
(butnth (+ -1 n) (cdr lst))))))
(define (mtrx_minor m i j)
(butnth (+ -1 i)
(map (lambda (x) (butnth (+ -1 j) x))
m)))
(define (transpose m)
(cond ((matrix? m)
(bunch_norm (apply map list m)))
((bunch? m) (map list m))
(else m)))
(define (elim_test)
(define a (sexp->var 'A))
(define x (sexp->var 'X))
(define y (sexp->var 'Y))
(test (list (list a 0 0 124 81 11 3 45))
poly_elim
(list (list y (list x (list a 0 0 2) (list a 0 1)) 1)
(list y (list x (list a 5 1) 0 -1) 0 1)
(list y (list x (list a -1 3) 5) -1))
(list x y)))
(define (mtrx_square? m)
(if (not (matrix? m))
#f
(let ((rows (length m)))
(every (lambda (r) (= (length r) rows)) m))))
(define (mtrx_genmatrix fun i2 j2 i1 j1)
(do ((i i2 (+ -1 i))
(mat '()
(cons
(do ((j j2 (+ -1 j))
(row '() (cons (app* fun i j) row)))
((< j j1) row))
mat)))
((< i i1)
(bunch_norm mat))))
;;;Create a matrix of indeterminates
(define (nmatrix var . rst)
(apply genmatrix
(lambda subs (apply make-subscripted-var var subs))
rst))
(define (dotproduct a b)
(let ((cols (matrix? a))
(colbs (matrix? b)))
(cond ((or (and cols (bunch? b)) (and colbs (bunch? a)))
(if (not (and (eqv? cols colbs) (eqv? (length a) (length b))))
(math-error "dotproduct of unequal size matrices: " a
" " b))
(reduce (lambda (x y) (app* _@1+@2 x y))
(reduce (lambda (x y) (app* _@1+@2 x y))
(map (lambda (x y) (app* _@1*@2 x y))
a b))))
((and (row? a) (row? b))
(reduce (lambda (x y) (app* _@1+@2 x y))
(map (lambda (x y) (app* _@1*@2 x y))
a b)))
((bunch? a) (map (lambda (x) (app* _@1*@2 x b))
a))
((bunch? b) (map (lambda (x) (app* _@1*@2 a x))
b))
(else (app* _@1*@2 a b)))))
(define (ncmult a b)
(let ((cols (matrix? a))
(colbs (matrix? b)))
(cond ((or (and cols (bunch? b)) (and colbs (bunch? a)))
(if (not colbs) (set! b (list b)))
(if (not (= (or cols (length a)) (length b)))
(math-error "matrix product of unequal size matrices: " a
" " b))
(or cols (set! a (matrix a)))
(bunch_norm
(map (lambda (arow)
(apply map
(lambda bcol
(dotproduct bcol arow))
b))
a)))
(else (deferop 'ncmult a b)))))
(define (mtrx_expt a pow)
(cond ((zero? pow) 1)
((one? pow) a)
((negative? pow) (mtrx_expt (mtrx_inverse a) (- pow)))
(else
(ipow-by-squaring a pow 1 dotproduct))))
(define (crossproduct a b)
(if (and (row? a) (row? b))
(cond ((not (= (length a) (length b)))
(math-error "Cross product of unequal length vectors: " a
" " b))
((= 2 (length a))
(app* _@1-@2
(app* _@1*@2 (car a) (cadr b))
(app* _@1*@2 (car b) (cadr a))))
((= 3 (length a))
(list
(app* _@1-@2
(app* _@1*@2 (cadr a) (caddr b))
(app* _@1*@2 (cadr b) (caddr a)))
(app* _@1-@2
(app* _@1*@2 (caddr a) (car b))
(app* _@1*@2 (caddr b) (car a)))
(app* _@1-@2
(app* _@1*@2 (car a) (cadr b))
(app* _@1*@2 (car b) (cadr a)))))
(else (math-error "Cross product of funny length vectors: " a
" " b)))
(dotproduct a b)))
(define (mtrx_scalarmatrix n x)
(do ((i 1 (+ 1 i))
(rows (list (nconc (make-list (+ -1 n) 0) (list x)))
(cons (append (cdar rows) (list 0)) rows)))
((>= i n) rows)))
(define (mtrx_diagmatrix diag)
(set! diag (reverse diag))
(do ((i (+ -1 (length diag)) (+ -1 i))
(j 0 (+ 1 j))
(diag diag (cdr diag))
(rows '()
(cons (nconc (make-list i 0) (list (car diag)) (make-list j 0))
rows)))
((null? diag) rows)))
;;; Something is subtly worng here.
(define (determinant m)
(if (not (mtrx_square? m)) (eval-error "determinant of non-square matrix" m))
(letrec
((x11 1)
(ds (lambda (m)
(cond ((null? m) 0)
((expr_0? (caar m))
(set! x11 (- x11))
(cons (cdar m) (ds (cdr m))))
(else
(set! x11
(if (negative? x11) (app* _-@1 (caar m)) (caar m)))
(map (lambda (ro)
(if (expr_0? (car ro))
(cdr ro)
(app* _@1-@2*@3
(cdr ro)
(cdar m)
(app* _@1/@2 (car ro) (caar m)))))
(cdr m)))))))
(cond ((null? (cdr m)) (caar m))
((null? (cddr m)) (crossproduct (car m) (cadr m)))
(else (let ((subdet (determinant (ds m))))
(app* _@1*@2 x11 subdet))))))
(define (charpoly m var)
(determinant (app* _@1-@2
m
(mtrx_scalarmatrix (length m) var))))
(define (cofactor
m i j)
(let ((det (determinant (mtrx_minor m i j))))
(if (odd? (+ i j))
(app* _-@1 det)
det)))
(define (coefmatrix eqns vars)
(map
(lambda (eq)
(map (lambda (var)
(let ((c (poly_coeff eq var 1)))
(set! eq (poly_coeff eq var 0))
c))
vars))
eqns))
(define (augcoefmatrix eqns vars)
(map
(lambda (eq)
(nconc
(map (lambda (var)
(let ((value (poly_coeff eq var 1)))
(set! eq (poly_coeff eq var 0))
value))
vars)
(list eq)))
eqns))
(define (echelon m) m)
;;; This is not correct
(define (triangularize m)
(let* ((cols (matrix? m))
(n (length m))
(c (make-list n -1)))
(cond ((null? cols) (set! m (list m))
(set! n 1))
((> cols n) (set! m (copymatrix m)))
(else (set! m (transpose m))
(set! n cols)))
(do ((k 0 (+ 1 k)))
((>= k n) m)
(let* ((j
(do ((j 0 (+ 1 j)))
((>= j n) #f)
(if (and (< (list-ref c j) 0)
(not (poly_0? (list-ref (list-ref m k) j))))
(return j))))
(rowj (and j (list-ref m j))))
(cond (j
(set! rowj (app* _-@1/@2
rowj
(list-ref rowj k)))
(set-car! (list-tail m j) rowj)
(do ((restrows m (cdr restrows)))
((null? restrows))
(if (not (eq? (car restrows) rowj))
(set-car! restrows
(app* _@1*@2+@3
rowj
(list-ref (car restrows) k)
(car restrows)))))
(set-car! (list-tail c j) k))
(else #f))))))
;;; This algorithm taken from:
;;; Knuth, D. E.,
;;; The Art Of Computer Programming, Vol. 2: Seminumerical Algorithms,
;;; Addison Wesley, Reading, MA 1969.
;;; pp 425-426
(define (rank m)
(let* ((cols (matrix? m))
(n (length m))
(c (make-list n -1))
(r 0))
(cond ((null? cols) (set! m (list m))
(set! n 1))
((> cols n) (set! m (copymatrix m)))
(else (set! m (transpose m))
(set! n cols)))
(do ((k 0 (+ 1 k)))
((>= k n) (- n r))
(let* ((j
(do ((j 0 (+ 1 j)))
((>= j n) #f)
(if (and (< (list-ref c j) 0)
(not (poly_0? (list-ref (list-ref m k) j))))
(return j))))
(rowj (and j (list-ref m j))))
(cond (j
(set! rowj (app* _-@1/@2
rowj
(list-ref rowj k)))
(set-car! (list-tail m j) rowj)
(do ((restrows m (cdr restrows)))
((null? restrows))
(if (not (eq? (car restrows) rowj))
(set-car! restrows
(app* _@1*@2+@3
rowj
(list-ref (car restrows) k)
(car restrows)))))
(set-car! (list-tail c j) k))
(else (set! r (+ 1 r))))))))
(define (matrix_test)
(let ((d2 (list (list 2
(list a 1 1)
(list b 0 -5))
(list (list a 0 1)
(list b 0 1)
(list c 0 1)))))
(test (list (list 1 (make-rat 1 2) (make-rat 1 3))
(list (make-rat 1 2) (make-rat 1 3) (make-rat 1 4))
(list (make-rat 1 3) (make-rat 1 4) (make-rat 1 5)))
mtrx_genmatrix
(make-rat 1 (list _@2 (list _@1 -1 1) -1)) 3 3 1 1)
(test d2
augcoefmatrix
(list (list y (list x (list b 0 5) -2) (list a -1 1))
(list y (list x (list c 0 1) (list a 0 1)) (list b 0 1)))
(list x y))
(test (list (list 1
(make-rat (list a 1 1) 2)
(make-rat (list b 0 5) 2))
(list 0
1
(make-rat (list c (list b 0 (list a 0 5)) -2)
(list b (list a 0 -1 -1) 2))))
echelon
d2)
(test (list (list 2
(list a 1 1)
(list b 0 -5))
(list 0
(list b (list a 0 -1 1) 2)
(list c (list b 0 (list a 0 5)) 2)))
triangularize
d2)
(test 2
rank
d2)
(test (list y 10 -7 1)
charpoly
(list (list 3 1)
(list 2 4))
(list y))))
;;; Copyright 1989, 1990, 1991, 1992 Aubrey Jaffer.