lizfcm/doc/software_manual.org

LIZFCM Software Manual (v0.2)

• Design
• Compilation
• The LIZFCM API
  • Simple Routines
    • smaceps
    • dmaceps
  • Derivative Routines
    • central_derivative_at
    • forward_derivative_at
    • backward_derivative_at
  • Vector Routines
    • Vector Arithmetic: add_v, minus_v
    • Norms: l1_norm, l2_norm, linf_norm
    • vector_distance
    • Distances: l1_distance, l2_distance, linf_distance
    • sum_v
    • scale_v
    • free_vector
    • copy_vector
    • format_vector_into
  • Matrix Routines
    • lu_decomp
    • bsubst
    • fsubst
    • solve_matrix
    • m_dot_v
    • put_identity_diagonal
    • copy_matrix
    • free_matrix
    • format_matrix_into
  • Root Finding Methods
    • find_ivt_range
    • bisect_find_root
    • bisect_find_root_with_error_assumption
  • Linear Routines
    • least_squares_lin_reg
  • Appendix / Miscellaneous
    • Data Types
      • Line
      • The Array_<type> and Matrix_<type>
    • Macros
      • c_max and c_min
      • InitArray
      • InitArrayWithSize
      • InitMatrixWithSize

Design

The LIZFCM static library (at [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 5, the testing routines are provided in test and utilize the utest microlibrary. They compile to a binary in ./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".

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

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".

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

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.

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

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.

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

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.

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

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

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

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

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

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

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

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.

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

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

double sum_v(Array_double *v) {
  double sum = 0;
  for (size_t i = 0; i < v->size; i++)
    sum += v->data[i];
  return sum;
}

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

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

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

void free_vector(Array_double *v) {
  free(v->data);
  free(v);
}

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

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

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

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

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

  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(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");
        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;
  }

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

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

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

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

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.

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

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

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

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

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

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

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

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.

void free_matrix(Matrix_double *m) {
  for (size_t y = 0; y < m->rows; ++y)
    free_vector(m->data[y]);
  free(m);
}

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

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

Root Finding Methods

find_ivt_range

• Author: Elizabeth Hunt
• Name: find_ivt_range
• Location: src/roots.c
• Input: a pointer to a oneary function taking a double and producing a double, the beginning point in $R$ to search for a range, a delta step that is taken, and a max_steps number of maximum iterations to perform.
• Output: a pair of double's representing a closed closed interval [beginning, end]

double *find_ivt_range(double (*f)(double), double start_x, double delta,
                       size_t max_steps) {
  double *range = malloc(sizeof(double) * 2);

  double a = start_x;

  while (f(a) * f(start_x) >= 0 && max_steps-- > 0)
    a += delta;

  if (max_steps == 0 && f(a) * f(start_x) > 0)
    return NULL;

  range[0] = start_x;
  range[1] = a + delta;
  return range;
}

bisect_find_root

• Author: Elizabeth Hunt
• Name(s): bisect_find_root
• Input: a one-ary function taking a double and producing a double, a closed interval represented by a and b: [a, b], a tolerance at which we return the estimated root, and a max_iterations to break us out of a loop if we can never reach the tolerance
• Output: a double representing the estimated root
• Description: recursively uses binary search to split the interval until we reach tolerance. We also assume the function f is continuous on [a, b].

double bisect_find_root(double (*f)(double), double a, double b,
                        double tolerance, size_t max_iterations) {
  assert(a <= b);
  // guarantee there's a root somewhere between a and b by IVT
  assert(f(a) * f(b) < 0);

  double c = (1.0 / 2) * (a + b);
  if (b - a < tolerance || max_iterations == 0)
    return c;
  if (f(a) * f(c) < 0)
    return bisect_find_root(f, a, c, tolerance, max_iterations - 1);
  return bisect_find_root(f, c, b, tolerance, max_iterations - 1);
}

bisect_find_root_with_error_assumption

• Author: Elizabeth Hunt
• Name: bisect_find_root_with_error_assumption
• Input: a one-ary function taking a double and producing a double, a closed interval represented by a and b: [a, b], and a tolerance at which we return the estimated root
• Output: a double representing the estimated root
• Description: using the bisection method we know that $e_k \le (\frac{1}{2})^k (b_0 - a_0)$. So we can calculate $k$ at the worst possible case (that the error is exactly the tolerance) to be $\frac{log(tolerance) - log(b_0 - a_0)}{log(\frac{1}{2})}$. We pass this value into the max_iterations of bisect_find_root as above.

double bisect_find_root_with_error_assumption(double (*f)(double), double a,
                                              double b, double tolerance) {
  assert(a <= b);

  uint64_t max_iterations =
      (uint64_t)ceil((log(tolerance) - log(b - a)) / log(1 / 2.0));
  return bisect_find_root(f, a, b, tolerance, max_iterations);
}

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

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

Appendix / Miscellaneous

Data Types

Line

• Author: Elizabeth Hunt
• Location: inc/types.h

typedef struct Line {
  double m;
  double a;
} Line;

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:

  typedef struct {
    int* data;
    size_t size;
  } Array_int

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.

  typedef struct {
    Array_int **data;
    size_t cols;
    size_t rows;
  } Matrix_int

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

#define c_max(x, y) (((x) >= (y)) ? (x) : (y))
#define c_min(x, y) (((x) <= (y)) ? (x) : (y))

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

#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;                                                                       \
  })

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

#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;                                                                       \
  })

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

#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;                                                                    \
  })