diff --git a/extension b/extension
deleted file mode 160000
index 2dc6ce7..0000000
--- a/extension
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit 2dc6ce7eabce9d86c0f5f44606d0cb18fc64d983
diff --git a/extension/src/Makefile b/extension/src/Makefile
new file mode 100644
index 0000000..3d728fc
--- /dev/null
+++ b/extension/src/Makefile
@@ -0,0 +1,18 @@
+CC = gcc
+CFLAGS = -std=c17 -g\
+	-D_POSIX_SOURCE -D_DEFAULT_SOURCE\
+	-Wall -Werror -pedantic\
+
+LDLIBS = -lm
+
+.SUFFIXES: .c .o
+
+.PHONY: all clean
+
+all: test/test_math test/test_ann
+
+test/test_math: test/test_math.o mymath.o
+test/test_ann: test/test_ann.o ann.o mymath.o
+
+clean:
+	$(RM) *.o test/*.o
diff --git a/extension/src/ann.c b/extension/src/ann.c
new file mode 100644
index 0000000..47ec042
--- /dev/null
+++ b/extension/src/ann.c
@@ -0,0 +1,81 @@
+#include <stdlib.h>
+#include "ann.h"
+
+// Helper function that simulate 'bias' term by setting a neuron that always outputs 1.0
+static void network_fillbias(Network network, unsigned int layerIndex) {
+
+    for (int j = 0; j < network->trainingWidth; j++) {
+        *(network->outputs[layerIndex]->data + network->layers[layerIndex].neuronCount * network->trainingWidth + j) = 1.0;
+    }
+}
+
+// Defines a sequence of random numbers sampled from a standard normal distribution
+double normalseq(unsigned int i) {
+    return stdnormal();
+}
+
+// Creates a new network instance given an array of layers and the size of the array
+Network network_create(Layer *layers, unsigned int noLayers, unsigned int trainingWidth) {
+
+    // Initialize network struct and its fields
+    Network network = malloc(sizeof(struct Network));
+    if (network == NULL) {
+        fprintf(stderr, "ERROR: Couldn't allocate sufficient memory for Network struct!\n");
+        abort();
+    }
+
+    network->layers = layers;
+    network->noLayers = noLayers;
+    network->trainingWidth = trainingWidth;
+
+    network->outputs = malloc(sizeof(struct Matrix) * noLayers);
+    if (network->outputs == NULL) {
+        fprintf(stderr, "ERROR: Couldn't allocate sufficient memory for Network output matrix array!\n");
+        abort();
+    }
+
+    network->weights = malloc(sizeof(struct Matrix) * (noLayers - 1));
+    if (network->weights == NULL) {
+        fprintf(stderr, "ERROR: Couldn't allocate sufficient memory for Network weight matrix array!\n");
+        abort();
+    }
+
+    // Initialize output and weight matrix array data
+    for (unsigned int i = 0; i < noLayers; i++) {
+        network->outputs[i] = matrix_create(layers[i].neuronCount + 1, trainingWidth);
+        
+        network_fillbias(network, i);
+
+        if (i < noLayers - 1) {
+            network->weights[i] = matrix_seqcreate(layers[i + 1].neuronCount + 1, layers[i].neuronCount + 1, &normalseq);
+        }
+    }
+
+    return network;
+}
+
+
+// Returns the output of the ANN given the specified input
+Matrix network_pass(Network network, Matrix features) {
+
+    Matrix currentOut = features;
+    for (int i = 0; i < network->noLayers; i++) {
+        
+        // Update current layer's outputs
+        matrix_free(network->outputs[i]);
+        network->outputs[i] = currentOut;
+        network_fillbias(network, i);
+
+        // Calculate outputs of next layer if not in last layer
+        if (i < network->noLayers - 1) {
+            currentOut = matrix_apply(matrix_multiply(network->weights[i], currentOut), network->layers[i].activation);
+        }
+    }
+
+    return currentOut;
+}
+
+
+Matrix log_loss(Matrix result, Matrix expected) {
+    
+}
diff --git a/extension/src/ann.h b/extension/src/ann.h
new file mode 100644
index 0000000..189ae3e
--- /dev/null
+++ b/extension/src/ann.h
@@ -0,0 +1,41 @@
+/* @file ann.h
+   @brief Data structures and functions for feedforward neural networks.
+
+   @author Themis Demetriades */
+
+#ifndef __ANN__
+#define __ANN__
+
+#include "mymath.h"
+
+// Defines function pointer type for activation functions
+typedef double (*activationFunc)(double);
+
+// Specifies structure of a given layer within a network
+typedef struct {
+   unsigned int neuronCount;
+   activationFunc activation;
+} Layer;
+
+// Defines a specific instance of an active network
+struct Network {
+   Layer *layers; // Array of layers specifying network structure
+   unsigned int noLayers;
+   Matrix *weights; // Array of matrices specifying weights of *outputs* to each layer
+   Matrix *outputs;
+
+   unsigned int trainingWidth; // Number of training examples to train per gradient descent step
+};
+
+typedef struct Network *Network;
+
+// Creates a new network instance given an array of layers and the size of the array
+Network network_create(Layer *layers, unsigned int noLayers, unsigned int trainingWidth);
+
+// Returns the output of the ANN given the specified input
+Matrix network_pass(Network network, Matrix features);
+
+// Calculates the log loss of a result wrt to the expected result
+Matrix log_loss(Matrix result, Matrix expected);
+
+#endif
diff --git a/extension/src/mymath.c b/extension/src/mymath.c
new file mode 100644
index 0000000..caf814a
--- /dev/null
+++ b/extension/src/mymath.c
@@ -0,0 +1,166 @@
+#include <stddef.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "mymath.h"
+
+#define LRELU_SFACTOR 0.001
+
+// Maximum number of characters that would ever be required to represent a valid matrix entry
+#define MAX_ENTRY_LENGTH 100
+
+// Defines the format used to print matrix entries
+#define VALUE_PRINT_FORMAT "%.2lf"
+
+// Pretty print matrix data to the given output stream
+void matrix_print(Matrix m, FILE *stream) {
+
+    double (*data)[m->cols] = (double (*)[m->cols]) m->data;
+    char valstr[MAX_ENTRY_LENGTH];
+    unsigned int maxlen = 1;
+    
+    // Get maximum length of a matrix entry
+    for (int i = 0; i < m->rows; i++) {
+        for (int j = 0; j < m->cols; j++) {
+            sprintf(valstr, VALUE_PRINT_FORMAT, data[i][j]);
+            unsigned int curlen = strlen(valstr);
+            if (curlen > maxlen) maxlen = curlen;
+        }
+    }
+
+    // Print matrix entries
+    for (int i = 0; i < m->rows; i++) {
+        fputc('[', stream);
+        for (int j = 0; j < m->cols; j++) {
+            sprintf(valstr, VALUE_PRINT_FORMAT, data[i][j]);
+            fprintf(stream, "%s", valstr);
+            for (int k = 0; k < maxlen - strlen(valstr); k++) fputc(' ', stream);
+            if (j < m->cols - 1) fputc(' ', stream);
+        }
+        fprintf(stream, "]\n");
+    }
+}
+
+// Return matrix with specified dimensions and random values as its entries
+Matrix matrix_create(unsigned int rows, unsigned int cols) {
+
+    // Allocate memory for matrix struct, aborting if couldn't allocate memory
+    Matrix m = malloc(sizeof(struct Matrix));
+    if (m == NULL) {
+        fprintf(stderr, "ERROR: Couldn't allocate sufficient memory for matrix struct!\n");
+        abort();
+    }
+
+    // Allocate memory for matrix data (values), aborting if couldn't allocate memory
+    m->data = malloc(rows * cols * sizeof(double));
+    if (m->data == NULL) {
+        fprintf(stderr, "ERROR: Couldn't allocate sufficient memory for matrix data!\n");
+        abort();
+    }
+
+    // Initialize matrix struct data
+    m->rows = rows;
+    m->cols = cols;
+    
+    return m;
+}
+
+// Free data associated with the specified matrix
+void matrix_free(Matrix m) {
+
+    free(m->data);
+    free(m);
+}
+
+// Return matrix with specified dimensions and entries whose values are the first values of the specified matrix sequence 
+Matrix matrix_seqcreate(unsigned int rows, unsigned int cols, matrix_sequence seq) {
+
+    // Create matrix with specified dimensions containing entries with random values
+    Matrix m = matrix_create(rows, cols);
+    
+    double (*data)[cols] = (double (*)[cols]) m->data;
+    // Fill matrix with correct entries
+    for (int i = 0; i < rows; i++) {
+        for (int j = 0; j < cols; j++) {
+           data[i][j] = seq(i * cols + j);
+        }
+    }
+    
+    return m;
+}
+
+// Multiply left matrix (l) with right matrix (r), returning a new instance representing the result
+Matrix matrix_multiply(Matrix l, Matrix r) {
+
+    // Ensure that dimensions of matrices are compatible for multiplication
+    if (l->cols != r->rows) {
+        fprintf(stderr, "ERROR: Attempting to multiply matrices with incompatible dimensions!\n");
+        abort();
+    }
+    
+    // Create matrix instance to store result of product
+    Matrix prod = matrix_create(l->rows, r->cols);
+
+    double (*prod_data)[prod->cols] = (double (*)[prod->cols]) prod->data;
+    double (*l_data)[l->cols] = (double (*)[l->cols]) l->data;
+    double (*r_data)[r->cols] = (double (*)[r->cols]) r->data;
+    // Perform matrix multiplication, storing result in prod
+    for (int i = 0; i < l->rows; i++) {
+        for (int j = 0; j < r->cols; j++) {
+            
+            prod_data[i][j] = 0.0;
+            for (int k = 0; k < l->cols; k++) {
+                prod_data[i][j] += l_data[i][k] * r_data[k][j];
+            }
+        }
+    }
+
+    // Return result
+    return prod; 
+}
+
+// Return result of applying given function to each entry in the matrix independently
+Matrix matrix_apply(Matrix m, matrix_entrytrans f) {
+     
+    Matrix result = matrix_create(m->rows, m->cols);
+    double (*result_data)[result->cols] = (double (*)[result->cols])result->data;
+    double (*m_data)[m->cols] = (double (*)[m->cols])m->data;
+
+    for (int i = 0; i < m->rows; i++) {
+        for (int j = 0; j < m->cols; j++) {
+           result_data[i][j] = f(m_data[i][j]);
+        }
+    }
+
+    return result;
+}
+
+// Draws a value from the standard normal distribution
+double stdnormal() {
+
+    // Generate 2 random numbers in the interval [0, 1)
+    double s1 = ((double)rand() / RAND_MAX);
+    double s2 = ((double)rand() / RAND_MAX);
+
+    // Peform box-muller transform
+    return sqrt(-2 * log(s1)) * cos(2 * PI * s2);
+}
+
+// Defines the leaky rectified linear unit activation function
+double lrelu(double x) {
+    if (x > 0) {
+        return x;
+    } else {
+        return LRELU_SFACTOR * x;
+    }
+}
+
+// Defines the identity activation function f(x) = x 
+double identity(double x) {
+    return x;
+}
+
+// Defines a simple counting sequence
+double countseq(unsigned int i) {
+    return i + 1;
+}
diff --git a/extension/src/mymath.h b/extension/src/mymath.h
new file mode 100644
index 0000000..4487d3b
--- /dev/null
+++ b/extension/src/mymath.h
@@ -0,0 +1,58 @@
+/* @file math.h
+   @brief Math functions and types for feedforward ANNs.
+   
+   @author Themis Demetriades */
+
+#ifndef __MATH__
+#define __MATH__
+
+#include <stdio.h>
+
+#define PI 3.1415926535
+
+// Definition of matrix types
+struct Matrix {
+    unsigned int rows;
+    unsigned int cols;
+    double *data;
+};
+
+typedef struct Matrix *Matrix;
+
+// Pretty print matrix data to the given output stream
+void matrix_print(Matrix m, FILE *stream);
+
+// Defines a sequence that maps matrix entry indices to floating point values
+typedef double (*matrix_sequence)(unsigned int);
+
+// Defines a function that transforms matrix entries
+typedef double (*matrix_entrytrans)(double);
+
+// Return matrix with specified dimensions and random values as its entries
+Matrix matrix_create(unsigned int rows, unsigned int cols);
+
+// Free data associated with the specified matrix
+void matrix_free(Matrix m);
+
+// Return matrix with specified dimensions and entries whose values are the first values of the specified matrix sequence 
+Matrix matrix_seqcreate(unsigned int rows, unsigned int cols, matrix_sequence seq);
+
+// Multiply left matrix (l) with right matrix (r), returning a new instance representing the result
+Matrix matrix_multiply(Matrix l, Matrix r);
+
+// Return result of applying given function to each entry in the matrix independently
+Matrix matrix_apply(Matrix m, matrix_entrytrans f);
+
+// Draws a value from the standard normal distribution
+double stdnormal();
+
+// Defines the identity activation function f(x) = x 
+double identity(double x);
+
+// Defines the rectified linear unit activation function
+double lrelu(double x);
+
+// Defines a simple counting sequence
+double countseq(unsigned int i);
+
+#endif
diff --git a/extension/src/test/test_ann b/extension/src/test/test_ann
new file mode 100755
index 0000000..80b301e
Binary files /dev/null and b/extension/src/test/test_ann differ
diff --git a/extension/src/test/test_ann.c b/extension/src/test/test_ann.c
new file mode 100644
index 0000000..9cb92f8
--- /dev/null
+++ b/extension/src/test/test_ann.c
@@ -0,0 +1,28 @@
+#include "../ann.h"
+
+#define TEST1_LAYERS 4
+
+int main(int argc, char **argv) {
+
+    Layer layers[] = {{1, &lrelu}, {3, &lrelu}, {3, &lrelu}, {2, &identity}};
+    Network network = network_create(layers, TEST1_LAYERS, 1);
+
+    for (int i = 0; i < TEST1_LAYERS; i++) {
+        matrix_print(network->outputs[i], stdout);
+        printf("\n");
+        if (i < TEST1_LAYERS - 1) {
+            matrix_print(network->weights[i], stdout);
+            printf("\n ------------- \n");
+        }
+    }
+    printf("===========================\n");
+
+    Matrix features = matrix_create(2, 1);
+    *(features->data + 0) = 1;
+
+    Matrix out = network_pass(network, features);
+
+    matrix_print(out, stdout);
+
+    return 0;
+}
\ No newline at end of file
diff --git a/extension/src/test/test_math b/extension/src/test/test_math
new file mode 100755
index 0000000..77f595a
Binary files /dev/null and b/extension/src/test/test_math differ
diff --git a/extension/src/test/test_math.c b/extension/src/test/test_math.c
new file mode 100644
index 0000000..9a6ec7e
--- /dev/null
+++ b/extension/src/test/test_math.c
@@ -0,0 +1,20 @@
+#include "../mymath.h"
+
+int main(int argc, char **argv) {
+
+    Matrix l = matrix_seqcreate(4, 3, &countseq);
+    Matrix r = matrix_seqcreate(3, 2, &countseq); 
+    
+    matrix_print(l, stdout);
+    printf("\n");
+    matrix_print(r, stdout);
+    printf("\n");
+
+    Matrix p = matrix_multiply(l, r);
+    matrix_print(p, stdout);
+
+    matrix_free(l);
+    matrix_free(r);
+    matrix_free(p);
+    return 0;
+}
\ No newline at end of file