diff --git a/include/J3ML/LinearAlgebra/Angle2D.h b/include/J3ML/LinearAlgebra/Angle2D.h
index 0a0c2e5..8f7ea94 100644
--- a/include/J3ML/LinearAlgebra/Angle2D.h
+++ b/include/J3ML/LinearAlgebra/Angle2D.h
@@ -1,4 +1,5 @@
 #pragma once
+#include <J3ML/J3ML.h>
 
 
 namespace J3ML::LinearAlgebra {
@@ -19,4 +20,4 @@ namespace J3ML::LinearAlgebra {
         }
     };
 
-}
\ No newline at end of file
+}
diff --git a/include/J3ML/LinearAlgebra/Matrices.inl b/include/J3ML/LinearAlgebra/Matrices.inl
index fe19c73..aee8b80 100644
--- a/include/J3ML/LinearAlgebra/Matrices.inl
+++ b/include/J3ML/LinearAlgebra/Matrices.inl
@@ -3,11 +3,143 @@
 /// Template Parameterized (Generic) Matrix Functions.
 
 #include <J3ML/LinearAlgebra/Quaternion.h>
-
+#include <J3ML/J3ML.h>
 
 
 namespace J3ML::LinearAlgebra {
 
+    template <typename Matrix>
+    bool InverseMatrix(Matrix &mat, float epsilon)
+    {
+        Matrix inversed = Matrix::Identity;
+
+        const int nc = std::min<int>(Matrix::Rows, Matrix::Cols);
+
+        for (int column = 0; column < nc; ++column)
+        {
+            // find the row i with i >= j such that M has the largest absolute value.
+            int greatest = column;
+            float greatestVal = std::abs(mat[greatest][column]);
+            for (int i = column+1; i < Matrix::Rows; i++)
+            {
+                float val = std::abs(mat[i][column]);
+                if (val > greatestVal) {
+                    greatest = i;
+                    greatestVal = val;
+                }
+            }
+
+            if (greatestVal < epsilon) {
+                mat = inversed;
+                return false;
+            }
+
+            // exchange rows
+            if (greatest != column) {
+                inversed.SwapRows(greatest, column);
+                mat.SwapRows(greatest, column);
+            }
+
+            // multiply rows
+            assert(!Math::EqualAbs(mat[column][column], 0.f, epsilon));
+            float scale = 1.f / mat[column][column];
+            inversed.ScaleRow(column, scale);
+            mat.ScaleRow(column, scale);
+
+            // add rows
+            for (int i = 0; i < column; i++) {
+                inversed.SetRow(i, inversed.Row(i) - inversed.Row(column) * mat[i][column]);
+                mat.SetRow(i, mat.Row(i) - mat.Row(column) * mat[i][column]);
+            }
+
+            for (int i = column + 1; i < Matrix::Rows; i++) {
+                inversed.SetRow(i, inversed.Row(i) - inversed.Row(column) * mat[i][column]);
+                mat.SetRow(i, mat.Row(i) - mat.Row(column) * mat[i][column]);
+            }
+        }
+        mat = inversed;
+        return true;
+    }
+
+
+    /// Computes the LU-decomposition on the given square matrix.
+    /// @return True if the composition was successful, false otherwise. If the return value is false, the contents of the output matrix are unspecified.
+    template <typename Matrix>
+    bool LUDecomposeMatrix(const Matrix &mat, Matrix &lower, Matrix &upper)
+    {
+        lower = Matrix::Identity;
+        upper = Matrix::Zero;
+
+        for (int i = 0; i < Matrix::Rows; ++i)
+        {
+            for (int col = i; col < Matrix::Cols; ++col)
+            {
+                upper[i][col] = mat[i][col];
+                for (int k = 0; k < i; ++k)
+                    upper[i][col] -= lower[i][k] * upper[k][col];
+            }
+            for (int row = i+1; row < Matrix::Rows; ++row)
+            {
+                lower[row][i] = mat[row][i];
+                for (int k = 0; k < i; ++k)
+                    lower[row][i] -= lower[row][k] * upper[k][i];
+                if (Math::EqualAbs(upper[i][i], 0.f))
+                    return false;
+                lower[row][i] /= upper[i][i];
+            }
+        }
+        return true;
+    }
+
+    /// Computes the Cholesky decomposition on the given square matrix *on the real domain*.
+    /// @return True if successful, false otherwise. If the return value is false, the contents of the output matrix are uspecified.
+    template <typename Matrix>
+    bool CholeskyDecomposeMatrix(const Matrix &mat, Matrix& lower)
+    {
+        lower = Matrix::Zero;
+        for (int i = 0; i < Matrix::Rows; ++i)
+        {
+            for (int j = 0; j < i; ++i)
+            {
+                lower[i][j] = mat[i][j];
+                for (int k = 0; k < j; ++k)
+                    lower[i][j] -= lower[i][j] * lower[j][k];
+                if (Math::EqualAbs(lower[j][j], 0.f))
+                    return false;
+                lower[i][j] /= lower[j][j];
+            }
+            lower[i][i] = mat[i][i];
+            if (lower[i][i])
+                return false;
+            for (int k = 0; k < i; ++k)
+                lower[i][i] -= lower[i][k] * lower[i][k];
+            lower[i][i] = std::sqrt(lower[i][i]);
+        }
+
+        return false;
+    }
+
+
+    template<typename Matrix>
+    void SetMatrixRotatePart(Matrix &m, const Quaternion &q) {
+        // See https://www.geometrictools.com/Documentation/LinearAlgebraicQuaternions.pdf .
+
+        assert(q.IsNormalized(1e-3f));
+        const float x = q.x;
+        const float y = q.y;
+        const float z = q.z;
+        const float w = q.w;
+        m[0][0] = 1 - 2 * (y * y + z * z);
+        m[0][1] = 2 * (x * y - z * w);
+        m[0][2] = 2 * (x * y + y * w);
+        m[1][0] = 2 * (x * y + z * w);
+        m[1][1] = 1 - 2 * (x * x + z * z);
+        m[1][2] = 2 * (y * z - x * w);
+        m[2][0] = 2 * (x * z - y * w);
+        m[2][1] = 2 * (y * z + x * w);
+        m[2][2] = 1 - 2 * (x * x + y * y);
+    }
+
     /** Sets the top-left 3x3 area of the matrix to the rotation matrix about the X-axis. Elements
         outside the top-left 3x3 area are ignored. This matrix rotates counterclockwise if multiplied
         in the order M*v, and clockwise if rotated in the order v*M.
diff --git a/include/J3ML/LinearAlgebra/Matrix2x2.h b/include/J3ML/LinearAlgebra/Matrix2x2.h
index 8d77231..4d23adb 100644
--- a/include/J3ML/LinearAlgebra/Matrix2x2.h
+++ b/include/J3ML/LinearAlgebra/Matrix2x2.h
@@ -1,9 +1,9 @@
 #pragma once
 
 
-#include <J3ML/LinearAlgebra/Vector2.h>
+
 #include <J3ML/LinearAlgebra/Common.h>
-#include <J3ML/LinearAlgebra/Matrices.inl>
+
 
 namespace J3ML::LinearAlgebra {
     class Matrix2x2 {
diff --git a/include/J3ML/LinearAlgebra/Matrix3x3.h b/include/J3ML/LinearAlgebra/Matrix3x3.h
index 502da64..1ff0a12 100644
--- a/include/J3ML/LinearAlgebra/Matrix3x3.h
+++ b/include/J3ML/LinearAlgebra/Matrix3x3.h
@@ -1,7 +1,7 @@
 #pragma once
 
 
-#include <J3ML/LinearAlgebra/Common.h>
+
 #include <J3ML/LinearAlgebra/Vector2.h>
 #include <J3ML/LinearAlgebra/Vector3.h>
 #include <J3ML/LinearAlgebra/Quaternion.h>
@@ -13,25 +13,7 @@ using namespace J3ML::Algorithm;
 
 namespace J3ML::LinearAlgebra {
 
-    template<typename Matrix>
-    void SetMatrixRotatePart(Matrix &m, const Quaternion &q) {
-        // See https://www.geometrictools.com/Documentation/LinearAlgebraicQuaternions.pdf .
 
-        assert(q.IsNormalized(1e-3f));
-        const float x = q.x;
-        const float y = q.y;
-        const float z = q.z;
-        const float w = q.w;
-        m[0][0] = 1 - 2 * (y * y + z * z);
-        m[0][1] = 2 * (x * y - z * w);
-        m[0][2] = 2 * (x * y + y * w);
-        m[1][0] = 2 * (x * y + z * w);
-        m[1][1] = 1 - 2 * (x * x + z * z);
-        m[1][2] = 2 * (y * z - x * w);
-        m[2][0] = 2 * (x * z - y * w);
-        m[2][1] = 2 * (y * z + x * w);
-        m[2][2] = 1 - 2 * (x * x + y * y);
-    }
 
     /// A 3-by-3 matrix for linear transformations of 3D geometry.
     /* This can represent any kind of linear transformations of 3D geometry, which include
diff --git a/include/J3ML/LinearAlgebra/Matrix4x4.h b/include/J3ML/LinearAlgebra/Matrix4x4.h
index 813e93b..b0e381c 100644
--- a/include/J3ML/LinearAlgebra/Matrix4x4.h
+++ b/include/J3ML/LinearAlgebra/Matrix4x4.h
@@ -1,11 +1,7 @@
 #pragma once
 
-#include <J3ML/LinearAlgebra/Common.h>
 
-#include <J3ML/LinearAlgebra/Matrix3x3.h>
-#include <J3ML/LinearAlgebra/Quaternion.h>
-#include <J3ML/LinearAlgebra/Vector4.h>
-#include <J3ML/LinearAlgebra/Matrices.inl>
+#include <J3ML/LinearAlgebra/Common.h>
 #include <J3ML/Algorithm/RNG.h>
 
 #include <algorithm>
@@ -14,114 +10,7 @@ using namespace J3ML::Algorithm;
 
 namespace J3ML::LinearAlgebra {
 
-    template <typename Matrix>
-    bool InverseMatrix(Matrix &mat, float epsilon)
-    {
-        Matrix inversed = Matrix::Identity;
 
-        const int nc = std::min<int>(Matrix::Rows, Matrix::Cols);
-
-        for (int column = 0; column < nc; ++column)
-        {
-            // find the row i with i >= j such that M has the largest absolute value.
-            int greatest = column;
-            float greatestVal = std::abs(mat[greatest][column]);
-            for (int i = column+1; i < Matrix::Rows; i++)
-            {
-                float val = std::abs(mat[i][column]);
-                if (val > greatestVal) {
-                    greatest = i;
-                    greatestVal = val;
-                }
-            }
-
-            if (greatestVal < epsilon) {
-                mat = inversed;
-                return false;
-            }
-
-            // exchange rows
-            if (greatest != column) {
-                inversed.SwapRows(greatest, column);
-                mat.SwapRows(greatest, column);
-            }
-
-            // multiply rows
-            assert(!Math::EqualAbs(mat[column][column], 0.f, epsilon));
-            float scale = 1.f / mat[column][column];
-            inversed.ScaleRow(column, scale);
-            mat.ScaleRow(column, scale);
-
-            // add rows
-            for (int i = 0; i < column; i++) {
-                inversed.SetRow(i, inversed.Row(i) - inversed.Row(column) * mat[i][column]);
-                mat.SetRow(i, mat.Row(i) - mat.Row(column) * mat[i][column]);
-            }
-
-            for (int i = column + 1; i < Matrix::Rows; i++) {
-                inversed.SetRow(i, inversed.Row(i) - inversed.Row(column) * mat[i][column]);
-                mat.SetRow(i, mat.Row(i) - mat.Row(column) * mat[i][column]);
-            }
-        }
-        mat = inversed;
-        return true;
-    }
-
-
-    /// Computes the LU-decomposition on the given square matrix.
-    /// @return True if the composition was successful, false otherwise. If the return value is false, the contents of the output matrix are unspecified.
-    template <typename Matrix>
-    bool LUDecomposeMatrix(const Matrix &mat, Matrix &lower, Matrix &upper)
-    {
-        lower = Matrix::Identity;
-        upper = Matrix::Zero;
-
-        for (int i = 0; i < Matrix::Rows; ++i)
-        {
-            for (int col = i; col < Matrix::Cols; ++col)
-            {
-                upper[i][col] = mat[i][col];
-                for (int k = 0; k < i; ++k)
-                    upper[i][col] -= lower[i][k] * upper[k][col];
-            }
-            for (int row = i+1; row < Matrix::Rows; ++row)
-            {
-                lower[row][i] = mat[row][i];
-                for (int k = 0; k < i; ++k)
-                    lower[row][i] -= lower[row][k] * upper[k][i];
-                if (Math::EqualAbs(upper[i][i], 0.f))
-                    return false;
-                lower[row][i] /= upper[i][i];
-            }
-        }
-        return true;
-    }
-
-    /// Computes the Cholesky decomposition on the given square matrix *on the real domain*.
-    /// @return True if successful, false otherwise. If the return value is false, the contents of the output matrix are uspecified.
-    template <typename Matrix>
-    bool CholeskyDecomposeMatrix(const Matrix &mat, Matrix& lower)
-    {
-        lower = Matrix::Zero;
-        for (int i = 0; i < Matrix::Rows; ++i)
-        {
-            for (int j = 0; j < i; ++i)
-            {
-                lower[i][j] = mat[i][j];
-                for (int k = 0; k < j; ++k)
-                    lower[i][j] -= lower[i][j] * lower[j][k];
-                if (Math::EqualAbs(lower[j][j], 0.f))
-                    return false;
-                lower[i][j] /= lower[j][j];
-            }
-            lower[i][i] = mat[i][i];
-            if (lower[i][i])
-                return false;
-            for (int k = 0; k < i; ++k)
-                lower[i][i] -= lower[i][k] * lower[i][k];
-            lower[i][i] = std::sqrt(lower[i][i]);
-        }
-    }
 
     /// @brief A 4-by-4 matrix for affine transformations and perspective projections of 3D geometry.
     /// This matrix can represent the most generic form of transformations for 3D objects,
@@ -416,7 +305,8 @@ namespace J3ML::LinearAlgebra {
 
         /// Returns only the first three elements of the given column.
         Vector3 GetColumn3(int index) const;
-        Vector3 Column3(int index) const { return GetColumn3(index);}
+        Vector3 Column3(int index) const;
+
         Vector3 Col3(int i) const;
 
         /// Returns the scaling performed by this matrix. This function assumes taht the last row is [0 0 0 1].
@@ -539,11 +429,7 @@ namespace J3ML::LinearAlgebra {
         /// @note This function assumes that this matrix does not contain projection (the fourth row of this matrix is [0 0 0 1]).
         /// @note This function assumes that this matrix has orthogonal basis vectors (row and column vector sets are orthogonal).
         /// @note This function does not remove reflection (-1 scale along some axis).
-        void RemoveScale()
-        {
-            float tx = Row3(0).Normalize();
-            float ty = Row3(1).Normalize();
-        }
+        void RemoveScale();
 
 
         /// Decomposes this matrix to translate, rotate, and scale parts.
diff --git a/include/J3ML/LinearAlgebra/Quaternion.h b/include/J3ML/LinearAlgebra/Quaternion.h
index c64f6cc..074652a 100644
--- a/include/J3ML/LinearAlgebra/Quaternion.h
+++ b/include/J3ML/LinearAlgebra/Quaternion.h
@@ -1,10 +1,6 @@
 #pragma once
 
 #include <J3ML/LinearAlgebra/Common.h>
-#include <J3ML/LinearAlgebra/Matrix3x3.h>
-#include <J3ML/LinearAlgebra/Matrix4x4.h>
-#include <J3ML/LinearAlgebra/Vector4.h>
-#include <J3ML/LinearAlgebra/AxisAngle.h>
 #include <J3ML/Algorithm/RNG.h>
 #include <cmath>
 
@@ -104,7 +100,7 @@ namespace J3ML::LinearAlgebra
 
 
         /// Returns a uniformly random unitary quaternion.
-        static Quaternion RandomRotation(J3ML::Algorithm::RNG &rng);
+        static Quaternion RandomRotation(RNG &rng);
     public:
         void SetFromAxisAngle(const Vector3 &vector3, float between);
 
@@ -216,8 +212,12 @@ namespace J3ML::LinearAlgebra
         // The rotation q2 is applied first before q1.
         Quaternion operator * (const Quaternion& rhs) const;
 
+        // Unsafe
         Quaternion operator * (float scalar) const;
 
+        // Unsafe
+        Quaternion operator / (float scalar) const;
+
         // Transforms the given vector by this Quaternion.
         Vector3 operator * (const Vector3& rhs) const;
         Vector4 operator * (const Vector4& rhs) const;
@@ -225,10 +225,9 @@ namespace J3ML::LinearAlgebra
         // Divides a quaternion by another. Divison "a / b" results in a quaternion that rotates the orientation b to coincide with orientation of
         Quaternion operator / (const Quaternion& rhs) const;
         Quaternion operator + (const Quaternion& rhs) const;
-        Quaternion operator / (float scalar) const
-        {
-            assert(!Math::EqualAbs(scalar, 0.f));
-        }
+
+
+
         Quaternion operator + () const;
         Quaternion operator - () const;
 
diff --git a/include/J3ML/LinearAlgebra/Vector2.h b/include/J3ML/LinearAlgebra/Vector2.h
index ce32ff7..214a757 100644
--- a/include/J3ML/LinearAlgebra/Vector2.h
+++ b/include/J3ML/LinearAlgebra/Vector2.h
@@ -1,5 +1,6 @@
 #pragma once
-#include <J3ML/J3ML.h>
+
+
 #include <cstddef>
 
 namespace J3ML::LinearAlgebra {
diff --git a/include/J3ML/LinearAlgebra/Vector4.h b/include/J3ML/LinearAlgebra/Vector4.h
index 817ca4b..bc46f5e 100644
--- a/include/J3ML/LinearAlgebra/Vector4.h
+++ b/include/J3ML/LinearAlgebra/Vector4.h
@@ -1,7 +1,9 @@
 #pragma once
 
+#include <cstddef>
+#include <cmath>
 
-#include <J3ML/LinearAlgebra/Vector3.h>
+#include <J3ML/LinearAlgebra/Common.h>
 
 namespace J3ML::LinearAlgebra {
     class Vector4 {
@@ -20,10 +22,7 @@ namespace J3ML::LinearAlgebra {
         {
             return &x;
         }
-        Vector3 XYZ() const
-        {
-            return {x, y, z};
-        }
+        Vector3 XYZ() const;
 
         float GetX() const { return x; }
         float GetY() const { return y; }
diff --git a/src/J3ML/LinearAlgebra/Matrix2x2.cpp b/src/J3ML/LinearAlgebra/Matrix2x2.cpp
index 1fd5323..6fcbaa6 100644
--- a/src/J3ML/LinearAlgebra/Matrix2x2.cpp
+++ b/src/J3ML/LinearAlgebra/Matrix2x2.cpp
@@ -1,4 +1,5 @@
 #include <J3ML/LinearAlgebra/Matrix2x2.h>
+#include <J3ML/LinearAlgebra/Vector2.h>
 
 namespace J3ML::LinearAlgebra {
 
diff --git a/src/J3ML/LinearAlgebra/Matrix3x3.cpp b/src/J3ML/LinearAlgebra/Matrix3x3.cpp
index e0856ff..124f1bd 100644
--- a/src/J3ML/LinearAlgebra/Matrix3x3.cpp
+++ b/src/J3ML/LinearAlgebra/Matrix3x3.cpp
@@ -1,5 +1,8 @@
 #include <J3ML/LinearAlgebra/Matrix3x3.h>
+#include <J3ML/LinearAlgebra/Matrix4x4.h>
 #include <J3ML/LinearAlgebra/Matrices.inl>
+#include <J3ML/LinearAlgebra/Vector3.h>
+#include <J3ML/LinearAlgebra/Vector4.h>
 #include <cmath>
 
 namespace J3ML::LinearAlgebra {
@@ -771,6 +774,10 @@ namespace J3ML::LinearAlgebra {
         x[2] = b[2] / v22;
         x[1] = b[1] - x[2] * v12;
         x[0] = b[0] - x[2] * v02 - x[1] * v01;
+
+
+
+        return true;
     }
 
     void Matrix3x3::ScaleCol(int col, float scalar) {
diff --git a/src/J3ML/LinearAlgebra/Matrix4x4.cpp b/src/J3ML/LinearAlgebra/Matrix4x4.cpp
index 26c467c..44de7cd 100644
--- a/src/J3ML/LinearAlgebra/Matrix4x4.cpp
+++ b/src/J3ML/LinearAlgebra/Matrix4x4.cpp
@@ -1,5 +1,9 @@
 #include <J3ML/LinearAlgebra/Matrix4x4.h>
 #include <J3ML/LinearAlgebra/Vector4.h>
+#include <J3ML/LinearAlgebra/Matrix3x3.h>
+#include <J3ML/LinearAlgebra/Quaternion.h>
+#include <J3ML/LinearAlgebra/Matrices.inl>
+
 
 namespace J3ML::LinearAlgebra {
 
@@ -400,6 +404,8 @@ namespace J3ML::LinearAlgebra {
         return Vector3{At(0, index), At(1, index), At(2, index)};
     }
 
+    Vector3 Matrix4x4::Column3(int index) const { return GetColumn3(index);}
+
     Vector2 Matrix4x4::operator*(const Vector2 &rhs) const { return this->Transform(rhs);}
 
     Vector3 Matrix4x4::operator*(const Vector3 &rhs) const { return this->Transform(rhs);}
@@ -617,6 +623,11 @@ namespace J3ML::LinearAlgebra {
         return {GetColumn3(0).Length(), GetColumn3(1).Length(), GetColumn3(2).Length()};
     }
 
+    void Matrix4x4::RemoveScale() {
+        float tx = Row3(0).Normalize();
+        float ty = Row3(1).Normalize();
+    }
+
     bool Matrix4x4::IsColOrthogonal(float epsilon) const {
         return IsColOrthogonal3(epsilon);
     }
@@ -804,6 +815,8 @@ namespace J3ML::LinearAlgebra {
         p[1][0] = 0;           p[1][1] = 2.f * n / v; p[1][2] = 0;         p[1][3] = 0.f;
         p[2][0] = 0;           p[2][1] = 0;           p[2][2] = f / (f-n); p[2][3] = n * f / (n-f);
         p[3][0] = 0;           p[3][1] = 0;           p[3][2] = 1.f;       p[3][3] = 0.f;
+
+        return p;
     }
 
     Matrix4x4 Matrix4x4::D3DPerspProjRH(float n, float f, float h, float v) {
diff --git a/src/J3ML/LinearAlgebra/Quaternion.cpp b/src/J3ML/LinearAlgebra/Quaternion.cpp
index 33d5539..667c7fa 100644
--- a/src/J3ML/LinearAlgebra/Quaternion.cpp
+++ b/src/J3ML/LinearAlgebra/Quaternion.cpp
@@ -3,6 +3,7 @@
 #include <J3ML/LinearAlgebra/Matrix3x3.h>
 #include <J3ML/LinearAlgebra/Matrix4x4.h>
 #include <J3ML/LinearAlgebra/Quaternion.h>
+#include <J3ML/LinearAlgebra/AxisAngle.h>
 
 namespace J3ML::LinearAlgebra {
 
@@ -69,6 +70,12 @@ namespace J3ML::LinearAlgebra {
         return Quaternion(x * scalar, y * scalar, z * scalar, w * scalar);
     }
 
+    Quaternion Quaternion::operator/(float scalar) const {
+        assert(!Math::EqualAbs(scalar, 0.f));
+
+        return *this * (1.f / scalar);
+    }
+
     Quaternion Quaternion::operator+() const { return *this; }
 
     Quaternion::Quaternion() {}
diff --git a/src/J3ML/LinearAlgebra/Vector4.cpp b/src/J3ML/LinearAlgebra/Vector4.cpp
index c90bd84..c4471f2 100644
--- a/src/J3ML/LinearAlgebra/Vector4.cpp
+++ b/src/J3ML/LinearAlgebra/Vector4.cpp
@@ -150,6 +150,10 @@ Vector4 Vector4::operator-(const Vector4& rhs) const
         return *this;
     }
 
+    Vector3 Vector4::XYZ() const {
+        return {x, y, z};
+    }
+
     bool Vector4::Equals(const Vector4 &rhs, float epsilon) const {
         return std::abs(x - rhs.x) < epsilon &&
                std::abs(y - rhs.y) < epsilon &&