3 files changed, 78 insertions, 47 deletions

diff --git a/src/lin.c b/src/lin.c
index 2df6f28..d531025 100644
--- a/src/lin.c
+++ b/src/lin.c
@@ -7,8 +7,8 @@ Line *least_squares_lin_reg(Array_double *x, Array_double *y) {
   uint64_t n = x->size;
   double sum_x = sum_v(x);
   double sum_y = sum_v(y);
-  double sum_xy = dot_v(x, y);
-  double sum_xx = dot_v(x, x);
+  double sum_xy = v_dot_v(x, y);
+  double sum_xx = v_dot_v(x, x);
   double denom = ((n * sum_xx) - (sum_x * sum_x));
 
   Line *line = malloc(sizeof(Line));
diff --git a/src/matrix.c b/src/matrix.c
index 438a468..51de22c 100644
--- a/src/matrix.c
+++ b/src/matrix.c
@@ -3,11 +3,23 @@
 #include <stdio.h>
 #include <string.h>
 
-void put_identity_diagonal(Matrix_double *m) {
-  assert(m->rows == m->cols);
+Array_double *m_dot_v(Matrix_double *m, Array_double *v) {
+  assert(v->size == m->cols);
+
+  Array_double *product = copy_vector(v);
+
+  for (size_t row = 0; row < v->size; ++row)
+    product->data[row] = v_dot_v(m->data[row], v);
+
+  return product;
+}
 
+Matrix_double *put_identity_diagonal(Matrix_double *m) {
+  assert(m->rows == m->cols);
+  Matrix_double *copy = copy_matrix(m);
   for (size_t y = 0; y < m->rows; ++y)
-    m->data[y]->data[y] = 1.0;
+    copy->data[y]->data[y] = 1.0;
+  return copy;
 }
 
 Matrix_double *copy_matrix(Matrix_double *m) {
@@ -89,6 +101,25 @@ Array_double *fsubst(Matrix_double *l, Array_double *b) {
   return x;
 }
 
+Array_double *solve_matrix(Matrix_double *m, Array_double *b) {
+  assert(b->size == m->rows);
+  assert(m->rows == m->cols);
+
+  Array_double *x = copy_vector(b);
+  Matrix_double **u_l = lu_decomp(m);
+  Matrix_double *u = u_l[0];
+  Matrix_double *l = u_l[1];
+
+  Array_double *b_fsub = fsubst(l, b);
+  x = bsubst(u, b_fsub);
+  free_vector(b_fsub);
+
+  free_matrix(u);
+  free_matrix(l);
+
+  return x;
+}
+
 void free_matrix(Matrix_double *m) {
   for (size_t y = 0; y < m->rows; ++y)
     free_vector(m->data[y]);
diff --git a/src/vector.c b/src/vector.c
index a7ff783..4a6250c 100644
--- a/src/vector.c
+++ b/src/vector.c
@@ -2,30 +2,18 @@
 #include <assert.h>
 #include <float.h>
 #include <math.h>
-#include <string.h>
 #include <stdio.h>
+#include <string.h>
 
-double l2_norm(Array_double *v) {
-  double norm = 0;
-  for (size_t i = 0; i < v->size; ++i)
-    norm += v->data[i] * v->data[i];
-  return sqrt(norm);
-}
+Array_double *add_v(Array_double *v1, Array_double *v2) {
+  assert(v1->size == v2->size);
 
-double l1_norm(Array_double *v) {
-  double sum = 0;
-  for (size_t i = 0; i < v->size; ++i)
-    sum += fabs(v->data[i]);
+  Array_double *sum = copy_vector(v1);
+  for (size_t i = 0; i < v1->size; i++)
+    sum->data[i] += v2->data[i];
   return sum;
 }
 
-double linf_norm(Array_double *v) {
-  double max = -DBL_MAX;
-  for (size_t i = 0; i < v->size; ++i)
-    max = c_max(v->data[i], max);
-  return max;
-}
-
 Array_double *minus_v(Array_double *v1, Array_double *v2) {
   assert(v1->size == v2->size);
 
@@ -35,13 +23,6 @@ Array_double *minus_v(Array_double *v1, Array_double *v2) {
   return sub;
 }
 
-double sum_v(Array_double *v) {
-  double sum = 0;
-  for (size_t i = 0; i < v->size; i++)
-    sum += v->data[i];
-  return sum;
-}
-
 Array_double *scale_v(Array_double *v, double m) {
   Array_double *copy = copy_vector(v);
   for (size_t i = 0; i < v->size; i++)
@@ -49,16 +30,29 @@ Array_double *scale_v(Array_double *v, double m) {
   return copy;
 }
 
-Array_double *add_v(Array_double *v1, Array_double *v2) {
-  assert(v1->size == v2->size);
-
-  Array_double *sum = copy_vector(v1);
-  for (size_t i = 0; i < v1->size; i++)
-    sum->data[i] += v2->data[i];
+double l1_norm(Array_double *v) {
+  double sum = 0;
+  for (size_t i = 0; i < v->size; ++i)
+    sum += fabs(v->data[i]);
   return sum;
 }
 
-double dot_v(Array_double *v1, Array_double *v2) {
+double l2_norm(Array_double *v) {
+  double norm = 0;
+  for (size_t i = 0; i < v->size; ++i)
+    norm += v->data[i] * v->data[i];
+  return sqrt(norm);
+}
+
+double linf_norm(Array_double *v) {
+  assert(v->size > 0);
+  double max = v->data[0];
+  for (size_t i = 0; i < v->size; ++i)
+    max = c_max(v->data[i], max);
+  return max;
+}
+
+double v_dot_v(Array_double *v1, Array_double *v2) {
   assert(v1->size == v2->size);
 
   double dot = 0;
@@ -67,25 +61,24 @@ double dot_v(Array_double *v1, Array_double *v2) {
   return dot;
 }
 
-double l2_distance(Array_double *v1, Array_double *v2) {
+double vector_distance(Array_double *v1, Array_double *v2,
+                       double (*norm)(Array_double *)) {
   Array_double *minus = minus_v(v1, v2);
-  double dist = l2_norm(minus);
+  double dist = (*norm)(minus);
   free(minus);
   return dist;
 }
 
 double l1_distance(Array_double *v1, Array_double *v2) {
-  Array_double *minus = minus_v(v1, v2);
-  double dist = l1_norm(minus);
-  free(minus);
-  return dist;
+  return vector_distance(v1, v2, &l1_norm);
+}
+
+double l2_distance(Array_double *v1, Array_double *v2) {
+  return vector_distance(v1, v2, &l2_norm);
 }
 
 double linf_distance(Array_double *v1, Array_double *v2) {
-  Array_double *minus = minus_v(v1, v2);
-  double dist = linf_norm(minus);
-  free(minus);
-  return dist;
+  return vector_distance(v1, v2, &linf_norm);
 }
 
 Array_double *copy_vector(Array_double *v) {
@@ -115,3 +108,10 @@ void format_vector_into(Array_double *v, char *s) {
   }
   strcat(s, "\n");
 }
+
+double sum_v(Array_double *v) {
+  double sum = 0;
+  for (size_t i = 0; i < v->size; i++)
+    sum += v->data[i];
+  return sum;
+}