diff --git a/graph.c b/graph.c
deleted file mode 100644
index 6110d04..0000000
--- a/graph.c
+++ /dev/null
@@ -1,481 +0,0 @@
-#include "graph.h"
-#include <stdint.h>
-#include <stdio.h>
-#include <string.h>
-#include <time.h>
-#include <stdlib.h>
-#include "matrix.h"
-
-void random_adjacency(const uint64_t vertex_count, uint64_t matrix[vertex_count][vertex_count]) {
-    srand(time(NULL));
-
-    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-            if (column_index == row_index) {
-                matrix[row_index][column_index] = 0;
-            } else {
-                matrix[row_index][column_index] = rand() % 2;
-            }
-        }
-    }
-}
-
-void calculate_distance_matrix(const uint64_t vertex_count, const uint64_t adjacency_matrix[vertex_count][vertex_count], uint64_t distance_matrix[vertex_count][vertex_count]) {
-    uint64_t power_matrix[vertex_count][vertex_count];
-    uint64_t temp_power_matrix[vertex_count][vertex_count];
-    memcpy(power_matrix, adjacency_matrix, vertex_count * vertex_count * sizeof(uint64_t));
-
-    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-            if (row_index == column_index) {
-                distance_matrix[row_index][column_index] = 0;
-            } else if (adjacency_matrix[row_index][column_index] == 1) {
-                distance_matrix[row_index][column_index] = 1;
-            } else {
-                distance_matrix[row_index][column_index] = UINT64_MAX;
-            }
-        }
-    }
-
-    for(uint64_t k = 2; k <= vertex_count; k++) {
-        memcpy(temp_power_matrix, power_matrix, vertex_count * vertex_count * sizeof(uint64_t));
-        gemm_basic(vertex_count, vertex_count, adjacency_matrix, vertex_count, vertex_count, temp_power_matrix, power_matrix);
-
-        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-            for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-                if (power_matrix[row_index][column_index] != 0 && distance_matrix[row_index][column_index] == UINT64_MAX) {
-                    distance_matrix[row_index][column_index] = k;
-                }
-            }
-        }
-    }
-}
-
-int get_eccentricities(const uint64_t vertex_count, const uint64_t distance_matrix[vertex_count][vertex_count], uint64_t eccentricities[vertex_count]) {
-    uint64_t eccentricity;
-
-    // set all eccentricities to infinity in case this is a disconnected graph
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        eccentricities[index] = UINT64_MAX;
-    }
-
-    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-        eccentricity = 0;
-
-        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-            if (distance_matrix[row_index][column_index] > eccentricity) {
-                eccentricity = distance_matrix[row_index][column_index];
-            }
-        }
-        // in case of a disconnected graph
-        if (eccentricity == 0) {
-            return 1;
-        }
-        eccentricities[row_index] = eccentricity;
-    }
-
-    return 0;
-}
-
-uint64_t get_radius(const uint64_t vertex_count, const uint64_t eccentricities[vertex_count]) {
-    uint64_t radius = UINT64_MAX;
-
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        if (eccentricities[index] < radius) {
-            radius = eccentricities[index];
-        }
-    }
-
-    return radius;
-}
-
-uint64_t get_diameter(const uint64_t vertex_count, const uint64_t eccentricities[vertex_count]) {
-    uint64_t diamter = 0;
-
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        if (eccentricities[index] > diamter) {
-            diamter = eccentricities[index];
-        }
-    }
-
-    return diamter;
-}
-
-void get_centre(const uint64_t vertex_count, const uint64_t eccentricities[vertex_count], const uint64_t radius, uint64_t centre[vertex_count]) {
-    memset(centre, 0, vertex_count * sizeof(uint64_t));
-
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        if (eccentricities[index] == radius) {
-            centre[index] = 1;
-        }
-    }
-}
-
-void calculate_path_matrix(const uint64_t vertex_count, const uint64_t adjacency_matrix[vertex_count][vertex_count], uint64_t path_matrix[vertex_count][vertex_count]) {
-    uint64_t power_matrix[vertex_count][vertex_count];
-    uint64_t temp_power_matrix[vertex_count][vertex_count];
-    memcpy(power_matrix, adjacency_matrix, vertex_count * vertex_count * sizeof(uint64_t));
-
-    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-            if (row_index == column_index || adjacency_matrix[row_index][column_index] == 1) {
-                path_matrix[row_index][column_index] = 1;
-            } else {
-                path_matrix[row_index][column_index] = 0;
-            }
-        }
-    }
-
-    for(uint64_t k = 2; k <= vertex_count; k++) {
-        memcpy(temp_power_matrix, power_matrix, vertex_count * vertex_count * sizeof(uint64_t));
-        gemm_basic(vertex_count, vertex_count, adjacency_matrix, vertex_count, vertex_count, temp_power_matrix, power_matrix);
-
-        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-            for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-                if (power_matrix[row_index][column_index] != 0) {
-                    path_matrix[row_index][column_index] = 1;
-                }
-            }
-        }
-    }
-}
-
-void find_components_basic(const uint64_t vertex_count, const uint64_t path_matrix[vertex_count][vertex_count], uint64_t components[vertex_count][vertex_count]) {
-    uint64_t component[vertex_count];
-    int contains_component;
-
-    memset(components, 0, vertex_count * vertex_count * sizeof(uint64_t));
-
-    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-        memset(component, 0, vertex_count * sizeof(uint64_t));
-        contains_component = 0;
-
-        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-            if (path_matrix[row_index][column_index] == 1) {
-                component[column_index] = column_index + 1;
-            }
-        }
-
-        for (uint64_t index = 0; index < vertex_count; index++) {
-            if (memcmp(components[index], component, vertex_count * sizeof(uint64_t)) == 0) {
-                contains_component = 1;
-            }
-        }
-
-        if (!contains_component) {
-            for (uint64_t index = 0; index < vertex_count; index++) {
-                components[row_index][index] = component[index];
-            }
-        }
-    }
-}
-
-void dfs(const uint64_t vertex_count, const uint64_t adjacency_matrix[vertex_count][vertex_count], const uint64_t vertex, uint64_t visited[vertex_count]) {
-    visited[vertex] = 1;
-
-    for (uint64_t neighbor_vertex = 0; neighbor_vertex < vertex_count; neighbor_vertex++) {
-        if (adjacency_matrix[vertex][neighbor_vertex] != 1) {
-            continue;
-        }
-        if (visited[neighbor_vertex]) {
-            continue;
-        }
-        dfs(vertex_count, adjacency_matrix, neighbor_vertex, visited);
-    }
-}
-
-void find_components_dfs(const uint64_t vertex_count, const uint64_t adjacency_matrix[vertex_count][vertex_count], uint64_t components[vertex_count][vertex_count]) {
-    uint64_t component[vertex_count];
-    int contains_component = 0;
-
-    memset(components, 0, vertex_count * vertex_count * sizeof(uint64_t));
-
-    for (uint64_t vertex = 0; vertex < vertex_count; vertex++) {
-        memset(component, 0, vertex_count * sizeof(uint64_t));
-        contains_component = 0;
-
-        dfs(vertex_count, adjacency_matrix, vertex, component);
-
-        for (uint64_t index = 0; index < vertex_count; index++) {
-            if (memcmp(components[index], component, vertex_count * sizeof(uint64_t)) == 0) {
-                contains_component = 1;
-            }
-        }
-
-        if (!contains_component) {
-            for (uint64_t index = 0; index < vertex_count; index++) {
-                components[vertex][index] = component[index];
-            }
-        }
-    }
-}
-
-uint64_t amount_of_components(const uint64_t vertex_count, const uint64_t components[vertex_count][vertex_count]) {
-    uint64_t amount_of_components = 0;
-
-    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-            if (components[row_index][column_index] != 0) {
-                amount_of_components++;
-                break;
-            }
-        }
-    }
-
-    return amount_of_components;
-}
-
-int contains_bridge(const uint64_t vertex_count, const uint64_t bridges[vertex_count][2], const uint64_t bridge[2]) {
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        if (memcmp(bridges[index], bridge, 2 * sizeof(uint64_t)) == 0) {
-            return 1;
-        }
-    }
-
-    return 0;
-}
-
-void find_bridges_basic(const uint64_t vertex_count, const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                        const uint64_t components[vertex_count][vertex_count], uint64_t bridges[vertex_count][2]) {
-    uint64_t path_matrix[vertex_count][vertex_count];
-    uint64_t temp_adjacency_matrix[vertex_count][vertex_count];
-    uint64_t temp_components[vertex_count][vertex_count];
-    uint64_t bridge[2];
-
-    memset(bridges, 0, vertex_count * 2 * sizeof(uint64_t));
-
-    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-            if (row_index == column_index) {
-                continue;
-            }
-
-            if (column_index < row_index) {
-                bridge[0] = column_index + 1;
-                bridge[1] = row_index + 1;
-            } else {
-                bridge[0] = row_index + 1;
-                bridge[1] = column_index + 1;
-            }
-
-            memcpy(temp_adjacency_matrix, adjacency_matrix, vertex_count * vertex_count * sizeof(uint64_t));
-
-            temp_adjacency_matrix[row_index][column_index] = 0;
-            temp_adjacency_matrix[column_index][row_index] = 0;
-
-            calculate_path_matrix(vertex_count, temp_adjacency_matrix, path_matrix);
-            find_components_basic(vertex_count, path_matrix, temp_components);
-
-            if (amount_of_components(vertex_count, temp_components) > amount_of_components(vertex_count, components) && !contains_bridge(vertex_count, bridges, bridge)) {
-                bridges[row_index][0] = bridge[0];
-                bridges[row_index][1] = bridge[1];
-            }
-        }
-    }
-    print_matrix(vertex_count, 2, bridges);
-}
-
-void find_bridges_dfs_v1(const uint64_t vertex_count, const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                        const uint64_t components[vertex_count][vertex_count], uint64_t bridges[vertex_count][2]) {
-    uint64_t temp_adjacency_matrix[vertex_count][vertex_count];
-    uint64_t temp_components[vertex_count][vertex_count];
-    uint64_t bridge[2];
-
-    memset(bridges, 0, vertex_count * 2 * sizeof(uint64_t));
-
-    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-            if (row_index == column_index) {
-                continue;
-            }
-
-            if (column_index < row_index) {
-                bridge[0] = column_index + 1;
-                bridge[1] = row_index + 1;
-            } else {
-                bridge[0] = row_index + 1;
-                bridge[1] = column_index + 1;
-            }
-
-            memcpy(temp_adjacency_matrix, adjacency_matrix, vertex_count * vertex_count * sizeof(uint64_t));
-
-            temp_adjacency_matrix[row_index][column_index] = 0;
-            temp_adjacency_matrix[column_index][row_index] = 0;
-
-            find_components_dfs(vertex_count, temp_adjacency_matrix, temp_components);
-
-            if (amount_of_components(vertex_count, temp_components) > amount_of_components(vertex_count, components) && !contains_bridge(vertex_count, bridges, bridge)) {
-                bridges[row_index][0] = bridge[0];
-                bridges[row_index][1] = bridge[1];
-            }
-        }
-    }
-}
-
-void dfs_bridges(const uint64_t vertex_count, const uint64_t adjacency_matrix[vertex_count][vertex_count], const uint64_t vertex, const uint64_t parent_vertex, uint64_t visited[vertex_count],
-                 uint64_t current_time, uint64_t discovery_time[vertex_count], uint64_t lowest_time[vertex_count], uint64_t bridges[vertex_count][2]) {
-    current_time++;
-    visited[vertex] = 1;
-    discovery_time[vertex] = current_time;
-    lowest_time[vertex] = current_time;
-
-    for (uint64_t neighbor_vertex = 0; neighbor_vertex < vertex_count; neighbor_vertex++) {
-        if (adjacency_matrix[vertex][neighbor_vertex] != 1) {
-            continue;
-        }
-        if (parent_vertex != neighbor_vertex && visited[neighbor_vertex]) {
-            if (lowest_time[vertex] > discovery_time[neighbor_vertex]) {
-                lowest_time[vertex] = discovery_time[neighbor_vertex];
-            }
-            continue;
-        }
-        if (visited[neighbor_vertex]) {
-            continue;
-        }
-
-        dfs_bridges(vertex_count, adjacency_matrix, neighbor_vertex, vertex, visited, current_time, discovery_time, lowest_time, bridges);
-
-        if (lowest_time[vertex] > lowest_time[neighbor_vertex]) {
-            lowest_time[vertex] = lowest_time[neighbor_vertex];
-        }
-        if (discovery_time[vertex] < lowest_time[neighbor_vertex]) {
-            bridges[neighbor_vertex][0] = vertex + 1;
-            bridges[neighbor_vertex][1] = neighbor_vertex + 1;
-        }
-    }
-}
-
-void find_bridges_dfs_v2(const uint64_t vertex_count, const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                         const uint64_t components[vertex_count][vertex_count], uint64_t bridges[vertex_count][2]) {
-    uint64_t visited[vertex_count];
-    uint64_t current_time = 0;
-    uint64_t discovery_time[vertex_count];
-    uint64_t lowest_time[vertex_count];
-
-    memset(bridges, 0, vertex_count * 2 * sizeof(uint64_t));
-    memset(visited, 0, vertex_count * sizeof(uint64_t));
-    memset(discovery_time, 0, vertex_count * sizeof(uint64_t));
-    memset(lowest_time, 0, vertex_count * sizeof(uint64_t));
-
-    for (uint64_t vertex = 0; vertex < vertex_count; vertex++) {
-        if (!visited[vertex]) {
-            dfs_bridges(vertex_count, adjacency_matrix, vertex, UINT64_MAX, visited, current_time, discovery_time, lowest_time, bridges);
-        }
-    }
-}
-
-void find_articulations_basic(const uint64_t vertex_count, const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                              const uint64_t components[vertex_count][vertex_count], uint64_t articulations[vertex_count]) {
-    uint64_t path_matrix[vertex_count][vertex_count];
-    uint64_t temp_adjacency_matrix[vertex_count][vertex_count];
-    uint64_t temp_components[vertex_count][vertex_count];
-
-    memset(articulations, 0, vertex_count * sizeof(uint64_t));
-
-    for (uint64_t i = 0; i < vertex_count; i++) {
-        memcpy(temp_adjacency_matrix, adjacency_matrix, vertex_count * vertex_count * sizeof(uint64_t));
-
-        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-            for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-                temp_adjacency_matrix[row_index][i] = 0;
-                temp_adjacency_matrix[i][column_index] = 0;
-            }
-        }
-
-        calculate_path_matrix(vertex_count, temp_adjacency_matrix, path_matrix);
-        find_components_basic(vertex_count, path_matrix, temp_components);
-
-        // the + 1 is needed because I am not removing the vertex, I am just removing all of its edges
-        // removing all of its edges, means it itself becomes a component, which needs to be accounted for
-        if (amount_of_components(vertex_count, temp_components) > amount_of_components(vertex_count, components) + 1) {
-            articulations[i] = i + 1;
-        }
-    }
-}
-
-void find_articulations_dfs_v1(const uint64_t vertex_count, const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                              const uint64_t components[vertex_count][vertex_count], uint64_t articulations[vertex_count]) {
-    uint64_t temp_adjacency_matrix[vertex_count][vertex_count];
-    uint64_t temp_components[vertex_count][vertex_count];
-
-    memset(articulations, 0, vertex_count * sizeof(uint64_t));
-
-    for (uint64_t i = 0; i < vertex_count; i++) {
-        memcpy(temp_adjacency_matrix, adjacency_matrix, vertex_count * vertex_count * sizeof(uint64_t));
-
-        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-            for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-                temp_adjacency_matrix[row_index][i] = 0;
-                temp_adjacency_matrix[i][column_index] = 0;
-            }
-        }
-
-        find_components_dfs(vertex_count, temp_adjacency_matrix, temp_components);
-
-        // the + 1 is needed because I am not removing the vertex, I am just removing all of its edges
-        // removing all of its edges, means it itself becomes a component, which needs to be accounted for
-        if (amount_of_components(vertex_count, temp_components) > amount_of_components(vertex_count, components) + 1) {
-            articulations[i] = i + 1;
-        }
-    }
-}
-
-void dfs_articulations(const uint64_t vertex_count, const uint64_t adjacency_matrix[vertex_count][vertex_count], const uint64_t vertex, const uint64_t parent_vertex,
-                       uint64_t visited[vertex_count], uint64_t current_time, uint64_t discovery_time[vertex_count], uint64_t lowest_time[vertex_count],
-                       uint64_t articulations[vertex_count]) {
-    current_time++;
-    visited[vertex] = 1;
-    discovery_time[vertex] = current_time;
-    lowest_time[vertex] = current_time;
-    uint64_t child_count = 0;
-    uint64_t is_articulation = 0;
-
-    for (uint64_t neighbor_vertex = 0; neighbor_vertex < vertex_count; neighbor_vertex++) {
-        if (adjacency_matrix[vertex][neighbor_vertex] != 1) {
-            continue;
-        }
-        if (visited[neighbor_vertex]) {
-            if (lowest_time[vertex] > discovery_time[neighbor_vertex]) {
-                lowest_time[vertex] = discovery_time[neighbor_vertex];
-            }
-            continue;
-        }
-
-        child_count++;
-
-        dfs_articulations(vertex_count, adjacency_matrix, neighbor_vertex, vertex, visited, current_time, discovery_time, lowest_time, articulations);
-
-        if (lowest_time[vertex] > lowest_time[neighbor_vertex]) {
-            lowest_time[vertex] = lowest_time[neighbor_vertex];
-        }
-
-        if (parent_vertex != UINT64_MAX && discovery_time[vertex] <= lowest_time[neighbor_vertex]) {
-            is_articulation = 1;
-        }
-    }
-    if (parent_vertex == UINT64_MAX && child_count > 1) {
-        is_articulation = 1;
-    }
-    if (is_articulation) {
-        articulations[vertex] = vertex + 1;
-    }
-}
-
-void find_articulations_dfs_v2(const uint64_t vertex_count, const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                         const uint64_t components[vertex_count][vertex_count], uint64_t articulations[vertex_count]) {
-    uint64_t visited[vertex_count];
-    uint64_t current_time = 0;
-    uint64_t discovery_time[vertex_count];
-    uint64_t lowest_time[vertex_count];
-
-    memset(articulations, 0, vertex_count * sizeof(uint64_t));
-    memset(visited, 0, vertex_count * sizeof(uint64_t));
-    memset(discovery_time, 0, vertex_count * sizeof(uint64_t));
-    memset(lowest_time, 0, vertex_count * sizeof(uint64_t));
-
-    for (uint64_t vertex = 0; vertex < vertex_count; vertex++) {
-        if (!visited[vertex]) {
-            dfs_articulations(vertex_count, adjacency_matrix, vertex, UINT64_MAX, visited, current_time, discovery_time, lowest_time, articulations);
-        }
-    }
-}
diff --git a/graph.cpp b/graph.cpp
new file mode 100644
index 0000000..e2e7629
--- /dev/null
+++ b/graph.cpp
@@ -0,0 +1,498 @@
+#include "matrix.hpp"
+#include <stdint.h>
+#include <string.h>
+#include <time.h>
+#include <stdlib.h>
+#include <cstdint>
+#include <vector>
+#include <algorithm>
+
+std::vector<std::vector<uint64_t>> random_adjacency(const uint64_t vertex_count) {
+    std::vector<std::vector<uint64_t>> adjacency_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
+    srand(time(NULL));
+
+    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
+        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
+            if (column_index != row_index) {
+                adjacency_matrix[row_index][column_index] = rand() % 2;
+            }
+        }
+    }
+
+    return adjacency_matrix;
+}
+
+std::vector<std::vector<uint64_t>> calculate_distance_matrix(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
+    uint64_t vertex_count = adjacency_matrix.size();
+    std::vector<std::vector<uint64_t>> distance_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
+    std::vector<std::vector<uint64_t>> power_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
+    std::copy(adjacency_matrix.begin(), adjacency_matrix.end(), power_matrix.begin());
+
+    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
+        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
+            if (row_index == column_index) {
+                continue;
+            } else if (adjacency_matrix[row_index][column_index] == 1) {
+                distance_matrix[row_index][column_index] = 1;
+            } else {
+                distance_matrix[row_index][column_index] = UINT64_MAX;
+            }
+        }
+    }
+
+    for(uint64_t k = 2; k <= vertex_count; k++) {
+        power_matrix = gemm_basic(adjacency_matrix, power_matrix);
+
+        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
+            for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
+                if (power_matrix[row_index][column_index] != 0 && distance_matrix[row_index][column_index] == UINT64_MAX) {
+                    distance_matrix[row_index][column_index] = k;
+                }
+            }
+        }
+    }
+
+    return distance_matrix;
+}
+
+std::vector<uint64_t> get_eccentricities(const std::vector<std::vector<uint64_t>>& distance_matrix) {
+    uint64_t vertex_count = distance_matrix.size();
+    std::vector<uint64_t> eccentricities(vertex_count, UINT64_MAX);
+    uint64_t eccentricity;
+
+    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
+        eccentricity = 0;
+
+        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
+            if (distance_matrix[row_index][column_index] > eccentricity) {
+                eccentricity = distance_matrix[row_index][column_index];
+            }
+        }
+
+        if (eccentricity == 0) {
+            break;
+        }
+        eccentricities[row_index] = eccentricity;
+    }
+
+    return eccentricities;
+}
+
+uint64_t get_radius(const std::vector<uint64_t>& eccentricities) {
+    uint64_t radius = UINT64_MAX;
+
+    for (uint64_t index = 0; index < eccentricities.size(); index++) {
+        if (eccentricities[index] < radius) {
+            radius = eccentricities[index];
+        }
+    }
+
+    return radius;
+}
+
+uint64_t get_diameter(const std::vector<uint64_t>& eccentricities) {
+    uint64_t diamter = 0;
+
+    for (uint64_t index = 0; index < eccentricities.size(); index++) {
+        if (eccentricities[index] > diamter) {
+            diamter = eccentricities[index];
+        }
+    }
+
+    return diamter;
+}
+
+std::vector<uint64_t> get_centre(const std::vector<uint64_t>& eccentricities) {
+    std::vector<uint64_t> centre;
+    uint64_t radius = get_radius(eccentricities);
+
+    for (uint64_t index = 0; index < eccentricities.size(); index++) {
+        if (eccentricities[index] == radius) {
+            centre.push_back(index + 1);
+        }
+    }
+
+    centre.shrink_to_fit();
+    return centre;
+}
+
+std::vector<std::vector<uint64_t>> calculate_path_matrix(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
+    uint64_t vertex_count = adjacency_matrix.size();
+    std::vector<std::vector<uint64_t>> path_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
+    std::vector<std::vector<uint64_t>> power_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
+    std::copy(adjacency_matrix.begin(), adjacency_matrix.end(), power_matrix.begin());
+
+    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
+        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
+            if (row_index == column_index || adjacency_matrix[row_index][column_index] == 1) {
+                path_matrix[row_index][column_index] = 1;
+            }
+        }
+    }
+
+    for(uint64_t k = 2; k <= vertex_count; k++) {
+        power_matrix = gemm_basic(adjacency_matrix, power_matrix);
+
+        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
+            for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
+                if (power_matrix[row_index][column_index] != 0) {
+                    path_matrix[row_index][column_index] = 1;
+                }
+            }
+        }
+    }
+
+    return path_matrix;
+}
+
+std::vector<std::vector<uint64_t>> find_components_basic(const std::vector<std::vector<uint64_t>>& path_matrix) {
+    uint64_t vertex_count = path_matrix.size();
+    std::vector<std::vector<uint64_t>> components;
+    std::vector<uint64_t> component;
+    components.reserve(0);
+    component.reserve(vertex_count);
+
+    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
+        component.clear();
+
+        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
+            if (path_matrix[row_index][column_index] == 1) {
+                component.push_back(column_index + 1);
+            }
+        }
+
+        component.shrink_to_fit();
+        std::sort(component.begin(), component.end());
+
+        auto iterator = find(components.begin(), components.end(), component);
+
+        if (iterator == components.end()) {
+            components.push_back(component);
+        }
+    }
+
+    components.shrink_to_fit();
+    return components;
+}
+
+void dfs(const std::vector<std::vector<uint64_t>>& adjacency_matrix, const uint64_t vertex, std::vector<bool>& visited) {
+    visited[vertex] = true;
+
+    for (uint64_t neighbor_vertex = 0; neighbor_vertex < adjacency_matrix.size(); neighbor_vertex++) {
+        if (adjacency_matrix[vertex][neighbor_vertex] != 1) {
+            continue;
+        }
+        if (visited[neighbor_vertex]) {
+            continue;
+        }
+        dfs(adjacency_matrix, neighbor_vertex, visited);
+    }
+}
+
+std::vector<std::vector<uint64_t>> find_components_dfs(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
+    uint64_t vertex_count = adjacency_matrix.size();
+    std::vector<std::vector<uint64_t>> components;
+    std::vector<uint64_t> component;
+    std::vector<bool> visited(vertex_count, false);
+    components.reserve(0);
+    component.reserve(vertex_count);
+
+    for (uint64_t vertex = 0; vertex < vertex_count; vertex++) {
+        component.clear();
+
+        dfs(adjacency_matrix, vertex, visited);
+
+        for (uint64_t index = 0; index < vertex_count; index++) {
+            if (visited[index]) {
+                component.push_back(index + 1);
+            }
+        }
+
+        component.shrink_to_fit();
+        std::sort(component.begin(), component.end());
+
+        auto iterator = find(components.begin(), components.end(), component);
+
+        if (iterator == components.end()) {
+            components.push_back(component);
+        }
+    }
+
+    components.shrink_to_fit();
+    return components;
+}
+
+std::vector<std::vector<uint64_t>> find_bridges_basic(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
+    uint64_t vertex_count = adjacency_matrix.size();
+    std::vector<std::vector<uint64_t>> path_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
+    std::vector<std::vector<uint64_t>> temp_adjacency_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
+    std::vector<std::vector<uint64_t>> components = find_components_basic(path_matrix);
+    std::vector<std::vector<uint64_t>> temp_components;
+    std::vector<std::vector<uint64_t>> bridges;
+    std::vector<uint64_t> bridge(2, 0);
+
+    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
+        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
+            if (row_index == column_index) {
+                continue;
+            }
+
+            std::copy(adjacency_matrix.begin(), adjacency_matrix.end(), temp_adjacency_matrix.begin());
+
+            temp_adjacency_matrix[row_index][column_index] = 0;
+            temp_adjacency_matrix[column_index][row_index] = 0;
+
+            path_matrix = calculate_path_matrix(temp_adjacency_matrix);
+            temp_components = find_components_basic(path_matrix);
+
+            if (temp_components.size() <= components.size()) {
+                continue;
+            }
+
+            if (column_index < row_index) {
+                bridge[0] = column_index + 1;
+                bridge[1] = row_index + 1;
+            } else {
+                bridge[0] = row_index + 1;
+                bridge[1] = column_index + 1;
+            }
+
+            auto iterator = find(bridges.begin(), bridges.end(), bridge);
+
+            if (iterator == bridges.end()) {
+                bridges.push_back(bridge);
+            }
+        }
+    }
+
+    bridges.shrink_to_fit();
+    return bridges;
+}
+
+std::vector<std::vector<uint64_t>> find_bridges_dfs_v1(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
+    uint64_t vertex_count = adjacency_matrix.size();
+    std::vector<std::vector<uint64_t>> temp_adjacency_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
+    std::vector<std::vector<uint64_t>> components = find_components_dfs(adjacency_matrix);
+    std::vector<std::vector<uint64_t>> temp_components;
+    std::vector<std::vector<uint64_t>> bridges;
+    std::vector<uint64_t> bridge(2, 0);
+
+    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
+        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
+            if (row_index == column_index) {
+                continue;
+            }
+
+            std::copy(adjacency_matrix.begin(), adjacency_matrix.end(), temp_adjacency_matrix.begin());
+
+            temp_adjacency_matrix[row_index][column_index] = 0;
+            temp_adjacency_matrix[column_index][row_index] = 0;
+
+            temp_components = find_components_dfs(temp_adjacency_matrix);
+
+            if (temp_components.size() <= components.size()) {
+                continue;
+            }
+
+            if (column_index < row_index) {
+                bridge[0] = column_index + 1;
+                bridge[1] = row_index + 1;
+            } else {
+                bridge[0] = row_index + 1;
+                bridge[1] = column_index + 1;
+            }
+
+            auto iterator = find(bridges.begin(), bridges.end(), bridge);
+
+            if (iterator == bridges.end()) {
+                bridges.push_back(bridge);
+            }
+        }
+    }
+
+    bridges.shrink_to_fit();
+    return bridges;
+}
+
+void dfs_bridges(const std::vector<std::vector<uint64_t>>& adjacency_matrix, const uint64_t vertex, const uint64_t parent_vertex, std::vector<bool>& visited,
+                 uint64_t current_time, std::vector<uint64_t>& discovery_time, std::vector<uint64_t>& lowest_time, std::vector<std::vector<uint64_t>>& bridges) {
+    current_time++;
+    visited[vertex] = true;
+    discovery_time[vertex] = current_time;
+    lowest_time[vertex] = current_time;
+    std::vector<uint64_t> bridge(2, 0);
+
+    for (uint64_t neighbor_vertex = 0; neighbor_vertex < adjacency_matrix.size(); neighbor_vertex++) {
+        if (adjacency_matrix[vertex][neighbor_vertex] != 1) {
+            continue;
+        }
+        if (parent_vertex != neighbor_vertex && visited[neighbor_vertex]) {
+            if (lowest_time[vertex] > discovery_time[neighbor_vertex]) {
+                lowest_time[vertex] = discovery_time[neighbor_vertex];
+            }
+            continue;
+        }
+        if (visited[neighbor_vertex]) {
+            continue;
+        }
+
+        dfs_bridges(adjacency_matrix, neighbor_vertex, vertex, visited, current_time, discovery_time, lowest_time, bridges);
+
+        if (lowest_time[vertex] > lowest_time[neighbor_vertex]) {
+            lowest_time[vertex] = lowest_time[neighbor_vertex];
+        }
+        if (discovery_time[vertex] >= lowest_time[neighbor_vertex]) {
+            continue;
+        }
+        bridge[0] = vertex + 1;
+        bridge[1] = neighbor_vertex + 1;
+
+        auto iterator = find(bridges.begin(), bridges.end(), bridge);
+
+        if (iterator == bridges.end()) {
+            bridges.push_back(bridge);
+        }
+    }
+}
+
+std::vector<std::vector<uint64_t>> find_bridges_dfs_v2(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
+    uint64_t vertex_count = adjacency_matrix.size();
+    std::vector<bool> visited(vertex_count, false);
+    uint64_t current_time = 0;
+    std::vector<uint64_t> discovery_time(vertex_count, 0);
+    std::vector<uint64_t> lowest_time(vertex_count, 0);
+    std::vector<std::vector<uint64_t>> bridges;
+
+    for (uint64_t vertex = 0; vertex < vertex_count; vertex++) {
+        if (!visited[vertex]) {
+            dfs_bridges(adjacency_matrix, vertex, UINT64_MAX, visited, current_time, discovery_time, lowest_time, bridges);
+        }
+    }
+
+    bridges.shrink_to_fit();
+    return bridges;
+}
+
+std::vector<uint64_t> find_articulations_basic(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
+    uint64_t vertex_count = adjacency_matrix.size();
+    std::vector<std::vector<uint64_t>> path_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
+    std::vector<std::vector<uint64_t>> temp_adjacency_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
+    std::vector<std::vector<uint64_t>> components = find_components_basic(path_matrix);
+    std::vector<std::vector<uint64_t>> temp_components;
+    std::vector<uint64_t> articulations;
+
+    for (uint64_t i = 0; i < vertex_count; i++) {
+        std::copy(adjacency_matrix.begin(), adjacency_matrix.end(), temp_adjacency_matrix.begin());
+
+        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
+            for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
+                temp_adjacency_matrix[row_index][i] = 0;
+                temp_adjacency_matrix[i][column_index] = 0;
+            }
+        }
+
+        path_matrix = calculate_path_matrix(temp_adjacency_matrix);
+        temp_components = find_components_basic(path_matrix);
+
+        // the + 1 is needed because I am not removing the vertex, I am just removing all of its edges
+        // removing all of its edges, means it itself becomes a component, which needs to be accounted for
+        if (temp_components.size() > components.size() + 1) {
+            articulations.push_back(i + 1);
+        }
+    }
+
+    articulations.shrink_to_fit();
+    std::sort(articulations.begin(), articulations.end());
+    return articulations;
+}
+
+std::vector<uint64_t> find_articulations_dfs_v1(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
+    uint64_t vertex_count = adjacency_matrix.size();
+    std::vector<std::vector<uint64_t>> path_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
+    std::vector<std::vector<uint64_t>> temp_adjacency_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
+    std::vector<std::vector<uint64_t>> components = find_components_dfs(adjacency_matrix);
+    std::vector<std::vector<uint64_t>> temp_components;
+    std::vector<uint64_t> articulations;
+
+    for (uint64_t i = 0; i < vertex_count; i++) {
+        std::copy(adjacency_matrix.begin(), adjacency_matrix.end(), temp_adjacency_matrix.begin());
+
+        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
+            for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
+                temp_adjacency_matrix[row_index][i] = 0;
+                temp_adjacency_matrix[i][column_index] = 0;
+            }
+        }
+
+        temp_components = find_components_dfs(temp_adjacency_matrix);
+
+        // the + 1 is needed because I am not removing the vertex, I am just removing all of its edges
+        // removing all of its edges, means it itself becomes a component, which needs to be accounted for
+        if (temp_components.size() > components.size() + 1) {
+            articulations.push_back(i + 1);
+        }
+    }
+
+    articulations.shrink_to_fit();
+    std::sort(articulations.begin(), articulations.end());
+    return articulations;
+}
+
+void dfs_articulations(const std::vector<std::vector<uint64_t>>& adjacency_matrix, const uint64_t vertex, const uint64_t parent_vertex, std::vector<bool>& visited,
+                 uint64_t current_time, std::vector<uint64_t>& discovery_time, std::vector<uint64_t>& lowest_time, std::vector<uint64_t>& articulations) {
+    current_time++;
+    visited[vertex] = true;
+    discovery_time[vertex] = current_time;
+    lowest_time[vertex] = current_time;
+    uint64_t child_count = 0;
+    bool is_articulation = false;
+
+    for (uint64_t neighbor_vertex = 0; neighbor_vertex < adjacency_matrix.size(); neighbor_vertex++) {
+        if (adjacency_matrix[vertex][neighbor_vertex] != 1) {
+            continue;
+        }
+        if (visited[neighbor_vertex]) {
+            if (lowest_time[vertex] > discovery_time[neighbor_vertex]) {
+                lowest_time[vertex] = discovery_time[neighbor_vertex];
+            }
+            continue;
+        }
+
+        child_count++;
+
+        dfs_articulations(adjacency_matrix, neighbor_vertex, vertex, visited, current_time, discovery_time, lowest_time, articulations);
+
+        if (lowest_time[vertex] > lowest_time[neighbor_vertex]) {
+            lowest_time[vertex] = lowest_time[neighbor_vertex];
+        }
+
+        if (parent_vertex != UINT64_MAX && discovery_time[vertex] <= lowest_time[neighbor_vertex]) {
+            is_articulation = true;
+        }
+    }
+    if (parent_vertex == UINT64_MAX && child_count > 1) {
+        is_articulation = true;
+    }
+    if (is_articulation) {
+        articulations.push_back(vertex + 1);
+    }
+}
+
+std::vector<uint64_t> find_articulations_dfs_v2(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
+    uint64_t vertex_count = adjacency_matrix.size();
+    std::vector<bool> visited(vertex_count, false);
+    uint64_t current_time = 0;
+    std::vector<uint64_t> discovery_time(vertex_count, 0);
+    std::vector<uint64_t> lowest_time(vertex_count, 0);
+    std::vector<uint64_t> articulations;
+
+    for (uint64_t vertex = 0; vertex < vertex_count; vertex++) {
+        if (!visited[vertex]) {
+            dfs_articulations(adjacency_matrix, vertex, UINT64_MAX, visited, current_time, discovery_time, lowest_time, articulations);
+        }
+    }
+
+    articulations.shrink_to_fit();
+    std::sort(articulations.begin(), articulations.end());
+    return articulations;
+}
diff --git a/graph.h b/graph.h
deleted file mode 100644
index a88aa53..0000000
--- a/graph.h
+++ /dev/null
@@ -1,71 +0,0 @@
-#ifndef GRAPH_H
-#define GRAPH_H
-
-#include <stdint.h>
-
-void random_adjacency(const uint64_t vertex_count,
-                      uint64_t matrix[vertex_count][vertex_count]);
-
-void calculate_distance_matrix(const uint64_t vertex_count,
-                               const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                               uint64_t distance_matrix[vertex_count][vertex_count]);
-
-// returns 1 if it is a disconnected graph and sets all values in eccentricities to UINT64_MAX
-int get_eccentricities(const uint64_t vertex_count,
-                       const uint64_t distance_matrix[vertex_count][vertex_count],
-                       uint64_t eccentricities[vertex_count]);
-
-uint64_t get_radius(const uint64_t vertex_count,
-                    const uint64_t eccentricities[vertex_count]);
-
-uint64_t get_diameter(const uint64_t vertex_count,
-                      const uint64_t eccentricities[vertex_count]);
-
-void get_centre(const uint64_t vertex_count,
-                const uint64_t eccentricities[vertex_count],
-                const uint64_t radius,
-                uint64_t centre[vertex_count]);
-
-void calculate_path_matrix(const uint64_t vertex_count,
-                           const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                           uint64_t path_matrix[vertex_count][vertex_count]);
-
-void find_components_basic(const uint64_t vertex_count,
-                           const uint64_t path_matrix[vertex_count][vertex_count],
-                           uint64_t components[vertex_count][vertex_count]);
-
-void find_components_dfs(const uint64_t vertex_count,
-                         const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                         uint64_t components[vertex_count][vertex_count]);
-
-void find_bridges_basic(const uint64_t vertex_count,
-                        const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                        const uint64_t components[vertex_count][vertex_count],
-                        uint64_t bridges[vertex_count][2]);
-
-void find_bridges_dfs_v1(const uint64_t vertex_count,
-                         const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                         const uint64_t components[vertex_count][vertex_count],
-                         uint64_t bridges[vertex_count][2]);
-
-void find_bridges_dfs_v2(const uint64_t vertex_count,
-                         const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                         const uint64_t components[vertex_count][vertex_count],
-                         uint64_t bridges[vertex_count][2]);
-
-void find_articulations_basic(const uint64_t vertex_count,
-                              const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                              const uint64_t components[vertex_count][vertex_count],
-                              uint64_t articulations[vertex_count]);
-
-void find_articulations_dfs_v1(const uint64_t vertex_count,
-                               const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                               const uint64_t components[vertex_count][vertex_count],
-                               uint64_t articulations[vertex_count]);
-
-void find_articulations_dfs_v2(const uint64_t vertex_count,
-                               const uint64_t adjacency_matrix[vertex_count][vertex_count],
-                               const uint64_t components[vertex_count][vertex_count],
-                               uint64_t articulations[vertex_count]);
-
-#endif
diff --git a/graph.hpp b/graph.hpp
new file mode 100644
index 0000000..a618fa2
--- /dev/null
+++ b/graph.hpp
@@ -0,0 +1,33 @@
+#pragma once
+#include <cstdint>
+#include <vector>
+
+std::vector<std::vector<uint64_t>> random_adjacency(const uint64_t vertex_count);
+
+std::vector<std::vector<uint64_t>> calculate_distance_matrix(const std::vector<std::vector<uint64_t>>& adjacency_matrix);
+
+std::vector<uint64_t> get_eccentricities(const std::vector<std::vector<uint64_t>>& distance_matrix);
+
+uint64_t get_radius(const std::vector<uint64_t>& eccentricities);
+
+uint64_t get_diameter(const std::vector<uint64_t>& eccentricities);
+
+std::vector<uint64_t> get_centre(const std::vector<uint64_t>& eccentricities);
+
+std::vector<std::vector<uint64_t>> calculate_path_matrix(const std::vector<std::vector<uint64_t>>& adjacency_matrix);
+
+std::vector<std::vector<uint64_t>> find_components_basic(const std::vector<std::vector<uint64_t>>& path_matrix);
+
+std::vector<std::vector<uint64_t>> find_components_dfs(const std::vector<std::vector<uint64_t>>& adjacency_matrix);
+
+std::vector<std::vector<uint64_t>> find_bridges_basic(const std::vector<std::vector<uint64_t>>& adjacency_matrix);
+
+std::vector<std::vector<uint64_t>> find_bridges_dfs_v1(const std::vector<std::vector<uint64_t>>& adjacency_matrix);
+
+std::vector<std::vector<uint64_t>> find_bridges_dfs_v2(const std::vector<std::vector<uint64_t>>& adjacency_matrix);
+
+std::vector<uint64_t> find_articulations_basic(const std::vector<std::vector<uint64_t>>& adjacency_matrix);
+
+std::vector<uint64_t> find_articulations_dfs_v1(const std::vector<std::vector<uint64_t>>& adjacency_matrix);
+
+std::vector<uint64_t> find_articulations_dfs_v2(const std::vector<std::vector<uint64_t>>& adjacency_matrix);
diff --git a/main.c b/main.c
deleted file mode 100644
index 0c0f8f9..0000000
--- a/main.c
+++ /dev/null
@@ -1,280 +0,0 @@
-#include "graph.h"
-#include "matrix.h"
-#include <limits.h>
-#include <stdint.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <time.h>
-
-void benchmark_gemm() {
-    const uint64_t vertex_count1 = 1024;
-    uint64_t (*matrix1)[vertex_count1] = malloc(vertex_count1 * vertex_count1 * sizeof(uint64_t));
-
-    const uint64_t vertex_count2 = 1024;
-    uint64_t (*matrix2)[vertex_count2] = malloc(vertex_count2 * vertex_count2 * sizeof(uint64_t));
-
-    uint64_t (*new_matrix)[vertex_count2] = malloc(vertex_count1 * vertex_count2 * sizeof(uint64_t));
-
-    double elapsed_time = 0.0;
-    const uint64_t iterations = 10;
-    clock_t start_time;
-
-    for (uint64_t i = 0; i < iterations; i++) {
-        random_adjacency(vertex_count1, matrix1);
-        random_adjacency(vertex_count2, matrix2);
-
-        start_time = clock();
-
-        gemm_basic(vertex_count1, vertex_count1, matrix1,
-                              vertex_count2, vertex_count2, matrix2,
-                              new_matrix);
-
-        elapsed_time += (double)(clock() - start_time) / CLOCKS_PER_SEC;
-    }
-
-    printf("%lu iterations of gemm_basic took roughly %f seconds\n", iterations, elapsed_time);
-    printf("An iteration of gemm_basic took on average roughly %f seconds\n", elapsed_time/iterations);
-
-    free(matrix1);
-    free(matrix2);
-    free(new_matrix);
-}
-
-void benchmark_find_components() {
-    const uint64_t vertex_count = 100;
-    uint64_t adjacency_matrix[vertex_count][vertex_count];
-    uint64_t components[vertex_count][vertex_count];
-    uint64_t path_matrix[vertex_count][vertex_count];
-
-    double elapsed_time = 0.0;
-    const uint64_t iterations = 100;
-    clock_t start_time;
-
-    for (uint64_t i = 0; i < iterations; i++) {
-        random_adjacency(vertex_count, adjacency_matrix);
-
-        start_time = clock();
-
-        calculate_path_matrix(vertex_count, adjacency_matrix, path_matrix);
-        find_components_basic(vertex_count, path_matrix, components);
-
-        elapsed_time += (double)(clock() - start_time) / CLOCKS_PER_SEC;
-    }
-
-    printf("%lu iterations of find_components_basic took roughly %f seconds\n", iterations, elapsed_time);
-    printf("An iteration of find_components_basic took on average roughly %f seconds\n", elapsed_time/iterations);
-
-    elapsed_time = 0.0;
-
-    for (uint64_t i = 0; i < iterations; i++) {
-        random_adjacency(vertex_count, adjacency_matrix);
-
-        start_time = clock();
-
-        find_components_dfs(vertex_count, adjacency_matrix, components);
-
-        elapsed_time += (double)(clock() - start_time) / CLOCKS_PER_SEC;
-    }
-
-    printf("%lu iterations of find_components_dfs took roughly %f seconds\n", iterations, elapsed_time);
-    printf("An iteration of find_components_dfs took on average roughly %f seconds\n", elapsed_time/iterations);
-}
-
-void test_with_basic() {
-    const uint64_t vertex_count = 24;
-    uint64_t adjacency_matrix[vertex_count][vertex_count];
-    uint64_t distance_matrix[vertex_count][vertex_count];
-    uint64_t path_matrix[vertex_count][vertex_count];
-    uint64_t eccentricities[vertex_count];
-    uint64_t radius, diameter, centre[vertex_count];
-    uint64_t components[vertex_count][vertex_count];
-    uint64_t bridges[vertex_count][2];
-    uint64_t articulations[vertex_count];
-
-    if (read_csv("csv/24n.csv", vertex_count, vertex_count, adjacency_matrix) == 1) {
-        return;
-    }
-
-    calculate_distance_matrix(vertex_count, adjacency_matrix, distance_matrix);
-    get_eccentricities(vertex_count, distance_matrix, eccentricities);
-    radius = get_radius(vertex_count, eccentricities);
-    diameter = get_diameter(vertex_count, eccentricities);
-    get_centre(vertex_count, eccentricities, radius, centre);
-    calculate_path_matrix(vertex_count, adjacency_matrix, path_matrix);
-    find_components_basic(vertex_count, path_matrix, components);
-    find_bridges_basic(vertex_count, adjacency_matrix, components, bridges);
-    find_articulations_basic(vertex_count, adjacency_matrix, components, articulations);
-
-    puts("adjacency_matrix:");
-    print_matrix(vertex_count, vertex_count, adjacency_matrix);
-
-    puts("\ndistance_matrix:");
-    print_matrix(vertex_count, vertex_count, distance_matrix);
-
-    puts("\neccentricities:");
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        printf("\tVertex %lu: %lu\n", index + 1,  eccentricities[index]);
-    }
-
-    printf("\nradius: %lu", radius);
-    printf("\ndiameter: %lu", diameter);
-    puts("\ncentre:");
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        if (centre[index] == 1) {
-            printf("\tVertex %lu\n", index + 1);
-        }
-    }
-
-    puts("\npath_matrix:");
-    print_matrix(vertex_count, vertex_count, path_matrix);
-
-    puts("\ncomponents:");
-    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-        int empty = 1;
-
-        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-            if (components[row_index][column_index] != 0) {
-                empty = 0;
-            }
-        }
-
-        if (empty) {
-            continue;
-        }
-
-        printf("\tComponent %lu: {", row_index + 1);
-        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-            if (components[row_index][column_index] != 0) {
-                printf("%lu, ", components[row_index][column_index]);
-            }
-        }
-        puts("}");
-    }
-
-    puts("\nbridges:");
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        if (bridges[index][0] != 0) {
-            printf("\tBridge %lu: {%lu, %lu}\n", index + 1, bridges[index][0], bridges[index][1]);
-        }
-    }
-
-    puts("\narticulations:");
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        if (articulations[index] != 0) {
-            printf("\tVertex %lu\n", articulations[index]);
-        }
-    }
-}
-
-void test_with_dfs() {
-    const uint64_t vertex_count = 24;
-    uint64_t adjacency_matrix[vertex_count][vertex_count];
-    uint64_t distance_matrix[vertex_count][vertex_count];
-    uint64_t eccentricities[vertex_count];
-    uint64_t radius, diameter;
-    uint64_t centre[vertex_count];
-    uint64_t components[vertex_count][vertex_count];
-    uint64_t bridges[vertex_count][2];
-    uint64_t articulations[vertex_count];
-
-    if (read_csv("csv/24n.csv", vertex_count, vertex_count, adjacency_matrix) == 1) {
-        return;
-    }
-
-    /*
-    const uint64_t vertex_count = 1500;
-    uint64_t (*adjacency_matrix)[vertex_count] = malloc(vertex_count * vertex_count * sizeof(uint64_t));
-    uint64_t (*distance_matrix)[vertex_count] = malloc(vertex_count * vertex_count * sizeof(uint64_t));
-    uint64_t *eccentricities = malloc(vertex_count * sizeof(uint64_t));
-    uint64_t radius, diameter;
-    uint64_t *centre = malloc(vertex_count * sizeof(uint64_t));
-    uint64_t (*components)[vertex_count] = malloc(vertex_count * vertex_count * sizeof(uint64_t));
-    uint64_t (*bridges)[vertex_count] = malloc(vertex_count * 2 * sizeof(uint64_t));
-    uint64_t *articulations = malloc(vertex_count * sizeof(uint64_t));
-
-    random_adjacency(vertex_count, adjacency_matrix);
-    */
-
-    calculate_distance_matrix(vertex_count, adjacency_matrix, distance_matrix);
-    get_eccentricities(vertex_count, distance_matrix, eccentricities);
-    radius = get_radius(vertex_count, eccentricities);
-    diameter = get_diameter(vertex_count, eccentricities);
-    get_centre(vertex_count, eccentricities, radius, centre);
-    find_components_dfs(vertex_count, adjacency_matrix, components);
-    find_bridges_dfs_v2(vertex_count, adjacency_matrix, components, bridges);
-    find_articulations_dfs_v2(vertex_count, adjacency_matrix, components, articulations);
-
-    puts("\nadjacency_matrix:");
-    print_matrix(vertex_count, vertex_count, adjacency_matrix);
-
-    puts("\ndistance_matrix:");
-    print_matrix(vertex_count, vertex_count, distance_matrix);
-
-    puts("\neccentricities:");
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        printf("\tVertex %lu: %lu\n", index + 1,  eccentricities[index]);
-    }
-
-    printf("\nradius: %lu", radius);
-    printf("\ndiameter: %lu", diameter);
-    puts("\ncentre:");
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        if (centre[index] == 1) {
-            printf("\tVertex %lu\n", index + 1);
-        }
-    }
-
-    puts("\ncomponents:");
-    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
-        int empty = 1;
-
-        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-            if (components[row_index][column_index] != 0) {
-                empty = 0;
-            }
-        }
-
-        if (empty) {
-            continue;
-        }
-
-        printf("\tComponent %lu: {", row_index + 1);
-        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
-            if (components[row_index][column_index] != 0) {
-                printf("%lu, ", column_index + 1);
-            }
-        }
-        puts("}");
-    }
-
-    uint64_t bridge_number = 1;
-    puts("\nbridges:");
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        if (bridges[index][0] != 0) {
-            printf("\tBridge %lu: {%lu, %lu}\n", bridge_number, bridges[index][0], bridges[index][1]);
-            bridge_number++;
-        }
-    }
-
-    puts("\narticulations:");
-    for (uint64_t index = 0; index < vertex_count; index++) {
-        if (articulations[index] != 0) {
-            printf("\tVertex %lu\n", articulations[index]);
-        }
-    }
-
-    // free(adjacency_matrix);
-    // free(distance_matrix);
-    // free(eccentricities);
-    // free(centre);
-    // free(components);
-    // free(bridges);
-    // free(articulations);
-}
-
-int main(void) {
-    // test_with_basic();
-    test_with_dfs();
-    // benchmark_gemm();
-    // benchmark_find_components();
-}
diff --git a/main.cpp b/main.cpp
new file mode 100644
index 0000000..4bdeee4
--- /dev/null
+++ b/main.cpp
@@ -0,0 +1,247 @@
+#include "graph.hpp"
+#include "matrix.hpp"
+#include <cstdint>
+#include <limits.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <vector>
+#include <iostream>
+
+
+void benchmark_gemm() {
+    const uint64_t vertex_count1 = 1024;
+    std::vector<std::vector<uint64_t>> matrix1;
+
+    const uint64_t vertex_count2 = 1024;
+    std::vector<std::vector<uint64_t>> matrix2;
+
+    std::vector<std::vector<uint64_t>> new_matrix;
+
+    const uint64_t iterations = 10;
+    double elapsed_time = 0.0;
+    // clock_t start_time;
+
+    for (uint64_t i = 0; i < iterations; i++) {
+        matrix1 = random_adjacency(vertex_count1);
+        matrix2 = random_adjacency(vertex_count2);
+
+        // start_time = clock();
+
+        new_matrix = gemm_basic(matrix1, matrix2);
+
+        // elapsed_time += (double)(clock() - start_time) / CLOCKS_PER_SEC;
+    }
+
+    printf("%lu iterations of gemm_basic took roughly %f seconds\n", iterations, elapsed_time);
+    printf("An iteration of gemm_basic took on average roughly %f seconds\n", elapsed_time/iterations);
+}
+
+
+
+void benchmark_find_components() {
+    const uint64_t vertex_count = 100;
+    std::vector<std::vector<uint64_t>> adjacency_matrix;
+    std::vector<std::vector<uint64_t>> components;
+    std::vector<std::vector<uint64_t>> path_matrix;
+
+    const uint64_t iterations = 100;
+    double elapsed_time = 0.0;
+    // clock_t start_time;
+
+    for (uint64_t i = 0; i < iterations; i++) {
+        adjacency_matrix = random_adjacency(vertex_count);
+
+        // start_time = clock();
+
+        path_matrix = calculate_path_matrix(adjacency_matrix);
+        components = find_components_basic(path_matrix);
+
+        // elapsed_time += (double)(clock() - start_time) / CLOCKS_PER_SEC;
+    }
+
+    printf("%lu iterations of find_components_basic took roughly %f seconds\n", iterations, elapsed_time);
+    printf("An iteration of find_components_basic took on average roughly %f seconds\n", elapsed_time/iterations);
+
+    elapsed_time = 0.0;
+
+    for (uint64_t i = 0; i < iterations; i++) {
+        adjacency_matrix = random_adjacency(vertex_count);
+
+        // start_time = clock();
+
+        components = find_components_dfs(adjacency_matrix);
+
+        // elapsed_time += (double)(clock() - start_time) / CLOCKS_PER_SEC;
+    }
+
+    printf("%lu iterations of find_components_dfs took roughly %f seconds\n", iterations, elapsed_time);
+    printf("An iteration of find_components_dfs took on average roughly %f seconds\n", elapsed_time/iterations);
+}
+
+
+void test_with_basic() {
+    std::vector<std::vector<uint64_t>> adjacency_matrix = read_csv("csv/24n.csv");
+    std::vector<std::vector<uint64_t>> distance_matrix = calculate_distance_matrix(adjacency_matrix);
+    std::vector<std::vector<uint64_t>> path_matrix = calculate_path_matrix(adjacency_matrix);
+    std::vector<uint64_t> eccentricities = get_eccentricities(distance_matrix);
+    uint64_t radius = get_radius(eccentricities);
+    uint64_t diameter = get_diameter(eccentricities);
+    std::vector<uint64_t> centre = get_centre(eccentricities);
+    std::vector<std::vector<uint64_t>> components = find_components_basic(path_matrix);
+    std::vector<std::vector<uint64_t>> bridges = find_bridges_basic(adjacency_matrix);
+    std::vector<uint64_t> articulations = find_articulations_basic(adjacency_matrix);
+
+    std::cout << "\nadjacency_matrix:\n";
+    print_matrix(adjacency_matrix);
+
+    std::cout << "\ndistance_matrix:\n";
+    print_matrix(distance_matrix);
+
+    std::cout << "\neccentricities: {";
+    for (uint64_t index = 0; index < eccentricities.size(); index++) {
+        std::cout << eccentricities[index];
+
+        if (index < eccentricities.size() - 1) {
+            std::cout << ", ";
+        }
+    }
+    std::cout << "}";
+
+    std::cout << "\nradius: " << radius;
+    std::cout << "\ndiameter: " << diameter;
+    std::cout << "\ncentre: {";
+    for (uint64_t index = 0; index < centre.size(); index++) {
+        std::cout << centre[index];
+
+        if (index < centre.size() - 1) {
+            std::cout << ", ";
+        }
+    }
+    std::cout << "}";
+
+    std::cout << "\ncomponents:";
+    for (uint64_t component_index = 0; component_index < components.size(); component_index++) {
+        std::cout << "\ncomponent " << component_index + 1 << ": {";
+
+        for (uint64_t index = 0; index < components[component_index].size(); index++) {
+            std::cout << components[component_index][index];
+
+            if (index < components[component_index].size() - 1) {
+                std::cout << ", ";
+            }
+        }
+        std::cout << "}";
+    }
+
+    std::cout << "\nbridges:";
+    for (uint64_t bridge_index = 0; bridge_index < bridges.size(); bridge_index++) {
+        std::cout << "\nbridge " << bridge_index + 1 << ": {";
+
+        for (uint64_t index = 0; index < bridges[bridge_index].size(); index++) {
+            std::cout << bridges[bridge_index][index];
+
+            if (index < bridges[bridge_index].size() - 1) {
+                std::cout << ", ";
+            }
+        }
+        std::cout << "}";
+    }
+
+    std::cout << "\narticulations: {";
+    for (uint64_t index = 0; index < articulations.size(); index++) {
+        std::cout << articulations[index];
+
+        if (index < articulations.size() - 1) {
+            std::cout << ", ";
+        }
+    }
+    std::cout << "}\n";
+}
+
+void test_with_dfs() {
+    // const uint64_t vertex_count = 100;
+    // std::vector<std::vector<uint64_t>> adjacency_matrix = random_adjacency(vertex_count);
+    std::vector<std::vector<uint64_t>> adjacency_matrix = read_csv("csv/24n.csv");
+    std::vector<std::vector<uint64_t>> distance_matrix = calculate_distance_matrix(adjacency_matrix);
+    std::vector<std::vector<uint64_t>> path_matrix = calculate_path_matrix(adjacency_matrix);
+    std::vector<uint64_t> eccentricities = get_eccentricities(distance_matrix);
+    uint64_t radius = get_radius(eccentricities);
+    uint64_t diameter = get_diameter(eccentricities);
+    std::vector<uint64_t> centre = get_centre(eccentricities);
+    std::vector<std::vector<uint64_t>> components = find_components_dfs(adjacency_matrix);
+    std::vector<std::vector<uint64_t>> bridges = find_bridges_dfs_v2(adjacency_matrix);
+    std::vector<uint64_t> articulations = find_articulations_dfs_v2(adjacency_matrix);
+
+    std::cout << "\nadjacency_matrix:\n";
+    print_matrix(adjacency_matrix);
+
+    std::cout << "\ndistance_matrix:\n";
+    print_matrix(distance_matrix);
+
+    std::cout << "\neccentricities: {";
+    for (uint64_t index = 0; index < eccentricities.size(); index++) {
+        std::cout << eccentricities[index];
+
+        if (index < eccentricities.size() - 1) {
+            std::cout << ", ";
+        }
+    }
+    std::cout << "}";
+
+    std::cout << "\nradius: " << radius;
+    std::cout << "\ndiameter: " << diameter;
+    std::cout << "\ncentre: {";
+    for (uint64_t index = 0; index < centre.size(); index++) {
+        std::cout << centre[index];
+
+        if (index < centre.size() - 1) {
+            std::cout << ", ";
+        }
+    }
+    std::cout << "}";
+
+    std::cout << "\ncomponents: ";
+    for (uint64_t component_index = 0; component_index < components.size(); component_index++) {
+        std::cout << "\n\tcomponent " << component_index + 1 << ": {";
+
+        for (uint64_t index = 0; index < components[component_index].size(); index++) {
+            std::cout << components[component_index][index];
+
+            if (index < components[component_index].size() - 1) {
+                std::cout << ", ";
+            }
+        }
+        std::cout << "}";
+    }
+
+    std::cout << "\nbridges:";
+    for (uint64_t bridge_index = 0; bridge_index < bridges.size(); bridge_index++) {
+        std::cout << "\n\tbridge " << bridge_index + 1 << ": {";
+
+        for (uint64_t index = 0; index < bridges[bridge_index].size(); index++) {
+            std::cout << bridges[bridge_index][index];
+
+            if (index < bridges[bridge_index].size() - 1) {
+                std::cout << ", ";
+            }
+        }
+        std::cout << "}";
+    }
+
+    std::cout << "\narticulations: {";
+    for (uint64_t index = 0; index < articulations.size(); index++) {
+        std::cout << articulations[index];
+
+        if (index < articulations.size() - 1) {
+            std::cout << ", ";
+        }
+    }
+    std::cout << "}\n";
+}
+
+int main(void) {
+    // test_with_basic();
+    test_with_dfs();
+    // benchmark_gemm();
+    // benchmark_find_components();
+}
diff --git a/matrix.c b/matrix.c
deleted file mode 100644
index 3fc459a..0000000
--- a/matrix.c
+++ /dev/null
@@ -1,71 +0,0 @@
-#include <stddef.h>
-#include <stdint.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <time.h>
-
-void print_matrix(const uint64_t row_length, const uint64_t column_length, const uint64_t matrix[row_length][column_length]) {
-    for (uint64_t column_index=0; column_index < column_length; column_index++) {
-        for (uint64_t row_index=0; row_index < row_length; row_index++) {
-            printf("%lu ", matrix[row_index][column_index]);
-        }
-        puts("");
-    }
-}
-
-void gemm_basic(const uint64_t row_length1, const uint64_t column_length1, const uint64_t matrix1[row_length1][column_length1],
-                          const uint64_t row_length2, const uint64_t column_length2, const uint64_t matrix2[row_length2][column_length2],
-                          uint64_t output_matrix[row_length1][column_length2]) {
-    uint64_t sum;
-
-    for (uint64_t i = 0; i < row_length1; i++) {
-        for (uint64_t j = 0; j < column_length2; j++) {
-            sum = 0;
-
-            for (uint64_t k = 0; k < row_length1; k++) {
-                sum += matrix1[i][k] * matrix2[k][j];
-            }
-
-            output_matrix[i][j] = sum;
-        }
-    }
-}
-
-int read_csv(const char *file_name, const uint64_t row_length, const uint64_t column_length, uint64_t output_matrix[row_length][column_length]) {
-    FILE *file_ptr;
-    uint64_t bufsize = row_length*2+1; // have to account for delimiters
-    char buffer[bufsize];
-    char *value, *file_line;
-    uint64_t row_index = 0, column_index = 0;
-
-    file_ptr = fopen(file_name, "r");
-    if (file_ptr == NULL) {
-        puts("Unable to open csv");
-        return 1;
-    }
-
-    while ((file_line = fgets(buffer, bufsize, file_ptr)) != NULL) {
-        // This shit is just needed and I dont know why
-        file_line[strcspn(file_line, "\n")] = 0;
-        value = strtok(file_line, ";,");
-
-        // for some reason there are two NULLs at the end of a line
-        // and I dont wanna increment the column_index
-        if (value == NULL) {
-            continue;
-        }
-
-        while (value != NULL) {
-            output_matrix[row_index++][column_index] = strtoul(value, NULL, 0);
-            value = strtok(NULL, ";,");
-        }
-
-        row_index = 0;
-        column_index++;
-    }
-
-    fclose(file_ptr);
-
-    return 0;
-}
diff --git a/matrix.cpp b/matrix.cpp
new file mode 100644
index 0000000..eae5414
--- /dev/null
+++ b/matrix.cpp
@@ -0,0 +1,70 @@
+#include <cstdint>
+#include <stdlib.h>
+#include <iostream>
+#include <vector>
+#include <fstream>
+#include <sstream>
+
+void print_matrix(const std::vector<std::vector<uint64_t>>& matrix) {
+    for (uint64_t column_index = 0; column_index < matrix[0].size(); column_index++) {
+        for (uint64_t row_index = 0; row_index < matrix.size(); row_index++) {
+            std::cout << matrix[row_index][column_index] << " ";
+        }
+        std::cout << "\n";
+    }
+}
+
+std::vector<std::vector<uint64_t>> gemm_basic(const std::vector<std::vector<uint64_t>>& matrix1,
+                                              const std::vector<std::vector<uint64_t>>& matrix2) {
+    std::vector<std::vector<uint64_t>> output_matrix(matrix1.size(), std::vector<uint64_t>(matrix2[0].size(), 0));
+    uint64_t sum;
+
+    for (uint64_t i = 0; i < matrix1.size(); i++) {
+        for (uint64_t j = 0; j < matrix2[0].size(); j++) {
+            sum = 0;
+
+            for (uint64_t k = 0; k < matrix1.size(); k++) {
+                sum += matrix1[i][k] * matrix2[k][j];
+            }
+
+            output_matrix[i][j] = sum;
+        }
+    }
+
+    return output_matrix;
+}
+
+std::vector<std::vector<uint64_t>> read_csv(const std::string file_name) {
+    std::vector<std::vector<uint64_t>> output_matrix;
+    std::vector<uint64_t> row;
+    std::string line;
+    uint64_t column_index = 0;
+
+    std::ifstream input_stream(file_name);
+
+    while (std::getline(input_stream, line)) {
+        uint64_t vertex_count = (line.length()+1)/2;
+        output_matrix.resize(vertex_count);
+        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
+            output_matrix[row_index].push_back(0);
+        }
+
+        row.clear();
+
+        std::stringstream ss(line);
+        std::string token;
+        char delimiter = ';';
+
+        while (getline(ss, token, delimiter)) {
+            row.push_back(strtoul(token.c_str(), NULL, 0));
+        }
+
+        for (uint64_t row_index = 0; row_index < row.size(); row_index++) {
+            output_matrix[row_index][column_index] = row[row_index];
+        }
+
+        column_index++;
+    }
+
+    return output_matrix;
+}
diff --git a/matrix.h b/matrix.h
deleted file mode 100644
index 2666d32..0000000
--- a/matrix.h
+++ /dev/null
@@ -1,28 +0,0 @@
-#ifndef MATRIX_H
-#define MATRIX_H
-
-#include <stdint.h>
-
-void print_matrix(const uint64_t row_length,
-                  const uint64_t column_length,
-                  const uint64_t matrix[row_length][column_length]);
-
-/*
-    First two matrices will be multiplied and
-    restult will be written to output_matrix.
-    Matrix size requirements are as specified in the parameters.
-*/
-void gemm_basic(const uint64_t row_length1,
-               const uint64_t column_length1,
-               const uint64_t matrix1[row_length1][column_length1],
-               const uint64_t row_length2,
-               const uint64_t column_length2,
-               const uint64_t matrix2[row_length2][column_length2],
-               uint64_t output_matrix[row_length1][column_length2]);
-
-int read_csv(char *file_name,
-             uint64_t row_length,
-             uint64_t column_length,
-             uint64_t output_matrix[row_length][column_length]);
-
-#endif
diff --git a/matrix.hpp b/matrix.hpp
new file mode 100644
index 0000000..91e572e
--- /dev/null
+++ b/matrix.hpp
@@ -0,0 +1,15 @@
+#pragma once
+#include <cstdint>
+#include <vector>
+#include <string>
+
+void print_matrix(const std::vector<std::vector<uint64_t>>& matrix);
+
+/*
+    First two matrices will be multiplied and
+    restult will be written to output_matrix.
+*/
+std::vector<std::vector<uint64_t>> gemm_basic(const std::vector<std::vector<uint64_t>>& matrix1,
+                                              const std::vector<std::vector<uint64_t>>& matrix2);
+
+std::vector<std::vector<uint64_t>> read_csv(const std::string file_name);