summaryrefslogtreecommitdiff
path: root/Homework/math4610/src
diff options
context:
space:
mode:
Diffstat (limited to 'Homework/math4610/src')
-rw-r--r--Homework/math4610/src/approx_derivative.c38
-rw-r--r--Homework/math4610/src/eigen.c116
-rw-r--r--Homework/math4610/src/lin.c19
-rw-r--r--Homework/math4610/src/maceps.c28
-rw-r--r--Homework/math4610/src/matrix.c346
-rw-r--r--Homework/math4610/src/rand.c5
-rw-r--r--Homework/math4610/src/roots.c127
-rw-r--r--Homework/math4610/src/vector.c143
8 files changed, 822 insertions, 0 deletions
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;
+}
generated by cgit v1.3 (git 2.53.0) at 2026-07-18 15:45:29 +0000