hw 4

.gitignore

@@ -1,2 +1,3 @@
 build
 dist
+lib

Makefile

@@ -21,7 +21,7 @@ all: $(TEST_EXE)
 $(TEST_EXE): $(BIN_DIR) | $(LIBRARY)
 	$(CC) $(CPPFLAGS) $(CFLAGS) $(LDFLAGS) $(TEST_SRC) $(LIBRARY) -o $@
 
-$(LIBRARY): $(OBJ)
+$(LIBRARY): $(OBJ) | $(LIB_DIR)
 	ar rcs $(LIBRARY) $(OBJ_DIR)/*.o
 	ranlib $(LIBRARY)
 
@@ -32,6 +32,6 @@ $(BIN_DIR) $(OBJ_DIR) $(LIB_DIR):
 	mkdir -p $@
 
 clean:
-	@$(RM) -r $(BIN_DIR) $(OBJ_DIR)
+	@$(RM) -r $(BIN_DIR) $(OBJ_DIR) $(LIB_DIR)
 
 -include $(OBJ:.o=.d)

doc/software_manual.org

@@ -0,0 +1,727 @@
+#+TITLE: LIZFCM Software Manual (v0.1)
+#+AUTHOR: Elizabeth Hunt
+#+LATEX_HEADER: \notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+#+LATEX: \setlength\parindent{0pt}
+
+* Design
+The LIZFCM static library (at [[https://github.com/Simponic/math-4610][[https://github.com/Simponic/math-4610]] is a successor to my
+attempt at writing codes for the Fundamentals of Computational Mathematics course in Common
+Lisp, but the effort required to meet the requirement of creating a static library became
+too difficult to integrate outside of the ~ASDF~ solution that Common Lisp already brings
+to the table.
+
+All of the work established in ~deprecated-cl~ has been painstakingly translated into
+the C programming language. I have a couple tenets for its design:
+
++ Implemntations of routines should all be done immutably in respect to arguments.
++ Functional programming is good (it's... rough in C though).
++ Routines are separated into "module" c files, and not individual files per function.
+
+* Compilation
+A provided ~Makefile~ is added for convencience. It has been tested on an M1 machine running MacOS as
+well as Arch Linux.
+
+1. ~cd~ into the root of the repo
+2. ~make~
+   
+Then, as of homework 4, the testing routine in ~test/main.c~ can be run via
+~./dist/lizfcm.test~.
+
+Execution of the Makefile will perform compilation of individual routines.
+
+But, in the requirement of manual intervention (should the little alien workers
+inside the computer fail to do their job), one can use the following command to
+produce an object file:
+
+\begin{verbatim}
+  gcc -Iinc/ -lm -Wall -c src/<the_routine>.c -o build/<the_routine>.o
+\end{verbatim}
+
+Which is then bundled into a static library in ~lib/lizfcm.a~ which can be linked
+in the standard method.
+
+* The LIZFCM API
+** Simple Routines
+*** ~smaceps~
++ Author: Elizabeth Hunt
++ Name: ~smaceps~
++ Location: ~src/maceps.c~
++ Input: none
++ Output: a ~float~ returning the specific "Machine Epsilon" of a machine on a
+  single precision floating point number at which it becomes "indistinguishable".
+
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+
+*** ~dmaceps~
++ Author: Elizabeth Hunt
++ Name: ~dmaceps~
++ Location: ~src/maceps.c~
++ Input: none
++ Output: a ~double~ returning the specific "Machine Epsilon" of a machine on a
+  double precision floating point number at which it becomes "indistinguishable".
+
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+
+** Derivative Routines
+*** ~central_derivative_at~
++ Author: Elizabeth Hunt
++ Name: ~central_derivative_at~
++ Location: ~src/approx_derivative.c~
++ Input:
+  - ~f~ is a pointer to a one-ary function that takes a double as input and produces
+    a double as output
+  - ~a~ is the domain value at which we approximate ~f'~
+  - ~h~ is the step size
++ Output: a ~double~ of the approximate value of ~f'(a)~ via the central difference
+  method.
+
+#+BEGIN_SRC c
+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);
+}
+#+END_SRC
+
+*** ~forward_derivative_at~
++ Author: Elizabeth Hunt
++ Name: ~forward_derivative_at~
++ Location: ~src/approx_derivative.c~
++ Input:
+  - ~f~ is a pointer to a one-ary function that takes a double as input and produces
+    a double as output
+  - ~a~ is the domain value at which we approximate ~f'~
+  - ~h~ is the step size
++ Output: a ~double~ of the approximate value of ~f'(a)~ via the forward difference
+  method.
+
+#+BEGIN_SRC c
+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);
+}
+#+END_SRC
+
+*** ~backward_derivative_at~
++ Author: Elizabeth Hunt
++ Name: ~backward_derivative_at~
++ Location: ~src/approx_derivative.c~
++ Input:
+  - ~f~ is a pointer to a one-ary function that takes a double as input and produces
+    a double as output
+  - ~a~ is the domain value at which we approximate ~f'~
+  - ~h~ is the step size
++ Output: a ~double~ of the approximate value of ~f'(a)~ via the backward difference
+  method.
+
+#+BEGIN_SRC c
+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);
+}
+#+END_SRC
+
+** Vector Routines
+*** Vector Arithmetic: ~add_v, minus_v~
++ Author: Elizabeth Hunt
++ Name(s): ~add_v~, ~minus_v~
++ Location: ~src/vector.c~
++ Input: two pointers to locations in memory wherein ~Array_double~'s lie
++ Output: a pointer to a new ~Array_double~ as the result of addition or subtraction
+  of the two input ~Array_double~'s
+
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+
+*** Norms: ~l1_norm~, ~l2_norm~, ~linf_norm~
++ Author: Elizabeth Hunt
++ Name(s): ~l1_norm~, ~l2_norm~, ~linf_norm~
++ Location: ~src/vector.c~
++ Input: a pointer to a location in memory wherein an ~Array_double~ lies
++ Output: a ~double~ representing the value of the norm the function applies
+
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+
+*** ~vector_distance~
++ Author: Elizabeth Hunt
++ Name: ~vector_distance~
++ Location: ~src/vector.c~
++ Input: two pointers to locations in memory wherein ~Array_double~'s lie, and a pointer to a
+  one-ary function ~norm~ taking as input a pointer to an ~Array_double~ and returning a double
+  representing the norm of that ~Array_double~
+
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+
+*** Distances: ~l1_distance~, ~l2_distance~, ~linf_distance~
++ Author: Elizabeth Hunt
++ Name(s): ~l1_distance~, ~l2_distance~, ~linf_distance~
++ Location: ~src/vector.c~
++ Input: two pointers to locations in memory wherein ~Array_double~'s lie, and the distance
+  via the corresponding ~l1~, ~l2~, or ~linf~ norms
++ Output: A ~double~ representing the distance between the two ~Array_doubles~'s by the given
+  norm.
+  
+#+BEGIN_SRC c
+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);
+}
+#+END_SRC
+
+*** ~sum_v~
++ Author: Elizabeth Hunt
++ Name: ~sum_v~
++ Location: ~src/vector.c~
++ Input: a pointer to an ~Array_double~
++ Output: a ~double~ representing the sum of all the elements of an ~Array_double~
+
+#+BEGIN_SRC c
+double sum_v(Array_double *v) {
+  double sum = 0;
+  for (size_t i = 0; i < v->size; i++)
+    sum += v->data[i];
+  return sum;
+}
+#+END_SRC
+
+
+*** ~scale_v~
++ Author: Elizabeth Hunt
++ Name: ~scale_v~
++ Location: ~src/vector.c~
++ Input: a pointer to an ~Array_double~ and a scalar ~double~ to scale the vector
++ Output: a pointer to a new ~Array_double~ of the scaled input ~Array_double~
+
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+
+*** ~free_vector~
++ Author: Elizabeth Hunt
++ Name: ~free_vector~
++ Location: ~src/vector.c~
++ Input: a pointer to an ~Array_double~
++ Output: nothing.
++ Side effect: free the memory of the reserved ~Array_double~ on the heap
+
+#+BEGIN_SRC c
+void free_vector(Array_double *v) {
+  free(v->data);
+  free(v);
+}
+#+END_SRC
+
+*** ~copy_vector~
++ Author: Elizabeth Hunt
++ Name: ~copy_vector~
++ Location: ~src/vector.c~
++ Input: a pointer to an ~Array_double~
++ Output: a pointer to a new ~Array_double~ whose ~data~ and ~size~ are copied from the input
+ ~Array_double~
+
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+
+*** ~format_vector_into~
++ Author: Elizabeth Hunt
++ Name: ~format_vector_into~
++ Location: ~src/vector.c~
++ Input: a pointer to an ~Array_double~ and a pointer to a c-string ~s~ to "print" the vector out
+  into
++ Output: nothing.
++ Side effect: overwritten memory into ~s~
+
+#+BEGIN_SRC c
+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");
+}
+#+END_SRC
+
+** Matrix Routines
+*** ~lu_decomp~
++ Author: Elizabeth Hunt
++ Name: ~lu_decomp~
++ Location: ~src/matrix.c~
++ Input: a pointer to a ~Matrix_double~ $m$ to decompose into a lower triangular and upper triangular
+  matrix $L$, $U$, respectively such that $LU = m$.
++ Output: a pointer to the location in memory in which two ~Matrix_double~'s reside: the first
+  representing $L$, the second, $U$.
++ Errors: Exits and throws a status code of ~-1~ when encountering a matrix that cannot be
+  decomposed
+
+#+BEGIN_SRC c
+Matrix_double **lu_decomp(Matrix_double *m) {
+  assert(m->cols == m->rows);
+
+  Matrix_double *u = copy_matrix(m);
+  Matrix_double *l = InitMatrixWithSize(double, m->rows, m->cols, 0.0);
+  put_identity_diagonal(l);
+
+  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");
+      exit(-1);
+    }
+  }
+
+  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");
+          exit(-1);
+        }
+
+        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;
+}
+#+END_SRC
+*** ~bsubst~
++ Author: Elizabeth Hunt
++ Name: ~bsubst~
++ Location: ~src/matrix.c~
++ Input: a pointer to an upper-triangular ~Matrix_double~ $u$ and a ~Array_double~
+  $b$
++ Output: a pointer to a new ~Array_double~ whose entries are given by performing
+  back substitution
+
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+*** ~fsubst~
++ Author: Elizabeth Hunt
++ Name: ~fsubst~
++ Location: ~src/matrix.c~
++ Input: a pointer to a lower-triangular ~Matrix_double~ $l$ and a ~Array_double~
+  $b$
++ Output: a pointer to a new ~Array_double~ whose entries are given by performing
+  forward substitution
+
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+
+*** ~solve_matrix~
++ Author: Elizabeth Hunt
++ Location: ~src/matrix.c~
++ Input: a pointer to a ~Matrix_double~ $m$ and a pointer to an ~Array_double~ $b$
++ Output: $x$ such that $mx = b$ if such a solution exists (else it's non LU-factorable as discussed
+  above)
+
+Here we make use of forward substitution to first solve $Ly = b$ given $L$ as the $L$ factor in
+~lu_decomp~. Then we use back substitution to solve $Ux = y$ for $x$ similarly given $U$.
+
+Then, $LUx = b$, thus $x$ is a solution.
+
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+
+*** ~m_dot_v~
++ Author: Elizabeth Hunt
++ Location: ~src/matrix.c~
++ Input: a pointer to a ~Matrix_double~ $m$ and ~Array_double~ $v$
++ Output: the dot product $mv$ as an ~Array_double~
+  
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+
+*** ~put_identity_diagonal~
++ Author: Elizabeth Hunt
++ Location: ~src/matrix.c~
++ Input: a pointer to a ~Matrix_double~
++ Output: a pointer to a copy to ~Matrix_double~ whose diagonal is full of 1's
+
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+
+*** ~copy_matrix~
++ Author: Elizabeth Hunt
++ Location: ~src/matrix.c~
++ Input: a pointer to a ~Matrix_double~
++ Output: a pointer to a copy of the given ~Matrix_double~
+
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+  
+*** ~free_matrix~
++ Author: Elizabeth Hunt
++ Location: ~src/matrix.c~
++ Input: a pointer to a ~Matrix_double~
++ Output: none.
++ Side Effects: frees memory reserved by a given ~Matrix_double~ and its member
+ ~Array_double~ vectors describing its rows.
+
+#+BEGIN_SRC c
+void free_matrix(Matrix_double *m) {
+  for (size_t y = 0; y < m->rows; ++y)
+    free_vector(m->data[y]);
+  free(m);
+}
+#+END_SRC
+
+*** ~format_matrix_into~
++ Author: Elizabeth Hunt
++ Name: ~format_matrix_into~
++ Location: ~src/matrix.c~
++ Input: a pointer to a ~Matrix_double~ and a pointer to a c-string ~s~ to "print" the vector out
+  into
++ Output: nothing.
++ Side effect: overwritten memory into ~s~
+
+#+BEGIN_SRC c
+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[256];
+    strcpy(row_s, "");
+
+    format_vector_into(m->data[y], row_s);
+    strcat(s, row_s);
+  }
+  strcat(s, "\n");
+}
+#+END_SRC
+** Linear Routines
+*** ~least_squares_lin_reg~
++ Author: Elizabeth Hunt
++ Name: ~least_squares_lin_reg~
++ Location: ~src/lin.c~
++ Input: two pointers to ~Array_double~'s whose entries correspond two ordered
+  pairs in R^2
++ Output: a linear model best representing the ordered pairs via least squares
+  regression
+  
+#+BEGIN_SRC c
+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;
+}
+#+END_SRC
+** Appendix / Miscellaneous
+*** Data Types
+**** ~Line~
++ Author: Elizabeth Hunt
++ Location: ~inc/types.h~
+
+#+BEGIN_SRC c
+typedef struct Line {
+  double m;
+  double a;
+} Line;
+#+END_SRC
+**** The ~Array_<type>~ and ~Matrix_<type>~
++ Author: Elizabeth Hunt
++ Location: ~inc/types.h~
+
+We define two Pre processor Macros ~DEFINE_ARRAY~ and ~DEFINE_MATRIX~ that take
+as input a type, and construct a struct definition for the given type for
+convenient access to the vector or matrices dimensions.
+
+Such that ~DEFINE_ARRAY(int)~ would expand to:
+
+#+BEGIN_SRC c
+  typedef struct {
+    int* data;
+    size_t size;
+  } Array_int
+#+END_SRC
+
+And ~DEFINE_MATRIX(int)~ would expand a to ~Matrix_int~; containing a pointer to
+a collection of pointers of ~Array_int~'s and its dimensions.
+
+#+BEGIN_SRC c
+  typedef struct {
+    Array_int **data;
+    size_t cols;
+    size_t rows;
+  } Matrix_int
+#+END_SRC
+
+*** Macros
+**** ~c_max~ and ~c_min~
++ Author: Elizabeth Hunt
++ Location: ~inc/macros.h~
++ Input: two structures that define an order measure
++ Output: either the larger or smaller of the two depending on the measure
+
+#+BEGIN_SRC c
+#define c_max(x, y) (((x) >= (y)) ? (x) : (y))
+#define c_min(x, y) (((x) <= (y)) ? (x) : (y))
+#+END_SRC
+
+**** ~InitArray~
++ Author: Elizabeth Hunt
++ Location: ~inc/macros.h~
++ Input: a type and array of values to initialze an array with such type
++ Output: a new ~Array_type~ with the size of the given array and its data
+
+#+BEGIN_SRC c
+#define InitArray(TYPE, ...)                                                   \
+  ({                                                                           \
+    TYPE temp[] = __VA_ARGS__;                                                 \
+    Array_##TYPE *arr = malloc(sizeof(Array_##TYPE));                          \
+    arr->size = sizeof(temp) / sizeof(temp[0]);                                \
+    arr->data = malloc(arr->size * sizeof(TYPE));                              \
+    memcpy(arr->data, temp, arr->size * sizeof(TYPE));                         \
+    arr;                                                                       \
+  })
+#+END_SRC
+
+**** ~InitArrayWithSize~
++ Author: Elizabeth Hunt
++ Location: ~inc/macros.h~
++ Input: a type, a size, and initial value
++ Output: a new ~Array_type~ with the given size filled with the initial value
+
+#+BEGIN_SRC c
+#define InitArrayWithSize(TYPE, SIZE, INIT_VALUE)                              \
+  ({                                                                           \
+    Array_##TYPE *arr = malloc(sizeof(Array_##TYPE));                          \
+    arr->size = SIZE;                                                          \
+    arr->data = malloc(arr->size * sizeof(TYPE));                              \
+    for (size_t i = 0; i < arr->size; i++)                                     \
+      arr->data[i] = INIT_VALUE;                                               \
+    arr;                                                                       \
+  })
+#+END_SRC
+
+**** ~InitMatrixWithSize~
++ Author: Elizabeth Hunt
++ Location: ~inc/macros.h~
++ Input: a type, number of rows, columns, and initial value
++ Output: a new ~Matrix_type~ of size ~rows x columns~ filled with the initial
+  value
+
+#+BEGIN_SRC c
+#define InitMatrixWithSize(TYPE, ROWS, COLS, INIT_VALUE)                       \
+  ({                                                                           \
+    Matrix_##TYPE *matrix = malloc(sizeof(Matrix_##TYPE));                     \
+    matrix->rows = ROWS;                                                       \
+    matrix->cols = COLS;                                                       \
+    matrix->data = malloc(matrix->rows * sizeof(Array_##TYPE *));              \
+    for (size_t y = 0; y < matrix->rows; y++)                                  \
+      matrix->data[y] = InitArrayWithSize(TYPE, COLS, INIT_VALUE);             \
+    matrix;                                                                    \
+  })
+#+END_SRc
+

doc/software_manual.tex

@@ -0,0 +1,867 @@
+% Created 2023-10-13 Fri 20:48
+% Intended LaTeX compiler: pdflatex
+\documentclass[11pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{graphicx}
+\usepackage{longtable}
+\usepackage{wrapfig}
+\usepackage{rotating}
+\usepackage[normalem]{ulem}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{capt-of}
+\usepackage{hyperref}
+\notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+\author{Elizabeth Hunt}
+\date{\today}
+\title{LIZFCM Software Manual (v0.1)}
+\hypersetup{
+ pdfauthor={Elizabeth Hunt},
+ pdftitle={LIZFCM Software Manual (v0.1)},
+ pdfkeywords={},
+ pdfsubject={},
+ pdfcreator={Emacs 28.2 (Org mode 9.7-pre)}, 
+ pdflang={English}}
+\begin{document}
+
+\maketitle
+\tableofcontents
+
+\setlength\parindent{0pt}
+
+\section{Design}
+\label{sec:org23cc15b}
+The LIZFCM static library is a successor to my attempt at writing codes for the
+Fundamentals of Computational Mathematics course in Common Lisp, but the effort required
+to meet the requirement of creating a static library became too difficult to integrate
+outside of the \texttt{ASDF} solution that Common Lisp already brings to the table.
+
+All of the work established in \texttt{deprecated-cl} has been painstakingly translated into
+the C programming language. I have a couple tenets for its design:
+
+\begin{itemize}
+\item Implemntations of routines should all be done immutably in respect to arguments.
+\item Functional programming is good (it's\ldots{} rough in this language though).
+\end{itemize}
+
+\section{Compilation}
+\label{sec:orgc704fb9}
+A provided \texttt{Makefile} is added for convencience. It has been tested on an M1 machine running MacOS as
+well as Arch Linux.
+
+\begin{enumerate}
+\item \texttt{cd} into the root of the repo
+\item \texttt{make}
+\end{enumerate}
+
+Then, as of homework 4, the testing routine in \texttt{test/main.c} can be run via
+\texttt{./dist/lizfcm.test}.
+
+Execution of the Makefile will perform compilation of individual routines.
+
+But, in the requirement of manual intervention (should the little alien workers
+inside the computer fail to do their job), one can use the following command to
+produce an object file:
+
+\begin{verbatim}
+  gcc -Iinc/ -lm -Wall -c src/<the_routine>.c -o build/<the_routine>.o
+\end{verbatim}
+
+Which is then bundled into a static library in \texttt{lib/lizfcm.a} which can be linked
+in the standard method.
+
+\section{The LIZFCM API}
+\label{sec:org01b31c2}
+\subsection{Simple Routines}
+\label{sec:org5355145}
+\subsubsection{\texttt{smaceps}}
+\label{sec:org4ed063f}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{smaceps}
+\item Location: \texttt{src/maceps.c}
+\item Input: none
+\item Output: a \texttt{float} returning the specific "Machine Epsilon" of a machine on a
+single precision floating point number at which it becomes "indistinguishable".
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+
+\subsubsection{\texttt{dmaceps}}
+\label{sec:orgcd0dfff}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{dmaceps}
+\item Location: \texttt{src/maceps.c}
+\item Input: none
+\item Output: a \texttt{double} returning the specific "Machine Epsilon" of a machine on a
+double precision floating point number at which it becomes "indistinguishable".
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+
+\subsection{Derivative Routines}
+\label{sec:org8f54012}
+\subsubsection{\texttt{central\_derivative\_at}}
+\label{sec:org2c81fc1}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{central\_derivative\_at}
+\item Location: \texttt{src/approx\_derivative.c}
+\item Input:
+\begin{itemize}
+\item \texttt{f} is a pointer to a one-ary function that takes a double as input and produces
+a double as output
+\item \texttt{a} is the domain value at which we approximate \texttt{f'}
+\item \texttt{h} is the step size
+\end{itemize}
+\item Output: a \texttt{double} of the approximate value of \texttt{f'(a)} via the central difference
+method.
+\end{itemize}
+
+\begin{verbatim}
+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);
+}
+\end{verbatim}
+
+\subsubsection{\texttt{forward\_derivative\_at}}
+\label{sec:org149b09e}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{forward\_derivative\_at}
+\item Location: \texttt{src/approx\_derivative.c}
+\item Input:
+\begin{itemize}
+\item \texttt{f} is a pointer to a one-ary function that takes a double as input and produces
+a double as output
+\item \texttt{a} is the domain value at which we approximate \texttt{f'}
+\item \texttt{h} is the step size
+\end{itemize}
+\item Output: a \texttt{double} of the approximate value of \texttt{f'(a)} via the forward difference
+method.
+\end{itemize}
+
+\begin{verbatim}
+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);
+}
+\end{verbatim}
+
+\subsubsection{\texttt{backward\_derivative\_at}}
+\label{sec:orgbaa2238}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{backward\_derivative\_at}
+\item Location: \texttt{src/approx\_derivative.c}
+\item Input:
+\begin{itemize}
+\item \texttt{f} is a pointer to a one-ary function that takes a double as input and produces
+a double as output
+\item \texttt{a} is the domain value at which we approximate \texttt{f'}
+\item \texttt{h} is the step size
+\end{itemize}
+\item Output: a \texttt{double} of the approximate value of \texttt{f'(a)} via the backward difference
+method.
+\end{itemize}
+
+\begin{verbatim}
+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);
+}
+\end{verbatim}
+
+\subsection{Vector Routines}
+\label{sec:orgf9dd708}
+\subsubsection{Vector Arithmetic: \texttt{add\_v, minus\_v}}
+\label{sec:orgbb91d9d}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name(s): \texttt{add\_v}, \texttt{minus\_v}
+\item Location: \texttt{src/vector.c}
+\item Input: two pointers to locations in memory wherein \texttt{Array\_double}'s lie
+\item Output: a pointer to a new \texttt{Array\_double} as the result of addition or subtraction
+of the two input \texttt{Array\_double}'s
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+
+\subsubsection{Norms: \texttt{l1\_norm}, \texttt{l2\_norm}, \texttt{linf\_norm}}
+\label{sec:org88bb4c5}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name(s): \texttt{l1\_norm}, \texttt{l2\_norm}, \texttt{linf\_norm}
+\item Location: \texttt{src/vector.c}
+\item Input: a pointer to a location in memory wherein an \texttt{Array\_double} lies
+\item Output: a \texttt{double} representing the value of the norm the function applies
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+
+\subsubsection{\texttt{vector\_distance}}
+\label{sec:org4499de6}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{vector\_distance}
+\item Location: \texttt{src/vector.c}
+\item Input: two pointers to locations in memory wherein \texttt{Array\_double}'s lie, and a pointer to a
+one-ary function \texttt{norm} taking as input a pointer to an \texttt{Array\_double} and returning a double
+representing the norm of that \texttt{Array\_double}
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+
+\subsubsection{Distances: \texttt{l1\_distance}, \texttt{l2\_distance}, \texttt{linf\_distance}}
+\label{sec:org1e61aea}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name(s): \texttt{l1\_distance}, \texttt{l2\_distance}, \texttt{linf\_distance}
+\item Location: \texttt{src/vector.c}
+\item Input: two pointers to locations in memory wherein \texttt{Array\_double}'s lie, and the distance
+via the corresponding \texttt{l1}, \texttt{l2}, or \texttt{linf} norms
+\item Output: A \texttt{double} representing the distance between the two \texttt{Array\_doubles}'s by the given
+norm.
+\end{itemize}
+
+\begin{verbatim}
+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);
+}
+\end{verbatim}
+
+\subsubsection{\texttt{sum\_v}}
+\label{sec:org57e2591}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{sum\_v}
+\item Location: \texttt{src/vector.c}
+\item Input: a pointer to an \texttt{Array\_double}
+\item Output: a \texttt{double} representing the sum of all the elements of an \texttt{Array\_double}
+\end{itemize}
+
+\begin{verbatim}
+double sum_v(Array_double *v) {
+  double sum = 0;
+  for (size_t i = 0; i < v->size; i++)
+    sum += v->data[i];
+  return sum;
+}
+\end{verbatim}
+
+
+\subsubsection{\texttt{scale\_v}}
+\label{sec:org61b466a}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{scale\_v}
+\item Location: \texttt{src/vector.c}
+\item Input: a pointer to an \texttt{Array\_double} and a scalar \texttt{double} to scale the vector
+\item Output: a pointer to a new \texttt{Array\_double} of the scaled input \texttt{Array\_double}
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+
+\subsubsection{\texttt{free\_vector}}
+\label{sec:org398d778}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{free\_vector}
+\item Location: \texttt{src/vector.c}
+\item Input: a pointer to an \texttt{Array\_double}
+\item Output: nothing.
+\item Side effect: free the memory of the reserved \texttt{Array\_double} on the heap
+\end{itemize}
+
+\begin{verbatim}
+void free_vector(Array_double *v) {
+  free(v->data);
+  free(v);
+}
+\end{verbatim}
+
+\subsubsection{\texttt{copy\_vector}}
+\label{sec:orgf6b116b}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{copy\_vector}
+\item Location: \texttt{src/vector.c}
+\item Input: a pointer to an \texttt{Array\_double}
+\item Output: a pointer to a new \texttt{Array\_double} whose \texttt{data} and \texttt{size} are copied from the input
+\texttt{Array\_double}
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+
+\subsubsection{\texttt{format\_vector\_into}}
+\label{sec:org595519d}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{format\_vector\_into}
+\item Location: \texttt{src/vector.c}
+\item Input: a pointer to an \texttt{Array\_double} and a pointer to a c-string \texttt{s} to "print" the vector out
+into
+\item Output: nothing.
+\item Side effect: overwritten memory into \texttt{s}
+\end{itemize}
+
+\begin{verbatim}
+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");
+}
+\end{verbatim}
+
+\subsection{Matrix Routines}
+\label{sec:org53505d6}
+\subsubsection{\texttt{lu\_decomp}}
+\label{sec:org22ad28d}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{lu\_decomp}
+\item Location: \texttt{src/matrix.c}
+\item Input: a pointer to a \texttt{Matrix\_double} \(m\) to decompose into a lower triangular and upper triangular
+matrix \(L\), \(U\), respectively such that \(LU = m\).
+\item Output: a pointer to the location in memory in which two \texttt{Matrix\_double}'s reside: the first
+representing \(L\), the second, \(U\).
+\item Errors: Exits and throws a status code of \texttt{-1} when encountering a matrix that cannot be
+decomposed
+\end{itemize}
+
+\begin{verbatim}
+Matrix_double **lu_decomp(Matrix_double *m) {
+  assert(m->cols == m->rows);
+
+  Matrix_double *u = copy_matrix(m);
+  Matrix_double *l = InitMatrixWithSize(double, m->rows, m->cols, 0.0);
+  put_identity_diagonal(l);
+
+  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");
+      exit(-1);
+    }
+  }
+
+  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");
+          exit(-1);
+        }
+
+        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;
+}
+\end{verbatim}
+\subsubsection{\texttt{bsubst}}
+\label{sec:org15fec98}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{bsubst}
+\item Location: \texttt{src/matrix.c}
+\item Input: a pointer to an upper-triangular \texttt{Matrix\_double} \(u\) and a \texttt{Array\_double}
+\(b\)
+\item Output: a pointer to a new \texttt{Array\_double} whose entries are given by performing
+back substitution
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+\subsubsection{\texttt{fsubst}}
+\label{sec:orgdeab27c}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{fsubst}
+\item Location: \texttt{src/matrix.c}
+\item Input: a pointer to a lower-triangular \texttt{Matrix\_double} \(l\) and a \texttt{Array\_double}
+\(b\)
+\item Output: a pointer to a new \texttt{Array\_double} whose entries are given by performing
+forward substitution
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+
+\subsubsection{\texttt{solve\_matrix}}
+\label{sec:orge57c26b}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Location: \texttt{src/matrix.c}
+\item Input: a pointer to a \texttt{Matrix\_double} \(m\) and a pointer to an \texttt{Array\_double} \(b\)
+\item Output: \(x\) such that \(mx = b\) if such a solution exists (else it's non LU-factorable as discussed
+above)
+\end{itemize}
+
+Here we make use of forward substitution to first solve \(Ly = b\) given \(L\) as the \(L\) factor in
+\texttt{lu\_decomp}. Then we use back substitution to solve \(Ux = y\) for \(x\) similarly given \(U\).
+
+Then, \(LUx = b\), thus \(x\) is a solution.
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+
+\subsubsection{\texttt{m\_dot\_v}}
+\label{sec:org6afa7d5}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Location: \texttt{src/matrix.c}
+\item Input: a pointer to a \texttt{Matrix\_double} \(m\) and \texttt{Array\_double} \(v\)
+\item Output: the dot product \(mv\) as an \texttt{Array\_double}
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+
+\subsubsection{\texttt{put\_identity\_diagonal}}
+\label{sec:orgdd1c373}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Location: \texttt{src/matrix.c}
+\item Input: a pointer to a \texttt{Matrix\_double}
+\item Output: a pointer to a copy to \texttt{Matrix\_double} whose diagonal is full of 1's
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+
+\subsubsection{\texttt{copy\_matrix}}
+\label{sec:org3d1b7b0}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Location: \texttt{src/matrix.c}
+\item Input: a pointer to a \texttt{Matrix\_double}
+\item Output: a pointer to a copy of the given \texttt{Matrix\_double}
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+
+\subsubsection{\texttt{free\_matrix}}
+\label{sec:org697f6cc}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Location: \texttt{src/matrix.c}
+\item Input: a pointer to a \texttt{Matrix\_double}
+\item Output: none.
+\item Side Effects: frees memory reserved by a given \texttt{Matrix\_double} and its member
+\texttt{Array\_double} vectors describing its rows.
+\end{itemize}
+
+\begin{verbatim}
+void free_matrix(Matrix_double *m) {
+  for (size_t y = 0; y < m->rows; ++y)
+    free_vector(m->data[y]);
+  free(m);
+}
+\end{verbatim}
+
+\subsubsection{\texttt{format\_matrix\_into}}
+\label{sec:orgc43bda3}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{format\_matrix\_into}
+\item Location: \texttt{src/matrix.c}
+\item Input: a pointer to a \texttt{Matrix\_double} and a pointer to a c-string \texttt{s} to "print" the vector out
+into
+\item Output: nothing.
+\item Side effect: overwritten memory into \texttt{s}
+\end{itemize}
+
+\begin{verbatim}
+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[256];
+    strcpy(row_s, "");
+
+    format_vector_into(m->data[y], row_s);
+    strcat(s, row_s);
+  }
+  strcat(s, "\n");
+}
+\end{verbatim}
+\subsection{Linear Routines}
+\label{sec:org1e850f2}
+\subsubsection{\texttt{least\_squares\_lin\_reg}}
+\label{sec:org02e6d37}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Name: \texttt{least\_squares\_lin\_reg}
+\item Location: \texttt{src/lin.c}
+\item Input: two pointers to \texttt{Array\_double}'s whose entries correspond two ordered
+pairs in R\textsuperscript{2}
+\item Output: a linear model best representing the ordered pairs via least squares
+regression
+\end{itemize}
+
+\begin{verbatim}
+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;
+}
+\end{verbatim}
+\subsection{Appendix / Miscellaneous}
+\label{sec:org83c0f8d}
+\subsubsection{Data Types}
+\label{sec:org22f30f4}
+\begin{enumerate}
+\item \texttt{Line}
+\label{sec:orgd014841}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Location: \texttt{inc/types.h}
+\end{itemize}
+
+\begin{verbatim}
+typedef struct Line {
+  double m;
+  double a;
+} Line;
+\end{verbatim}
+\item The \texttt{Array\_<type>} and \texttt{Matrix\_<type>}
+\label{sec:org3f90e03}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Location: \texttt{inc/types.h}
+\end{itemize}
+
+We define two Pre processor Macros \texttt{DEFINE\_ARRAY} and \texttt{DEFINE\_MATRIX} that take
+as input a type, and construct a struct definition for the given type for
+convenient access to the vector or matrices dimensions.
+
+Such that \texttt{DEFINE\_ARRAY(int)} would expand to:
+
+\begin{verbatim}
+typedef struct {
+  int* data;
+  size_t size;
+} Array_int
+\end{verbatim}
+
+And \texttt{DEFINE\_MATRIX(int)} would expand a to \texttt{Matrix\_int}; containing a pointer to
+a collection of pointers of \texttt{Array\_int}'s and its dimensions.
+
+\begin{verbatim}
+typedef struct {
+  Array_int **data;
+  size_t cols;
+  size_t rows;
+} Matrix_int
+\end{verbatim}
+\end{enumerate}
+
+\subsubsection{Macros}
+\label{sec:org60b549e}
+\begin{enumerate}
+\item \texttt{c\_max} and \texttt{c\_min}
+\label{sec:org04ff2db}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Location: \texttt{inc/macros.h}
+\item Input: two structures that define an order measure
+\item Output: either the larger or smaller of the two depending on the measure
+\end{itemize}
+
+\begin{verbatim}
+#define c_max(x, y) (((x) >= (y)) ? (x) : (y))
+#define c_min(x, y) (((x) <= (y)) ? (x) : (y))
+\end{verbatim}
+
+\item \texttt{InitArray}
+\label{sec:orgf67f153}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Location: \texttt{inc/macros.h}
+\item Input: a type and array of values to initialze an array with such type
+\item Output: a new \texttt{Array\_type} with the size of the given array and its data
+\end{itemize}
+
+\begin{verbatim}
+#define InitArray(TYPE, ...)                                                   \
+  ({                                                                           \
+    TYPE temp[] = __VA_ARGS__;                                                 \
+    Array_##TYPE *arr = malloc(sizeof(Array_##TYPE));                          \
+    arr->size = sizeof(temp) / sizeof(temp[0]);                                \
+    arr->data = malloc(arr->size * sizeof(TYPE));                              \
+    memcpy(arr->data, temp, arr->size * sizeof(TYPE));                         \
+    arr;                                                                       \
+  })
+\end{verbatim}
+
+\item \texttt{InitArrayWithSize}
+\label{sec:org47e5e66}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Location: \texttt{inc/macros.h}
+\item Input: a type, a size, and initial value
+\item Output: a new \texttt{Array\_type} with the given size filled with the initial value
+\end{itemize}
+
+\begin{verbatim}
+#define InitArrayWithSize(TYPE, SIZE, INIT_VALUE)                              \
+  ({                                                                           \
+    Array_##TYPE *arr = malloc(sizeof(Array_##TYPE));                          \
+    arr->size = SIZE;                                                          \
+    arr->data = malloc(arr->size * sizeof(TYPE));                              \
+    for (size_t i = 0; i < arr->size; i++)                                     \
+      arr->data[i] = INIT_VALUE;                                               \
+    arr;                                                                       \
+  })
+\end{verbatim}
+
+\item \texttt{InitMatrixWithSize}
+\label{sec:org3b96b75}
+\begin{itemize}
+\item Author: Elizabeth Hunt
+\item Location: \texttt{inc/macros.h}
+\item Input: a type, number of rows, columns, and initial value
+\item Output: a new \texttt{Matrix\_type} of size \texttt{rows x columns} filled with the initial
+value
+\end{itemize}
+
+\begin{verbatim}
+#define InitMatrixWithSize(TYPE, ROWS, COLS, INIT_VALUE)                       \
+  ({                                                                           \
+    Matrix_##TYPE *matrix = malloc(sizeof(Matrix_##TYPE));                     \
+    matrix->rows = ROWS;                                                       \
+    matrix->cols = COLS;                                                       \
+    matrix->data = malloc(matrix->rows * sizeof(Array_##TYPE *));              \
+    for (size_t y = 0; y < matrix->rows; y++)                                  \
+      matrix->data[y] = InitArrayWithSize(TYPE, COLS, INIT_VALUE);             \
+    matrix;                                                                    \
+  })
+\end{verbatim}
+\end{enumerate}
+\end{document}

homeworks/hw-4.org

@@ -1,10 +1,34 @@
-* My Basic Routines
-* Linear Solution Routines
-** Gaussian Elimination with Backwards Substitution
-** LU factorization
-* Matrix Routines
-** dotproduct
-** crossproduct
-** matrixvector
-** matrixmatrix
+#+TITLE: Homework 4
+#+AUTHOR: Elizabeth Hunt
+#+LATEX_HEADER: \notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+#+LATEX: \setlength\parindent{0pt}
+#+OPTIONS: toc:nil
 
+* Question 1
+See the attached LIZFCM Software Manual.
+
+* Question 2, 3, 4
+#+attr_latex: :width 350px
+[[./img/make_run.png]]
+
+* Question 5
+#+attr_latex: :width 350px
+[[./img/test_routine_1.png]]
+
+#+attr_latex: :width 350px
+[[./img/test_routine_2.png]]
+
+* Question 6
+See the LIZFCM Software Manual.
+
+* Question 7
+See ~src/matrix.c -> lu_decomp, fsubst, bsubst, solve_matrix~
+
+* Question 8
+See ~test/main.c -> lines 109 - 113~ in correspondence to the run in Question 5
+
+* Question 9
+See ~test/main.c -> lines 118 - 121~ in correspondence to the run in Question 5
+
+* Question 10
+See the TOC on the first page of the LIZFCM Software Manual.

homeworks/hw-4.tex

@@ -0,0 +1,70 @@
+% Created 2023-10-13 Fri 20:54
+% Intended LaTeX compiler: pdflatex
+\documentclass[11pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{graphicx}
+\usepackage{longtable}
+\usepackage{wrapfig}
+\usepackage{rotating}
+\usepackage[normalem]{ulem}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{capt-of}
+\usepackage{hyperref}
+\notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+\author{Elizabeth Hunt}
+\date{\today}
+\title{Homework 4}
+\hypersetup{
+ pdfauthor={Elizabeth Hunt},
+ pdftitle={Homework 4},
+ pdfkeywords={},
+ pdfsubject={},
+ pdfcreator={Emacs 28.2 (Org mode 9.7-pre)}, 
+ pdflang={English}}
+\begin{document}
+
+\maketitle
+\setlength\parindent{0pt}
+
+\section{Question 1}
+\label{sec:org0ce46da}
+See the attached LIZFCM Software Manual.
+
+\section{Question 2, 3, 4}
+\label{sec:org0794d50}
+\begin{center}
+\includegraphics[width=350px]{./img/make_run.png}
+\end{center}
+
+\section{Question 5}
+\label{sec:orgf7b84cc}
+\begin{center}
+\includegraphics[width=350px]{./img/test_routine_1.png}
+\end{center}
+
+\begin{center}
+\includegraphics[width=350px]{./img/test_routine_2.png}
+\end{center}
+
+\section{Question 6}
+\label{sec:orgee1884d}
+See the LIZFCM Software Manual.
+
+\section{Question 7}
+\label{sec:org6f86e09}
+See \texttt{src/matrix.c -> lu\_decomp, fsubst, bsubst, solve\_matrix}
+
+\section{Question 8}
+\label{sec:orga79c31e}
+See \texttt{test/main.c -> lines 109 - 113} in correspondence to the run in Question 5
+
+\section{Question 9}
+\label{sec:org5b1b74d}
+See \texttt{test/main.c -> lines 118 - 121} in correspondence to the run in Question 5
+
+\section{Question 10}
+\label{sec:org7409683}
+See the TOC on the first page of the LIZFCM Software Manual.
+\end{document}

inc/lizfcm.h

@@ -15,25 +15,29 @@ extern double sum_v(Array_double *v);
 extern Array_double *scale_v(Array_double *v1, double m);
 extern Array_double *add_v(Array_double *v1, Array_double *v2);
 extern Array_double *minus_v(Array_double *v1, Array_double *v2);
-extern double dot_v(Array_double *v1, Array_double *v2);
+extern double v_dot_v(Array_double *v1, Array_double *v2);
 extern double l2_norm(Array_double *v);
 extern double l1_norm(Array_double *v);
 extern double linf_norm(Array_double *v);
+extern double vector_distance(Array_double *v1, Array_double *v2,
+                              double (*norm)(Array_double *));
 extern double l2_distance(Array_double *v1, Array_double *v2);
 extern double l1_distance(Array_double *v1, Array_double *v2);
 extern double linf_distance(Array_double *v1, Array_double *v2);
 extern Array_double *copy_vector(Array_double *v1);
 extern void free_vector(Array_double *v);
 extern void format_vector_into(Array_double *v, char *s);
-extern Line *least_squares_lin_reg(Array_double *x, Array_double *y);
 
-extern void put_identity_diagonal(Matrix_double *m);
-extern Matrix_double *copy_matrix(Matrix_double *m);
-extern void free_matrix(Matrix_double *m);
+extern Matrix_double *put_identity_diagonal(Matrix_double *m);
 extern Matrix_double **lu_decomp(Matrix_double *m);
 extern Array_double *bsubst(Matrix_double *u, Array_double *b);
 extern Array_double *fsubst(Matrix_double *l, Array_double *b);
-
+extern Array_double *solve_matrix(Matrix_double *m, Array_double *b);
+extern Array_double *m_dot_v(Matrix_double *m, Array_double *v);
+extern Matrix_double *copy_matrix(Matrix_double *m);
+extern void free_matrix(Matrix_double *m);
 extern void format_matrix_into(Matrix_double *m, char *s);
 
+extern Line *least_squares_lin_reg(Array_double *x, Array_double *y);
+
 #endif // LIZFCM_H

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));

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]);

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;
+}

test/main.c

@@ -14,12 +14,12 @@ int main() {
   printf("========\n");
 
   printf("Norm, Distance\n");
-  Array_double *v = InitArray(double, {3, 1, -4, 1, 5, -9, 3});
+  Array_double *v = InitArray(double, {3, 1, -4});
   strcpy(s, "");
   format_vector_into(v, s);
   printf("v: %s", s);
 
-  Array_double *w = InitArray(double, {-2, 7, 1, -8, -2, 8, 5});
+  Array_double *w = InitArray(double, {-2, 7, 1});
   strcpy(s, "");
   format_vector_into(w, s);
   printf("w: %s", s);
@@ -51,7 +51,6 @@ int main() {
   printf("w: %s", s);
 
   Line *line = least_squares_lin_reg(v, w);
-
   printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a);
   v = InitArray(double, {1, 2, 3, 4, 5, 6, 7});
   strcpy(s, "");
@@ -64,6 +63,7 @@ int main() {
 
   line = least_squares_lin_reg(v, w);
   printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a);
+  free(line);
   printf("========\n");
 
   printf("LU Decomp\n");
@@ -77,8 +77,9 @@ int main() {
   format_matrix_into(a, s);
   printf("a = %s", s);
 
-  uint32_t solution = 100;
-  Array_double *b = InitArrayWithSize(double, n, (double)solution);
+  uint32_t solution = 1;
+  Array_double *y = InitArrayWithSize(double, n, (double)solution);
+  Array_double *b = m_dot_v(a, y); // q8 steps
   Matrix_double **u_l = lu_decomp(a);
   Matrix_double *u = u_l[0];
   Matrix_double *l = u_l[1];
@@ -91,12 +92,11 @@ int main() {
   printf("l = %s", s);
   strcpy(s, "");
   format_vector_into(b, s);
-  printf("b = %s", s);
+  printf("(after following q8) b = %s", s);
   printf("========\n");
   printf("Forward / Backward Substitution Solution to ax=b\n");
 
   Array_double *b_fsub = fsubst(l, b);
-  free_vector(b);
   strcpy(s, "");
   format_vector_into(b_fsub, s);
   printf("b_fsub: %s", s);
@@ -106,16 +106,19 @@ int main() {
   format_vector_into(x_bsub, s);
   printf("x_bsub: %s", s);
 
+  Array_double *x_solve_matrix = solve_matrix(a, b);
+  free_vector(b);
+  strcpy(s, "");
+  format_vector_into(x_solve_matrix, s);
+  printf("\\--> == x_solve_matrix: %s", s);
+
   free_vector(b_fsub);
+  free_vector(x_solve_matrix);
 
   printf("Verifications\n");
-  for (size_t row = 0; row < a->rows; row++) {
-    double curr = 0;
-    for (size_t col = 0; col < a->cols; col++)
-      curr += a->data[row]->data[col] * x_bsub->data[col];
-    printf("Substituions for values in row %zu = %f, true value err=%.10e\n",
-           row, curr, fabs(curr - solution));
-  }
+  for (size_t i = 0; i < x_bsub->size; i++)
+    printf("in row %zu, solution = %f, true value err=%.10e\n", i,
+           x_bsub->data[i], fabs(x_bsub->data[i] - solution));
 
   return 0;
 }