tb_matrix.hpp

#ifndef TurtleBrains_Matrix_hpp
#define TurtleBrains_Matrix_hpp

#include <turtle_brains/math/tb_vector.hpp>
#include <turtle_brains/math/tb_angle.hpp>
#include <turtle_brains/core/tb_error.hpp>
#include <turtle_brains/core/tb_defines.hpp>

#include <cmath>
#include <cstring> //memcpy / memset

// TODO: TurtleBrains: Make Identity() a constexpr static member function. Readup on constexpr

//Pre-Multiply (vector on left), Bases Contiguous:
//  Xx, Xy, Xz, 0
//  Yx, Yy, Yz, 0
//  Zx, Zy, Zz, 0
//  Tx, Ty, Tz, 1
//  Xx, Xy, Xz, 0, Yx, Yy, Yz, 0, Zx, Zy, Zz, 0, Tx, Ty, Tz, 1

//
//  The above is close to what I desire, but the data in matrix test kit is showing that we have something much different
//  probably looking like (basis/rotation 3x3 unconfirmed): Xx, Xy, Xz, Tx, Yx, Yy, Yz, Ty, Zx, Zy Zz, Tz, 0, 0, 0, 1
//  as that is how our test data is being shoved into the matrix object... 2020-11-25
//
//  This was found because FastInverse() did not work in doing physics; FastInverse() was using the final COLUMN as the
//  placement of translation/position and the above shows the final row being used. Changing FastInverse() breaks the
//  previous tests, so what gives?
//
//  To make this matter even more interesting, the CreateTranslation() function does put the position components in the
//  final ROW as it shows in the above 'image'. The mat(1, 3) operator is specified by: column, row and memory is accessed
//  by components[column + (row * 4)]; CreateTranslation() then puts translation at mat(0, 3), mat(1, 3), mat(2, 3) aka: Final row.
//
//  Intrets says pre-multiply, vector*matrix, should have the translation on the last row.
//
//  There was only a single instance of FastInverse() called in Rally of Rockets, and that instance didn't even use
//  the positional data. This tells me that the test likely contains bad data and we just got really lucky in Rally of Rockets
//  to use TransformNormal() instead of TransformCoordinate(). There were no occurances of FastInverse() in Turbo Boom!,
//  Track Builder, and LudumDare47.
//
//  2020-11-26 - FastInverse() was fixed and the test data was also fixed. Essentially the 16 skeletal animation matrices
//  provided by shakesoda had translation in the final column while our matrix actually does expect it in the final ROW.
//  While I was adjust the translation data (only moving translation from final column to final row), shakesoda said this
//  was suspect and I was doing it wrong, but never gave more information on why that was the case. Since every other test
//  was passing except for the known broken FastInverse, moving translation data in those matrices fixed everything once
//  FastInverse() was also fixed. (FastInverse had to swap final column for final row and also pre/post mulitply direction.
//

namespace TurtleBrains
{
    namespace Math
    {
        namespace Unstable
        {
            class Quaternion;
        };

        enum ColumnMajorMatrix
        {
            ColumnMajor = 1,
        };

        //Forward Declarations used in Matrix4 class.
        class Matrix3;
        class Matrix4;
        Vector3* Vector3MatrixMultiply(Vector3* result, const Vector3* inputVector, const Matrix3* inputMatrix);
        Vector3* MatrixVector3Multiply(Vector3* result, const Matrix3* inputMatrix, const Vector3* inputVector);
        Matrix3* MatrixMultiply(Matrix3* result, const Matrix3* left, const Matrix3* right);
        Matrix4* MatrixMultiply(Matrix4* result, const Matrix4* left, const Matrix4* right);
        Matrix3* MatrixMultiply(Matrix3* result, const Matrix3* input, const float scalar);
        Matrix4* MatrixMultiply(Matrix4* result, const Matrix4* input, const float scalar);

        Matrix3* MatrixComputeInverse(Matrix3* result, const Matrix3* input);
        Matrix4* MatrixComputeInverse(Matrix4* result, const Matrix4* input);
        Matrix3* MatrixFastInverse(Matrix3* result, const Matrix3* input);
        Matrix4* MatrixFastInverse(Matrix4* result, const Matrix4* input);
        Vector3* Vector3TransformCoordinate(Vector3* result, const Vector3* inputVector, const Matrix4* inputMatrix);
        Vector3* Vector3TransformNormal(Vector3* result, const Vector3* inputVector, const Matrix4* inputMatrix);
        Vector4* Vector4MatrixMultiply(Vector4* result, const Vector4* inputVector, const Matrix4* inputMatrix);
        Vector4* MatrixVector4Multiply(Vector4* result, const Matrix4* inputMatrix, const Vector4* inputVector);

        Matrix4* MatrixCreateTranslation(Matrix4* result, const Vector3* translation);
        Matrix4* MatrixCreateTranslation(Matrix4* result, const float translationX, const float translationY, const float translationZ);
        Matrix3* MatrixCreateScale(Matrix3* result, const Vector3* scale);
        Matrix4* MatrixCreateScale(Matrix4* result, const Vector3* scale);
        Matrix3* MatrixCreateScale(Matrix3* result, const float scaleX, const float scaleY, const float scaleZ);
        Matrix4* MatrixCreateScale(Matrix4* result, const float scaleX, const float scaleY, const float scaleZ);
        Matrix3* MatrixCreateRotationX(Matrix3* result, const Angle rotation);
        Matrix4* MatrixCreateRotationX(Matrix4* result, const Angle rotation);
        Matrix3* MatrixCreateRotationY(Matrix3* result, const Angle rotation);
        Matrix4* MatrixCreateRotationY(Matrix4* result, const Angle rotation);
        Matrix3* MatrixCreateRotationZ(Matrix3* result, const Angle rotation);
        Matrix4* MatrixCreateRotationZ(Matrix4* result, const Angle rotation);
        Matrix3* MatrixCreateRotationA(Matrix3* result, const Vector3* rotationAxis, const Angle rotation);
        Matrix4* MatrixCreateRotationA(Matrix4* result, const Vector3* rotationAxis, const Angle rotation);

        class Matrix3
        {
        public:
            static Matrix3 Zero(void)
            {
                return Matrix3(
                    0.0f, 0.0f, 0.0f,
                    0.0f, 0.0f, 0.0f,
                    0.0f, 0.0f, 0.0f);
            }

            static Matrix3 Identity(void)
            {
                return Matrix3(
                    1.0f, 0.0f, 0.0f,
                    0.0f, 1.0f, 0.0f,
                    0.0f, 0.0f, 1.0f);
            }

            static Matrix3 Scale(const tbMath::Vector3& scale)
            {
                Matrix3 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateScale(&matrix, &scale);
                return matrix;
            }

            static Matrix3 Scale(const float scaleX, const float scaleY, const float scaleZ)
            {
                Matrix3 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateScale(&matrix, scaleX, scaleY, scaleZ);
                return matrix;
            }

            static Matrix3 RotationX(const Angle& rotation)
            {
                Matrix3 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateRotationX(&matrix, rotation);
                return matrix;
            }

            static Matrix3 RotationY(const Angle& rotation)
            {
                Matrix3 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateRotationY(&matrix, rotation);
                return matrix;
            }

            static Matrix3 RotationZ(const Angle& rotation)
            {
                Matrix3 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateRotationZ(&matrix, rotation);
                return matrix;
            }

            static Matrix3 RotationA(const tbMath::Vector3& rotationAxis, const Angle& rotation)
            {
                Matrix3 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateRotationA(&matrix, &rotationAxis, rotation);
                return matrix;
            }

            //Defined in tb_matrix_quaternion.h
            static Matrix3 FromQuaternion(const tbMath::Unstable::Quaternion& quaternion);

            union
            {
                float mComponents[9];

#if defined(tb_visual_cpp)
#pragma warning(push)
#pragma warning(disable: 4201)
                struct {    //May be deprecating!?
                    float m_f11; float m_f21; float m_f31;
                    float m_f12; float m_f22; float m_f32;
                    float m_f13; float m_f23; float m_f33;
                };
#pragma warning(pop)
#else
                struct {    //May be deprecating!?
                    float m_f11; float m_f21; float m_f31;
                    float m_f12; float m_f22; float m_f32;
                    float m_f13; float m_f23; float m_f33;
                };
#endif
            };

            inline Matrix3(void) :
                mComponents{
                1.0f, 0.0f, 0.0f,
                0.0f, 1.0f, 0.0f,
                0.0f, 0.0f, 1.0f }
            {
            }

            inline explicit Matrix3(float f11, float f21, float f31,
                float f12, float f22, float f32,
                float f13, float f23, float f33) :
                mComponents{
                f11, f21, f31,
                f12, f22, f32,
                f13, f23, f33 }
            {
            }

            inline explicit Matrix3(const float * const componentArray) :
                mComponents{
                    componentArray[0], componentArray[1], componentArray[2],
                    componentArray[3], componentArray[4], componentArray[5],
                    componentArray[6], componentArray[7], componentArray[8] }
            {
            }

            inline explicit Matrix3(const SkipInitialization& fastAndStupid)
            {
                tb_unused(fastAndStupid);
            }

            inline Matrix3(const Matrix3& other) :
                mComponents{
                    other.mComponents[0], other.mComponents[1], other.mComponents[2],
                    other.mComponents[3], other.mComponents[4], other.mComponents[5],
                    other.mComponents[6], other.mComponents[7], other.mComponents[8]}
            {
            }

