oct 16, 18 notes. add unit tests with utest, and bisection root finding methods

3 files changed, 69 insertions, 1 deletions

diff --git a/src/matrix.c b/src/matrix.c
index ac4d138..5766d94 100644
--- a/src/matrix.c
+++ b/src/matrix.c
@@ -37,7 +37,7 @@ Matrix_double **lu_decomp(Matrix_double *m) {
   Matrix_double *u = copy_matrix(m);
   Matrix_double *l_empt = InitMatrixWithSize(double, m->rows, m->cols, 0.0);
   Matrix_double *l = put_identity_diagonal(l_empt);
-  free(l_empt);
+  free_matrix(l_empt);
 
   Matrix_double **u_l = malloc(sizeof(Matrix_double *) * 2);
 
@@ -140,3 +140,14 @@ void format_matrix_into(Matrix_double *m, char *s) {
   }
   strcat(s, "\n");
 }
+
+int matrix_equal(Matrix_double *a, Matrix_double *b) {
+  if (a->cols != b->cols || a->rows != b->rows)
+    return false;
+
+  for (size_t y = 0; y < a->rows; ++y)
+    if (!vector_equal(a->data[y], b->data[y])) {
+      return false;
+    }
+  return true;
+}
diff --git a/src/roots.c b/src/roots.c
new file mode 100644
index 0000000..de16adf
--- /dev/null
+++ b/src/roots.c
@@ -0,0 +1,46 @@
+#include "lizfcm.h"
+#include <assert.h>
+#include <math.h>
+
+// f is well defined at start_x + delta*n for all n on the integer range [0,
+// max_steps]
+double *find_ivt_range(double (*f)(double), double start_x, double delta,
+                       size_t max_steps) {
+  double *range = malloc(sizeof(double) * 2);
+
+  double a = start_x;
+
+  while (f(a) * f(start_x) >= 0 && max_steps-- > 0)
+    a += delta;
+
+  if (max_steps == 0 && f(a) * f(start_x) > 0)
+    return NULL;
+
+  range[0] = start_x;
+  range[1] = a + delta;
+  return range;
+}
+
+// f is continuous on [a, b]
+double bisect_find_root(double (*f)(double), double a, double b,
+                        double tolerance, size_t max_iterations) {
+  assert(a <= b);
+  // guarantee there's a root somewhere between a and b by IVT
+  assert(f(a) * f(b) < 0);
+
+  double c = (1.0 / 2) * (a + b);
+  if (b - a < tolerance || max_iterations == 0)
+    return c;
+  if (f(a) * f(c) < 0)
+    return bisect_find_root(f, a, c, tolerance, max_iterations - 1);
+  return bisect_find_root(f, c, b, tolerance, max_iterations - 1);
+}
+
+double bisect_find_root_with_error_assumption(double (*f)(double), double a,
+                                              double b, double tolerance) {
+  assert(a <= b);
+
+  uint64_t max_iterations =
+      (uint64_t)ceil((log(tolerance) - log(b - a)) / log(1 / 2.0));
+  return bisect_find_root(f, a, b, tolerance, max_iterations);
+}
diff --git a/src/vector.c b/src/vector.c
index 4a6250c..3e4f62d 100644
--- a/src/vector.c
+++ b/src/vector.c
@@ -115,3 +115,14 @@ double sum_v(Array_double *v) {
     sum += v->data[i];
   return sum;
 }
+
+int vector_equal(Array_double *a, Array_double *b) {
+  if (a->size != b->size)
+    return false;
+
+  for (size_t i = 0; i < a->size; ++i) {
+    if (a->data[i] != b->data[i])
+      return false;
+  }
+  return true;
+}