diff --git a/setup.nqp b/setup.nqp
index a2689eb..82cabe2 100644
--- a/setup.nqp
+++ b/setup.nqp
@@ -35,7 +35,8 @@ sub MAIN(@argv) {
         run_test_harness(%PLA);
         pir::exit__vI(0);
     }
-    setup_dynpmc_sources(%PLA);
+    setup_c_library_files(%PLA);
+    setup_dynpmc(%PLA);
     setup_testlib(%PLA);
     setup_nqp_bootstrapper(%PLA);
 
@@ -49,14 +50,14 @@ sub get_args() {
 }
 
 sub probe_for_cblas(%PLA) {
-    if probe_include("cblas.h", :verbose(0)) {
+    if probe_include("cblas.h", :verbose(1)) {
         pir::say("Cannot find cblas.h\nPlease install libatlas-base-dev");
         pir::exit__vI(1);
     }
 }
 
 sub system_linker_settings(%PLA) {
-    %PLA{'dynpmc_cflags'} := '-g';
+    %PLA{'dynpmc_cflags'} := '-g -Isrc/include/';
     my %config := get_config();
     my $osname := %config{'osname'};
     if $osname eq 'linux' {
@@ -67,7 +68,10 @@ sub system_linker_settings(%PLA) {
             my $searchloc := $_;
             my $test_ldd := pir::spawnw__IS('ldd ' ~ $searchloc);
             if $test_ldd == 0 {
-                %PLA{'dynpmc_ldflags'} := %searches{$searchloc};
+                my $flags := %PLA{'dynpmc_ldflags'};
+                my $libflags := %searches{$searchloc};
+                $flags := ~$flags ~ $libflags;
+                %PLA{'dynpmc_ldflags'} := $flags;
                 return;
             }
         }
@@ -84,7 +88,60 @@ sub run_test_harness(%PLA) {
     pir::exit(+$result);
 }
 
-sub setup_dynpmc_sources(%PLA) {
+sub setup_c_library_files(%PLA) {
+    register_step_before('build', compile_c_library_files);
+    register_step_after('clean', clean_c_library_files);
+
+    my @cfiles := <
+        src/lib/matrix_common
+        src/lib/math_common
+    >;
+
+    %PLA{'cc_dynpmc'}{'linalg_group'} := @cfiles;
+
+    my @files := %PLA{'cc_dynpmc'}{'linalg_group'};
+    my $obj := get_obj();
+    my $ldflags := "";
+    for @files {
+        my $file := $_;
+        my $objfile := $file ~ $obj;
+        $ldflags := $ldflags ~ " ";
+        $ldflags := $ldflags ~ $objfile;
+    }
+    my $flags := %PLA{'dynpmc_ldflags'};
+    #$flags := $flags ~ $ldflags;
+    $flags := ~$flags ~ " " ~ $ldflags;
+    %PLA{'dynpmc_ldflags'} := $flags;
+}
+
+sub compile_c_library_files(*%args) {
+    my @files := %args{'cc_dynpmc'}{'linalg_group'};
+    my $obj := get_obj();
+    my $cflags := get_cflags() ~ " -g -Isrc/include";
+    my $ldflags := "";
+    my $libheader := "src/include/pla_matrix_library.h";
+    for @files {
+        my $file := $_;
+        my $objfile := $file ~ $obj;
+        my $cfile := $file ~ ".c";
+        unless newer($objfile, [$cfile]) {
+            __compile_cc($objfile, $cfile, $cflags);
+        }
+        $ldflags := $ldflags ~ " " ~ $objfile;
+    }
+}
+
+sub clean_c_library_files(*%kv) {
+    my @files := %kv{'cc_dynpmc'}{'linalg_group'};
+    my $obj := get_obj();
+    for @files {
+        my $file := $_;
+        my $objfile := $file ~ $obj;
+        unlink($objfile, :verbose(1));
+    }
+}
+
+sub setup_dynpmc(%PLA) {
     %PLA{'dynpmc'}{'linalg_group'} := <
         src/pmc/nummatrix2d.pmc
         src/pmc/pmcmatrix2d.pmc
diff --git a/src/include/pla.h b/src/include/pla.h
new file mode 100644
index 0000000..3ea05f7
--- /dev/null
+++ b/src/include/pla.h
@@ -0,0 +1,10 @@
+#ifndef _PLA_H_
+#define _PLA_H_
+
+#include <cblas.h>
+#include <math.h>
+#include <parrot/parrot.h>
+#include "pla_matrix_types.h"
+#include "pla_matrix_library.h"
+
+#endif /* _PLA_H_ */
diff --git a/src/include/pla_matrix_library.h b/src/include/pla_matrix_library.h
new file mode 100644
index 0000000..43bcfb9
--- /dev/null
+++ b/src/include/pla_matrix_library.h
@@ -0,0 +1,10 @@
+#ifndef _PLA_MATRIX_LIBRARY_H_
+#define _PLA_MATRIX_LIBRARY_H_
+
+INTVAL floats_are_equal(FLOATVAL a, FLOATVAL b);
+void get_complex_value_from_pmc(PARROT_INTERP, PMC * value, FLOATVAL * real, FLOATVAL * imag);
+
+void intkey_to_coords(PARROT_INTERP, const INTVAL rows, const INTVAL cols, const INTVAL key, INTVAL * row, INTVAL * col);
+
+
+#endif /* _PLA_MATRIX_LIBRARY_H */
diff --git a/src/pmc/pla_matrix_types.h b/src/include/pla_matrix_types.h
similarity index 99%
rename from src/pmc/pla_matrix_types.h
rename to src/include/pla_matrix_types.h
index 49d844b..2089b06 100644
--- a/src/pmc/pla_matrix_types.h
+++ b/src/include/pla_matrix_types.h
@@ -1,8 +1,6 @@
 #ifndef _PLA_MATRIX_TYPES_H_
 #define _PLA_MATRIX_TYPES_H_
 
-#include <cblas.h>
-
 extern INTVAL __PLA_Type_NumMatrix2D;
 extern INTVAL __PLA_Type_ComplexMatrix2D;
 extern INTVAL __PLA_Type_PMCMatrix2D;
diff --git a/src/pmc/charmatrix2d.pmc b/src/pmc/charmatrix2d.pmc
index 761f89e..463d6eb 100644
--- a/src/pmc/charmatrix2d.pmc
+++ b/src/pmc/charmatrix2d.pmc
@@ -1,4 +1,4 @@
-#include "pla_matrix_types.h"
+#include "pla.h"
 
 #define ALLOCATE_STORAGE(s) (char *)mem_sys_allocate_zeroed(s * sizeof (char))
 
diff --git a/src/pmc/complexmatrix2d.pmc b/src/pmc/complexmatrix2d.pmc
index b43cad1..2daa97d 100644
--- a/src/pmc/complexmatrix2d.pmc
+++ b/src/pmc/complexmatrix2d.pmc
@@ -1,5 +1,4 @@
-#include "pla_matrix_types.h"
-#include <cblas.h>
+#include "pla.h"
 
 #define ALLOCATE_STORAGE(s) (FLOATVAL *)mem_sys_allocate_zeroed(s * sizeof (FLOATVAL) * 2)
 #define PLATYPENAME "ComplexMatrix2D"
@@ -8,34 +7,6 @@
 
 INTVAL __PLA_Type_ComplexMatrix2D;
 
-static void
-get_complex_value_from_pmc(PARROT_INTERP, PMC * value, FLOATVAL * real, FLOATVAL * imag)
-{
-    switch (value->vtable->base_type) {
-        case enum_class_String:
-            value = Parrot_pmc_new_init(interp, enum_class_Complex, value);
-        case enum_class_Complex:
-            *real = VTABLE_get_number_keyed_int(interp, value, 0);
-            *imag = VTABLE_get_number_keyed_int(interp, value, 1);
-            break;
-        case enum_class_Float:
-            *real = VTABLE_get_number(interp, value);
-            *imag = 0.0;
-            break;
-        case enum_class_Integer:
-            {
-                const INTVAL _r = VTABLE_get_integer(interp, value);
-                *real = (FLOATVAL)_r;
-            }
-            *imag = 0.0;
-            break;
-        default:
-            /* TODO: We should support HLL-mapped types here. */
-            Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS,
-                PLATYPENAME ": cannot set unknown PMC type");
-    }
-}
-
 /* If the matrix is lazily transposed, actually transpose the physical memory
    layout. This is necessary for calculations, especially BLAS calculations,
    which aren't lazy-transpose-aware. */
@@ -114,6 +85,7 @@ resize_matrix(PARROT_INTERP, PMC * self, INTVAL row, INTVAL col)
     FLOATVAL * new_s = ALLOCATE_STORAGE(newsize);
 
     INTVAL i, j;
+
     for (i = 0; i < old_rows; i++) {
         for (j = 0; j < old_cols; j++) {
             R_ITEM_XY_ROWMAJOR(new_s, new_rows, new_cols, i, j) =
@@ -375,26 +347,6 @@ add_matrices(PARROT_INTERP, PMC * A, PMC * B, INTVAL sub)
     }
 }
 
-
-static void
-intkey_to_coords(PARROT_INTERP, const INTVAL rows, const INTVAL cols,
-    const INTVAL key, INTVAL * row, INTVAL * col)
-{
-    if (rows == 0) {
-        Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS,
-            PLATYPENAME ": Index out of bounds.");
-    } else {
-        const INTVAL r = key / rows;
-        const INTVAL c = key % rows;
-        if (key < 0 || r > rows || c > cols) {
-            Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS,
-                PLATYPENAME ": Index out of bounds.");
-        }
-        *row = r;
-        *col = c;
-    }
-}
-
 pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix {
     ATTR FLOATVAL * storage;
     ATTR INTVAL rows;
diff --git a/src/pmc/matrixproxy.pmc b/src/pmc/matrixproxy.pmc
index f83d3d5..3c33408 100644
--- a/src/pmc/matrixproxy.pmc
+++ b/src/pmc/matrixproxy.pmc
@@ -1,4 +1,4 @@
-#include "pla_matrix_types.h"
+#include "pla.h"
 
 #define DECLATTRS(p, a) Parrot_MatrixProxy_attributes * const (a) = \
     (Parrot_MatrixProxy_attributes *)((p)->data)
diff --git a/src/pmc/nummatrix2d.pmc b/src/pmc/nummatrix2d.pmc
index f6f1b00..d54ab40 100644
--- a/src/pmc/nummatrix2d.pmc
+++ b/src/pmc/nummatrix2d.pmc
@@ -1,21 +1,10 @@
-#include "pla_matrix_types.h"
-#include <math.h>
-#include <cblas.h>
+#include "pla.h"
 
 #define ALLOCATE_STORAGE(s) (FLOATVAL *)mem_sys_allocate_zeroed(s * sizeof (FLOATVAL))
 #define PLATYPENAME "NumMatrix2D"
 #define DECLATTRS(p, a) Parrot_NumMatrix2D_attributes * const (a) = \
     (Parrot_NumMatrix2D_attributes *)((p)->data)
 
-/* Relative error between two large numbers */
-#define FLOATVAL_EPSILON_REL 0.0001
-
-/* Absolute difference between two small numbers */
-#define FLOATVAL_EPSILON_ABS 0.00001
-
-/* Above we do relative error, below we use abs */
-#define FLOATVAL_REL_THRESHOLD 10
-
 INTVAL __PLA_Type_NumMatrix2D;
 
 /* Resize the matrix internal storage to be able to hold a point at position
@@ -236,30 +225,6 @@ call_gemm(PARROT_INTERP, FLOATVAL alpha, PMC * A, PMC *B, FLOATVAL beta, PMC *C)
     );
 }
 
-/* Check if two floats are equalish. Above a certain threshold we want to do a
-   relative error calculation because precision between subsequent large
-   numbers can be much larger than standard epsilon values. Below the threshold
-   we want to do an absolute calculation because as the values approach zero the
-   error value quickly begins to diverge, even for very similar values. The
-   values for absolute difference, relative error, and threshold need to be
-   tuned. */
-static INTVAL
-floats_are_equal(FLOATVAL a, FLOATVAL b)
-{
-    FLOATVAL diff = fabs(a - b);
-    if (a == b)
-        return 1;
-    if (diff > FLOATVAL_EPSILON_ABS)
-        return 0;
-    if (a > FLOATVAL_REL_THRESHOLD) {
-        FLOATVAL err = diff / b;
-        if (err > FLOATVAL_EPSILON_REL)
-            return 0;
-    }
-    return 1;
-}
-
-
 pmclass NumMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix {
     ATTR FLOATVAL * storage;
     ATTR INTVAL rows;
diff --git a/src/pmc/pmcmatrix2d.pmc b/src/pmc/pmcmatrix2d.pmc
index 69e0911..2306697 100644
--- a/src/pmc/pmcmatrix2d.pmc
+++ b/src/pmc/pmcmatrix2d.pmc
@@ -1,4 +1,4 @@
-#include "pla_matrix_types.h"
+#include "pla.h"
 
 #define ALLOCATE_STORAGE(s) (PMC **)mem_sys_allocate_zeroed(s * sizeof (PMC *))
 #define PLATYPENAME "PMCMatrix2D"
@@ -55,8 +55,6 @@ resize_matrix(PARROT_INTERP, PMC * self, INTVAL row, INTVAL col)
         mem_sys_free(old_s);
 }
 
-
-
 static void
 init_from_pmc_array(PARROT_INTERP, PMC * self, INTVAL rows_size, INTVAL cols_size, PMC * values) {
     DECLATTRS(self, attrs);