            inline Matrix3& operator=(const Matrix3& other)
            {
                memcpy(mComponents, other.mComponents, sizeof(float) * 9);
                return (*this);
            }

            inline bool operator==(const Matrix3& other) const
            {
                for (size_t index(0); index < 9; ++index)
                {
                    if (false == IsEqual(mComponents[index], other[index]))
                    {
                        return false;
                    }
                }

                return true;
            }

            inline bool operator!=(const Matrix3& other) const
            {
                for (size_t index(0); index < 9; ++index)
                {
                    if (false == IsEqual(mComponents[index], other[index]))
                    {
                        return true;
                    }
                }

                return false;
            }

            inline const float& operator[](const int& index) const { return mComponents[index]; }

            inline float& operator[](const int& index) { return mComponents[index]; }

            inline const float& operator()(const size_t& column, const size_t& row) const { return mComponents[column + (row * 3)]; }

            inline float& operator()(const size_t& column, const size_t& row) { return mComponents[column + (row * 3)]; }

            inline const float& operator()(int column, int row) const { return mComponents[column + (row * 3)]; }

            inline float& operator()(int column, int row) { return mComponents[column + (row * 3)]; }

            inline explicit operator float*(void) { return mComponents; }

            inline explicit operator const float*(void) const { return mComponents; }

            inline const Vector3 GetBasis(const size_t& basisIndex) const { return Vector3(&mComponents[basisIndex * 3]); }

            inline void SetBasis(const size_t& basisIndex, const Vector3 &basis) { *((Vector3*)&mComponents[basisIndex * 3]) = basis; }

            inline void SetBasis(const size_t& basisIndex, float basisX, float basisY, float basisZ)
            {
                mComponents[(basisIndex * 3) + 0] = basisX;
                mComponents[(basisIndex * 3) + 1] = basisY;
                mComponents[(basisIndex * 3) + 2] = basisZ;
            }

            inline Vector3 Transform(const Vector3& inputVector) const
            {
                Vector3 result;
                Vector3MatrixMultiply(&result, &inputVector, this);
                return result;
            }

#if defined(tb_with_math_operators)
            inline Matrix3 operator+(const Matrix3& rightSide) const
            {
                const Matrix3& leftSide(*this);
                return Matrix3(
                    leftSide(0, 0) + rightSide(0, 0), leftSide(1, 0) + rightSide(1, 0), leftSide(2, 0) + rightSide(2, 0),
                    leftSide(0, 1) + rightSide(0, 1), leftSide(1, 1) + rightSide(1, 1), leftSide(2, 1) + rightSide(2, 1),
                    leftSide(0, 2) + rightSide(0, 2), leftSide(1, 2) + rightSide(1, 2), leftSide(2, 2) + rightSide(2, 2));
            }

            inline Matrix3& operator+=(const Matrix3& rightSide)
            {
                for (size_t index = 0; index < 9; ++index)
                {
                    (*this)[index] += rightSide[index];
                }
                return *this;
            }

            inline Matrix3 operator-(const Matrix3& rightSide) const
            {
                const Matrix3& leftSide(*this);
                return Matrix3(
                    leftSide(0, 0) - rightSide(0, 0), leftSide(1, 0) - rightSide(1, 0), leftSide(2, 0) - rightSide(2, 0),
                    leftSide(0, 1) - rightSide(0, 1), leftSide(1, 1) - rightSide(1, 1), leftSide(2, 1) - rightSide(2, 1),
                    leftSide(0, 2) - rightSide(0, 2), leftSide(1, 2) - rightSide(1, 2), leftSide(2, 2) - rightSide(2, 2));
            }

            inline Matrix3& operator-=(const Matrix3& rightSide)
            {
                for (size_t index = 0; index < 9; ++index)
                {
                    (*this)[index] -= rightSide[index];
                }
                return *this;
            }


            inline Matrix3 operator*(const Matrix3& rightSide) const
            {
                Matrix3 result;
                MatrixMultiply(&result, this, &rightSide);
                return result;
            }

            inline Matrix3 operator*=(const Matrix3& rightSide)
            {
                Matrix3 leftSide = *this;
                MatrixMultiply(this, &leftSide, &rightSide);
                return *this;
            }

            inline Matrix3 GetTransposed(void) const
            {
                const Matrix3& self(*this);
                return Matrix3(
                    self(0, 0), self(0, 1), self(0, 2),
                    self(1, 0), self(1, 1), self(1, 2),
                    self(2, 0), self(2, 1), self(2, 2));
            }

            inline Matrix3 ComputeInverse(void) const
            {
                Matrix3 result;
                MatrixComputeInverse(&result, this);
                return result;
            }

            inline Matrix3 FastInverse(void) const
            {
                Matrix3 result;
                MatrixFastInverse(&result, this);
                return result;
            }
#endif
        };


        class Matrix4
        {
        public:
            static Matrix4 Zero(void)
            {
                return Matrix4(
                    0.0f, 0.0f, 0.0f, 0.0f,
                    0.0f, 0.0f, 0.0f, 0.0f,
                    0.0f, 0.0f, 0.0f, 0.0f,
                    0.0f, 0.0f, 0.0f, 0.0f);
            }

            static Matrix4 Identity(void)
            {
                return Matrix4(
                    1.0f, 0.0f, 0.0f, 0.0f,
                    0.0f, 1.0f, 0.0f, 0.0f,
                    0.0f, 0.0f, 1.0f, 0.0f,
                    0.0f, 0.0f, 0.0f, 1.0f);
            }

            static Matrix4 Translation(const tbMath::Vector3& translation)
            {
                Matrix4 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateTranslation(&matrix, &translation);
                return matrix;
            }

            static Matrix4 Translation(const float translationX, const float translationY, const float translationZ)
            {
                Matrix4 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateTranslation(&matrix, translationX, translationY, translationZ);
                return matrix;
            }

            static Matrix4 Scale(const tbMath::Vector3& scale)
            {
                Matrix4 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateScale(&matrix, &scale);
                return matrix;
            }

            static Matrix4 Scale(const float scaleX, const float scaleY, const float scaleZ)
            {
                Matrix4 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateScale(&matrix, scaleX, scaleY, scaleZ);
                return matrix;
            }

            static Matrix4 RotationX(const Angle& rotation)
            {
                Matrix4 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateRotationX(&matrix, rotation);
                return matrix;
            }

            static Matrix4 RotationY(const Angle& rotation)
            {
                Matrix4 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateRotationY(&matrix, rotation);
                return matrix;
            }

            static Matrix4 RotationZ(const Angle& rotation)
            {
                Matrix4 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateRotationZ(&matrix, rotation);
                return matrix;
            }

            static Matrix4 RotationA(const tbMath::Vector3& rotationAxis, const Angle& rotation)
            {
                Matrix4 matrix(tbMath::kSkipInitialization);
                tbMath::MatrixCreateRotationA(&matrix, &rotationAxis, rotation);
                return matrix;
            }

            //Defined in tb_matrix_quaternion.h
            static Matrix4 FromQuaternion(const tbMath::Unstable::Quaternion& quaternion, const tbMath::Vector3& translation);

            union
            {
                float mComponents[16];

#if defined(tb_visual_cpp)
#pragma warning(push)
#pragma warning(disable: 4201)
#endif
                struct {
                    float xx; float xy; float xz; float hx;
                    float yx; float yy; float yz; float hy;
                    float zx; float zy; float zz; float hz;
                    float tx; float ty; float tz; float ht;
                };
#if defined(tb_visual_cpp)
#pragma warning(pop)
#endif
            };

            inline Matrix4(void) :
                mComponents{
                    1.0f, 0.0f, 0.0f, 0.0f,
                    0.0f, 1.0f, 0.0f, 0.0f,
                    0.0f, 0.0f, 1.0f, 0.0f,
                    0.0f, 0.0f, 0.0f, 1.0f}
            {
            }

            inline explicit Matrix4(
                float f11, float f21, float f31, float f41,
                float f12, float f22, float f32, float f42,
                float f13, float f23, float f33, float f43,
                float f14, float f24, float f34, float f44) :
                mComponents{
                    f11, f21, f31, f41,
                    f12, f22, f32, f42,
                    f13, f23, f33, f43,
                    f14, f24, f34, f44}
            {
            }

            inline explicit Matrix4(
                float f11, float f21, float f31, float f41,
                float f12, float f22, float f32, float f42,
                float f13, float f23, float f33, float f43,
                float f14, float f24, float f34, float f44, const ColumnMajorMatrix& columnMajor) :
                mComponents{
                    f11, f12, f13, f14,
                    f21, f22, f23, f24,
                    f31, f32, f33, f34,
                    f41, f42, f43, f44 }
            {
                tb_unused(columnMajor);
            }

            inline explicit Matrix4(const float* const componentArray) :
                mComponents{
                    componentArray[0], componentArray[1], componentArray[2], componentArray[3],
                    componentArray[4], componentArray[5], componentArray[6], componentArray[7],
                    componentArray[8], componentArray[9], componentArray[10], componentArray[11],
                    componentArray[12], componentArray[13], componentArray[14], componentArray[15] }
            {
            }

            inline explicit Matrix4(const float* const componentArray, const ColumnMajorMatrix& columnMajor) :
                mComponents{
                    componentArray[0], componentArray[4], componentArray[8], componentArray[12],
                    componentArray[1], componentArray[5], componentArray[9], componentArray[13],
                    componentArray[2], componentArray[6], componentArray[10], componentArray[14],
                    componentArray[3], componentArray[7], componentArray[11], componentArray[15] }
            {
                tb_unused(columnMajor);
            }

