home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #31 / NN_1992_31.iso / spool / comp / lang / scheme / 2801 < prev    next >
Encoding:
Text File  |  1992-12-21  |  4.6 KB  |  160 lines
  1. Newsgroups: comp.lang.scheme
  2. Path: sparky!uunet!brunix!brunix!mj
  3. From: mj@cs.brown.edu (Mark Johnson)
  4. Subject: Re: Numerical analysis routines, anyone?
  5. Message-ID: <1992Dec21.154653.26697@cs.brown.edu>
  6. Sender: news@cs.brown.edu
  7. Organization: Brown University Department of Computer Science
  8. References: <1gtos5INN66o@master.cs.rose-hulman.edu>
  9. Date: Mon, 21 Dec 1992 15:46:53 GMT
  10. Lines: 148
  11.  
  12. This is certainly not complete, but it is at least a start.
  13. I wrote this package for some statistical language processing
  14. that I'm doing.  The matrix routines all use a generic
  15. arithmetic interface because I hope at some stage to be
  16. working with more complex numerical types (specifically,
  17. a value and its first derivative).
  18.  
  19. Anyway, I hope this is useful,
  20.  
  21. Mark
  22.  
  23. ;;; Matrix manipulation package
  24. ;;;
  25. ;;; requires a DOTIMES macro a la CommonLisp
  26.  
  27. ;;; generic arithmetic interface
  28.  
  29. (define multiply *)
  30. (define add +)
  31. (define subtract -)
  32. (define divide /)
  33. (define zero 0)
  34. (define =zero? zero?)
  35. (define one 1)
  36. (define (inverse x) (/ 1.0 x))
  37.  
  38. (define (sum n fn)
  39.   (define (sum-helper i partial-sum)
  40.     (if (zero? i)
  41.         (add (fn 0) partial-sum)
  42.         (sum-helper (1- i) (add (fn i) partial-sum))))
  43.   (sum-helper (1- n) zero))
  44.     
  45. ;;; generic matrix interface.  This version uses vectors of vectors,
  46. ;;; and uses the generic arithemtic above
  47.  
  48. (define (make-zero-matrix dim)
  49.   (let ((mat (make-vector dim)))
  50.     (dotimes (i dim)
  51.       (vector-set! mat i (make-vector dim zero)))
  52.     mat))
  53.  
  54. (define (make-matrix dim)
  55.   (make-zero-matrix dim))
  56.  
  57. (define (matrix-ref mat i j)
  58.   (vector-ref (vector-ref mat i) j))
  59.  
  60. (define (matrix-set! mat i j value)
  61.   (vector-set! (vector-ref mat i) j value)
  62.   mat)
  63.  
  64. (define (matrix-inc! mat i j inc)
  65.   (matrix-set! mat i j (add (matrix-ref mat i j) inc)))
  66.  
  67. (define (matrix-dec! mat i j inc)
  68.   (matrix-set! mat i j (subtract (matrix-ref mat i j) inc)))
  69.  
  70. (define (matrix-mult! mat i j inc)
  71.   (matrix-set! mat i j (multiply (matrix-ref mat i j) inc)))
  72.  
  73. (define (matrix-div! mat i j inc)
  74.   (matrix-set! mat i j (divide (matrix-ref mat i j)  inc)))
  75.  
  76. (define (matrix-dimension mat)
  77.   (vector-length mat))
  78.  
  79. (define (make-unit-matrix dim)
  80.   (let ((mat (make-zero-matrix dim)))
  81.     (dotimes (i dim)
  82.       (matrix-set! mat i i one))
  83.     mat))
  84.  
  85. ;(define (matrix-swap-elts mat i1 j1 i2 j2)
  86. ;  (let ((temp (matrix-ref mat i1 j1)))
  87. ;    (matrix-set! mat i1 j1 (matrix-ref mat i2 j2))
  88. ;    (matrix-set! mat i2 j2 temp)))
  89.  
  90. (define (matrix-swap-rows! mat row1 row2)
  91.       (let ((temp (vector-ref mat row1)))
  92.      (vector-set! mat row1 (vector-ref mat row2))
  93.      (vector-set! mat row2 temp)))
  94.  
  95. (define (copy-matrix mat)
  96.   (let* ((n (matrix-dimension mat))
  97.          (copy (make-matrix n)))
  98.     (dotimes (i n)
  99.       (dotimes (j n)
  100.         (matrix-set! copy i j (matrix-ref mat i j))))
  101.     copy))
  102.  
  103. ;;; matrix multiplication
  104.  
  105. (define (matrix-multiply mat1 mat2)
  106.   (let* ((n (matrix-dimension mat1))
  107.          (ans (make-zero-matrix n)))
  108.     (dotimes (i n)
  109.       (dotimes (j n)
  110.         (matrix-set! ans i j
  111.                      (sum n (lambda (k)
  112.                               (multiply (matrix-ref mat1 i k)
  113.                                         (matrix-ref mat2 k j)))))))
  114.     ans))
  115.  
  116. (define (matrix-subtract mat1 mat2)
  117.   (let* ((n (matrix-dimension mat1))
  118.          (ans (make-matrix n)))
  119.     (dotimes (i n)
  120.       (dotimes (j n)
  121.         (matrix-set! ans i j (subtract (matrix-ref mat1 i j)
  122.                                        (matrix-ref mat2 i j)))))
  123.     ans))
  124.  
  125. ;;; gauss-jordan does matrix inversion.  Note that this destructively changes 
  126. ;;; the argument matrix mat
  127.  
  128. (define (gauss-jordan mat)
  129.   (let* ((n (matrix-dimension mat))
  130.          (inv (make-unit-matrix n)))
  131.     (dotimes (i n)
  132.          (do ((pi i (1+ pi)))
  133.           ((or (= pi n)
  134.                (not (=zero? (matrix-ref mat pi i))))
  135.            (cond ((= pi n) (error "GAUSS-JORDAN: matrix to invert is singular"))
  136.                  ((not (= i pi)) (matrix-swap-rows! mat i pi)
  137.                                  (matrix-swap-rows! inv i pi)))))
  138.       (let ((pivot (matrix-ref mat i i)))
  139.         (dotimes (j n)
  140.           (if (not (= j i))
  141.             (let ((factor (divide (matrix-ref mat j i) pivot)))
  142.               (matrix-set! mat j i zero)
  143.               (do ((k (1+ i) (1+ k)))
  144.                         ((= k n))
  145.                    (matrix-dec! mat j k (multiply (matrix-ref mat i k) factor)))
  146.               (dotimes (k n)
  147.                 (matrix-dec! inv j k (multiply (matrix-ref inv i k) factor))))))
  148.         (matrix-set! mat i i one)
  149.         (do ((j (1+ i) (1+ j)))
  150.             ((= j n))
  151.           (matrix-div! mat i j pivot))
  152.         (dotimes (j n)
  153.           (matrix-div! inv i j pivot))))
  154.     inv))
  155.                   
  156.  
  157. Mark Johnson
  158. Cognitive Science, Box 1978
  159. Brown University
  160.