From 6bf4b90c90f15f4ab60833bddf5b5756d1a6b1f6 Mon Sep 17 00:00:00 2001
From: Elizabeth Alexander Hunt <me@liz.coffee>
Date: Thu, 2 Jul 2026 11:55:17 -0700
Subject: Init

---
 Homework/math4610/src/approx_derivative.c |  38 ++++
 Homework/math4610/src/eigen.c             | 116 ++++++++++
 Homework/math4610/src/lin.c               |  19 ++
 Homework/math4610/src/maceps.c            |  28 +++
 Homework/math4610/src/matrix.c            | 346 ++++++++++++++++++++++++++++++
 Homework/math4610/src/rand.c              |   5 +
 Homework/math4610/src/roots.c             | 127 +++++++++++
 Homework/math4610/src/vector.c            | 143 ++++++++++++
 8 files changed, 822 insertions(+)
 create mode 100644 Homework/math4610/src/approx_derivative.c
 create mode 100644 Homework/math4610/src/eigen.c
 create mode 100644 Homework/math4610/src/lin.c
 create mode 100644 Homework/math4610/src/maceps.c
 create mode 100644 Homework/math4610/src/matrix.c
 create mode 100644 Homework/math4610/src/rand.c
 create mode 100644 Homework/math4610/src/roots.c
 create mode 100644 Homework/math4610/src/vector.c

(limited to 'Homework/math4610/src')