            inline explicit Matrix4(const SkipInitialization& fastAndStupid)
            {
                tb_unused(fastAndStupid);
            }

            inline Matrix4(const Matrix4& other) :
                mComponents{
                    other.mComponents[0], other.mComponents[1], other.mComponents[2], other.mComponents[3],
                    other.mComponents[4], other.mComponents[5], other.mComponents[6], other.mComponents[7],
                    other.mComponents[8], other.mComponents[9], other.mComponents[10], other.mComponents[11],
                    other.mComponents[12], other.mComponents[13], other.mComponents[14], other.mComponents[15]
                }
            {
            }

            inline Matrix4& operator=(const Matrix4& other)
            {
                memcpy(mComponents, other.mComponents, sizeof(float) << 4);  //Same as (sizeof(float) * 16)
                return (*this);
            }

            inline bool operator==(const Matrix4& other) const
            {
                for (size_t index(0); index < 16; ++index)
                {
                    if (false == IsEqual(mComponents[index], other[index]))
                    {
                        return false;
                    }
                }

                return true;
            }

            inline bool operator!=(const Matrix4& other) const
            {
                for (size_t index(0); index < 16; ++index)
                {
                    if (false == IsEqual(mComponents[index], other[index]))
                    {
                        return true;
                    }
                }

                return false;
            }

            inline const float& operator[](const size_t& index) const { return mComponents[index]; }

            inline float& operator[](const size_t& index) { return mComponents[index]; }

            inline const float& operator()(const size_t& column, const size_t& row) const { return mComponents[column + (row << 2)]; }

            inline float& operator()(const size_t& column, const size_t& row) { return mComponents[column + (row << 2)]; }

            inline const float& operator()(int column, int row) const { return mComponents[column + (row << 2)]; }

            inline float& operator()(int column, int row) { return mComponents[column + (row << 2)]; }


            inline explicit operator float*(void) { return mComponents; }

            inline explicit operator const float*(void) const { return mComponents; }

            inline Vector3 GetBasis(const size_t& basisIndex) const { return Vector3(&mComponents[basisIndex<<2]); }

            inline void SetBasis(const size_t& basisIndex, const Vector3 &basis) { *((Vector3*)&mComponents[basisIndex<<2]) = basis; }

            inline void SetBasis(const size_t& basisIndex, float basisX, float basisY, float basisZ)
            {
                mComponents[(basisIndex<<2) + 0] = basisX;
                mComponents[(basisIndex<<2) + 1] = basisY;
                mComponents[(basisIndex<<2) + 2] = basisZ;
            }

            inline Vector3 GetPosition(void) const
            {
                return Vector3(&mComponents[12]);
            }

            inline void SetPosition(const Vector3& position)
            {
                (*this)(0, 3) = position.x;
                (*this)(1, 3) = position.y;
                (*this)(2, 3) = position.z;
            }

            inline void SetPosition(const float x, const float y, const float z)
            {
                (*this)(0, 3) = x;
                (*this)(1, 3) = y;
                (*this)(2, 3) = z;
            }

#if defined(tb_with_math_operators)
            inline Matrix4 operator+(const Matrix4& rightSide) const
            {
                const Matrix4& leftSide(*this);
                return Matrix4(
                    leftSide(0, 0) + rightSide(0, 0), leftSide(1, 0) + rightSide(1, 0), leftSide(2, 0) + rightSide(2, 0), leftSide(3, 0) + rightSide(3, 0),
                    leftSide(0, 1) + rightSide(0, 1), leftSide(1, 1) + rightSide(1, 1), leftSide(2, 1) + rightSide(2, 1), leftSide(3, 1) + rightSide(3, 1),
                    leftSide(0, 2) + rightSide(0, 2), leftSide(1, 2) + rightSide(1, 2), leftSide(2, 2) + rightSide(2, 2), leftSide(3, 2) + rightSide(3, 2),
                    leftSide(0, 3) + rightSide(0, 3), leftSide(1, 3) + rightSide(1, 3), leftSide(2, 3) + rightSide(2, 3), leftSide(3, 3) + rightSide(3, 3));
            }

            inline Matrix4& operator+=(const Matrix4& rightSide)
            {
                for (size_t index = 0; index < 16; ++index)
                {
                    (*this)[index] += rightSide[index];
                }
                return *this;
            }

            inline Matrix4 operator-(const Matrix4& rightSide) const
            {
                const Matrix4& leftSide(*this);
                return Matrix4(
                    leftSide(0, 0) - rightSide(0, 0), leftSide(1, 0) - rightSide(1, 0), leftSide(2, 0) - rightSide(2, 0), leftSide(3, 0) - rightSide(3, 0),
                    leftSide(0, 1) - rightSide(0, 1), leftSide(1, 1) - rightSide(1, 1), leftSide(2, 1) - rightSide(2, 1), leftSide(3, 1) - rightSide(3, 1),
                    leftSide(0, 2) - rightSide(0, 2), leftSide(1, 2) - rightSide(1, 2), leftSide(2, 2) - rightSide(2, 2), leftSide(3, 2) - rightSide(3, 2),
                    leftSide(0, 3) - rightSide(0, 3), leftSide(1, 3) - rightSide(1, 3), leftSide(2, 3) - rightSide(2, 3), leftSide(3, 3) - rightSide(3, 3));
            }

            inline Matrix4& operator-=(const Matrix4& rightSide)
            {
                for (size_t index = 0; index < 16; ++index)
                {
                    (*this)[index] -= rightSide[index];
                }
                return *this;
            }

            inline Matrix4 operator*(const Matrix4& rightSide) const
            {
                Matrix4 result;
                MatrixMultiply(&result, this, &rightSide);
                return result;
            }

            inline Matrix4 operator*=(const Matrix4& rightSide)
            {
                Matrix4 leftSide = *this;
                MatrixMultiply(this, &leftSide, &rightSide);
                return *this;
            }

            inline Matrix4 GetTransposed(void) const
            {
                Matrix4 temp;
                const Matrix4& self(*this);
                temp(0, 0) = self(0, 0);    temp(0, 1) = self(1, 0);    temp(0, 2) = self(2, 0);    temp(0, 3) = self(3, 0);
                temp(1, 0) = self(0, 1);    temp(1, 1) = self(1, 1);    temp(1, 2) = self(2, 1);    temp(1, 3) = self(3, 1);
                temp(2, 0) = self(0, 2);    temp(2, 1) = self(1, 2);    temp(2, 2) = self(2, 2);    temp(2, 3) = self(3, 2);
                temp(3, 0) = self(0, 3);    temp(3, 1) = self(1, 3);    temp(3, 2) = self(2, 3);    temp(3, 3) = self(3, 3);
                return temp;
            }

            inline Matrix4 ComputeInverse(void) const
            {
                Matrix4 result;
                MatrixComputeInverse(&result, this);
                return result;
            }

            inline Matrix4 FastInverse(void) const
            {
                Matrix4 result;
                MatrixFastInverse(&result, this);
                return result;
            }

            inline Vector3 TransformCoordinate(const Vector3& inputVector) const
            {
                Vector3 result;
                Vector3TransformCoordinate(&result, &inputVector, this);
                return result;
            }

            inline Vector3 TransformNormal(const Vector3& inputVector) const
            {
                Vector3 result;
                Vector3TransformNormal(&result, &inputVector, this);
                return result;
            }

            inline Vector4 Transform(const Vector4& inputVector) const
            {
                Vector4 result;
                Vector4MatrixMultiply(&result, &inputVector, this);
                return result;
            }

#endif /* tb_with_math_operators */
        };

#if defined(tb_with_math_operators)

        inline Matrix3 operator*(const float leftSide, const Matrix3& rightSide)
        {
            tbMath::Matrix3 result;
            MatrixMultiply(&result, &rightSide, leftSide);
            return result;
        }

        inline Matrix3 operator*(const Matrix3& leftSide, const float rightSide)
        {
            tbMath::Matrix3 result;
            MatrixMultiply(&result, &leftSide, rightSide);
            return result;
        }

        inline Matrix4 operator*(const float leftSide, const Matrix4& rightSide)
        {
            tbMath::Matrix4 result;
            MatrixMultiply(&result, &rightSide, leftSide);
            return result;
        }

        inline Matrix4 operator*(const Matrix4& leftSide, const float rightSide)
        {
            tbMath::Matrix4 result;
            MatrixMultiply(&result, &leftSide, rightSide);
            return result;
        }

        inline Vector3 operator*(const Vector3& leftSide, const Matrix3& rightSide)
        {
            tbMath::Vector3 result;
            Vector3MatrixMultiply(&result, &leftSide, &rightSide);
            return result;
        }

        inline Vector3 operator*(const Matrix3& leftSide, const Vector3& rightSide)
        {
            tbMath::Vector3 result;
            MatrixVector3Multiply(&result, &leftSide, &rightSide);
            return result;
        }

        inline Vector4 operator*(const Vector4& leftSide, const Matrix4& rightSide)
        {
            tbMath::Vector4 result;
            Vector4MatrixMultiply(&result, &leftSide, &rightSide);
            return result;
        }

        inline Vector4 operator*(const Matrix4& leftSide, const Vector4& rightSide)
        {
            tbMath::Vector4 result;
            MatrixVector4Multiply(&result, &leftSide, &rightSide);
            return result;
        }

#endif /* tb_with_math_operators */



