diff --git a/src/matrix/3x3matrix.c b/src/matrix/3x3matrix.c
new file mode 100644
index 0000000000000000000000000000000000000000..63208317472e6184b773f31475d4c637360c7e68
--- /dev/null
+++ b/src/matrix/3x3matrix.c
@@ -0,0 +1,120 @@
+#include "3x3matrix.h"
+
+_3x3Matrix add_3x3(_3x3Matrix m1, _3x3Matrix m2){
+    _3x3Matrix ma = {
+        m1.a11 + m2.a11,m1.a12 + m2.a12,m1.a13 + m2.a13,
+        m1.a21 + m2.a21,m1.a22 + m2.a22,m1.a23 + m2.a23,
+        m1.a31 + m2.a31,m1.a32 + m2.a32,m1.a33 + m2.a33
+    };
+
+    return ma;
+}
+
+_3x3Matrix minus_3x3(_3x3Matrix m1, _3x3Matrix m2) {
+    _3x3Matrix ma = {
+        m1.a11 - m2.a11,m1.a12 - m2.a12,m1.a13 - m2.a13,
+        m1.a21 - m2.a21,m1.a22 - m2.a22,m1.a23 - m2.a23,
+        m1.a31 - m2.a31,m1.a32 - m2.a32,m1.a33 - m2.a33
+    };
+
+    return ma;
+}
+_3x3Matrix multiply_3x3(_3x3Matrix m1, _3x3Matrix m2) {
+    _3x3Matrix ma = {
+            (m1.a11*m2.a11)+(m1.a12*m2.a21)+(m1.a13*m2.a31),
+            (m1.a11*m2.a12)+(m1.a12*m2.a22)+(m1.a13*m2.a32),
+            (m1.a11*m2.a13)+(m1.a12*m2.a23)+(m1.a13*m2.a33),
+            (m1.a21*m2.a11)+(m1.a22*m2.a21)+(m1.a23*m2.a31),
+            (m1.a21*m2.a12)+(m1.a22*m2.a22)+(m1.a23*m2.a32),
+            (m1.a21*m2.a13)+(m1.a22*m2.a23)+(m1.a23*m2.a33),
+            (m1.a31*m2.a11)+(m1.a32*m2.a21)+(m1.a33*m2.a31),
+            (m1.a31*m2.a12)+(m1.a32*m2.a22)+(m1.a33*m2.a32),
+            (m1.a31*m2.a13)+(m1.a32*m2.a23)+(m1.a33*m2.a32)
+    };
+    return ma;
+}
+
+_3x3Matrix scale_3x3(_3x3Matrix m1, float scale){
+    _3x3Matrix ma = {
+        m1.a11 * scale,m1.a12 * scale,m1.a13 * scale,
+        m1.a21 * scale,m1.a22 * scale,m1.a23 * scale,
+        m1.a31 * scale,m1.a32 * scale,m1.a33 * scale
+    };
+    return ma;
+}
+float det_3x3(_3x3Matrix m){
+    return (m.a12*m.a23*m.a31) + (m.a13*m.a21*m.a32) + (m.a33 * m.a11 * m.a22)
+            -(m.a33*m.a12*m.a21) - (m.a11*m.a23*m.a32) - (m.a13*m.a31*m.a22);
+}
+
+_3x3Matrix inverse_3x3(_3x3Matrix m) {
+    float detM = det_3x3(m);
+    if(detM){
+        _3x3Matrix adjM = {
+                (m._mI*m._mE) - (m._mF*m._mH),
+                (m._mC*m._mH) - (m._mI*m._mB),
+                (m._mB*m._mF) - (m._mC*m._mE),
+                (m._mF*m._mG) - (m._mI*m._mD),
+                (m._mA*m._mI) - (m._mC*m._mG),
+                (m._mC*m._mD) - (m._mA*m._mF),
+                (m._mD*m._mH) - (m._mG*m._mE),
+                (m._mB*m._mG) - (m._mA*m._mH),
+                (m._mA*m._mE) - (m._mB*m._mD)
+        };
+        return scale_3x3(adjM, 1/detM);
+    }
+    return m;
+}
+_3x1Matrix transpose_3x3(_3x3Matrix m1, _3x1Matrix m2){
+    _3x1Matrix ma = {
+            (m1.a11 * m2.a11) + (m1.a12 * m2.a21) + (m1.a13 * m2.a31),
+            (m1.a21 * m2.a11) + (m1.a22 * m2.a21) + (m1.a23 * m2.a31),
+            (m1.a31 * m2.a11) + (m1.a32 * m2.a21) + (m1.a33 * m2.a31)
+
+    };
+    return ma;
+}
+unsigned char equals(_3x3Matrix m1, _3x3Matrix m2){
+    if(m1.a11 == m2.a11 &&
+        m1.a12 == m2.a12 &&
+        m1.a13 == m2.a13 &&
+        m1.a21 == m2.a21 &&
+        m1.a22 == m2.a22 &&
+        m1.a23 == m2.a23 &&
+        m1.a31 == m2.a31 &&
+        m1.a32 == m2.a32 &&
+        m1.a33 == m2.a33)
+        return 42;
+    return 0;
+}
+
+_3x3Matrix combine(_3x1Matrix m1, _3x1Matrix m2, _3x1Matrix m3){
+    _3x3Matrix m = {m1.a11, m2.a11, m3.a11,
+                    m1.a21, m2.a21, m3.a21,
+                    m1.a31, m2.a31, m3.a31};
+    return m;
+}
+
+_3x1Matrix add_3x1(_3x1Matrix m1, _3x1Matrix m2){
+    _3x1Matrix mT = {m1.a11 + m2.a11,
+                     m1.a21 + m2.a21,
+                     m1.a31 + m2.a31};
+
+    return mT;
+}
+
+_3x1Matrix minus_3x1(_3x1Matrix m1, _3x1Matrix m2){
+    _3x1Matrix mT = {m1.a11 - m2.a11,
+                     m1.a21 - m2.a21,
+                     m1.a31 - m2.a31};
+
+    return mT;
+}
+
+_3x1Matrix crossProduct(_3x1Matrix m1, _3x1Matrix m2){
+    _3x1Matrix m = {(m1.a21*m2.a31) - (m1.a31*m2.a21),
+                    (m1.a31*m2.a11) - (m1.a11*m2.a31),
+                    (m1.a11*m2.a21) - (m1.a21*m2.a11)};
+
+    return m;
+}