; matrix operations

(define (make-matrix rows columns . value)
  (do ((m (make-vector rows)) (i 0 (+ i 1)))
      ((= i rows) m)
    (if (null? value)
        (vector-set! m i (make-vector columns))
        (vector-set! m i (make-vector columns (car value))))))

(define (matrix-rows x) (vector-length x))

(define (matrix-cols x) (vector-length (vector-ref x 0)))

(define (matrix-ref m i j) (vector-ref (vector-ref m i) j))

(define (matrix-set! m i j x) (vector-set! (vector-ref m i) j x))

(define-syntax for
  (syntax-rules ()
    ((for (var first past step) body ...)
      (let ((ge? (if (< first past) >= <=)))
        (do ((var first (+ var step)))
            ((ge? var past))
          body ...)))
    ((for (var first past) body ...)
      (let* ((f first) (p past) (s (if (< first past) 1 -1)))
        (for (var f p s) body ...)))
    ((for (var past) body ...)
      (let* ((p past)) (for (var 0 p) body ...)))))

(define (matrix-add a b)
  (let ((ar (matrix-rows a)) (ac (matrix-cols a))
        (br (matrix-rows b)) (bc (matrix-cols b)))
    (if (or (not (= ar br)) (not (= ac bc)))
        (error 'matrix-add "incompatible matrices")
        (let ((c (make-matrix ar ac)))
          (for (i ar)
            (for (j ac)
              (matrix-set! c i j
                (+ (matrix-ref a i j)
                   (matrix-ref b i j)))))
          c))))

(define (matrix-scalar-multiply n a)
  (let* ((ar (matrix-rows a))
         (ac (matrix-cols a))
         (c (make-matrix ar ac)))
    (for (i ar)
      (for (j ac)
        (matrix-set! c i j
          (* n (matrix-ref a i j)))))
    c))


(define (matrix-multiply a b)
  (let ((ar (matrix-rows a)) (ac (matrix-cols a))
        (br (matrix-rows b)) (bc (matrix-cols b)))
    (if (not (= ac br))
        (error 'matrix-multiply "incompatible matrices")
        (let ((c (make-matrix ar bc 0)))
          (for (i ar)
            (for (j bc)
              (for (k ac)
                (matrix-set! c i j
                  (+ (matrix-ref c i j)
                     (* (matrix-ref a i k)
                        (matrix-ref b k j)))))))
          c))))

(define (matrix-transpose a)
  (let* ((ar (matrix-rows a))
         (ac (matrix-cols a))
         (c (make-matrix ac ar)))
    (for (i ar)
      (for (j ac)
        (matrix-set! c j i
          (matrix-ref a i j))))
    c))

(define a #(#(1)
            #(2)
            #(3)))

(define b #(#(1 2 3)
            #(4 5 6)))

(define c #(#(2 3 4)
            #(3 4 5)))

(define d #(#(1 2 3 4)
            #(2 3 4 5)
            #(3 4 5 6)))

(display (matrix-add b c)) (newline)

(display (matrix-scalar-multiply 2 b)) (newline)

(display (matrix-multiply b d)) (newline)

(display (matrix-transpose b)) (newline)

(display (matrix-multiply a d)) (newline)

#(#(3 5 7) #(7 9 11))
#(#(2 4 6) #(8 10 12))
#(#(14 20 26 32) #(32 47 62 77))
#(#(1 4) #(2 5) #(3 6))
matrix-multiply: incompatible matrices