        inline Matrix3* MatrixCreateIdentity(Matrix3* result)
        {
            *result = Matrix3::Identity();
            return result;
        }

        inline Matrix4* MatrixCreateIdentity(Matrix4* result)
        {
            *result = Matrix4::Identity();
            return result;
        }

        inline Matrix4* MatrixCreateLookAt(Matrix4* result, const Vector3& eyePosition, const Vector3& targetPosition, const Vector3& unitUp)
        {
            tb_error_if(nullptr == result, "tbExternalError: Invalid parameter; Expected result to be valid pointer.");

            Vector3 right, finalUp, forward;
            Vector3Subtract(&forward, &eyePosition, &targetPosition);
            Vector3Normalize(&forward, &forward);

#if defined(tb_debug_build) //Going to refer to these as the Math Idiot proofing tests.
            tb_error_if(eyePosition == targetPosition, "tbExternalError: Invalid parameter; This would result in looking directly at eye position.");
            tb_error_if(fabsf(Vector3Magnitude(&unitUp) - 1.0f) > 0.00001f, "tbExternalError: Invalid unitUp, this is not a unit length vector.");
            tb_error_if(fabsf(Vector3DotProduct(&forward, &unitUp)) > 0.9999f, "tbExternalError: Invalid looking direction, parallel to unitUp.");
#endif

            Vector3CrossProduct(&right, &unitUp, &forward);
            Vector3Normalize(&right, &right);

            Vector3CrossProduct(&finalUp, &forward, &right);

            //result->SetBasis(0, right);
            //result->SetBasis(1, finalUp);
            //result->SetBasis(2, forward);
            result->SetBasis(0, right.x, finalUp.x, forward.x);
            result->SetBasis(1, right.y, finalUp.y, forward.y);
            result->SetBasis(2, right.z, finalUp.z, forward.z);

            (*result)(0, 3) = -(eyePosition.x * (*result)(0, 0) + eyePosition.y * (*result)(0, 1) + eyePosition.z * (*result)(0, 2));
            (*result)(1, 3) = -(eyePosition.x * (*result)(1, 0) + eyePosition.y * (*result)(1, 1) + eyePosition.z * (*result)(1, 2));
            (*result)(2, 3) = -(eyePosition.x * (*result)(2, 0) + eyePosition.y * (*result)(2, 1) + eyePosition.z * (*result)(2, 2));

            return result;
        }

        inline Matrix4* MatrixCreateTranslation(Matrix4* result, const Vector3* translation)
        {   //Changed September 19, 2010, changed back March 22nd, 2018, changed again (to bottom row) May 9th 2018
            MatrixCreateIdentity(result);
            (*result)(0, 3) = translation->x;
            (*result)(1, 3) = translation->y;
            (*result)(2, 3) = translation->z;
            return result;
        }

        inline Matrix4* MatrixCreateTranslation(Matrix4* result, const float translationX, const float translationY, const float translationZ)
        {   //Changed September 19, 2010, changed back March 22nd, 2018, changed again (to bottom row) May 9th 2018
            MatrixCreateIdentity(result);
            (*result)(0, 3) = translationX;
            (*result)(1, 3) = translationY;
            (*result)(2, 3) = translationZ;
            return result;
        }

        inline Matrix3* MatrixCreateScale(Matrix3* result, const Vector3* scale)
        {
            memset(result->mComponents, 0, sizeof(float) * 9);
            (*result)(0, 0) = scale->x;
            (*result)(1, 1) = scale->y;
            (*result)(2, 2) = scale->z;
            return result;
        }

        inline Matrix4* MatrixCreateScale(Matrix4* result, const Vector3* scale)
        {
            memset(result->mComponents, 0, sizeof(float) << 4);
            (*result)(0, 0) = scale->x;
            (*result)(1, 1) = scale->y;
            (*result)(2, 2) = scale->z;
            (*result)(3, 3) = 1.0f;
            return result;
        }

        inline Matrix3* MatrixCreateScale(Matrix3* result, const float scaleX, const float scaleY, const float scaleZ)
        {
            memset(result->mComponents, 0, sizeof(float) * 9);
            (*result)(0, 0) = scaleX;
            (*result)(1, 1) = scaleY;
            (*result)(2, 2) = scaleZ;
            return result;
        }

        inline Matrix4* MatrixCreateScale(Matrix4* result, const float scaleX, const float scaleY, const float scaleZ)
        {
            memset(result->mComponents, 0, sizeof(float) << 4);
            (*result)(0, 0) = scaleX;
            (*result)(1, 1) = scaleY;
            (*result)(2, 2) = scaleZ;
            (*result)(3, 3) = 1.0f;
            return result;
        }

        inline Matrix3* MatrixCreateRotationX(Matrix3* result, const Angle rotation)
        {
            MatrixCreateIdentity(result);
            const float sinOfRotation = sinf(rotation.AsRadians());
            const float cosOfRotation = cosf(rotation.AsRadians());
            (*result)(1, 1) = cosOfRotation;
            (*result)(2, 1) = -sinOfRotation;
            (*result)(1, 2) = sinOfRotation;
            (*result)(2, 2) = cosOfRotation;
            return result;
        }

        inline Matrix4* MatrixCreateRotationX(Matrix4* result, const Angle rotation)
        {
            MatrixCreateIdentity(result);
            const float sinOfRotation = sinf(rotation.AsRadians());
            const float cosOfRotation = cosf(rotation.AsRadians());
            (*result)(1, 1) = cosOfRotation;
            (*result)(2, 1) = -sinOfRotation;
            (*result)(1, 2) = sinOfRotation;
            (*result)(2, 2) = cosOfRotation;
            return result;
        }

        inline Matrix3* MatrixCreateRotationY(Matrix3* result, const Angle rotation)
        {
            MatrixCreateIdentity(result);
            const float sinOfRotation = sinf(rotation.AsRadians());
            const float cosOfRotation = cosf(rotation.AsRadians());
            (*result)(0, 0) = cosOfRotation;
            (*result)(2, 0) = sinOfRotation;
            (*result)(0, 2) = -sinOfRotation;
            (*result)(2, 2) = cosOfRotation;
            return result;
        }

        inline Matrix4* MatrixCreateRotationY(Matrix4* result, const Angle rotation)
        {
            MatrixCreateIdentity(result);
            const float sinOfRotation = sinf(rotation.AsRadians());
            const float cosOfRotation = cosf(rotation.AsRadians());
            (*result)(0, 0) = cosOfRotation;
            (*result)(2, 0) = sinOfRotation;
            (*result)(0, 2) = -sinOfRotation;
            (*result)(2, 2) = cosOfRotation;
            return result;
        }

        inline Matrix3* MatrixCreateRotationZ(Matrix3* result, const Angle rotation)
        {
            MatrixCreateIdentity(result);
            const float sinOfRotation = sinf(rotation.AsRadians());
            const float cosOfRotation = cosf(rotation.AsRadians());
            (*result)(0, 0) = cosOfRotation;
            (*result)(1, 0) = -sinOfRotation;
            (*result)(0, 1) = sinOfRotation;
            (*result)(1, 1) = cosOfRotation;
            return result;
        }

        inline Matrix4* MatrixCreateRotationZ(Matrix4* result, const Angle rotation)
        {
            MatrixCreateIdentity(result);
            const float sinOfRotation = sinf(rotation.AsRadians());
            const float cosOfRotation = cosf(rotation.AsRadians());
            (*result)(0, 0) = cosOfRotation;
            (*result)(1, 0) = -sinOfRotation;
            (*result)(0, 1) = sinOfRotation;
            (*result)(1, 1) = cosOfRotation;
            return result;
        }

        inline Matrix3* MatrixCreateRotationA(Matrix3* result, const Vector3* rotationAxis, const Angle rotation)
        {
            const float s = sinf(rotation.AsRadians());
            const float c = cosf(rotation.AsRadians());
            const float t = 1.0f - c;

            const auto av = rotationAxis->GetNormalized();
            const auto tv = tbMath::Vector3(av * t);
            const auto sv = tbMath::Vector3(av * s);

            (*result)(0, 0) = c + tv.x * av.x;
            (*result)(0, 1) = tv.x * av.y + sv.z;
            (*result)(0, 2) = tv.x * av.z - sv.y;

            (*result)(1, 0) = tv.y * av.x - sv.z;
            (*result)(1, 1) = c + tv.y * av.y;
            (*result)(1, 2) = tv.y * av.z + sv.x;

            (*result)(2, 0) = tv.z * av.x + sv.y;
            (*result)(2, 1) = tv.z * av.y - sv.x;
            (*result)(2, 2) = c + tv.z * av.z;

            return result;
        }

        inline Matrix4* MatrixCreateRotationA(Matrix4* result, const Vector3* rotationAxis, const Angle rotation)
        {
            const float s = sinf(rotation.AsRadians());
            const float c = cosf(rotation.AsRadians());
            const float t = 1.0f - c;

            const auto av = rotationAxis->GetNormalized();
            const auto tv = tbMath::Vector3(av * t);
            const auto sv = tbMath::Vector3(av * s);

            (*result)(0, 0) = c + tv.x * av.x;
            (*result)(0, 1) = tv.x * av.y + sv.z;
            (*result)(0, 2) = tv.x * av.z - sv.y;
            (*result)(0, 3) = 0.0f;

            (*result)(1, 0) = tv.y * av.x - sv.z;
            (*result)(1, 1) = c + tv.y * av.y;
            (*result)(1, 2) = tv.y * av.z + sv.x;
            (*result)(1, 3) = 0.0f;

            (*result)(2, 0) = tv.z * av.x + sv.y;
            (*result)(2, 1) = tv.z * av.y - sv.x;
            (*result)(2, 2) = c + tv.z * av.z;
            (*result)(2, 3) = 0.0f;

            (*result)(3, 0) = 0.0f;
            (*result)(3, 1) = 0.0f;
            (*result)(3, 2) = 0.0f;
            (*result)(3, 3) = 1.0f;

            return result;
        }

