# HG changeset patch
# User Nolan Prescott <mail@nprescott.com>
# Date 1584676153 14400
#      Thu Mar 19 23:49:13 2020 -0400
# Node ID 972c81dcde9b073d30818260c8cc6fb7552c7ba1
# Parent  3bd215111329ad2a23773da355978d0b2c82c9f3
add a few basic matrix operations

diff --git a/tracer.lisp b/tracer.lisp
--- a/tracer.lisp
+++ b/tracer.lisp
@@ -90,7 +90,7 @@
       (push (subseq text offset previous) lines)
       (setq offset (1+ previous)))))
 
-(flatten (ls)
+(defun flatten (ls)
   (labels ((mklist (x) (if (listp x) x (list x))))
     (mapcan #'(lambda (x) (if (atom x) (mklist x) (flatten x))) ls)))
 
@@ -114,6 +114,145 @@
   (with-open-file (stream path :direction :output :if-exists :supersede)
     (format stream ppm-string)))
 
+(defun matrix*vector (matrix vec)
+  (let* ((m (array-dimension matrix 0))
+         (n (length vec))
+         (result (make-array m :initial-element 0)))
+    (dotimes (i m result)
+      (dotimes (j n)
+        (incf (aref result i)
+              (* (aref matrix i j)
+                 (aref vec j)))))))
+
+(defun matrix*matrix (a b)
+  (let* ((m (array-dimension a 0))
+         (n (array-dimension b 1))
+         (common (array-dimension b 0))
+         (result (make-array (list m n) :initial-element 0)))
+    (dotimes (i m result)
+      (dotimes (j n)
+        (dotimes (k common)
+          (incf (aref result i j)
+                (* (aref a i k)
+                   (aref b k j))))))))
+
+(defun determinant (m)
+  ;; only operates on 2x2 matrices
+  (- (* (aref m 0 0) (aref m 1 1))
+     (* (aref m 0 1) (aref m 1 0))))
+
+(defun submatrix (m row column)
+  (let ((result (make-array (mapcar #'1- (array-dimensions m)))))
+    (loop for i below (array-dimension m 0)
+       unless (= i row)
+       do (loop for j below (array-dimension m 1)
+             unless (= j column)
+             do (let ((ii i)
+                      (jj j))
+                  (if (> i row) (setf ii (1- i)))
+                  (if (> j column) (setf jj (1- j)))
+                  (setf (aref result ii jj)
+                        (aref m i j)))))
+    result))
+
+;;; since we're never taking the identity but for a 4x4 matrix...
+(defvar identity-matrix #2A((1 0 0 0)
+                            (0 1 0 0)
+                            (0 0 1 0)
+                            (0 0 0 1)))
+
+(defun transpose-matrix (a)
+  (let ((result (make-array (array-dimensions a))))
+    (dotimes (i (array-dimension a 0) result)
+      (dotimes (j (array-dimension a 1))
+        (setf (aref result i j)
+              (aref   a    j i))))))
+
+(deftest multiply-4x4-matrices
+  (let ((a #2A((1 2 3 4)
+               (5 6 7 8)
+               (9 8 7 6)
+               (5 4 3 2)))
+        (b #2A((-2 1 2  3)
+               ( 3 2 1 -1)
+               ( 4 3 6  5)
+               ( 1 2 7  8))))
+    (check (equalp (matrix*matrix a b)
+                   #2A((20 22  50  48)
+                       (44 54 114 108)
+                       (40 58 110 102)
+                       (16 26  46  42))))))
+
+(deftest multiply-matrix-with-vector
+  (let ((a #2A((1 2 3 4)
+               (2 4 4 2)
+               (8 6 4 1)
+               (0 0 0 1)))
+        (b (vector 1 2 3 1)))
+    (check (equalp (matrix*vector a b) (vector 18 24 33 1)))))
+
+(deftest identity-matrix-returns-original-matrix
+  (check (equalp
+          (matrix*matrix identity-matrix #2A((1 2 3 4)
+                                             (2 4 4 2)
+                                             (8 6 4 1)
+                                             (0 0 0 1)))
+          #2A((1 2 3 4)
+              (2 4 4 2)
+              (8 6 4 1)
+              (0 0 0 1)))))
+
+(deftest identity-matrix-times-tuple
+  (check (equalp
+          (matrix*vector identity-matrix (vector 1 2 3 4))
+          (vector 1 2 3 4))))
+
+(deftest transpose-a-matrix
+  (let ((a #2A((0 9 3 0)
+               (9 8 0 8)
+               (1 8 5 3)
+               (0 0 5 6))))
+    (check (equalp (transpose-matrix a) #2A((0 9 1 0)
+                                            (9 8 8 0)
+                                            (3 0 5 5)
+                                            (0 8 3 6))))))
+
+(deftest transpose-identity-matrix
+  (check (equalp (transpose-matrix identity-matrix)
+                 identity-matrix)))
+
+(deftest verify-2x2-determinant
+  (let ((a #2A(( 1 5)
+               (-3 2))))
+    (check (equalp (determinant a) 17))))
+
+(deftest submatrix-of-3x3-is-2x2
+  (let ((m #2A(( 1 5  0)
+               (-3 2  7)
+               ( 0 6 -3))))
+    (check (equalp (submatrix m 0 2) #2A((-3 2)
+                                         ( 0 6))))))
+
+(deftest submatrix-of-4x4-is-3x3
+  (let ((m #2A((-6 1  1 6)
+               (-8 5  8 6)
+               (-1 0  8 2)
+               (-7 1 -1 1))))
+    (check (equalp (submatrix m 2 1) #2A((-6  1 6)
+                                         (-8  8 6)
+                                         (-7 -1 1))))))
+
+(deftest matrix-basics
+  (multiply-4x4-matrices)
+  (multiply-matrix-with-vector)
+  (identity-matrix-returns-original-matrix)
+  (identity-matrix-times-tuple)
+  (transpose-a-matrix)
+  (transpose-identity-matrix)
+  (verify-2x2-determinant)
+  (submatrix-of-3x3-is-2x2)
+  (submatrix-of-4x4-is-3x3))
+
 ;;; tests:
 (deftest adding-vec3
   (let ((a (vector 3 -2 5))
@@ -275,7 +414,7 @@
   (vector-basics)
   (canvas-and-visuals))
 
-;;; projectile tracking
+;;; sample
 (defstruct projectile position velocity)
 (defstruct environment gravity wind)
 
@@ -285,27 +424,28 @@
                                                     (environment-wind env)))))
     (make-projectile :position p :velocity v)))
 
