fork download
copy 
  1. ; matrix determinant and inverse
  2.  
  3. (define (make-matrix rows columns . value)
  4.   (do ((m (make-vector rows)) (i 0 (+ i 1)))
  5.       ((= i rows) m)
  6.     (if (null? value)
  7.         (vector-set! m i (make-vector columns))
  8.         (vector-set! m i (make-vector columns (car value))))))
  9.  
  10. (define (matrix-rows x) (vector-length x))
  11.  
  12. (define (matrix-cols x) (vector-length (vector-ref x 0)))
  13.  
  14. (define (matrix-ref m i j) (vector-ref (vector-ref m i) j))
  15.  
  16. (define (matrix-set! m i j x) (vector-set! (vector-ref m i) j x))
  17.  
  18. (define-syntax for
  19.   (syntax-rules ()
  20.     ((for (var first past step) body ...)
  21.       (let ((ge? (if (< first past) >= <=)))
  22.         (do ((var first (+ var step)))
  23.             ((ge? var past))
  24.           body ...)))
  25.     ((for (var first past) body ...)
  26.       (let* ((f first) (p past) (s (if (< first past) 1 -1)))
  27.         (for (var f p s) body ...)))
  28.     ((for (var past) body ...)
  29.       (let* ((p past)) (for (var 0 p) body ...)))))
  30.  
  31. (define (matrix-add a b)
  32.   (let ((ar (matrix-rows a)) (ac (matrix-cols a))
  33.         (br (matrix-rows b)) (bc (matrix-cols b)))
  34.     (if (or (not (= ar br)) (not (= ac bc)))
  35.         (error 'matrix-add "incompatible matrices")
  36.         (let ((c (make-matrix ar ac)))
  37.           (for (i ar)
  38.             (for (j ac)
  39.               (matrix-set! c i j
  40.                 (+ (matrix-ref a i j)
  41.                    (matrix-ref b i j)))))
  42.           c))))
  43.  
  44. (define (matrix-scalar-multiply n a)
  45.   (let* ((ar (matrix-rows a))
  46.          (ac (matrix-cols a))
  47.          (c (make-matrix ar ac)))
  48.     (for (i ar)
  49.       (for (j ac)
  50.         (matrix-set! c i j
  51.           (* n (matrix-ref a i j)))))
  52.     c))
  53.  
  54. (define (matrix-multiply a b)
  55.   (let ((ar (matrix-rows a)) (ac (matrix-cols a))
  56.         (br (matrix-rows b)) (bc (matrix-cols b)))
  57.     (if (not (= ac br))
  58.         (error 'matrix-multiply "incompatible matrices")
  59.         (let ((c (make-matrix ar bc 0)))
  60.           (for (i ar)
  61.             (for (j bc)
  62.               (for (k ac)
  63.                 (matrix-set! c i j
  64.                   (+ (matrix-ref c i j)
  65.                      (* (matrix-ref a i k)
  66.                         (matrix-ref b k j)))))))
  67.           c))))
  68.  
  69. (define (matrix-transpose a)
  70.   (let* ((ar (matrix-rows a))
  71.          (ac (matrix-cols a))
  72.          (c (make-matrix ac ar)))
  73.     (for (i ar)
  74.       (for (j ac)
  75.         (matrix-set! c j i
  76.           (matrix-ref a i j))))
  77.     c))
  78.  
  79. (define (sub-matrix a i j)
  80.   (let ((r (matrix-rows a)) (c (matrix-cols a)))
  81.     (let ((m (make-matrix (- r 1) (- c 1))) (new-i -1))
  82.       (for (old-i c)
  83.         (when (not (= old-i i))
  84.           (set! new-i (+ new-i 1))
  85.           (let ((new-j -1))
  86.             (for (old-j r)
  87.               (when (not (= old-j j))
  88.                 (set! new-j (+ new-j 1))
  89.                 (matrix-set! m new-i new-j
  90.                   (matrix-ref a old-i old-j)))))))
  91.       m)))
  92.  
  93. (define (matrix-determinant a) ; assume a is square
  94.   (let ((n (matrix-rows a)))
  95.     (if (= n 2)
  96.         (- (* (matrix-ref a 0 0) (matrix-ref a 1 1))
  97.            (* (matrix-ref a 1 0) (matrix-ref a 0 1)))
  98.         (let loop ((j 0) (k 1) (d 0))
  99.           (if (= j n) d
  100.             (loop (+ j 1) (* k -1)
  101.                   (+ d (* k (matrix-ref a 0 j)
  102.                           (matrix-determinant
  103.                             (sub-matrix a 0 j))))))))))
  104.  
  105. (define (matrix-cofactors a) ; assume a is square
  106.   (let* ((n (matrix-rows a)) (cof (make-matrix n n)))
  107.     (if (= n 2)
  108.         (for (i n)
  109.           (for (j n)
  110.             (matrix-set! cof i j
  111.               (* (expt -1 (+ i j))
  112.                  (matrix-ref a (- 1 i) (- 1 j))))))
  113.         (for (i n)
  114.           (for (j n)
  115.             (matrix-set! cof i j
  116.               (* (expt -1 (+ i j))
  117.                  (matrix-determinant (sub-matrix a i j)))))))
  118.     cof))
  119.  
  120. (define (matrix-adjugate a)
  121.   (matrix-transpose (matrix-cofactors a)))
  122.  
  123. (define (matrix-inverse a)
  124.   (matrix-scalar-multiply
  125.     (/ (matrix-determinant a))
  126.     (matrix-adjugate a)))
  127.  
  128. (define a '#( #(6 24 1) #(13 16 10) #(20 17 15)))
  129. (display (matrix-multiply a (matrix-inverse a))) (newline)
Success #stdin #stdout 0.07s 43712KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
#(#(1 0 0) #(0 1 0) #(0 0 1))
https://ideone.com/nZGfo3
language:
Scheme (guile 2.2.4)
created:
6 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!