        inline Matrix4* MatrixCreatePerspectiveRH(Matrix4* result, const float fieldOfView, const float aspectRatio,
            const float nearPlane, const float farPlane)
        {
            memset(result->mComponents, 0, sizeof(float) << 4);

            //From GLM, https://pastebin.com/Gen40yLY except flipped for row-major.
            //Remember that GLM accesses Y, X (row, column)...
            const float t = tan(fieldOfView / 2.0f);
            (*result)(0, 0) = 1.0f / (aspectRatio * t);
            (*result)(1, 1) = 1.0f / t;
            (*result)(2, 2) = -(farPlane + nearPlane) / (farPlane - nearPlane);
            (*result)(3, 2) = -1.0f;
            (*result)(2, 3) = -(2.0f * farPlane * nearPlane) / (farPlane - nearPlane);

            return result;
        }

        inline Matrix4* MatrixCreatePerspectiveLH(Matrix4* result, const float fieldOfView, const float aspectRatio,
            const float nearPlane, const float farPlane)
        {
            memset(result->mComponents, 0, sizeof(float) << 4);

            //From GLM, https://pastebin.com/Gen40yLY except flipped for row-major.
            //Remember that GLM accesses Y, X (row, column)...
            const float t = tan(fieldOfView / 2.0f);
            (*result)(0, 0) = 1.0f / (aspectRatio * t);
            (*result)(1, 1) = 1.0f / t;
            (*result)(2, 2) = (farPlane + nearPlane) / (farPlane - nearPlane);
            (*result)(3, 2) = 1.0f;
            (*result)(2, 3) = -(2.0f * farPlane * nearPlane) / (farPlane - nearPlane);

            return result;
        }

        inline Matrix4* MatrixCreateOrthoRH(Matrix4* result, const float left, const float right,
            const float top, const float bottom, const float nearPlane, const float farPlane)
        {
            tbMath::Matrix4& output(*result);
            memset(result->mComponents, 0, sizeof(float) << 4);

            output(0, 0) = 2.0f / (right - left);
            output(1, 1) = 2.0f / (top - bottom);
            output(2, 2) = -2.0f / (farPlane - nearPlane);
            output(3, 3) = 1.0f;
            output(0, 3) = -(right + left) / (right - left);
            output(1, 3) = -(top + bottom) / (top - bottom);
            output(2, 3) = -((farPlane + nearPlane) / (farPlane - nearPlane));

            return result;
        }

        inline Matrix4* MatrixCreateOrthoLH(Matrix4* result, const float left, const float right,
            const float top, const float bottom, const float nearPlane, const float farPlane)
        {
            tbMath::Matrix4& output(*result);
            memset(result->mComponents, 0, sizeof(float) << 4);

            output(0, 0) = 2.0f / (right - left);
            output(1, 1) = 2.0f / (top - bottom);
            output(2, 2) = 2.0f / (farPlane - nearPlane);
            output(3, 3) = 1.0f;
            output(0, 3) = -(right + left) / (right - left);
            output(1, 3) = -(top + bottom) / (top - bottom);
            output(2, 3) = -(farPlane + nearPlane) / (farPlane - nearPlane);

            return result;
        }

        inline Matrix3* MatrixAdd(Matrix3* result, const Matrix3* leftSide, const Matrix3* rightSide)
        {
            float *r = result->mComponents;
            const float *m1 = leftSide->mComponents;
            const float *m2 = rightSide->mComponents;

            r[0] = m1[0] + m2[0];       r[1] = m1[1] + m2[1];       r[2] = m1[2] + m2[2];
            r[3] = m1[3] + m2[3];       r[4] = m1[4] + m2[4];       r[5] = m1[5] + m2[5];
            r[6] = m1[6] + m2[6];       r[7] = m1[7] + m2[7];       r[8] = m1[8] + m2[8];
            return result;
        }

        inline Matrix4* MatrixAdd(Matrix4* result, const Matrix4* leftSide, const Matrix4* rightSide)
        {
            float *r = result->mComponents;
            const float *m1 = leftSide->mComponents;
            const float *m2 = rightSide->mComponents;

            r[0] = m1[0] + m2[0];       r[1] = m1[1] + m2[1];       r[2] = m1[2] + m2[2];       r[3] = m1[3] + m2[3];
            r[4] = m1[4] + m2[4];       r[5] = m1[5] + m2[5];       r[6] = m1[6] + m2[6];       r[7] = m1[7] + m2[7];
            r[8] = m1[8] + m2[8];       r[9] = m1[9] + m2[9];       r[10] = m1[10] + m2[10];    r[11] = m1[11] + m2[11];
            r[12] = m1[12] + m2[12];    r[13] = m1[13] + m2[13];    r[14] = m1[14] + m2[14];    r[15] = m1[15] + m2[15];
            return result;
        }


        inline Matrix3* MatrixSubtract(Matrix3* result, const Matrix3* leftSide, const Matrix3* rightSide)
        {
            float *r = result->mComponents;
            const float *m1 = leftSide->mComponents;
            const float *m2 = rightSide->mComponents;

            r[0] = m1[0] - m2[0];       r[1] = m1[1] - m2[1];        r[2] = m1[2] - m2[2];
            r[3] = m1[3] - m2[3];       r[4] = m1[4] - m2[4];        r[5] = m1[5] - m2[5];
            r[6] = m1[6] - m2[6];       r[7] = m1[7] - m2[7];        r[8] = m1[8] - m2[8];
            return result;
        }

        inline Matrix4* MatrixSubtract(Matrix4* result, const Matrix4* leftSide, const Matrix4* rightSide)
        {
            float *r = result->mComponents;
            const float *m1 = leftSide->mComponents;
            const float *m2 = rightSide->mComponents;

            r[0] = m1[0] - m2[0];       r[1] = m1[1] - m2[1];        r[2] = m1[2] - m2[2];       r[3] = m1[3] - m2[3];
            r[4] = m1[4] - m2[4];       r[5] = m1[5] - m2[5];        r[6] = m1[6] - m2[6];       r[7] = m1[7] - m2[7];
            r[8] = m1[8] - m2[8];       r[9] = m1[9] - m2[9];        r[10] = m1[10] - m2[10];    r[11] = m1[11] - m2[11];
            r[12] = m1[12] - m2[12];    r[13] = m1[13] - m2[13];     r[14] = m1[14] - m2[14];    r[15] = m1[15] - m2[15];
            return result;
        }

        inline Matrix3* MatrixTranspose(Matrix3* result, const Matrix3* input)
        {
            tb_error_if(result == input, "tbExternalError: Invalid parameter; expected result to be different than input");

            Matrix3& a(*result);
            const Matrix3& b(*input);
            a(0, 0) = b(0, 0);    a(0, 1) = b(1, 0);    a(0, 2) = b(2, 0);
            a(1, 0) = b(0, 1);    a(1, 1) = b(1, 1);    a(1, 2) = b(2, 1);
            a(2, 0) = b(0, 2);    a(2, 1) = b(1, 2);    a(2, 2) = b(2, 2);
            return result;
        }

        inline Matrix4* MatrixTranspose(Matrix4* result, const Matrix4* input)
        {
            tb_error_if(result == input, "tbExternalError: Invalid parameter; expected result to be different than input");

            Matrix4& a(*result);
            const Matrix4& b(*input);
            a(0, 0) = b(0, 0);    a(0, 1) = b(1, 0);    a(0, 2) = b(2, 0);    a(0, 3) = b(3, 0);
            a(1, 0) = b(0, 1);    a(1, 1) = b(1, 1);    a(1, 2) = b(2, 1);    a(1, 3) = b(3, 1);
            a(2, 0) = b(0, 2);    a(2, 1) = b(1, 2);    a(2, 2) = b(2, 2);    a(2, 3) = b(3, 2);
            a(3, 0) = b(0, 3);    a(3, 1) = b(1, 3);    a(3, 2) = b(2, 3);    a(3, 3) = b(3, 3);
            return result;
        }

        inline Matrix3* MatrixMultiply(Matrix3* result, const Matrix3* leftSide, const Matrix3* rightSide)
        {
            tb_error_if(result == leftSide || result == rightSide, "tbExternalError: Invalid parameter; expected result to be different than leftSide and rightSide");

            const tbMath::Matrix3& a = *leftSide;
            const tbMath::Matrix3& b = *rightSide;

            (*result)(0, 0) = a(0, 0) * b(0, 0) + a(1, 0) * b(0, 1) + a(2, 0) * b(0, 2);
            (*result)(1, 0) = a(0, 0) * b(1, 0) + a(1, 0) * b(1, 1) + a(2, 0) * b(1, 2);
            (*result)(2, 0) = a(0, 0) * b(2, 0) + a(1, 0) * b(2, 1) + a(2, 0) * b(2, 2);
            (*result)(0, 1) = a(0, 1) * b(0, 0) + a(1, 1) * b(0, 1) + a(2, 1) * b(0, 2);
            (*result)(1, 1) = a(0, 1) * b(1, 0) + a(1, 1) * b(1, 1) + a(2, 1) * b(1, 2);
            (*result)(2, 1) = a(0, 1) * b(2, 0) + a(1, 1) * b(2, 1) + a(2, 1) * b(2, 2);
            (*result)(0, 2) = a(0, 2) * b(0, 0) + a(1, 2) * b(0, 1) + a(2, 2) * b(0, 2);
            (*result)(1, 2) = a(0, 2) * b(1, 0) + a(1, 2) * b(1, 1) + a(2, 2) * b(1, 2);
            (*result)(2, 2) = a(0, 2) * b(2, 0) + a(1, 2) * b(2, 1) + a(2, 2) * b(2, 2);

            return result;
        }