diff --git a/Homework/math4610/src/approx_derivative.c b/Homework/math4610/src/approx_derivative.c
new file mode 100644
index 0000000..63d0b05
--- /dev/null
+++ b/Homework/math4610/src/approx_derivative.c
@@ -0,0 +1,38 @@
+#include "lizfcm.h"
+#include <assert.h>
+
+double central_derivative_at(double (*f)(double), double a, double h) {
+  assert(h > 0);
+
+  double x2 = a + h;
+  double x1 = a - h;
+
+  double y2 = f(x2);
+  double y1 = f(x1);
+
+  return (y2 - y1) / (x2 - x1);
+}
+
+double forward_derivative_at(double (*f)(double), double a, double h) {
+  assert(h > 0);
+
+  double x2 = a + h;
+  double x1 = a;
+
+  double y2 = f(x2);
+  double y1 = f(x1);
+
+  return (y2 - y1) / (x2 - x1);
+}
+
+double backward_derivative_at(double (*f)(double), double a, double h) {
+  assert(h > 0);
+
+  double x2 = a;
+  double x1 = a - h;
+
+  double y2 = f(x2);
+  double y1 = f(x1);
+
+  return (y2 - y1) / (x2 - x1);
+}
diff --git a/Homework/math4610/src/eigen.c b/Homework/math4610/src/eigen.c
new file mode 100644
index 0000000..49cc0e4
--- /dev/null
+++ b/Homework/math4610/src/eigen.c
@@ -0,0 +1,116 @@
+#include "lizfcm.h"
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+
+double least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
+                                 double tolerance, size_t max_iterations) {
+  return shift_inverse_power_eigenvalue(m, v, 0.0, tolerance, max_iterations);
+}
+
+double dominant_eigenvalue(Matrix_double *m, Array_double *v, double tolerance,
+                           size_t max_iterations) {
+  assert(m->rows == m->cols);
+  assert(m->rows == v->size);
+
+  double error = tolerance;
+  size_t iter = max_iterations;
+  double lambda = 0.0;
+  Array_double *eigenvector_1 = copy_vector(v);
+
+  while (error >= tolerance && (--iter) > 0) {
+    Array_double *eigenvector_2 = m_dot_v(m, eigenvector_1);
+    Array_double *normalized_eigenvector_2 =
+        scale_v(eigenvector_2, 1.0 / linf_norm(eigenvector_2));
+    free_vector(eigenvector_2);
+    eigenvector_2 = normalized_eigenvector_2;
+
+    Array_double *mx = m_dot_v(m, eigenvector_2);
+    double new_lambda =
+        v_dot_v(mx, eigenvector_2) / v_dot_v(eigenvector_2, eigenvector_2);
+
+    error = fabs(new_lambda - lambda);
+    lambda = new_lambda;
+    free_vector(eigenvector_1);
+    eigenvector_1 = eigenvector_2;
+  }
+
+  return lambda;
+}
+
+double shift_inverse_power_eigenvalue(Matrix_double *m, Array_double *v,
+                                      double shift, double tolerance,
+                                      size_t max_iterations) {
+  assert(m->rows == m->cols);
+  assert(m->rows == v->size);
+
+  Matrix_double *m_c = copy_matrix(m);
+  for (size_t y = 0; y < m_c->rows; ++y)
+    m_c->data[y]->data[y] = m_c->data[y]->data[y] - shift;
+
+  double error = tolerance;
+  size_t iter = max_iterations;
+  double lambda = shift;
+  Array_double *eigenvector_1 = copy_vector(v);
+
+  while (error >= tolerance && (--iter) > 0) {
+    Array_double *eigenvector_2 = solve_matrix_lu_bsubst(m_c, eigenvector_1);
+    Array_double *normalized_eigenvector_2 =
+        scale_v(eigenvector_2, 1.0 / linf_norm(eigenvector_2));
+    free_vector(eigenvector_2);
+
+    Array_double *mx = m_dot_v(m, normalized_eigenvector_2);
+    double new_lambda =
+        v_dot_v(mx, normalized_eigenvector_2) /
+        v_dot_v(normalized_eigenvector_2, normalized_eigenvector_2);
+
+    error = fabs(new_lambda - lambda);
+    lambda = new_lambda;
+    free_vector(eigenvector_1);
+    eigenvector_1 = normalized_eigenvector_2;
+  }
+
+  return lambda;
+}
+
+Array_double *partition_find_eigenvalues(Matrix_double *m,
+                                         Matrix_double *guesses,
+                                         double tolerance,
+                                         size_t max_iterations) {
+  assert(guesses->rows >=
+         2); // we need at least, the most and least dominant eigenvalues
+
+  double end = dominant_eigenvalue(m, guesses->data[guesses->rows - 1],
+                                   tolerance, max_iterations);
+  double begin =
+      least_dominant_eigenvalue(m, guesses->data[0], tolerance, max_iterations);
+
+  double delta = (end - begin) / guesses->rows;
+  Array_double *eigenvalues = InitArrayWithSize(double, guesses->rows, 0.0);
+  for (size_t i = 0; i < guesses->rows; i++) {
+    double box_midpoint = ((delta * i) + (delta * (i + 1))) / 2;
+
+    double nearest_eigenvalue = shift_inverse_power_eigenvalue(
+        m, guesses->data[i], box_midpoint, tolerance, max_iterations);
+
+    eigenvalues->data[i] = nearest_eigenvalue;
+  }
+
+  return eigenvalues;
+}
+
+Matrix_double *leslie_matrix(Array_double *age_class_surivor_ratio,
+                             Array_double *age_class_offspring) {
+  assert(age_class_surivor_ratio->size + 1 == age_class_offspring->size);
+
+  Matrix_double *leslie = InitMatrixWithSize(double, age_class_offspring->size,
+                                             age_class_offspring->size, 0.0);
+
+  free_vector(leslie->data[0]);
+  leslie->data[0] = copy_vector(age_class_offspring);
+
+  for (size_t i = 0; i < age_class_surivor_ratio->size; i++)
+    leslie->data[i + 1]->data[i] = age_class_surivor_ratio->data[i];
+  return leslie;
+}
diff --git a/Homework/math4610/src/lin.c b/Homework/math4610/src/lin.c
new file mode 100644
index 0000000..d531025
--- /dev/null
+++ b/Homework/math4610/src/lin.c
@@ -0,0 +1,19 @@
+#include "lizfcm.h"
+#include <assert.h>
+
+Line *least_squares_lin_reg(Array_double *x, Array_double *y) {
+  assert(x->size == y->size);
+
+  uint64_t n = x->size;
+  double sum_x = sum_v(x);
+  double sum_y = sum_v(y);
+  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));
+  line->m = ((sum_xy * n) - (sum_x * sum_y)) / denom;
+  line->a = ((sum_y * sum_xx) - (sum_x * sum_xy)) / denom;
+
+  return line;
+}
diff --git a/Homework/math4610/src/maceps.c b/Homework/math4610/src/maceps.c
new file mode 100644
index 0000000..23bc9db
--- /dev/null
+++ b/Homework/math4610/src/maceps.c
@@ -0,0 +1,28 @@
+#include "lizfcm.h"
+#include <math.h>
+
+float smaceps() {
+  float one = 1.0;
+  float machine_epsilon = 1.0;
+  float one_approx = one + machine_epsilon;
+
+  while (fabsf(one_approx - one) > 0) {
+    machine_epsilon /= 2;
+    one_approx = one + machine_epsilon;
+  }
+
+  return machine_epsilon;
+}
+
+double dmaceps() {
+  double one = 1.0;
+  double machine_epsilon = 1.0;
+  double one_approx = one + machine_epsilon;
+
+  while (fabs(one_approx - one) > 0) {
+    machine_epsilon /= 2;
+    one_approx = one + machine_epsilon;
+  }
+
+  return machine_epsilon;
+}
diff --git a/Homework/math4610/src/matrix.c b/Homework/math4610/src/matrix.c
new file mode 100644
index 0000000..901a426
--- /dev/null
+++ b/Homework/math4610/src/matrix.c
@@ -0,0 +1,346 @@
+#include "lizfcm.h"
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+
+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;
+}
+
+Array_double *col_v(Matrix_double *m, size_t x) {
+  assert(x < m->cols);
+
+  Array_double *col = InitArrayWithSize(double, m->rows, 0.0);
+  for (size_t y = 0; y < m->rows; y++)
+    col->data[y] = m->data[y]->data[x];
+
+  return col;
+}
+
+Matrix_double *m_dot_m(Matrix_double *a, Matrix_double *b) {
+  assert(a->cols == b->rows);
+
+  Matrix_double *prod = InitMatrixWithSize(double, a->rows, b->cols, 0.0);
+
+  Array_double *curr_col;
+  for (size_t y = 0; y < a->rows; y++) {
+    for (size_t x = 0; x < b->cols; x++) {
+      curr_col = col_v(b, x);
+      prod->data[y]->data[x] = v_dot_v(curr_col, a->data[y]);
+      free_vector(curr_col);
+    }
+  }
+
+  return prod;
+}
+
+Matrix_double *transpose(Matrix_double *m) {
+  Matrix_double *transposed = InitMatrixWithSize(double, m->cols, m->rows, 0.0);
+
+  for (size_t x = 0; x < m->rows; x++)
+    for (size_t y = 0; y < m->cols; y++)
+      transposed->data[y]->data[x] = m->data[x]->data[y];
+
+  return transposed;
+}
+
+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)
+    copy->data[y]->data[y] = 1.0;
+  return copy;
+}
+
+Matrix_double *copy_matrix(Matrix_double *m) {
+  Matrix_double *copy = InitMatrixWithSize(double, m->rows, m->cols, 0.0);
+  for (size_t y = 0; y < copy->rows; y++) {
+    free_vector(copy->data[y]);
+    copy->data[y] = copy_vector(m->data[y]);
+  }
+  return copy;
+}
+
+Matrix_double **lu_decomp(Matrix_double *m) {
+  assert(m->cols == m->rows);
+
+  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_matrix(l_empt);
+
+  Matrix_double **u_l = malloc(sizeof(Matrix_double *) * 2);
+
+  for (size_t y = 0; y < m->rows; y++) {
+    if (u->data[y]->data[y] == 0) {
+      printf("ERROR: a pivot is zero in given matrix\n");
+      assert(false);
+    }
+  }
+
+  if (u && l) {
+    for (size_t x = 0; x < m->cols; x++) {
+      for (size_t y = x + 1; y < m->rows; y++) {
+        double denom = u->data[x]->data[x];
+
+        if (denom == 0) {
+          printf("ERROR: non-factorable matrix\n");
+          assert(false);
+        }
+
+        double factor = -(u->data[y]->data[x] / denom);
+
+        Array_double *scaled = scale_v(u->data[x], factor);
+        Array_double *added = add_v(scaled, u->data[y]);
+        free_vector(scaled);
+        free_vector(u->data[y]);
+
+        u->data[y] = added;
+        l->data[y]->data[x] = -factor;
+      }
+    }
+  }
+
+  u_l[0] = u;
+  u_l[1] = l;
+  return u_l;
+}
+
+Array_double *bsubst(Matrix_double *u, Array_double *b) {
+  assert(u->rows == b->size && u->cols == u->rows);
+
+  Array_double *x = copy_vector(b);
+  for (int64_t row = b->size - 1; row >= 0; row--) {
+    for (size_t col = b->size - 1; col > row; col--)
+      x->data[row] -= x->data[col] * u->data[row]->data[col];
+    x->data[row] /= u->data[row]->data[row];
+  }
+  return x;
+}
+
+Array_double *fsubst(Matrix_double *l, Array_double *b) {
+  assert(l->rows == b->size && l->cols == l->rows);
+
+  Array_double *x = copy_vector(b);
+
+  for (size_t row = 0; row < b->size; row++) {
+    for (size_t col = 0; col < row; col++)
+      x->data[row] -= x->data[col] * l->data[row]->data[col];
+    x->data[row] /= l->data[row]->data[row];
+  }
+
+  return x;
+}
+
+Array_double *solve_matrix_lu_bsubst(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);
+  free(u_l);
+
+  return x;
+}
+
+Matrix_double *gaussian_elimination(Matrix_double *m) {
+  uint64_t h = 0, k = 0;
+
+  Matrix_double *m_cp = copy_matrix(m);
+
+  while (h < m_cp->rows && k < m_cp->cols) {
+    uint64_t max_row = h;
+    double max_val = 0.0;
+
+    for (uint64_t row = h; row < m_cp->rows; row++) {
+      double val = fabs(m_cp->data[row]->data[k]);
+      if (val > max_val) {
+        max_val = val;
+        max_row = row;
+      }
+    }
+
+    if (max_val == 0.0) {
+      k++;
+      continue;
+    }
+
+    if (max_row != h) {
+      Array_double *swp = m_cp->data[max_row];
+      m_cp->data[max_row] = m_cp->data[h];
+      m_cp->data[h] = swp;
+    }
+
+    for (uint64_t row = h + 1; row < m_cp->rows; row++) {
+      double factor = m_cp->data[row]->data[k] / m_cp->data[h]->data[k];
+      m_cp->data[row]->data[k] = 0.0;
+
+      for (uint64_t col = k + 1; col < m_cp->cols; col++) {
+        m_cp->data[row]->data[col] -= m_cp->data[h]->data[col] * factor;
+      }
+    }
+
+    h++;
+    k++;
+  }
+
+  return m_cp;
+}
+
+Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b) {
+  assert(b->size == m->rows);
+  assert(m->rows == m->cols);
+
+  Matrix_double *m_augment_b = add_column(m, b);
+  Matrix_double *eliminated = gaussian_elimination(m_augment_b);
+
+  Array_double *b_gauss = col_v(eliminated, m->cols);
+  Matrix_double *u = slice_column(eliminated, m->rows);
+
+  Array_double *solution = bsubst(u, b_gauss);
+
+  free_matrix(m_augment_b);
+  free_matrix(eliminated);
+  free_matrix(u);
+  free_vector(b_gauss);
+
+  return solution;
+}
+
+Array_double *jacobi_solve(Matrix_double *m, Array_double *b,
+                           double l2_convergence_tolerance,
+                           size_t max_iterations) {
+  assert(m->rows == m->cols);
+  assert(b->size == m->cols);
+  size_t iter = max_iterations;
+
+  Array_double *x_k = InitArrayWithSize(double, b->size, 0.0);
+  Array_double *x_k_1 =
+      InitArrayWithSize(double, b->size, rand_from(0.1, 10.0));
+
+  while ((--iter) > 0 && l2_distance(x_k_1, x_k) > l2_convergence_tolerance) {
+    for (size_t i = 0; i < m->rows; i++) {
+      double delta = 0.0;
+      for (size_t j = 0; j < m->cols; j++) {
+        if (i == j)
+          continue;
+        delta += m->data[i]->data[j] * x_k->data[j];
+      }
+      x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i];
+    }
+
+    Array_double *tmp = x_k;
+    x_k = x_k_1;
+    x_k_1 = tmp;
+  }
+
+  free_vector(x_k);
+  return x_k_1;
+}
+
+Array_double *gauss_siedel_solve(Matrix_double *m, Array_double *b,
+                                 double l2_convergence_tolerance,
+                                 size_t max_iterations) {
+  assert(m->rows == m->cols);
+  assert(b->size == m->cols);
+  size_t iter = max_iterations;
+
+  Array_double *x_k = InitArrayWithSize(double, b->size, 0.0);
+  Array_double *x_k_1 =
+      InitArrayWithSize(double, b->size, rand_from(0.1, 10.0));
+
+  while ((--iter) > 0) {
+    for (size_t i = 0; i < x_k->size; i++)
+      x_k->data[i] = x_k_1->data[i];
+
+    for (size_t i = 0; i < m->rows; i++) {
+      double delta = 0.0;
+      for (size_t j = 0; j < m->cols; j++) {
+        if (i == j)
+          continue;
+        delta += m->data[i]->data[j] * x_k_1->data[j];
+      }
+      x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i];
+    }
+
+    if (l2_distance(x_k_1, x_k) <= l2_convergence_tolerance)
+      break;
+  }
+
+  free_vector(x_k);
+  return x_k_1;
+}
+
+Matrix_double *slice_column(Matrix_double *m, size_t x) {
+  Matrix_double *sliced = copy_matrix(m);
+
+  for (size_t row = 0; row < m->rows; row++) {
+    Array_double *old_row = sliced->data[row];
+    sliced->data[row] = slice_element(old_row, x);
+    free_vector(old_row);
+  }
+  sliced->cols--;
+
+  return sliced;
+}
+
+Matrix_double *add_column(Matrix_double *m, Array_double *v) {
+  Matrix_double *pushed = copy_matrix(m);
+
+  for (size_t row = 0; row < m->rows; row++) {
+    Array_double *old_row = pushed->data[row];
+    pushed->data[row] = add_element(old_row, v->data[row]);
+    free_vector(old_row);
+  }
+
+  pushed->cols++;
+  return pushed;
+}
+
+void free_matrix(Matrix_double *m) {
+  for (size_t y = 0; y < m->rows; ++y)
+    free_vector(m->data[y]);
+  free(m);
+}
+
+void format_matrix_into(Matrix_double *m, char *s) {
+  if (m->rows == 0)
+    strcpy(s, "empty");
+
+  for (size_t y = 0; y < m->rows; ++y) {
+    char row_s[5192];
+    strcpy(row_s, "");
+
+    format_vector_into(m->data[y], row_s);
+    strcat(s, row_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/Homework/math4610/src/rand.c b/Homework/math4610/src/rand.c
new file mode 100644
index 0000000..574a955
--- /dev/null
+++ b/Homework/math4610/src/rand.c
@@ -0,0 +1,5 @@
+#include "lizfcm.h"
+
+double rand_from(double min, double max) {
+  return min + (rand() / (RAND_MAX / (max - min)));
+}
diff --git a/Homework/math4610/src/roots.c b/Homework/math4610/src/roots.c
new file mode 100644
index 0000000..d6b22af
--- /dev/null
+++ b/Homework/math4610/src/roots.c
@@ -0,0 +1,127 @@
+#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_iterations]
+Array_double *find_ivt_range(double (*f)(double), double start_x, double delta,
+                             size_t max_iterations) {
+  double a = start_x;
+
+  while (f(a) * f(a + delta) >= 0 && max_iterations > 0) {
+    max_iterations--;
+    a += delta;
+  }
+
+  double end = a + delta;
+  double begin = a - delta;
+
+  if (max_iterations == 0 && f(begin) * f(end) >= 0)
+    return NULL;
+  return InitArray(double, {begin, end});
+}
+
+// f is continuous on [a, b]
+Array_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 InitArray(double, {a, b, 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));
+
+  Array_double *a_b_root = bisect_find_root(f, a, b, tolerance, max_iterations);
+  double root = a_b_root->data[2];
+  free_vector(a_b_root);
+
+  return root;
+}
+
+double fixed_point_iteration_method(double (*f)(double), double (*g)(double),
+                                    double x_0, double tolerance,
+                                    size_t max_iterations) {
+  if (max_iterations <= 0)
+    return x_0;
+
+  double root = g(x_0);
+  if (tolerance >= fabs(f(root)))
+    return root;
+
+  return fixed_point_iteration_method(f, g, root, tolerance,
+                                      max_iterations - 1);
+}
+
+double fixed_point_newton_method(double (*f)(double), double (*fprime)(double),
+                                 double x_0, double tolerance,
+                                 size_t max_iterations) {
+  if (max_iterations <= 0)
+    return x_0;
+
+  double root = x_0 - f(x_0) / fprime(x_0);
+  if (tolerance >= fabs(f(root)))
+    return root;
+
+  return fixed_point_newton_method(f, fprime, root, tolerance,
+                                   max_iterations - 1);
+}
+
+double fixed_point_secant_method(double (*f)(double), double x_0, double x_1,
+                                 double tolerance, size_t max_iterations) {
+  if (max_iterations == 0)
+    return x_1;
+
+  double root = x_1 - f(x_1) * ((x_1 - x_0) / (f(x_1) - f(x_0)));
+
+  if (tolerance >= fabs(f(root)))
+    return root;
+
+  return fixed_point_secant_method(f, x_1, root, tolerance, max_iterations - 1);
+}
+
+double fixed_point_secant_bisection_method(double (*f)(double), double x_0,
+                                           double x_1, double tolerance,
+                                           size_t max_iterations) {
+  double begin = x_0;
+  double end = x_1;
+  double root = x_0;
+
+  while (tolerance < fabs(f(root)) && max_iterations > 0) {
+    max_iterations--;
+
+    double secant_root = fixed_point_secant_method(f, begin, end, tolerance, 1);
+
+    if (secant_root < begin || secant_root > end) {
+      Array_double *range_root = bisect_find_root(f, begin, end, tolerance, 1);
+
+      begin = range_root->data[0];
+      end = range_root->data[1];
+      root = range_root->data[2];
+
+      free_vector(range_root);
+      continue;
+    }
+
+    root = secant_root;
+
+    if (f(root) * f(begin) < 0)
+      end = secant_root; // the root exists in [begin, secant_root]
+    else
+      begin = secant_root;
+  }
+
+  return root;
+}
diff --git a/Homework/math4610/src/vector.c b/Homework/math4610/src/vector.c
new file mode 100644
index 0000000..1b3e0b0
--- /dev/null
+++ b/Homework/math4610/src/vector.c
@@ -0,0 +1,143 @@
+#include "lizfcm.h"
+#include <assert.h>
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+
+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];
+  return sum;
+}
+
+Array_double *minus_v(Array_double *v1, Array_double *v2) {
+  assert(v1->size == v2->size);
+
+  Array_double *sub = InitArrayWithSize(double, v1->size, 0);
+  for (size_t i = 0; i < v1->size; i++)
+    sub->data[i] = v1->data[i] - v2->data[i];
+  return sub;
+}
+
+Array_double *scale_v(Array_double *v, double m) {
+  Array_double *copy = copy_vector(v);
+  for (size_t i = 0; i < v->size; i++)
+    copy->data[i] *= m;
+  return copy;
+}
+
+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 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;
+  for (size_t i = 0; i < v1->size; i++)
+    dot += v1->data[i] * v2->data[i];
+  return dot;
+}
+
+double vector_distance(Array_double *v1, Array_double *v2,
+                       double (*norm)(Array_double *)) {
+  Array_double *minus = minus_v(v1, v2);
+  double dist = (*norm)(minus);
+  free(minus);
+  return dist;
+}
+
+double l1_distance(Array_double *v1, Array_double *v2) {
+  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) {
+  return vector_distance(v1, v2, &linf_norm);
+}
+
+Array_double *copy_vector(Array_double *v) {
+  Array_double *copy = InitArrayWithSize(double, v->size, 0.0);
+  for (size_t i = 0; i < copy->size; ++i)
+    copy->data[i] = v->data[i];
+  return copy;
+}
+
+Array_double *add_element(Array_double *v, double x) {
+  Array_double *pushed = InitArrayWithSize(double, v->size + 1, 0.0);
+  for (size_t i = 0; i < v->size; ++i)
+    pushed->data[i] = v->data[i];
+  pushed->data[v->size] = x;
+  return pushed;
+}
+
+Array_double *slice_element(Array_double *v, size_t x) {
+  Array_double *sliced = InitArrayWithSize(double, v->size - 1, 0.0);
+  for (size_t i = 0; i < v->size - 1; ++i)
+    sliced->data[i] = i >= x ? v->data[i + 1] : v->data[i];
+  return sliced;
+}
+
+void free_vector(Array_double *v) {
+  free(v->data);
+  free(v);
+}
+
+void format_vector_into(Array_double *v, char *s) {
+  if (v->size == 0) {
+    strcat(s, "empty");
+    return;
+  }
+
+  for (size_t i = 0; i < v->size; ++i) {
+    char num[64];
+    strcpy(num, "");
+
+    sprintf(num, "%f,", v->data[i]);
+    strcat(s, num);
+  }
+  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;
+}
+
+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;
+}
-- 
cgit v1.3