add gaussian-elim

src/matrix/auto.clj

@@ -29,6 +29,24 @@
 
 (defn gaussian-column
   [m col]
+  (if (zero? (at m col col))
+    []
     (cons {:mul-row col :by (/ 1 (at m col col))}
           (for [y (range (inc col) (matrix-height m))]
-          {:add-row col :times (- (at m col y)) :to-row y})))
+            {:add-row col :times (- (at m col y)) :to-row y}))))
+
+(defn gaussian-elim
+  "Takes a matrix and returns the steps to put in the form
+   [[1 a ... | sum1]
+    [0 1 ... | sum2]
+    ...]
+
+   Assumes a lot of things lol"
+  ([m]
+   (gaussian-elim m 0))
+  ([m col]
+   (if (= col (matrix-height m))
+     []
+     (let [steps (gaussian-column m col)
+           new   (apply-all-at-once m steps)]
+       (concat steps (gaussian-elim new (inc col)))))))

src/matrix/base.clj

@@ -72,8 +72,8 @@
   "Apply a spec-compliant row-op to a matrix, curriable"
   ([op]
    (fn [m]
-     (apply-row-op op m)))
+     (apply-row-op m op)))
-  ([op m]
+  ([m op]
    (if-let [[type inner] (opspec/conform op)]
      (case type
        :swap (swap-rows
@@ -84,12 +84,23 @@
               m
               (:times   inner)
               (:add-row inner)
-              (:to      inner))
+              (:to-row  inner))
        :mul  (mul-row
               m
               (:by      inner)
               (:mul-row inner))))))
 
+(defn apply-steps
+  [m steps]
+  (if (empty? steps)
+    [m]
+    (let [new (apply-row-op m (first steps))]
+      (cons new (apply-steps new (rest steps))))))
+
+(defn apply-all-at-once
+  [m steps]
+  (last (apply-steps m steps)))
+
 (defn at
   "Returns the number at x y in a matrix"
   [m x y]