        inline Matrix4* MatrixMultiply(Matrix4* result, const Matrix4* leftSide, const Matrix4* rightSide)
        {
            tb_error_if(result == leftSide || result == rightSide, "tbExternalError: Invalid parameter; expected result to be different than leftSide and rightSide");

            const tbMath::Matrix4& a = *leftSide;
            const tbMath::Matrix4& b = *rightSide;

            (*result)(0, 0) = a(0, 0) * b(0, 0) + a(1, 0) * b(0, 1) + a(2, 0) * b(0, 2) + a(3, 0) * b(0, 3);
            (*result)(1, 0) = a(0, 0) * b(1, 0) + a(1, 0) * b(1, 1) + a(2, 0) * b(1, 2) + a(3, 0) * b(1, 3);
            (*result)(2, 0) = a(0, 0) * b(2, 0) + a(1, 0) * b(2, 1) + a(2, 0) * b(2, 2) + a(3, 0) * b(2, 3);
            (*result)(3, 0) = a(0, 0) * b(3, 0) + a(1, 0) * b(3, 1) + a(2, 0) * b(3, 2) + a(3, 0) * b(3, 3);
            (*result)(0, 1) = a(0, 1) * b(0, 0) + a(1, 1) * b(0, 1) + a(2, 1) * b(0, 2) + a(3, 1) * b(0, 3);
            (*result)(1, 1) = a(0, 1) * b(1, 0) + a(1, 1) * b(1, 1) + a(2, 1) * b(1, 2) + a(3, 1) * b(1, 3);
            (*result)(2, 1) = a(0, 1) * b(2, 0) + a(1, 1) * b(2, 1) + a(2, 1) * b(2, 2) + a(3, 1) * b(2, 3);
            (*result)(3, 1) = a(0, 1) * b(3, 0) + a(1, 1) * b(3, 1) + a(2, 1) * b(3, 2) + a(3, 1) * b(3, 3);
            (*result)(0, 2) = a(0, 2) * b(0, 0) + a(1, 2) * b(0, 1) + a(2, 2) * b(0, 2) + a(3, 2) * b(0, 3);
            (*result)(1, 2) = a(0, 2) * b(1, 0) + a(1, 2) * b(1, 1) + a(2, 2) * b(1, 2) + a(3, 2) * b(1, 3);
            (*result)(2, 2) = a(0, 2) * b(2, 0) + a(1, 2) * b(2, 1) + a(2, 2) * b(2, 2) + a(3, 2) * b(2, 3);
            (*result)(3, 2) = a(0, 2) * b(3, 0) + a(1, 2) * b(3, 1) + a(2, 2) * b(3, 2) + a(3, 2) * b(3, 3);
            (*result)(0, 3) = a(0, 3) * b(0, 0) + a(1, 3) * b(0, 1) + a(2, 3) * b(0, 2) + a(3, 3) * b(0, 3);
            (*result)(1, 3) = a(0, 3) * b(1, 0) + a(1, 3) * b(1, 1) + a(2, 3) * b(1, 2) + a(3, 3) * b(1, 3);
            (*result)(2, 3) = a(0, 3) * b(2, 0) + a(1, 3) * b(2, 1) + a(2, 3) * b(2, 2) + a(3, 3) * b(2, 3);
            (*result)(3, 3) = a(0, 3) * b(3, 0) + a(1, 3) * b(3, 1) + a(2, 3) * b(3, 2) + a(3, 3) * b(3, 3);

            return result;
        }

        inline Matrix3* MatrixMultiply(Matrix3* result, const Matrix3* input, const float scalar)
        {
            float* resultComponents = result->mComponents;
            const float* inputComponents = input->mComponents;
            for (size_t index = 0; index < 9; ++index)
            {
                resultComponents[index] = inputComponents[index] * scalar;
            }
            return result;
        }

        inline Matrix4* MatrixMultiply(Matrix4* result, const Matrix4* input, const float scalar)
        {
            float* resultComponents = result->mComponents;
            const float* inputComponents = input->mComponents;
            for (size_t index = 0; index < 16; ++index)
            {
                resultComponents[index] = inputComponents[index] * scalar;
            }
            return result;
        }

        //inline Vector3* Vector3MatrixMultiply(Vector3* result, const Vector3* inputVector, const Matrix3* inputMatrix)
        //{
        //  tb_error_if(result == inputVector, "tbExternalError: Invalid parameter; expected result to be different than inputVector.");
        //  result->x = (inputVector->x * inputMatrix->m_f11) + (inputVector->y * inputMatrix->m_f12) + (inputVector->z * inputMatrix->m_f13);
        //  result->y = (inputVector->x * inputMatrix->m_f21) + (inputVector->y * inputMatrix->m_f22) + (inputVector->z * inputMatrix->m_f23);
        //  result->z = (inputVector->x * inputMatrix->m_f31) + (inputVector->y * inputMatrix->m_f32) + (inputVector->z * inputMatrix->m_f33);
        //  return result;
        //}




        inline Vector3* MatrixVector3Multiply(Vector3* result, const Matrix3* inputMatrix, const Vector3* inputVector)
        {
            tb_error_if(result == inputVector, "tbExternalError: Invalid parameter; expected result to be different than inputVector.");
            result->x = (inputVector->x * (*inputMatrix)(0, 0)) + (inputVector->y * (*inputMatrix)(1, 0)) + (inputVector->z * (*inputMatrix)(2, 0));
            result->y = (inputVector->x * (*inputMatrix)(0, 1)) + (inputVector->y * (*inputMatrix)(1, 1)) + (inputVector->z * (*inputMatrix)(2, 1));
            result->z = (inputVector->x * (*inputMatrix)(0, 2)) + (inputVector->y * (*inputMatrix)(1, 2)) + (inputVector->z * (*inputMatrix)(2, 2));
            return result;
        }

        inline Vector3* Vector3MatrixMultiply(Vector3* result, const Vector3* inputVector, const Matrix3* inputMatrix)
        {
            tb_error_if(result == inputVector, "tbExternalError: Invalid parameter; expected result to be different than inputVector.");
            result->x = (inputVector->x * (*inputMatrix)(0, 0)) + (inputVector->y * (*inputMatrix)(0, 1)) + (inputVector->z * (*inputMatrix)(0, 2));
            result->y = (inputVector->x * (*inputMatrix)(1, 0)) + (inputVector->y * (*inputMatrix)(1, 1)) + (inputVector->z * (*inputMatrix)(1, 2));
            result->z = (inputVector->x * (*inputMatrix)(2, 0)) + (inputVector->y * (*inputMatrix)(2, 1)) + (inputVector->z * (*inputMatrix)(2, 2));
            return result;
        }

        inline Vector4* MatrixVector4Multiply(Vector4* result, const Matrix4* inputMatrix, const Vector4* inputVector)
        {
            tb_error_if(result == inputVector, "tbExternalError: Invalid parameter; expected result to be different than inputVector.");
            result->x = (inputVector->x * (*inputMatrix)(0, 0)) + (inputVector->y * (*inputMatrix)(1, 0)) + (inputVector->z * (*inputMatrix)(2, 0)) + (inputVector->w * (*inputMatrix)(3, 0));
            result->y = (inputVector->x * (*inputMatrix)(0, 1)) + (inputVector->y * (*inputMatrix)(1, 1)) + (inputVector->z * (*inputMatrix)(2, 1)) + (inputVector->w * (*inputMatrix)(3, 1));
            result->z = (inputVector->x * (*inputMatrix)(0, 2)) + (inputVector->y * (*inputMatrix)(1, 2)) + (inputVector->z * (*inputMatrix)(2, 2)) + (inputVector->w * (*inputMatrix)(3, 2));
            result->w = (inputVector->x * (*inputMatrix)(0, 3)) + (inputVector->y * (*inputMatrix)(1, 3)) + (inputVector->z * (*inputMatrix)(2, 3)) + (inputVector->w * (*inputMatrix)(3, 3));
            return result;
        }

        inline Vector4* Vector4MatrixMultiply(Vector4* result, const Vector4* inputVector, const Matrix4* inputMatrix)
        {
            tb_error_if(result == inputVector , "tbExternalError: Invalid parameter; expected result to be different than inputVector.");
            result->x = (inputVector->x * (*inputMatrix)(0, 0)) + (inputVector->y * (*inputMatrix)(0, 1)) + (inputVector->z * (*inputMatrix)(0, 2)) + (inputVector->w * (*inputMatrix)(0, 3));
            result->y = (inputVector->x * (*inputMatrix)(1, 0)) + (inputVector->y * (*inputMatrix)(1, 1)) + (inputVector->z * (*inputMatrix)(1, 2)) + (inputVector->w * (*inputMatrix)(1, 3));
            result->z = (inputVector->x * (*inputMatrix)(2, 0)) + (inputVector->y * (*inputMatrix)(2, 1)) + (inputVector->z * (*inputMatrix)(2, 2)) + (inputVector->w * (*inputMatrix)(2, 3));
            result->w = (inputVector->x * (*inputMatrix)(3, 0)) + (inputVector->y * (*inputMatrix)(3, 1)) + (inputVector->z * (*inputMatrix)(3, 2)) + (inputVector->w * (*inputMatrix)(3, 3));
            return result;
        }