-(do* ((env (make-environment :gravity (vector 0 -0.1 0)
-                             :wind (vector -0.01 0 0)))
-      (p (make-projectile :position (vector 0 1 0)
-                          :velocity (normalize (vector 1 1 0)))
-         (tick env p)))
-     ((if (<= (Y (projectile-position p)) 0) (return p)))
-  (format t "~a~%" (projectile-position p)))
+(defun projectile-tracking ()
+  (do* ((env (make-environment :gravity (vector 0 -0.1 0)
+                               :wind (vector -0.01 0 0)))
+        (p (make-projectile :position (vector 0 1 0)
+                            :velocity (normalize (vector 1 1 0)))
+           (tick env p)))
+       ((if (<= (Y (projectile-position p)) 0) (return p)))
+    (format t "~a~%" (projectile-position p))))
 
-;;; projectile visualization
-(do* ((env (make-environment :gravity (vector 0 -0.09 0)
-                             :wind (vector -0.01 0 0)))
-      (p (make-projectile :position (vector 0 1 0)
-                          :velocity (scale-vec3 (normalize (vector 1 1.5 0)) 10))
-         (tick env p))
-      (c (make-canvas :width 900 :height 550)))
-     ((<= (Y (projectile-position p)) 0)
-      (ppm->file (canvas->ppm c) "/home/nolan/test-output.ppm"))
-  (format t "~s ~s~%"
-          (max 0 (floor (X (projectile-position p))))
-          (max 0 (floor (- (canvas-height c) 100))))
-  (write-pixel c
-               (max 0 (floor (X (projectile-position p))))
-               (max 0 (floor (- (canvas-height c) (Y (projectile-position p)))))
-               (vector 255 255 255)))
+(defun projectile-visualization ()
+  (do* ((env (make-environment :gravity (vector 0 -0.093 0)
+                               :wind (vector -0.01 0 0)))
+        (p (make-projectile :position (vector 0 1 0)
+                            :velocity (scale-vec3 (normalize (vector 1 1.5 0)) 10))
+           (tick env p))
+        (c (make-canvas :width 900 :height 550)))
+       ((<= (Y (projectile-position p)) 0)
+        (ppm->file (canvas->ppm c) "/home/nolan/test-output.ppm"))
+    (format t "~s ~s~%"
+            (max 0 (floor (X (projectile-position p))))
+            (max 0 (floor (- (canvas-height c) 100))))
+    (write-pixel c
+                 (max 0 (floor (X (projectile-position p))))
+                 (max 0 (floor (- (canvas-height c) (Y (projectile-position p)))))
+                 (vector 255 255 255))))