        inline Vector3* Vector3TransformCoordinate(Vector3* result, const Vector3* inputVector, const Matrix4* inputMatrix)
        {
            tb_error_if(result == inputVector, "tbExternalError: Invalid parameter; expected result to be different than inputVector");

            const Matrix4& matrix(*inputMatrix);
            result->x = (inputVector->x * matrix(0, 0)) + (inputVector->y * matrix(0, 1)) + (inputVector->z * matrix(0, 2)) + matrix(0, 3);
            result->y = (inputVector->x * matrix(1, 0)) + (inputVector->y * matrix(1, 1)) + (inputVector->z * matrix(1, 2)) + matrix(1, 3);
            result->z = (inputVector->x * matrix(2, 0)) + (inputVector->y * matrix(2, 1)) + (inputVector->z * matrix(2, 2)) + matrix(2, 3);
            return result;
        }

        inline Vector3* Vector3TransformNormal(Vector3* result, const Vector3* inputVector, const Matrix4* inputMatrix)
        {
            tb_error_if(result == inputVector, "tbExternalError: Invalid parameter; expected result to be different than inputVector");

            const Matrix4& matrix(*inputMatrix);
            result->x = (inputVector->x * matrix(0, 0)) + (inputVector->y * matrix(0, 1)) + (inputVector->z * matrix(0, 2));
            result->y = (inputVector->x * matrix(1, 0)) + (inputVector->y * matrix(1, 1)) + (inputVector->z * matrix(1, 2));
            result->z = (inputVector->x * matrix(2, 0)) + (inputVector->y * matrix(2, 1)) + (inputVector->z * matrix(2, 2));
            return result;
        }

        inline float Matrix2x2Determinant(const float f11, const float f12, const float f21, const float f22)
        {
            return f11 * f22 - f12 * f21;
        }

        inline float Matrix3x3Determinant(const float f11, const float f12, const float f13, const float f21, const float f22, const float f23, const float f31, const float f32, const float f33)
        {
            return f11 * (f22 * f33 - f32 * f23) - f12 * (f21 * f33 - f31 * f23) + f13 * (f21 * f32 - f31 * f22);
        }

        inline float MatrixDeterminant(const Matrix3* input)
        {
            const Matrix3& m(*input);
            return m(0, 0) * (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) - m(0, 1) * (m(1, 0) * m(2, 2) - m(2, 0) * m(1, 2)) + m(0, 2) * (m(1, 0) * m(2, 1) - m(2, 0) * m(1, 1));
        }

        inline float MatrixDeterminant(const Matrix4* input)
        {
            const Matrix4& m(*input);
            return  (m(0, 0) * Matrix3x3Determinant(m(1, 1), m(2, 1), m(3, 1), m(1, 2), m(2, 2), m(3, 2), m(1, 3), m(2, 3), m(3, 3))) -
                    (m(1, 0) * Matrix3x3Determinant(m(0, 1), m(2, 1), m(3, 1), m(0, 2), m(2, 2), m(3, 2), m(0, 3), m(2, 3), m(3, 3))) +
                    (m(2, 0) * Matrix3x3Determinant(m(0, 1), m(1, 1), m(3, 1), m(0, 2), m(1, 2), m(3, 2), m(0, 3), m(1, 3), m(3, 3))) -
                    (m(3, 0) * Matrix3x3Determinant(m(0, 1), m(1, 1), m(2, 1), m(0, 2), m(1, 2), m(2, 2), m(0, 3), m(1, 3), m(2, 3)));
        }

        inline Matrix3* MatrixComputeInverse(Matrix3* result, const Matrix3* input)
        {
            //tb_error("This is not yet implemented in full, use the FastInverse which requires certain ... matrices.");

            tb_error_if(result == input, "tbExternalError: Invalid parameter, expected result to be different than input.");
            const float determinant = MatrixDeterminant(input);
            if (true == IsZero(determinant))
            {
                *result = *input;
                return result;
            }

            const float inverseDeterminant = 1.0f / determinant;

            const Matrix3& m(*input);
            (*result)(0, 0) = Matrix2x2Determinant(m(1, 1), m(1, 2), m(2, 1), m(2, 2)) * inverseDeterminant;
            (*result)(1, 0) = -Matrix2x2Determinant(m(1, 0), m(1, 2), m(2, 0), m(2, 2)) * inverseDeterminant;
            (*result)(2, 0) = Matrix2x2Determinant(m(1, 0), m(1, 1), m(2, 0), m(2, 1)) * inverseDeterminant;
            (*result)(0, 1) = -Matrix2x2Determinant(m(0, 1), m(0, 2), m(2, 1), m(2, 2)) * inverseDeterminant;
            (*result)(1, 1) = Matrix2x2Determinant(m(0, 0), m(0, 2), m(2, 0), m(2, 2)) * inverseDeterminant;
            (*result)(2, 1) = -Matrix2x2Determinant(m(0, 0), m(0, 1), m(2, 0), m(2, 1)) * inverseDeterminant;
            (*result)(0, 2) = Matrix2x2Determinant(m(0, 1), m(0, 2), m(1, 1), m(1, 2)) * inverseDeterminant;
            (*result)(1, 2) = -Matrix2x2Determinant(m(0, 0), m(0, 2), m(1, 0), m(1, 2)) * inverseDeterminant;
            (*result)(2, 2) = Matrix2x2Determinant(m(0, 0), m(0, 1), m(1, 0), m(1, 1)) * inverseDeterminant;
            return result;
        }

        inline Matrix4* MatrixComputeInverse(Matrix4* result, const Matrix4* input)
        {
            tb_error_if(result == input, "tbExternalError: Invalid parameter, expected result to be different than input.");
            const float determinant = MatrixDeterminant(input);
            if (true == IsZero(determinant))
            {
                *result = *input;
                return result;
            }

            const float inverseDeterminant = 1.0f / determinant;
            const Matrix4& m(*input);
            (*result)(0, 0) = Matrix3x3Determinant (m(1,1), m(2,1), m(3,1), m(1,2), m(2,2), m(3,2), m(1,3), m(2,3), m(3,3)) * inverseDeterminant;
            (*result)(1, 0) = -Matrix3x3Determinant(m(1,0), m(2,0), m(3,0), m(1,2), m(2,2), m(3,2), m(1,3), m(2,3), m(3,3)) * inverseDeterminant;
            (*result)(2, 0) = Matrix3x3Determinant (m(1,0), m(2,0), m(3,0), m(1,1), m(2,1), m(3,1), m(1,3), m(2,3), m(3,3)) * inverseDeterminant;
            (*result)(3, 0) = -Matrix3x3Determinant(m(1,0), m(2,0), m(3,0), m(1,1), m(2,1), m(3,1), m(1,2), m(2,2), m(3,2)) * inverseDeterminant;
            (*result)(0, 1) = -Matrix3x3Determinant(m(0,1), m(2,1), m(3,1), m(0,2), m(2,2), m(3,2), m(0,3), m(2,3), m(3,3)) * inverseDeterminant;
            (*result)(1, 1) = Matrix3x3Determinant (m(0,0), m(2,0), m(3,0), m(0,2), m(2,2), m(3,2), m(0,3), m(2,3), m(3,3)) * inverseDeterminant;
            (*result)(2, 1) = -Matrix3x3Determinant(m(0,0), m(2,0), m(3,0), m(0,1), m(2,1), m(3,1), m(0,3), m(2,3), m(3,3)) * inverseDeterminant;
            (*result)(3, 1) = Matrix3x3Determinant (m(0,0), m(2,0), m(3,0), m(0,1), m(2,1), m(3,1), m(0,2), m(2,2), m(3,2)) * inverseDeterminant;
            (*result)(0, 2) = Matrix3x3Determinant (m(0,1), m(1,1), m(3,1), m(0,2), m(1,2), m(3,2), m(0,3), m(1,3), m(3,3)) * inverseDeterminant;
            (*result)(1, 2) = -Matrix3x3Determinant(m(0,0), m(1,0), m(3,0), m(0,2), m(1,2), m(3,2), m(0,3), m(1,3), m(3,3)) * inverseDeterminant;
            (*result)(2, 2) = Matrix3x3Determinant (m(0,0), m(1,0), m(3,0), m(0,1), m(1,1), m(3,1), m(0,3), m(1,3), m(3,3)) * inverseDeterminant;
            (*result)(3, 2) = -Matrix3x3Determinant(m(0,0), m(1,0), m(3,0), m(0,1), m(1,1), m(3,1), m(0,2), m(1,2), m(3,2)) * inverseDeterminant;
            (*result)(0, 3) = -Matrix3x3Determinant(m(0,1), m(1,1), m(2,1), m(0,2), m(1,2), m(2,2), m(0,3), m(1,3), m(2,3)) * inverseDeterminant;
            (*result)(1, 3) = Matrix3x3Determinant (m(0,0), m(1,0), m(2,0), m(0,2), m(1,2), m(2,2), m(0,3), m(1,3), m(2,3)) * inverseDeterminant;
            (*result)(2, 3) = -Matrix3x3Determinant(m(0,0), m(1,0), m(2,0), m(0,1), m(1,1), m(2,1), m(0,3), m(1,3), m(2,3)) * inverseDeterminant;
            (*result)(3, 3) = Matrix3x3Determinant (m(0,0), m(1,0), m(2,0), m(0,1), m(1,1), m(2,1), m(0,2), m(1,2), m(2,2)) * inverseDeterminant;
            return result;
        }

        inline Matrix3* MatrixFastInverse(Matrix3* result, const Matrix3* input)
        {
            tb_error_if(result == input, "tbExternalError: Invalid parameter, expected result to be different than input.");

            (*result)(0, 0) = (*input)(0, 0);   //Diagonal is equal
            (*result)(1, 1) = (*input)(1, 1);
            (*result)(2, 2) = (*input)(2, 2);

            (*result)(1, 0) = (*input)(0, 1);   //Rotation is transposed
            (*result)(2, 0) = (*input)(0, 2);
            (*result)(0, 1) = (*input)(1, 0);
            (*result)(2, 1) = (*input)(1, 2);
            (*result)(0, 2) = (*input)(2, 0);
            (*result)(1, 2) = (*input)(2, 1);

            return result;
        }

        inline Matrix4* MatrixFastInverse(Matrix4* result, const Matrix4* input)
        {
            tb_error_if(result == input, "tbExternalError: Invalid parameter, expected result to be different than input.");

            (*result)(0, 0) = (*input)(0, 0);   //Diagonal is equal
            (*result)(1, 1) = (*input)(1, 1);
            (*result)(2, 2) = (*input)(2, 2);
            (*result)(3, 3) = (*input)(3, 3);

            (*result)(1, 0) = (*input)(0, 1);   //Rotation is transposed
            (*result)(2, 0) = (*input)(0, 2);
            (*result)(0, 1) = (*input)(1, 0);
            (*result)(2, 1) = (*input)(1, 2);
            (*result)(0, 2) = (*input)(2, 0);
            (*result)(1, 2) = (*input)(2, 1);

            //This was the 'old' way that worked with translation in final COLUMN; it was a discrepency with where it
            //should have been; the test data changed 2020-11-26 and this will now succeed in the new way below. See
            //the comment at the very top of file for more information.
            //
            //const tbMath::Matrix4& inversed(*result); //Translation is multiplied by inversed 3x3, then negated.
            //(*result)(3, 0) = -((*input)(3, 0) * inversed(0, 0) + (*input)(3, 1) * inversed(1, 0) + (*input)(3, 2) * inversed(2, 0));
            //(*result)(3, 1) = -((*input)(3, 0) * inversed(0, 1) + (*input)(3, 1) * inversed(1, 1) + (*input)(3, 2) * inversed(2, 1));
            //(*result)(3, 2) = -((*input)(3, 0) * inversed(0, 2) + (*input)(3, 1) * inversed(1, 2) + (*input)(3, 2) * inversed(2, 2));

            //(*result)(0, 3) = 0.0f;
            //(*result)(1, 3) = 0.0f;
            //(*result)(2, 3) = 0.0f;

            const tbMath::Matrix4& inversed(*result);   //Translation is multiplied by inversed 3x3, then negated.
            (*result)(0, 3) = -((*input)(0, 3) * inversed(0, 0) + (*input)(1, 3) * inversed(0, 1) + (*input)(2, 3) * inversed(0, 2));
            (*result)(1, 3) = -((*input)(0, 3) * inversed(1, 0) + (*input)(1, 3) * inversed(1, 1) + (*input)(2, 3) * inversed(1, 2));
            (*result)(2, 3) = -((*input)(0, 3) * inversed(2, 0) + (*input)(1, 3) * inversed(2, 1) + (*input)(2, 3) * inversed(2, 2));

            (*result)(3, 0) = 0.0f;
            (*result)(3, 1) = 0.0f;
            (*result)(3, 2) = 0.0f;

            return result;
        }

        namespace Unstable
        {
            inline Matrix4 Translate(const Matrix4& originalMatrix, const Vector3& translation)
            {
                Matrix4 outputMatrix(kSkipInitialization);
                Matrix4 translationMatrix(kSkipInitialization);

                MatrixCreateTranslation(&translationMatrix, &translation);
                MatrixMultiply(&outputMatrix, &translationMatrix, &originalMatrix);

                return outputMatrix;
            }

            inline Matrix4 Translate(const Matrix4& originalMatrix, float x, float y, float z)
            {
                return Translate(originalMatrix, tbMath::Vector3(x, y, z));
            }

            inline Matrix4 TranslateLocal(const Matrix4& originalMatrix, const Vector3& translation)
            {
                Matrix4 outputMatrix(kSkipInitialization);
                Matrix4 translationMatrix(kSkipInitialization);

                MatrixCreateTranslation(&translationMatrix, &translation);
                MatrixMultiply(&outputMatrix, &originalMatrix, &translationMatrix);

                return outputMatrix;
            }

            inline Matrix4 TranslateLocal(const Matrix4& originalMatrix, float x, float y, float z)
            {
                return TranslateLocal(originalMatrix, tbMath::Vector3(x, y, z));
            }

            inline Matrix3 Rotate(const Matrix3& originalMatrix, const Vector3& aroundAxis, const tbMath::Angle angle)
            {
                Matrix3 outputMatrix(kSkipInitialization);
                Matrix3 rotationMatrix(kSkipInitialization);

                //TODO: TurtleBrains: MathKit: Consistency check, do we negate here. CW / CCW seems to depend on the
                //handedness of the world rather than the negation here.
                MatrixCreateRotationA(&rotationMatrix, &aroundAxis, angle);
                MatrixMultiply(&outputMatrix, &rotationMatrix, &originalMatrix);

                return outputMatrix;
            }

            inline Matrix3 Rotate(const size_t& axisIndex, const Matrix3& originalMatrix, const tbMath::Angle angle)
            {
                return Rotate(originalMatrix, originalMatrix.GetBasis(axisIndex), angle);
            }

            inline Matrix4 Rotate(const Matrix4& originalMatrix, const Vector3& aroundAxis, const tbMath::Angle angle)
            {
                Matrix4 outputMatrix(kSkipInitialization);
                Matrix4 rotationMatrix(kSkipInitialization);

                //TODO: TurtleBrains: MathKit: Consistency check, do we negate here. CW / CCW seems to depend on the
                //handedness of the world rather than the negation here.
                MatrixCreateRotationA(&rotationMatrix, &aroundAxis, angle);
                MatrixMultiply(&outputMatrix, &rotationMatrix, &originalMatrix);

                return outputMatrix;
            }

            inline Matrix4 Rotate(const size_t& axisIndex, const Matrix4& originalMatrix, const tbMath::Angle angle)
            {
                return Rotate(originalMatrix, originalMatrix.GetBasis(axisIndex), angle);
            }

            //TODO: TurtleBrains: MathKit: Does RotateWorldSpace make sense, it is a little odd because of: Translate vs TranslateLocal
            //   and here is using Rotate vs RotateWorld. Consistency concerns.
            inline Matrix4 RotateWorldSpace(const size_t& axisIndex, const Matrix4& originalMatrix, const tbMath::Angle angle)
            {
                Matrix4 outputMatrix(kSkipInitialization);
                Matrix4 rotationMatrix(kSkipInitialization);

                //Negating here for positive rotationInDegrees to spin clockwise around the axis.
                //TODO: TurtleBrains: MathKit: Consistency check, do we negate here. CW / CCW seems to depend on the
                //handedness of the world rather than the negation here.
                const tbMath::Vector3 worldAxis(Matrix4::Identity().GetBasis(axisIndex));
                MatrixCreateRotationA(&rotationMatrix, &worldAxis, -angle);
                const Matrix4 rotationTransposed(rotationMatrix.GetTransposed());
                MatrixMultiply(&outputMatrix, &originalMatrix, &rotationTransposed);

                return outputMatrix;
            }

            inline Matrix3 Scale(const Matrix3& originalMatrix, const Vector3& scale)
            {
                Matrix3 outputMatrix(kSkipInitialization);
                Matrix3 scalingMatrix(kSkipInitialization);

                MatrixCreateScale(&scalingMatrix, &scale);
                MatrixMultiply(&outputMatrix, &scalingMatrix, &originalMatrix);

                return outputMatrix;
            }

            inline Matrix3 Scale(const Matrix3 originalMatrix, float x, float y, float z)
            {
                return Scale(originalMatrix, Vector3(x, y, z));
            }

            inline Matrix4 Scale(const Matrix4& originalMatrix, const Vector3& scale)
            {
                Matrix4 outputMatrix(kSkipInitialization);
                Matrix4 scalingMatrix(kSkipInitialization);

                MatrixCreateScale(&scalingMatrix, &scale);
                MatrixMultiply(&outputMatrix, &scalingMatrix, &originalMatrix);

                return outputMatrix;
            }

            inline Matrix4 Scale(const Matrix4& originalMatrix, float x, float y, float z)
            {
                return Scale(originalMatrix, Vector3(x, y, z));
            }

        }; /* namespace Unstable */

    }; /* namespace Math */
}; /* namespace TurtleBrains */

namespace tbMath = TurtleBrains::Math;

#endif /* TurtleBrains_Matrix_hpp */