feat: sync with Whoa implementation

tempest/CMakeLists.txt

@@ -1,7 +1,9 @@
 file(GLOB TEMPEST_SOURCES
     "*.cpp"
+    "math/*.cpp"
     "matrix/*.cpp"
     "quaternion/*.cpp"
+    "random/*.cpp"
     "rect/*.cpp"
     "vector/*.cpp"
 )

tempest/Random.hpp

@@ -0,0 +1,6 @@
+#ifndef TEMPEST_RANDOM_HPP
+#define TEMPEST_RANDOM_HPP
+
+#include "tempest/random/CRandom.hpp"
+
+#endif

tempest/Vector.hpp

@@ -4,6 +4,7 @@
 #include "tempest/vector/C2Vector.hpp"
 #include "tempest/vector/C2iVector.hpp"
 #include "tempest/vector/C3Vector.hpp"
+#include "tempest/vector/C3iVector.hpp"
 #include "tempest/vector/C4Vector.hpp"
 #include "tempest/vector/CImVector.hpp"
 

tempest/math/CMath.cpp

@@ -0,0 +1,7 @@
+#include "tempest/math/CMath.hpp"
+
+const float CMath::PI = 3.1415927f;
+const float CMath::TWO_PI = 6.2831855f;
+const float CMath::OO_TWO_PI = 1.0f / CMath::TWO_PI;
+const float CMath::EPSILON = 0.00000023841858f;
+const float CMath::DEG2RAD = CMath::PI / 180.0f;

tempest/math/CMath.hpp

@@ -8,9 +8,11 @@
 class CMath {
     public:
     // Static variables
-    static constexpr float PI = 3.1415927f;
+    static const float PI;
-    static constexpr float TWO_PI = 6.2831855f;
+    static const float TWO_PI;
-    static constexpr float OO_TWO_PI = 1.0f / TWO_PI;
+    static const float OO_TWO_PI;
+    static const float EPSILON;
+    static const float DEG2RAD;
 
     // Static functions
     static float acos(float x) {
@@ -18,25 +20,21 @@ class CMath {
     }
 
     static float cos(float x) {
-        return ::cosf(x);
+        return std::cos(x);
-    }
-
-    static float sin(float x) {
-        return ::sinf(x);
     }
 
     static float fabs(float x) {
-        return ::fabs(x);
+        return std::fabs(x);
+    }
+
+    static bool fequal(float a, float b) {
+        return CMath::fequalz(a, b, EPSILON);
     }
 
     static bool fequalz(float a, float b, float z) {
         return z > CMath::fabs(a - b);
     }
 
-    static bool fequal(float a, float b) {
-        return CMath::fequalz(a, b, 0.00000023841858f);
-    }
-
     static int32_t fint(float n) {
         return static_cast<int32_t>(n);
     }
@@ -45,6 +43,14 @@ class CMath {
         return n <= 0.0f ? static_cast<int32_t>(n - 0.5f) : static_cast<int32_t>(n + 0.5f);
     }
 
+    static bool fnotequal(float a, float b) {
+        return !CMath::fequal(a, b);
+    }
+
+    static bool fnotequalz(float a, float b, float z) {
+        return !CMath::fequalz(a, b, z);
+    }
+
     static uint32_t fuint(float n) {
         return static_cast<uint32_t>(n);
     }
@@ -57,26 +63,67 @@ class CMath {
         return static_cast<uint32_t>(n + 0.99994999);
     }
 
+    static float hypotinv(float x, float y) {
+        float s = (x * x) + (y * y);
+        STORM_ASSERT(s >= 0.0f);
+
+        return CMath::sqrtinv(s);
+    }
+
+    static float hypotinv(float x, float y, float z) {
+        float s = (x * x) + (y * y) + (z * z);
+        STORM_ASSERT(s >= 0.0f);
+
+        return CMath::sqrtinv(s);
+    }
+
+    static uint32_t mulhwu(uint32_t x, uint32_t y) {
+        return (y * static_cast<uint64_t>(x)) >> 32;
+    }
+
+    static void normalize(float& x, float& y) {
+        auto hi = CMath::hypotinv(x, y);
+
+        x *= hi;
+        y *= hi;
+    }
+
+    static void normalize(float& x, float& y, float& z) {
+        auto hi = CMath::hypotinv(x, y, z);
+
+        x *= hi;
+        y *= hi;
+        z *= hi;
+    }
+
+    static float normalizeangle0to2pi(float angle) {
+        angle = fmodf(angle, TWO_PI);
+        return angle < 0.0f ? angle + TWO_PI : angle;
+    }
+
+    static uint32_t rotl3(uint32_t v) {
+        return (v << 3) | (v >> 29);
+    }
+
+    static uint32_t rotl2(uint32_t v) {
+        return (v << 2) | (v >> 30);
+    }
+
+    static uint32_t rotl1(uint32_t v) {
+        return (v << 1) | (v >> 31);
+    }
+
+    static float sin(float x) {
+        return std::sin(x);
+    }
+
     static float sqrt(float x) {
         STORM_ASSERT(x >= 0.0f);
         return ::sqrt(x);
     }
 
-    static void normalize(float& x, float& y) {
+    static float sqrtinv(float x) {
-        float m = x * x + y * y;
+        return 1.0f / ::sqrt(x);
-        STORM_ASSERT(m >= 0.0f);
-        m = 1.0f / CMath::sqrt(m);
-        x *= m;
-        y *= m;
-    }
-
-    static void normalize(float& x, float& y, float& z) {
-        float m = x * x + y * y + z * z;
-        STORM_ASSERT(m >= 0.0f);
-        m = 1.0f / CMath::sqrt(m);
-        x *= m;
-        y *= m;
-        z *= m;
     }
 };
 

tempest/matrix/C33Matrix.cpp

@@ -13,6 +13,45 @@ float C33Matrix::Det(float a, float b, float c, float d) {
     return (a * d) - (b * c);
 }
 
+C33Matrix C33Matrix::Rotation(float angle) {
+    return C33Matrix::RotationAroundZ(angle);
+}
+
+C33Matrix C33Matrix::Rotation(float angle, const C3Vector& axis, bool unit) {
+    C3Vector axis_ = axis;
+    if (!unit) {
+        axis_.Normalize();
+    }
+    STORM_ASSERT(axis_.IsUnit());
+
+    float sina = CMath::sin(angle);
+    float cosa = CMath::cos(angle);
+
+    float xs = axis_.x * sina;
+    float ys = axis_.y * sina;
+    float zs = axis_.z * sina;
+
+    float one_c = 1.0f - cosa;
+
+    // https://en.wikipedia.org/wiki/Rotation_matrix
+    // Rotation matrix from axis and angle
+    C33Matrix result;
+
+    result.a0 = axis_.x * axis_.x * one_c + cosa;
+    result.a1 = axis_.x * axis_.y * one_c + zs;
+    result.a2 = axis_.x * axis_.z * one_c - ys;
+
+    result.b0 = axis_.x * axis_.y * one_c - zs;
+    result.b1 = axis_.y * axis_.y * one_c + cosa;
+    result.b2 = axis_.y * axis_.z * one_c + xs;
+
+    result.c0 = axis_.x * axis_.z * one_c + ys;
+    result.c1 = axis_.y * axis_.z * one_c - xs;
+    result.c2 = axis_.z * axis_.z * one_c + cosa;
+
+    return result;
+}
+
 C33Matrix C33Matrix::RotationAroundX(float angle) {
     float cosa = CMath::cos(angle);
     float sina = CMath::sin(angle);
@@ -56,63 +95,22 @@ C33Matrix C33Matrix::RotationAroundY(float angle) {
 }
 
 C33Matrix C33Matrix::RotationAroundZ(float angle) {
-    float cosa = CMath::cos(angle);
+    float cosAngle = cos(angle);
-    float sina = CMath::sin(angle);
+    float sinAngle = sin(angle);
 
-    C33Matrix result;
+    float a0 = cosAngle;
+    float a1 = sinAngle;
+    float a2 = 0.0f;
 
-    result.a0 = cosa;
+    float b0 = -sinAngle;
-    result.a1 = sina;
+    float b1 = cosAngle;
-    result.a2 = 0.0f;
+    float b2 = 0.0f;
 
-    result.b0 = -sina;
+    float c0 = 0.0f;
-    result.b1 = cosa;
+    float c1 = 0.0f;
-    result.b2 = 0.0f;
+    float c2 = 1.0f;
 
-    result.c0 = 0.0f;
+    return { a0, a1, a2, b0, b1, b2, c0, c1, c2 };
-    result.c1 = 0.0f;
-    result.c2 = 1.0f;
-
-    return result;
-}
-
-C33Matrix C33Matrix::Rotation(float angle) {
-    return C33Matrix::RotationAroundZ(angle);
-}
-
-C33Matrix C33Matrix::Rotation(float angle, const C3Vector& axis, bool unit) {
-    C3Vector axis_ = axis;
-    if (!unit) {
-        axis_.Normalize();
-    }
-    STORM_ASSERT(axis_.IsUnit());
-
-    float sina = CMath::sin(angle);
-    float cosa = CMath::cos(angle);
-
-    float xs = axis_.x * sina;
-    float ys = axis_.y * sina;
-    float zs = axis_.z * sina;
-
-    float one_c = 1.0f - cosa;
-
-    // https://en.wikipedia.org/wiki/Rotation_matrix
-    // Rotation matrix from axis and angle
-    C33Matrix result;
-
-    result.a0 = axis_.x * axis_.x * one_c + cosa;
-    result.a1 = axis_.x * axis_.y * one_c + zs;
-    result.a2 = axis_.x * axis_.z * one_c - ys;
-
-    result.b0 = axis_.x * axis_.y * one_c - zs;
-    result.b1 = axis_.y * axis_.y * one_c + cosa;
-    result.b2 = axis_.y * axis_.z * one_c + xs;
-
-    result.c0 = axis_.x * axis_.z * one_c + ys;
-    result.c1 = axis_.y * axis_.z * one_c - xs;
-    result.c2 = axis_.z * axis_.z * one_c + cosa;
-
-    return result;
 }
 
 C33Matrix::C33Matrix(const C3Vector& r0, const C3Vector& r1, const C3Vector& r2)
@@ -127,6 +125,18 @@ C33Matrix::C33Matrix(const C3Vector& r0, const C3Vector& r1, const C3Vector& r2)
     , c2(r2.z) {
 }
 
+C33Matrix::C33Matrix(float a0, float a1, float a2, float b0, float b1, float b2, float c0, float c1, float c2)
+    : a0(a0)
+    , a1(a1)
+    , a2(a2)
+    , b0(b0)
+    , b1(b1)
+    , b2(b2)
+    , c0(c0)
+    , c1(c1)
+    , c2(c2) {
+}
+
 C33Matrix::C33Matrix(const C34Matrix& m)
     : a0(m.a0)
     , a1(m.a1)
@@ -169,18 +179,6 @@ C33Matrix::C33Matrix(const C4Quaternion& rotation) {
     this->c2 = 1.0f - ((twox * rotation.x) + (twoy * rotation.y));
 }
 
-C33Matrix::C33Matrix(float a0, float a1, float a2, float b0, float b1, float b2, float c0, float c1, float c2)
-    : a0(a0)
-    , a1(a1)
-    , a2(a2)
-    , b0(b0)
-    , b1(b1)
-    , b2(b2)
-    , c0(c0)
-    , c1(c1)
-    , c2(c2) {
-}
-
 C33Matrix::C33Matrix(float a)
     : a0(a)
     , a1(a)
@@ -193,6 +191,216 @@ C33Matrix::C33Matrix(float a)
     , c2(a) {
 }
 
+C33Matrix C33Matrix::Adjoint() const {
+    float a0 = C33Matrix::Det(this->b1, this->b2, this->c1, this->c2);
+    float a1 = -C33Matrix::Det(this->a1, this->a2, this->c1, this->c2);
+    float a2 = C33Matrix::Det(this->a1, this->a2, this->b1, this->b2);
+
+    float b0 = -C33Matrix::Det(this->b0, this->b2, this->c0, this->c2);
+    float b1 = C33Matrix::Det(this->a0, this->a2, this->c0, this->c2);
+    float b2 = -C33Matrix::Det(this->a0, this->a2, this->b0, this->b2);
+
+    float c0 = C33Matrix::Det(this->b0, this->b1, this->c0, this->c1);
+    float c1 = -C33Matrix::Det(this->a0, this->a1, this->c0, this->c1);
+    float c2 = C33Matrix::Det(this->a0, this->a1, this->b0, this->b1);
+
+    return { a0, a1, a2, b0, b1, b2, c0, c1, c2 };
+}
+
+C33Matrix C33Matrix::AffineInverse(const C3Vector& v) const {
+    C33Matrix rotationScale = *this;
+    C3Vector s = { 1.0f / v.x, 1.0f / v.y, 1.0f / v.z };
+    rotationScale.Scale(s);
+    C33Matrix matrix = rotationScale.Transpose();
+    matrix.Scale(s);
+    return matrix;
+}
+
+C33Matrix C33Matrix::AffineInverse(float uniformScale) const {
+    if (CMath::fequalz(uniformScale, 1.0f, 0.00000095367432f)) {
+        return this->Transpose();
+    }
+
+    C33Matrix matrix = this->Transpose();
+    matrix.Scale(1.0f / (uniformScale * uniformScale));
+    return matrix;
+}
+
+C33Matrix C33Matrix::Cofactors() const {
+    float a0 = C33Matrix::Det(this->b1, this->b2, this->c1, this->c2);
+    float a1 = -C33Matrix::Det(this->b0, this->b2, this->c0, this->c2);
+    float a2 = C33Matrix::Det(this->b0, this->b1, this->c0, this->c1);
+
+    float b0 = -C33Matrix::Det(this->a1, this->a2, this->c1, this->c2);
+    float b1 = C33Matrix::Det(this->a0, this->a2, this->c0, this->c2);
+    float b2 = -C33Matrix::Det(this->a0, this->a1, this->c0, this->c1);
+
+    float c0 = C33Matrix::Det(this->a1, this->a2, this->b1, this->b2);
+    float c1 = -C33Matrix::Det(this->a0, this->a2, this->b0, this->b2);
+    float c2 = C33Matrix::Det(this->a0, this->a1, this->b0, this->b1);
+
+    return { a0, a1, a2, b0, b1, b2, c0, c1, c2 };
+}
+
+float C33Matrix::Determinant() const {
+    return this->c0 * this->a1 * this->b2 + this->a2 * this->b0 * this->c1 + this->c2 * this->a0 * this->b1 - this->c0 * this->a2 * this->b1 - this->c2 * this->a1 * this->b0 - this->a0 * this->c1 * this->b2;
+}
+
+void C33Matrix::FromEulerAnglesZYX(float yaw, float pitch, float roll) {
+    float cos_ = CMath::cos(yaw);
+    float sin_ = CMath::sin(yaw);
+
+    auto z = C33Matrix(
+        cos_, -sin_, 0.0f,
+        sin_,  cos_, 0.0f,
+        0.0f,  0.0f, 1.0f
+    );
+
+    cos_ = CMath::cos(pitch);
+    sin_ = CMath::sin(pitch);
+
+    auto y = C33Matrix(
+         cos_, 0.0f, sin_,
+         0.0f, 1.0f, 0.0f,
+        -sin_, 0.0f, cos_
+    );
+
+    cos_ = CMath::cos(roll);
+    sin_ = CMath::sin(roll);
+
+    auto x = C33Matrix(
+        1.0f, 0.0f,  0.0f,
+        0.0f, cos_, -sin_,
+        0.0f, sin_,  cos_
+    );
+
+    *this = (z * (y * x)).Transpose();
+}
+
+void C33Matrix::Identity() {
+    this->a0 = 1.0f;
+    this->a1 = 0.0f;
+    this->a2 = 0.0f;
+
+    this->b0 = 0.0f;
+    this->b1 = 1.0f;
+    this->b2 = 0.0f;
+
+    this->c0 = 0.0f;
+    this->c1 = 0.0f;
+    this->c2 = 1.0f;
+}
+
+C33Matrix C33Matrix::Inverse() const {
+    return this->Inverse(this->Determinant());
+}
+
+C33Matrix C33Matrix::Inverse(float det) const {
+    STORM_ASSERT(CMath::fequal(det, 0.0f) == false);
+    return this->Adjoint() * (1.0f / det);
+}
+
+void C33Matrix::Rotate(const C4Quaternion& rotation) {
+    float twox = 2.0f * rotation.x;
+    float twoy = 2.0f * rotation.y;
+    float twoz = 2.0f * rotation.z;
+
+    C33Matrix left;
+
+    left.a0 = 1.0f - ((twoy * rotation.y) + (twoz * rotation.z));
+    left.a1 = (twoy * rotation.x) + (twoz * rotation.w);
+    left.a2 = (twoz * rotation.x) - (twoy * rotation.w);
+
+    left.b0 = (twoy * rotation.x) - (twoz * rotation.w);
+    left.b1 = 1.0f - ((twox * rotation.x) + (twoz * rotation.z));
+    left.b2 = (twoz * rotation.y) + (twox * rotation.w);
+
+    left.c0 = (twoz * rotation.x) + (twoy * rotation.w);
+    left.c1 = (twoz * rotation.y) - (twox * rotation.w);
+    left.c2 = 1.0f - ((twox * rotation.x) + (twoy * rotation.y));
+
+    *this = left * (*this);
+}
+
+void C33Matrix::Rotate(float angle) {
+    *this = C33Matrix::Rotation(angle) * (*this);
+}
+
+void C33Matrix::Rotate(float angle, const C3Vector& axis, bool unit) {
+    *this = C33Matrix::Rotation(angle, axis, unit) * (*this);
+}
+
+void C33Matrix::Scale(const C2Vector& scale) {
+    Scale(scale.x, scale.y);
+}
+
+void C33Matrix::Scale(const C3Vector& scale) {
+    Scale(scale.x, scale.y, scale.z);
+}
+
+void C33Matrix::Scale(float scale) {
+    Scale(scale, scale, scale);
+}
+
+void C33Matrix::Scale(float x, float y) {
+    this->a0 *= x;
+    this->a1 *= x;
+
+    this->b0 *= y;
+    this->b1 *= y;
+}
+
+void C33Matrix::Scale(float x, float y, float z) {
+    this->a0 *= x;
+    this->a1 *= x;
+    this->a2 *= x;
+
+    this->b0 *= y;
+    this->b1 *= y;
+    this->b2 *= y;
+
+    this->c0 *= z;
+    this->c1 *= z;
+    this->c2 *= z;
+}
+
+float C33Matrix::Trace() {
+    return this->a0 + this->b1 + this->c2;
+}
+
+void C33Matrix::Translate(const C2Vector& move) {
+    this->c0 += this->a0 * move.x + this->b0 * move.y;
+    this->c1 += this->a1 * move.x + this->b1 * move.y;
+}
+
+C33Matrix C33Matrix::Transpose() const {
+    return {
+        this->a0,
+        this->b0,
+        this->c0,
+        this->a1,
+        this->b1,
+        this->c1,
+        this->a2,
+        this->b2,
+        this->c2
+    };
+}
+
+void C33Matrix::Zero() {
+    this->a0 = 0.0f;
+    this->a1 = 0.0f;
+    this->a2 = 0.0f;
+
+    this->b0 = 0.0f;
+    this->b1 = 0.0f;
+    this->b2 = 0.0f;
+
+    this->c0 = 0.0f;
+    this->c1 = 0.0f;
+    this->c2 = 0.0f;
+}
+
 C33Matrix& C33Matrix::operator+=(const C33Matrix& a) {
     this->a0 += a.a0;
     this->a1 += a.a1;
@@ -262,282 +470,6 @@ C33Matrix& C33Matrix::operator/=(float a) {
     return *this;
 }
 
-void C33Matrix::Zero() {
-    this->a0 = 0.0f;
-    this->a1 = 0.0f;
-    this->a2 = 0.0f;
-
-    this->b0 = 0.0f;
-    this->b1 = 0.0f;
-    this->b2 = 0.0f;
-
-    this->c0 = 0.0f;
-    this->c1 = 0.0f;
-    this->c2 = 0.0f;
-}
-
-void C33Matrix::Identity() {
-    this->a0 = 1.0f;
-    this->a1 = 0.0f;
-    this->a2 = 0.0f;
-
-    this->b0 = 0.0f;
-    this->b1 = 1.0f;
-    this->b2 = 0.0f;
-
-    this->c0 = 0.0f;
-    this->c1 = 0.0f;
-    this->c2 = 1.0f;
-}
-
-float C33Matrix::Trace() {
-    return this->a0 + this->b1 + this->c2;
-}
-
-void C33Matrix::Scale(float scale) {
-    Scale(scale, scale, scale);
-}
-
-void C33Matrix::Scale(float x, float y) {
-    this->a0 *= x;
-    this->a1 *= x;
-
-    this->b0 *= y;
-    this->b1 *= y;
-}
-
-void C33Matrix::Scale(const C2Vector& scale) {
-    Scale(scale.x, scale.y);
-}
-
-void C33Matrix::Scale(float x, float y, float z) {
-    this->a0 *= x;
-    this->a1 *= x;
-    this->a2 *= x;
-
-    this->b0 *= y;
-    this->b1 *= y;
-    this->b2 *= y;
-
-    this->c0 *= z;
-    this->c1 *= z;
-    this->c2 *= z;
-}
-
-void C33Matrix::Scale(const C3Vector& scale) {
-    Scale(scale.x, scale.y, scale.z);
-}
-
-void C33Matrix::Rotate(float angle) {
-    *this = C33Matrix::Rotation(angle) * (*this);
-}
-
-void C33Matrix::Rotate(const C4Quaternion& rotation) {
-    float twox = 2.0f * rotation.x;
-    float twoy = 2.0f * rotation.y;
-    float twoz = 2.0f * rotation.z;
-
-    C33Matrix left;
-
-    left.a0 = 1.0f - ((twoy * rotation.y) + (twoz * rotation.z));
-    left.a1 = (twoy * rotation.x) + (twoz * rotation.w);
-    left.a2 = (twoz * rotation.x) - (twoy * rotation.w);
-
-    left.b0 = (twoy * rotation.x) - (twoz * rotation.w);
-    left.b1 = 1.0f - ((twox * rotation.x) + (twoz * rotation.z));
-    left.b2 = (twoz * rotation.y) + (twox * rotation.w);
-
-    left.c0 = (twoz * rotation.x) + (twoy * rotation.w);
-    left.c1 = (twoz * rotation.y) - (twox * rotation.w);
-    left.c2 = 1.0f - ((twox * rotation.x) + (twoy * rotation.y));
-
-    *this = left * (*this);
-}
-
-void C33Matrix::Rotate(float angle, const C3Vector& axis, bool unit) {
-    *this = C33Matrix::Rotation(angle, axis, unit) * (*this);
-}
-
-void C33Matrix::Translate(const C2Vector& move) {
-    this->c0 += this->a0 * move.x + this->b0 * move.y;
-    this->c1 += this->a1 * move.x + this->b1 * move.y;
-}
-
-C33Matrix C33Matrix::Transpose() const {
-    return {
-        this->a0,
-        this->b0,
-        this->c0,
-        this->a1,
-        this->b1,
-        this->c1,
-        this->a2,
-        this->b2,
-        this->c2
-    };
-}
-
-float C33Matrix::Determinant() const {
- return this->c0 * this->a1 * this->b2
-    + this->a2 * this->b0 * this->c1
-    + this->c2 * this->a0 * this->b1
-    - this->c0 * this->a2 * this->b1
-    - this->c2 * this->a1 * this->b0
-    - this->a0 * this->c1 * this->b2;
-}
-
-C33Matrix C33Matrix::Cofactors() const {
-    float a0 = C33Matrix::Det(this->b1, this->b2, this->c1, this->c2);
-    float a1 = -C33Matrix::Det(this->b0, this->b2, this->c0, this->c2);
-    float a2 = C33Matrix::Det(this->b0, this->b1, this->c0, this->c1);
-
-    float b0 = -C33Matrix::Det(this->a1, this->a2, this->c1, this->c2);
-    float b1 = C33Matrix::Det(this->a0, this->a2, this->c0, this->c2);
-    float b2 = -C33Matrix::Det(this->a0, this->a1, this->c0, this->c1);
-
-    float c0 = C33Matrix::Det(this->a1, this->a2, this->b1, this->b2);
-    float c1 = -C33Matrix::Det(this->a0, this->a2, this->b0, this->b2);
-    float c2 = C33Matrix::Det(this->a0, this->a1, this->b0, this->b1);
-
-    return { a0, a1, a2, b0, b1, b2, c0, c1, c2 };
-}
-
-C33Matrix C33Matrix::Adjoint() const {
-    float a0 = C33Matrix::Det(this->b1, this->b2, this->c1, this->c2);
-    float a1 = -C33Matrix::Det(this->a1, this->a2, this->c1, this->c2);
-    float a2 = C33Matrix::Det(this->a1, this->a2, this->b1, this->b2);
-
-    float b0 = -C33Matrix::Det(this->b0, this->b2, this->c0, this->c2);
-    float b1 = C33Matrix::Det(this->a0, this->a2, this->c0, this->c2);
-    float b2 = -C33Matrix::Det(this->a0, this->a2, this->b0, this->b2);
-
-    float c0 = C33Matrix::Det(this->b0, this->b1, this->c0, this->c1);
-    float c1 = -C33Matrix::Det(this->a0, this->a1, this->c0, this->c1);
-    float c2 = C33Matrix::Det(this->a0, this->a1, this->b0, this->b1);
-
-    return { a0, a1, a2, b0, b1, b2, c0, c1, c2 };
-}
-
-C33Matrix C33Matrix::Inverse() const {
-    return this->Inverse(this->Determinant());
-}
-
-C33Matrix C33Matrix::Inverse(float det) const {
-    STORM_ASSERT(CMath::fequal(det, 0.0f) == false);
-    return this->Adjoint() * (1.0f / det);
-}
-
-C33Matrix C33Matrix::AffineInverse(const C3Vector& v) const {
-    C33Matrix rotationScale = *this;
-    C3Vector s = { 1.0f / v.x, 1.0f / v.y, 1.0f / v.z };
-    rotationScale.Scale(s);
-    C33Matrix matrix = rotationScale.Transpose();
-    matrix.Scale(s);
-    return matrix;
-}
-
-C33Matrix C33Matrix::AffineInverse(float a) const {
-    if (CMath::fequalz(a, 1.0f, 0.00000095367432f)) {
-        return this->Transpose();
-    }
-
-    C33Matrix matrix = this->Transpose();
-    matrix.Scale(1.0f / (a * a));
-    return matrix;
-}
-
-bool C33Matrix::ToEulerAnglesXYZ(float& xa_, float& ya_, float& za_) {
-    // TODO
-    throw;
-    return false;
-}
-
-bool C33Matrix::ToEulerAnglesXZY(float& xa_, float& za_, float& ya_) {
-    // TODO
-    throw;
-    return false;
-}
-
-bool C33Matrix::ToEulerAnglesYXZ(float& ya_, float& xa_, float& za_) {
-    // TODO
-    throw;
-    return false;
-}
-
-bool C33Matrix::ToEulerAnglesYZX(float& ya_, float& za_, float& xa_) {
-    // TODO
-    throw;
-    return false;
-}
-
-bool C33Matrix::ToEulerAnglesZXY(float& za_, float& xa_, float& ya_) {
-    // TODO
-    throw;
-    return false;
-}
-
-bool C33Matrix::ToEulerAnglesZYX(float& za_, float& ya_, float& xa_) {
-    // TODO
-    throw;
-    return false;
-}
-
-void C33Matrix::FromEulerAnglesXYZ(float yaw, float pitch, float roll) {
-    // TODO
-    throw;
-}
-
-void C33Matrix::FromEulerAnglesXZY(float yaw, float pitch, float roll) {
-    // TODO
-    throw;
-}
-
-void C33Matrix::FromEulerAnglesYXZ(float yaw, float pitch, float roll) {
-    // TODO
-    throw;
-}
-
-void C33Matrix::FromEulerAnglesYZX(float yaw, float pitch, float roll) {
-    // TODO
-    throw;
-}
-
-void C33Matrix::FromEulerAnglesZXY(float yaw, float pitch, float roll) {
-    // TODO
-    throw;
-}
-
-void C33Matrix::FromEulerAnglesZYX(float yaw, float pitch, float roll) {
-    float siny = CMath::sin(yaw);
-    float cosy = CMath::cos(yaw);
-
-    float sinp = CMath::sin(pitch);
-    float cosp = CMath::cos(pitch);
-
-    float sinr = CMath::sin(roll);
-    float cosr = CMath::cos(roll);
-
-    C33Matrix x_ = {
-        1.0f, 0.0f, 0.0f,
-        0.0f, cosr, -sinr,
-        0.0f, sinr, cosr
-    };
-
-    C33Matrix y_ = {
-        cosp, 0.0f, sinp,
-        0.0f, 1.0f, 0.0f,
-        -sinp, 0.0f, cosp
-    };
-
-    C33Matrix z_ = {
-        cosy, -siny, 0.0f,
-        siny, cosy, 0.0f,
-        0.0f, 0.0f, 1.0f
-    };
-
-    *this = (z_ * y_ * x_).Transpose();
-}
-
 C33Matrix operator*(const C33Matrix& l, const C33Matrix& r) {
     float a0 = l.a0 * r.a0 + l.a1 * r.b0 + l.a2 * r.c0;
     float a1 = l.a0 * r.a1 + l.a1 * r.b1 + l.a2 * r.c1;

tempest/matrix/C33Matrix.hpp

@@ -17,12 +17,12 @@ class C33Matrix {
 
     // Static functions
     static float Det(float a, float b, float c, float d);
+    static C33Matrix Rotation(float angle);
+    static C33Matrix Rotation(float angle, const C3Vector& axis, bool unit);
     static C33Matrix RotationAroundX(float angle);
     static C33Matrix RotationAroundY(float angle);
     static C33Matrix RotationAroundZ(float angle);
 
-    static C33Matrix Rotation(float angle);
-    static C33Matrix Rotation(float angle, const C3Vector& axis, bool unit);
 
     // Member variables
     float a0 = 1.0f;
@@ -38,51 +38,41 @@ class C33Matrix {
     // Member functions
     C33Matrix() = default;
     C33Matrix(const C3Vector& r0, const C3Vector& r1, const C3Vector& r2);
+    C33Matrix(float a0, float a1, float a2, float b0, float b1, float b2, float c0, float c1, float c2);
     explicit C33Matrix(const C34Matrix& m);
     explicit C33Matrix(const C44Matrix& m);
     explicit C33Matrix(const C4Quaternion& rotation);
-    C33Matrix(float a0, float a1, float a2, float b0, float b1, float b2, float c0, float c1, float c2);
     explicit C33Matrix(float a);
 
+    C33Matrix Adjoint() const;
+    C33Matrix AffineInverse(const C3Vector& v) const;
+    C33Matrix AffineInverse(float uniformScale) const;
+    C33Matrix Cofactors() const;
+
+    float Determinant() const;
+    void FromEulerAnglesZYX(float yaw, float pitch, float roll);
+    void Identity();
+    C33Matrix Inverse() const;
+    C33Matrix Inverse(float det) const;
+    void Rotate(const C4Quaternion& rotation);
+    void Rotate(float angle);
+    void Rotate(float angle, const C3Vector& axis, bool unit);
+    void Scale(const C2Vector& scale);
+    void Scale(const C3Vector& scale);
+    void Scale(float scale);
+    void Scale(float x, float y);
+    void Scale(float x, float y, float z);
+    float Trace();
+    void Translate(const C2Vector& move);
+    C33Matrix Transpose() const;
+
+    void Zero();
+
     C33Matrix& operator+=(const C33Matrix& a);
     C33Matrix& operator-=(const C33Matrix& a);
     C33Matrix& operator*=(float a);
     C33Matrix& operator*=(const C33Matrix& a);
     C33Matrix& operator/=(float a);
-
-    void Zero();
-    void Identity();
-    float Trace();
-    void Scale(float scale);
-    void Scale(float x, float y);
-    void Scale(const C2Vector& scale);
-    void Scale(float x, float y, float z);
-    void Scale(const C3Vector& scale);
-    void Rotate(float angle);
-    void Rotate(const C4Quaternion& rotation);
-    void Rotate(float angle, const C3Vector& axis, bool unit);
-    void Translate(const C2Vector& move);
-    C33Matrix Transpose() const;
-    float Determinant() const;
-    C33Matrix Cofactors() const;
-    C33Matrix Adjoint() const;
-    C33Matrix Inverse() const;
-    C33Matrix Inverse(float det) const;
-    C33Matrix AffineInverse(const C3Vector& v) const;
-    C33Matrix AffineInverse(float a) const;
-
-    bool ToEulerAnglesXYZ(float& xa_, float& ya_, float& za_);
-    bool ToEulerAnglesXZY(float& xa_, float& za_, float& ya_);
-    bool ToEulerAnglesYXZ(float& ya_, float& xa_, float& za_);
-    bool ToEulerAnglesYZX(float& ya_, float& za_, float& xa_);
-    bool ToEulerAnglesZXY(float& za_, float& xa_, float& ya_);
-    bool ToEulerAnglesZYX(float& za_, float& ya_, float& xa_);
-    void FromEulerAnglesXYZ(float yaw, float pitch, float roll);
-    void FromEulerAnglesXZY(float yaw, float pitch, float roll);
-    void FromEulerAnglesYXZ(float yaw, float pitch, float roll);
-    void FromEulerAnglesYZX(float yaw, float pitch, float roll);
-    void FromEulerAnglesZXY(float yaw, float pitch, float roll);
-    void FromEulerAnglesZYX(float yaw, float pitch, float roll);
 };
 
 C33Matrix operator*(const C33Matrix& l, const C33Matrix& r);

tempest/matrix/C44Matrix.cpp

@@ -1,7 +1,6 @@
 #include "tempest/matrix/C44Matrix.hpp"
 #include "tempest/math/CMath.hpp"
 #include "tempest/vector/C3Vector.hpp"
-#include "tempest/vector/C4Vector.hpp"
 #include "tempest/matrix/C33Matrix.hpp"
 #include "tempest/quaternion/C4Quaternion.hpp"
 
@@ -12,35 +11,6 @@ float C44Matrix::Det(float a, float b, float c, float d, float e, float f, float
     return (b * f * g) + (c * d * h) + (a * e * i) - (c * e * g) - (b * d * i) - (a * f * h);
 }
 
-C44Matrix C44Matrix::RotationAroundZ(float angle) {
-    float cosa = CMath::cos(angle);
-    float sina = CMath::sin(angle);
-
-    C44Matrix result;
-
-    result.a0 = cosa;
-    result.a1 = sina;
-    result.a2 = 0.0f;
-    result.a3 = 0.0f;
-
-    result.b0 = -sina;
-    result.b1 = cosa;
-    result.b2 = 0.0f;
-    result.b3 = 0.0f;
-
-    result.c0 = 0.0f;
-    result.c1 = 0.0f;
-    result.c2 = 1.0f;
-    result.c3 = 0.0f;
-
-    result.d0 = 0.0f;
-    result.d1 = 0.0f;
-    result.d2 = 0.0f;
-    result.d3 = 1.0f;
-
-    return result;
-}
-
 C44Matrix C44Matrix::Rotation(float angle, const C3Vector& axis, bool unit) {
     C3Vector axis_ = axis;
     if (!unit) {
@@ -82,6 +52,33 @@ C44Matrix C44Matrix::Rotation(float angle, const C3Vector& axis, bool unit) {
     return result;
 }
 
+C44Matrix C44Matrix::RotationAroundZ(float angle) {
+    float cosAngle = cos(angle);
+    float sinAngle = sin(angle);
+
+    float a0 = cosAngle;
+    float a1 = sinAngle;
+    float a2 = 0.0f;
+    float a3 = 0.0f;
+
+    float b0 = -sinAngle;
+    float b1 = cosAngle;
+    float b2 = 0.0f;
+    float b3 = 0.0f;
+
+    float c0 = 0.0f;
+    float c1 = 0.0f;
+    float c2 = 1.0f;
+    float c3 = 0.0f;
+
+    float d0 = 0.0f;
+    float d1 = 0.0f;
+    float d2 = 0.0f;
+    float d3 = 1.0f;
+
+    return { a0, a1, a2, a3, b0, b1, b2, b3, c0, c1, c2, c3, d0, d1, d2, d3 };
+}
+
 C44Matrix::C44Matrix(const C4Vector& r0, const C4Vector& r1, const C4Vector& r2, const C4Vector& r3)
     : a0(r0.x)
     , a1(r0.y)
@@ -101,6 +98,25 @@ C44Matrix::C44Matrix(const C4Vector& r0, const C4Vector& r1, const C4Vector& r2,
     , d3(r3.w) {
 }
 
+C44Matrix::C44Matrix(float a0, float a1, float a2, float a3, float b0, float b1, float b2, float b3, float c0, float c1, float c2, float c3, float d0, float d1, float d2, float d3)
+    : a0(a0)
+    , a1(a1)
+    , a2(a2)
+    , a3(a3)
+    , b0(b0)
+    , b1(b1)
+    , b2(b2)
+    , b3(b3)
+    , c0(c0)
+    , c1(c1)
+    , c2(c2)
+    , c3(c3)
+    , d0(d0)
+    , d1(d1)
+    , d2(d2)
+    , d3(d3) {
+}
+
 C44Matrix::C44Matrix(const C33Matrix& m)
     : a0(m.a0)
     , a1(m.a1)
@@ -157,25 +173,6 @@ C44Matrix::C44Matrix(const C4Quaternion& rotation) {
     this->c2 = 1.0f - (txx + tyy);
 }
 
-C44Matrix::C44Matrix(float a0, float a1, float a2, float a3, float b0, float b1, float b2, float b3, float c0, float c1, float c2, float c3, float d0, float d1, float d2, float d3)
-    : a0(a0)
-    , a1(a1)
-    , a2(a2)
-    , a3(a3)
-    , b0(b0)
-    , b1(b1)
-    , b2(b2)
-    , b3(b3)
-    , c0(c0)
-    , c1(c1)
-    , c2(c2)
-    , c3(c3)
-    , d0(d0)
-    , d1(d1)
-    , d2(d2)
-    , d3(d3) {
-}
-
 C44Matrix::C44Matrix(float a)
     : a0(a)
     , a1(a)
@@ -195,270 +192,73 @@ C44Matrix::C44Matrix(float a)
     , d3(a) {
 }
 
-C44Matrix& C44Matrix::operator+=(const C44Matrix& a) {
+C44Matrix C44Matrix::Adjoint() const {
-    this->a0 += a.a0;
+    float a0 = this->b2 * this->c3 * this->d1 + this->b3 * this->c1 * this->d2 + this->b1 * this->c2 * this->d3 - this->b3 * this->c2 * this->d1 - this->b2 * this->c1 * this->d3 - this->b1 * this->c3 * this->d2;
-    this->a1 += a.a1;
+    float a1 = -(this->a2 * this->c3 * this->d1 + this->a3 * this->c1 * this->d2 + this->a1 * this->c2 * this->d3 - this->a3 * this->c2 * this->d1 - this->a2 * this->c1 * this->d3 - this->a1 * this->c3 * this->d2);
-    this->a2 += a.a2;
+    float a2 = this->a1 * this->b2 * this->d3 + this->a2 * this->b3 * this->d1 + this->a3 * this->b1 * this->d2 - this->a3 * this->b2 * this->d1 - this->a2 * this->b1 * this->d3 - this->a1 * this->b3 * this->d2;
-    this->a3 += a.a3;
+    float a3 = -(this->a2 * this->b3 * this->c1 + this->a3 * this->b1 * this->c2 + this->a1 * this->b2 * this->c3 - this->a3 * this->b2 * this->c1 - this->a2 * this->b1 * this->c3 - this->a1 * this->b3 * this->c2);
 
-    this->b0 += a.b0;
+    float b0 = -(this->b3 * this->c0 * this->d2 + this->b0 * this->c2 * this->d3 + this->b2 * this->c3 * this->d0 - this->b3 * this->c2 * this->d0 - this->b2 * this->c0 * this->d3 - this->b0 * this->c3 * this->d2);
-    this->b1 += a.b1;
+    float b1 = this->a2 * this->c3 * this->d0 + this->a3 * this->c0 * this->d2 + this->a0 * this->c2 * this->d3 - this->a3 * this->c2 * this->d0 - this->a2 * this->c0 * this->d3 - this->a0 * this->c3 * this->d2;
-    this->b2 += a.b2;
+    float b2 = -(this->a2 * this->b3 * this->d0 + this->a3 * this->b0 * this->d2 + this->a0 * this->b2 * this->d3 - this->a3 * this->b2 * this->d0 - this->a2 * this->b0 * this->d3 - this->a0 * this->b3 * this->d2);
-    this->b3 += a.b3;
+    float b3 = this->a2 * this->b3 * this->c0 + this->a3 * this->b0 * this->c2 + this->a0 * this->b2 * this->c3 - this->a3 * this->b2 * this->c0 - this->a2 * this->b0 * this->c3 - this->a0 * this->b3 * this->c2;
 
-    this->c0 += a.c0;
+    float c0 = this->b1 * this->c3 * this->d0 + this->b3 * this->c0 * this->d1 + this->b0 * this->c1 * this->d3 - this->b3 * this->c1 * this->d0 - this->b1 * this->c0 * this->d3 - this->b0 * this->c3 * this->d1;
-    this->c1 += a.c1;
+    float c1 = -(this->a1 * this->c3 * this->d0 + this->a3 * this->c0 * this->d1 + this->a0 * this->c1 * this->d3 - this->a3 * this->c1 * this->d0 - this->a1 * this->c0 * this->d3 - this->a0 * this->c3 * this->d1);
-    this->c2 += a.c2;
+    float c2 = this->a1 * this->b3 * this->d0 + this->a3 * this->b0 * this->d1 + this->a0 * this->b1 * this->d3 - this->a3 * this->b1 * this->d0 - this->a1 * this->b0 * this->d3 - this->a0 * this->b3 * this->d1;
-    this->c3 += a.c3;
+    float c3 = -(this->a1 * this->b3 * this->c0 + this->a3 * this->b0 * this->c1 + this->a0 * this->b1 * this->c3 - this->a3 * this->b1 * this->c0 - this->a1 * this->b0 * this->c3 - this->a0 * this->b3 * this->c1);
 
-    this->d0 += a.d0;
+    float d0 = -(this->b2 * this->c0 * this->d1 + this->b0 * this->c1 * this->d2 + this->b1 * this->c2 * this->d0 - this->b2 * this->c1 * this->d0 - this->b1 * this->c0 * this->d2 - this->b0 * this->c2 * this->d1);
-    this->d1 += a.d1;
+    float d1 = this->a1 * this->c2 * this->d0 + this->a2 * this->c0 * this->d1 + this->a0 * this->c1 * this->d2 - this->a2 * this->c1 * this->d0 - this->a1 * this->c0 * this->d2 - this->a0 * this->c2 * this->d1;
-    this->d2 += a.d2;
+    float d2 = -(this->a1 * this->b2 * this->d0 + this->a2 * this->b0 * this->d1 + this->a0 * this->b1 * this->d2 - this->a2 * this->b1 * this->d0 - this->a1 * this->b0 * this->d2 - this->a0 * this->b2 * this->d1);
-    this->d3 += a.d3;
+    float d3 = this->a0 * this->b1 * this->c2 + this->a1 * this->b2 * this->c0 + this->a2 * this->b0 * this->c1 - this->b1 * this->a2 * this->c0 - this->c2 * this->b0 * this->a1 - this->a0 * this->b2 * this->c1;
 
-    return *this;
+    return { a0, a1, a2, a3, b0, b1, b2, b3, c0, c1, c2, c3, d0, d1, d2, d3 };
 }
 
-C44Matrix& C44Matrix::operator-=(const C44Matrix& a) {
+/**
-    this->a0 -= a.a0;
+ * Calculate and return the inverse of this affine transform matrix. This function assumes no
-    this->a1 -= a.a1;
+ * scaling factor is present in the transform.
-    this->a2 -= a.a2;
+ *
-    this->a3 -= a.a3;
+ * @return Inverse of matrix
+ */
+C44Matrix C44Matrix::AffineInverse() const {
+    auto matrix = C44Matrix(C33Matrix(*this).Transpose());
+    matrix.Translate(C3Vector(-this->d0, -this->d1, -this->d2));
 
-    this->b0 -= a.b0;
+    return matrix;
-    this->b1 -= a.b1;
-    this->b2 -= a.b2;
-    this->b3 -= a.b3;
-
-    this->c0 -= a.c0;
-    this->c1 -= a.c1;
-    this->c2 -= a.c2;
-    this->c3 -= a.c3;
-
-    this->d0 -= a.d0;
-    this->d1 -= a.d1;
-    this->d2 -= a.d2;
-    this->d3 -= a.d3;
-
-    return *this;
 }
 
-C44Matrix& C44Matrix::operator*=(float a) {
+C44Matrix C44Matrix::AffineInverse(const C3Vector& v) const {
-    this->a0 *= a;
+    C33Matrix rotationScale(*this);
-    this->a1 *= a;
+    C3Vector s = { 1.0f / v.x, 1.0f / v.y, 1.0f / v.z };
-    this->a2 *= a;
+    rotationScale.Scale(s);
-    this->a3 *= a;
+    C44Matrix matrix(rotationScale.Transpose());
-
+    matrix.Scale(s);
-    this->b0 *= a;
+    C3Vector move(-this->d0, -this->d1, -this->d2);
-    this->b1 *= a;
+    matrix.Translate(move);
-    this->b2 *= a;
+    return matrix;
-    this->b3 *= a;
-
-    this->c0 *= a;
-    this->c1 *= a;
-    this->c2 *= a;
-    this->c3 *= a;
-
-    this->d0 *= a;
-    this->d1 *= a;
-    this->d2 *= a;
-    this->d3 *= a;
-
-    return *this;
 }
 
-C44Matrix& C44Matrix::operator*=(const C44Matrix& a) {
+/**
-    *this = *this * a;
+ * Calculate and return the inverse of this affine transform matrix. This function assumes a
-    return *this;
+ * uniform scaling factor may be present in the transform.
+ *
+ * @param uniformScale Uniform scaling factor
+ * @return Inverse of matrix
+ */
+C44Matrix C44Matrix::AffineInverse(float uniformScale) const {
+    // No scaling factor present
+    if (CMath::fequal(uniformScale, 1.0f)) {
+        return this->AffineInverse();
     }
 
-C44Matrix& C44Matrix::operator/=(float a) {
+    // Uniform scaling factor present
-    this->a0 /= a;
+    auto matrix = C44Matrix(C33Matrix(*this).Transpose());
-    this->a1 /= a;
+    matrix.Scale(1.0f / (uniformScale * uniformScale));
-    this->a2 /= a;
+    matrix.Translate(C3Vector(-this->d0, -this->d1, -this->d2));
-    this->a3 /= a;
 
-    this->b0 /= a;
+    return matrix;
-    this->b1 /= a;
-    this->b2 /= a;
-    this->b3 /= a;
-
-    this->c0 /= a;
-    this->c1 /= a;
-    this->c2 /= a;
-    this->c3 /= a;
-
-    this->d0 /= a;
-    this->d1 /= a;
-    this->d2 /= a;
-    this->d3 /= a;
-
-    return *this;
-}
-
-void C44Matrix::Zero() {
-    this->a0 = 0.0f;
-    this->a1 = 0.0f;
-    this->a2 = 0.0f;
-    this->a3 = 0.0f;
-
-    this->b0 = 0.0f;
-    this->b1 = 0.0f;
-    this->b2 = 0.0f;
-    this->b3 = 0.0f;
-
-    this->c0 = 0.0f;
-    this->c1 = 0.0f;
-    this->c2 = 0.0f;
-    this->c3 = 0.0f;
-
-    this->d0 = 0.0f;
-    this->d1 = 0.0f;
-    this->d2 = 0.0f;
-    this->d3 = 0.0f;
-}
-
-void C44Matrix::Identity() {
-    this->a0 = 1.0f;
-    this->a1 = 0.0f;
-    this->a2 = 0.0f;
-    this->a3 = 0.0f;
-
-    this->b0 = 0.0f;
-    this->b1 = 1.0f;
-    this->b2 = 0.0f;
-    this->b3 = 0.0f;
-
-    this->c0 = 0.0f;
-    this->c1 = 0.0f;
-    this->c2 = 1.0f;
-    this->c3 = 0.0f;
-
-    this->d0 = 0.0f;
-    this->d1 = 0.0f;
-    this->d2 = 0.0f;
-    this->d3 = 1.0f;
-}
-
-float C44Matrix::Trace() {
-    return this->a0 + this->b1 + this->c2 + this->d3;
-}
-
-void C44Matrix::Translate(const C3Vector& move) {
-    this->d0 = this->a0 * move.x + this->b0 * move.y + this->c0 * move.z + this->d0;
-    this->d1 = this->a1 * move.x + this->b1 * move.y + this->c1 * move.z + this->d1;
-    this->d2 = this->a2 * move.x + this->b2 * move.y + this->c2 * move.z + this->d2;
-}
-
-void C44Matrix::Scale(float scale) {
-    this->a0 *= scale;
-    this->a1 *= scale;
-    this->a2 *= scale;
-
-    this->b0 *= scale;
-    this->b1 *= scale;
-    this->b2 *= scale;
-
-    this->c0 *= scale;
-    this->c1 *= scale;
-    this->c2 *= scale;
-}
-
-void C44Matrix::Scale(const C3Vector& scale) {
-    this->a0 *= scale.x;
-    this->a1 *= scale.x;
-    this->a2 *= scale.x;
-
-    this->b0 *= scale.y;
-    this->b1 *= scale.y;
-    this->b2 *= scale.y;
-
-    this->c0 *= scale.z;
-    this->c1 *= scale.z;
-    this->c2 *= scale.z;
-}
-
-void C44Matrix::RotateAroundZ(float angle) {
-    *this = C44Matrix::RotationAroundZ(angle) * (*this);
-}
-
-void C44Matrix::Rotate(const C4Quaternion& rotation) {
-    *this = C44Matrix(rotation) * (*this);
-}
-
-void C44Matrix::Rotate(float angle, const C3Vector& axis, bool unit) {
-    *this = C44Matrix::Rotation(angle, axis, unit) * (*this);
-}
-
-C44Matrix C44Matrix::Transpose() const {
-    return {
-        this->a0,
-        this->b0,
-        this->c0,
-        this->d0,
-        this->a1,
-        this->b1,
-        this->c1,
-        this->d1,
-        this->a2,
-        this->b2,
-        this->c2,
-        this->d2,
-        this->a3,
-        this->b3,
-        this->c3,
-        this->d3
-    };
-}
-
-float C44Matrix::Determinant() const {
-    float v2;
-    float v3;
-    float v4;
-
-    v2 = this->a0 * C44Matrix::Det(
-        this->b1,
-        this->b2,
-        this->b3,
-        this->c1,
-        this->c2,
-        this->c3,
-        this->d1,
-        this->d2,
-        this->d3);
-
-    v3 = v2 - this->a1 * C44Matrix::Det(
-        this->b0,
-        this->b2,
-        this->b3,
-        this->c0,
-        this->c2,
-        this->c3,
-        this->d0,
-        this->d2,
-        this->d3);
-
-    v4 = v3 + this->a2 * C44Matrix::Det(
-        this->b0,
-        this->b1,
-        this->b3,
-        this->c0,
-        this->c1,
-        this->c3,
-        this->d0,
-        this->d1,
-        this->d3);
-
-    return v4 - this->a3 * C44Matrix::Det(
-        this->b0,
-        this->b1,
-        this->b2,
-        this->c0,
-        this->c1,
-        this->c2,
-        this->d0,
-        this->d1,
-        this->d2);
 }
 
 C44Matrix C44Matrix::Cofactors() const {
@@ -625,28 +425,46 @@ C44Matrix C44Matrix::Cofactors() const {
     return { a0, a1, a2, a3, b0, b1, b2, b3, c0, c1, c2, c3, d0, d1, d2, d3 };
 }
 
-C44Matrix C44Matrix::Adjoint() const {
+C4Vector C44Matrix::Col0() const {
-    float a0 = this->b2 * this->c3 * this->d1 + this->b3 * this->c1 * this->d2 + this->b1 * this->c2 * this->d3 - this->b3 * this->c2 * this->d1 - this->b2 * this->c1 * this->d3 - this->b1 * this->c3 * this->d2;
+    return { this->a0, this->b0, this->c0, this->d0 };
-    float a1 = -(this->a2 * this->c3 * this->d1 + this->a3 * this->c1 * this->d2 + this->a1 * this->c2 * this->d3 - this->a3 * this->c2 * this->d1 - this->a2 * this->c1 * this->d3 - this->a1 * this->c3 * this->d2);
+}
-    float a2 = this->a1 * this->b2 * this->d3 + this->a2 * this->b3 * this->d1 + this->a3 * this->b1 * this->d2 - this->a3 * this->b2 * this->d1 - this->a2 * this->b1 * this->d3 - this->a1 * this->b3 * this->d2;
-    float a3 = -(this->a2 * this->b3 * this->c1 + this->a3 * this->b1 * this->c2 + this->a1 * this->b2 * this->c3 - this->a3 * this->b2 * this->c1 - this->a2 * this->b1 * this->c3 - this->a1 * this->b3 * this->c2);
 
-    float b0 = -(this->b3 * this->c0 * this->d2 + this->b0 * this->c2 * this->d3 + this->b2 * this->c3 * this->d0 - this->b3 * this->c2 * this->d0 - this->b2 * this->c0 * this->d3 - this->b0 * this->c3 * this->d2);
+C4Vector C44Matrix::Col1() const {
-    float b1 = this->a2 * this->c3 * this->d0 + this->a3 * this->c0 * this->d2 + this->a0 * this->c2 * this->d3 - this->a3 * this->c2 * this->d0 - this->a2 * this->c0 * this->d3 - this->a0 * this->c3 * this->d2;
+    return { this->a1, this->b1, this->c1, this->d1 };
-    float b2 = -(this->a2 * this->b3 * this->d0 + this->a3 * this->b0 * this->d2 + this->a0 * this->b2 * this->d3 - this->a3 * this->b2 * this->d0 - this->a2 * this->b0 * this->d3 - this->a0 * this->b3 * this->d2);
+}
-    float b3 = this->a2 * this->b3 * this->c0 + this->a3 * this->b0 * this->c2 + this->a0 * this->b2 * this->c3 - this->a3 * this->b2 * this->c0 - this->a2 * this->b0 * this->c3 - this->a0 * this->b3 * this->c2;
 
-    float c0 = this->b1 * this->c3 * this->d0 + this->b3 * this->c0 * this->d1 + this->b0 * this->c1 * this->d3 - this->b3 * this->c1 * this->d0 - this->b1 * this->c0 * this->d3 - this->b0 * this->c3 * this->d1;
+C4Vector C44Matrix::Col2() const {
-    float c1 = -(this->a1 * this->c3 * this->d0 + this->a3 * this->c0 * this->d1 + this->a0 * this->c1 * this->d3 - this->a3 * this->c1 * this->d0 - this->a1 * this->c0 * this->d3 - this->a0 * this->c3 * this->d1);
+    return { this->a2, this->b2, this->c2, this->d2 };
-    float c2 = this->a1 * this->b3 * this->d0 + this->a3 * this->b0 * this->d1 + this->a0 * this->b1 * this->d3 - this->a3 * this->b1 * this->d0 - this->a1 * this->b0 * this->d3 - this->a0 * this->b3 * this->d1;
+}
-    float c3 = -(this->a1 * this->b3 * this->c0 + this->a3 * this->b0 * this->c1 + this->a0 * this->b1 * this->c3 - this->a3 * this->b1 * this->c0 - this->a1 * this->b0 * this->c3 - this->a0 * this->b3 * this->c1);
 
-    float d0 = -(this->b2 * this->c0 * this->d1 + this->b0 * this->c1 * this->d2 + this->b1 * this->c2 * this->d0 - this->b2 * this->c1 * this->d0 - this->b1 * this->c0 * this->d2 - this->b0 * this->c2 * this->d1);
+C4Vector C44Matrix::Col3() const {
-    float d1 = this->a1 * this->c2 * this->d0 + this->a2 * this->c0 * this->d1 + this->a0 * this->c1 * this->d2 - this->a2 * this->c1 * this->d0 - this->a1 * this->c0 * this->d2 - this->a0 * this->c2 * this->d1;
+    return { this->a3, this->b3, this->c3, this->d3 };
-    float d2 = -(this->a1 * this->b2 * this->d0 + this->a2 * this->b0 * this->d1 + this->a0 * this->b1 * this->d2 - this->a2 * this->b1 * this->d0 - this->a1 * this->b0 * this->d2 - this->a0 * this->b2 * this->d1);
+}
-    float d3 = this->a0 * this->b1 * this->c2 + this->a1 * this->b2 * this->c0 + this->a2 * this->b0 * this->c1 - this->b1 * this->a2 * this->c0 - this->c2 * this->b0 * this->a1 - this->a0 * this->b2 * this->c1;
 
-    return { a0, a1, a2, a3, b0, b1, b2, b3, c0, c1, c2, c3, d0, d1, d2, d3 };
+float C44Matrix::Determinant() const {
+    return (this->b1 * this->c2 * this->d3 + this->c3 * this->b2 * this->d1 + this->b3 * this->c1 * this->d2 - this->c2 * this->b3 * this->d1 - this->d3 * (this->c1 * this->b2) - this->b1 * this->c3 * this->d2) * this->a0 - (this->c2 * this->b0 * this->d3 + this->c3 * this->b2 * this->d0 + this->b3 * this->c0 * this->d2 - this->b3 * this->c2 * this->d0 - this->d3 * (this->c0 * this->b2) - this->b0 * this->c3 * this->d2) * this->a1 + (this->c1 * this->b0 * this->d3 + this->c3 * this->b1 * this->d0 + this->b3 * this->c0 * this->d1 - this->b3 * this->c1 * this->d0 - this->d3 * (this->c0 * this->b1) - this->b0 * this->c3 * this->d1) * this->a2 - (this->c1 * this->b0 * this->d2 + this->c2 * this->b1 * this->d0 + this->b2 * this->c0 * this->d1 - this->b2 * this->c1 * this->d0 - this->d2 * (this->c0 * this->b1) - this->b0 * this->c2 * this->d1) * this->a3;
+}
+
+void C44Matrix::Identity() {
+    this->a0 = 1.0f;
+    this->a1 = 0.0f;
+    this->a2 = 0.0f;
+    this->a3 = 0.0f;
+
+    this->b0 = 0.0f;
+    this->b1 = 1.0f;
+    this->b2 = 0.0f;
+    this->b3 = 0.0f;
+
+    this->c0 = 0.0f;
+    this->c1 = 0.0f;
+    this->c2 = 1.0f;
+    this->c3 = 0.0f;
+
+    this->d0 = 0.0f;
+    this->d1 = 0.0f;
+    this->d2 = 0.0f;
+    this->d3 = 1.0f;
 }
 
 C44Matrix C44Matrix::Inverse() const {
@@ -658,34 +476,233 @@ C44Matrix C44Matrix::Inverse(float det) const {
     return this->Adjoint() * (1.0f / det);
 }
 
-C44Matrix C44Matrix::AffineInverse(const C3Vector& v) const {
+void C44Matrix::Rotate(const C4Quaternion& rotation) {
-    C33Matrix rotationScale(*this);
+    *this = C44Matrix(rotation) * (*this);
-    C3Vector s = { 1.0f / v.x, 1.0f / v.y, 1.0f / v.z };
-    rotationScale.Scale(s);
-    C44Matrix matrix(rotationScale.Transpose());
-    matrix.Scale(s);
-    C3Vector move(-this->d0, -this->d1, -this->d2);
-    matrix.Translate(move);
-    return matrix;
 }
 
-C44Matrix C44Matrix::AffineInverse(float a) const {
+void C44Matrix::Rotate(float angle, const C3Vector& axis, bool unit) {
-    if (CMath::fequalz(a, 1.0f, 0.00000095367432f)) {
+    *this = C44Matrix::Rotation(angle, axis, unit) * (*this);
-        return this->AffineInverse();
 }
 
-    C44Matrix matrix(C33Matrix(*this).Transpose());
+void C44Matrix::RotateAroundZ(float angle) {
-    matrix.Scale(1.0f / (a * a));
+    *this = C44Matrix::RotationAroundZ(angle) * (*this);
-    C3Vector move(-this->d0, -this->d1, -this->d2);
-    matrix.Translate(move);
-    return matrix;
 }
 
-C44Matrix C44Matrix::AffineInverse() const {
+C4Vector C44Matrix::Row0() const {
-    C44Matrix matrix(C33Matrix(*this).Transpose());
+    return { this->a0, this->a1, this->a2, this->a3 };
-    C3Vector move(-this->d0, -this->d1, -this->d2);
+}
-    matrix.Translate(move);
+
-    return matrix;
+const C3Vector* C44Matrix::Row0AsVec3() const {
+    return reinterpret_cast<const C3Vector*>(&this->a0);
+}
+
+C4Vector C44Matrix::Row1() const {
+    return { this->b0, this->b1, this->b2, this->b3 };
+}
+
+const C3Vector* C44Matrix::Row1AsVec3() const {
+    return reinterpret_cast<const C3Vector*>(&this->b0);
+}
+
+C4Vector C44Matrix::Row2() const {
+    return { this->c0, this->c1, this->c2, this->c3 };
+}
+
+const C3Vector* C44Matrix::Row2AsVec3() const {
+    return reinterpret_cast<const C3Vector*>(&this->c0);
+}
+
+C4Vector C44Matrix::Row3() const {
+    return { this->d0, this->d1, this->d2, this->d3 };
+}
+
+const C3Vector* C44Matrix::Row3AsVec3() const {
+    return reinterpret_cast<const C3Vector*>(&this->d0);
+}
+
+void C44Matrix::Scale(const C3Vector& scale) {
+    this->a0 *= scale.x;
+    this->a1 *= scale.x;
+    this->a2 *= scale.x;
+
+    this->b0 *= scale.y;
+    this->b1 *= scale.y;
+    this->b2 *= scale.y;
+
+    this->c0 *= scale.z;
+    this->c1 *= scale.z;
+    this->c2 *= scale.z;
+}
+
+void C44Matrix::Scale(float scale) {
+    this->a0 *= scale;
+    this->a1 *= scale;
+    this->a2 *= scale;
+
+    this->b0 *= scale;
+    this->b1 *= scale;
+    this->b2 *= scale;
+
+    this->c0 *= scale;
+    this->c1 *= scale;
+    this->c2 *= scale;
+}
+
+float C44Matrix::Trace() {
+    return this->a0 + this->b1 + this->c2 + this->d3;
+}
+
+void C44Matrix::Translate(const C3Vector& move) {
+    this->d0 = this->a0 * move.x + this->b0 * move.y + this->c0 * move.z + this->d0;
+    this->d1 = this->a1 * move.x + this->b1 * move.y + this->c1 * move.z + this->d1;
+    this->d2 = this->a2 * move.x + this->b2 * move.y + this->c2 * move.z + this->d2;
+}
+
+C44Matrix C44Matrix::Transpose() const {
+    float a0 = this->a0;
+    float a1 = this->a1;
+    float a2 = this->a2;
+    float a3 = this->a3;
+
+    float b0 = this->b0;
+    float b1 = this->b1;
+    float b2 = this->b2;
+    float b3 = this->b3;
+
+    float c0 = this->c0;
+    float c1 = this->c1;
+    float c2 = this->c2;
+    float c3 = this->c3;
+
+    float d0 = this->d0;
+    float d1 = this->d1;
+    float d2 = this->d2;
+    float d3 = this->d3;
+
+    return { a0, b0, c0, d0, a1, b1, c1, d1, a2, b2, c2, d2, a3, b3, c3, d3 };
+}
+
+void C44Matrix::Zero() {
+    this->a0 = 0.0f;
+    this->a1 = 0.0f;
+    this->a2 = 0.0f;
+    this->a3 = 0.0f;
+
+    this->b0 = 0.0f;
+    this->b1 = 0.0f;
+    this->b2 = 0.0f;
+    this->b3 = 0.0f;
+
+    this->c0 = 0.0f;
+    this->c1 = 0.0f;
+    this->c2 = 0.0f;
+    this->c3 = 0.0f;
+
+    this->d0 = 0.0f;
+    this->d1 = 0.0f;
+    this->d2 = 0.0f;
+    this->d3 = 0.0f;
+}
+
+C44Matrix& C44Matrix::operator+=(const C44Matrix& a) {
+    this->a0 += a.a0;
+    this->a1 += a.a1;
+    this->a2 += a.a2;
+    this->a3 += a.a3;
+
+    this->b0 += a.b0;
+    this->b1 += a.b1;
+    this->b2 += a.b2;
+    this->b3 += a.b3;
+
+    this->c0 += a.c0;
+    this->c1 += a.c1;
+    this->c2 += a.c2;
+    this->c3 += a.c3;
+
+    this->d0 += a.d0;
+    this->d1 += a.d1;
+    this->d2 += a.d2;
+    this->d3 += a.d3;
+
+    return *this;
+}
+
+C44Matrix& C44Matrix::operator-=(const C44Matrix& a) {
+    this->a0 -= a.a0;
+    this->a1 -= a.a1;
+    this->a2 -= a.a2;
+    this->a3 -= a.a3;
+
+    this->b0 -= a.b0;
+    this->b1 -= a.b1;
+    this->b2 -= a.b2;
+    this->b3 -= a.b3;
+
+    this->c0 -= a.c0;
+    this->c1 -= a.c1;
+    this->c2 -= a.c2;
+    this->c3 -= a.c3;
+
+    this->d0 -= a.d0;
+    this->d1 -= a.d1;
+    this->d2 -= a.d2;
+    this->d3 -= a.d3;
+
+    return *this;
+}
+
+C44Matrix& C44Matrix::operator*=(float a) {
+    this->a0 *= a;
+    this->a1 *= a;
+    this->a2 *= a;
+    this->a3 *= a;
+
+    this->b0 *= a;
+    this->b1 *= a;
+    this->b2 *= a;
+    this->b3 *= a;
+
+    this->c0 *= a;
+    this->c1 *= a;
+    this->c2 *= a;
+    this->c3 *= a;
+
+    this->d0 *= a;
+    this->d1 *= a;
+    this->d2 *= a;
+    this->d3 *= a;
+
+    return *this;
+}
+
+C44Matrix& C44Matrix::operator*=(const C44Matrix& a) {
+    *this = *this * a;
+    return *this;
+}
+
+C44Matrix& C44Matrix::operator/=(float a) {
+    this->a0 /= a;
+    this->a1 /= a;
+    this->a2 /= a;
+    this->a3 /= a;
+
+    this->b0 /= a;
+    this->b1 /= a;
+    this->b2 /= a;
+    this->b3 /= a;
+
+    this->c0 /= a;
+    this->c1 /= a;
+    this->c2 /= a;
+    this->c3 /= a;
+
+    this->d0 /= a;
+    this->d1 /= a;
+    this->d2 /= a;
+    this->d3 /= a;
+
+    return *this;
 }
 
 C44Matrix operator+(const C44Matrix& l, const C44Matrix& r) {

tempest/matrix/C44Matrix.hpp

@@ -3,6 +3,8 @@
 
 #include <cstdint>
 
+#include "tempest/vector/C4Vector.hpp"
+
 class C3Vector;
 class C4Vector;
 class C33Matrix;
@@ -16,8 +18,8 @@ class C44Matrix {
 
     // Static functions
     static float Det(float a, float b, float c, float d, float e, float f, float g, float h, float i);
-    static C44Matrix RotationAroundZ(float angle);
     static C44Matrix Rotation(float angle, const C3Vector& axis, bool unit);
+    static C44Matrix RotationAroundZ(float angle);
 
     // Member variables
     float a0 = 1.0f;
@@ -40,35 +42,51 @@ class C44Matrix {
     // Member functions
     C44Matrix() = default;
     C44Matrix(const C4Vector& r0, const C4Vector& r1, const C4Vector& r2, const C4Vector& r3);
+    C44Matrix(float a0, float a1, float a2, float a3, float b0, float b1, float b2, float b3, float c0, float c1, float c2, float c3, float d0, float d1, float d2, float d3);
     explicit C44Matrix(const C33Matrix& m);
     explicit C44Matrix(const C4Quaternion& rotation);
-    C44Matrix(float a0, float a1, float a2, float a3, float b0, float b1, float b2, float b3, float c0, float c1, float c2, float c3, float d0, float d1, float d2, float d3);
     explicit C44Matrix(float a);
 
+    C44Matrix Adjoint() const;
+    C44Matrix AffineInverse() const;
+    C44Matrix AffineInverse(const C3Vector& v) const;
+    C44Matrix AffineInverse(float uniformScale) const;
+    C44Matrix Cofactors() const;
+
+    C4Vector Col0() const;
+    C4Vector Col1() const;
+    C4Vector Col2() const;
+    C4Vector Col3() const;
+
+    float Determinant() const;
+    void Identity();
+    C44Matrix Inverse() const;
+    C44Matrix Inverse(float det) const;
+    void Rotate(const C4Quaternion& rotation);
+    void Rotate(float angle, const C3Vector& axis, bool unit);
+    void RotateAroundZ(float angle);
+
+    C4Vector Row0() const;
+    const C3Vector* Row0AsVec3() const;
+    C4Vector Row1() const;
+    const C3Vector* Row1AsVec3() const;
+    C4Vector Row2() const;
+    const C3Vector* Row2AsVec3() const;
+    C4Vector Row3() const;
+    const C3Vector* Row3AsVec3() const;
+
+    void Scale(const C3Vector& scale);
+    void Scale(float scale);
+    float Trace();
+    void Translate(const C3Vector& move);
+    C44Matrix Transpose() const;
+    void Zero();
+
     C44Matrix& operator+=(const C44Matrix& a);
     C44Matrix& operator-=(const C44Matrix& a);
     C44Matrix& operator*=(float a);
     C44Matrix& operator*=(const C44Matrix& a);
     C44Matrix& operator/=(float a);
-
-    void Zero();
-    void Identity();
-    float Trace();
-    void Translate(const C3Vector& move);
-    void Scale(float scale);
-    void Scale(const C3Vector& scale);
-    void RotateAroundZ(float angle);
-    void Rotate(const C4Quaternion& rotation);
-    void Rotate(float angle, const C3Vector& axis, bool unit);
-    C44Matrix Transpose() const;
-    float Determinant() const;
-    C44Matrix Cofactors() const;
-    C44Matrix Adjoint() const;
-    C44Matrix Inverse() const;
-    C44Matrix Inverse(float det) const;
-    C44Matrix AffineInverse(const C3Vector& v) const;
-    C44Matrix AffineInverse(float a) const;
-    C44Matrix AffineInverse() const;
 };
 
 C44Matrix operator+(const C44Matrix& l, const C44Matrix& r);

tempest/random/CRandom.cpp

@@ -0,0 +1,54 @@
+#include "tempest/random/CRandom.hpp"
+#include "tempest/math/CMath.hpp"
+
+const uint32_t gnoise32[] = {
+    0x9927148E, 0x08C7AAFD, 0x1F3EE6D5, 0xDA55BBF6,
+    0x6A4AA075, 0xFF97BDE8, 0x9FBC9BDE, 0x46A18A81,
+    0x63E30B6E, 0x5D6C7A76, 0xCA69D388, 0x25B947C3,
+    0x3FA2AB83, 0xBA7C41A6, 0x0195ACE5, 0xC109CF7E,
+    0x717062D9, 0x0205DB8D, 0x54EF8724, 0x3037D4C6,
+    0x7BCB1BD0, 0xECD8E4B8, 0xDCADCE49, 0xC494A913,
+    0x0DAE398F, 0x0EDD5218, 0x85F5FA78, 0x6DAFD258,
+    0x3B53B2A4, 0xBE50A551, 0x11F42DFC, 0xF1169848,
+    0x663DDF86, 0x2F2E445E, 0x176B0736, 0xB64C298B,
+    0xE75F89E2, 0xE121A7CD, 0xED65C94D, 0x239CEEFE,
+    0x04B77D33, 0x402A9A9E, 0xF35B10B3, 0x921C7782,
+    0x571E4E20, 0x8C067222, 0xFB732C67, 0xBF0AC259,
+    0x0CF95C79, 0x68121A28, 0x42193474, 0xF884C0B1,
+    0x9D15F038, 0x6F3AF260, 0x91EB90B4, 0x61357F1D,
+    0x5603325A, 0x932BC5A3, 0x434B0F80, 0x3CE0A8F7,
+    0x2664D196, 0x4FCC45D7, 0xB5E9B0C8, 0xEA31D600,
+};
+
+uint32_t CRandom::dice(uint32_t sides, CRndSeed& seed) {
+    return CMath::mulhwu(sides, CRandom::uint32(seed));
+}
+
+uint32_t CRandom::uint32(CRndSeed& seed) {
+    auto acc = seed.rndacc;
+    auto vls = seed.rndvls;
+
+    uint8_t b1 = (vls >> 24) & 0xFF;
+    uint8_t b2 = (vls >> 16) & 0xFF;
+    uint8_t b3 = (vls >>  8) & 0xFF;
+    uint8_t b4 = (vls >>  0) & 0xFF;
+
+    int32_t r1 = b1 - 0x04 + (b1 - 0x04 < 0 ? 0xBC : 0);
+    int32_t r2 = b2 - 0x0C + (b2 - 0x0C < 0 ? 0xD4 : 0);
+    int32_t r3 = b3 - 0x18 + (b3 - 0x18 < 0 ? 0xEC : 0);
+    int32_t r4 = b4 - 0x1C + (b4 - 0x1C < 0 ? 0xF4 : 0);
+
+    auto noise = reinterpret_cast<const char*>(gnoise32);
+
+    auto n1 = *reinterpret_cast<const uint32_t*>(&noise[r1]);
+    auto n2 = *reinterpret_cast<const uint32_t*>(&noise[r2]);
+    auto n3 = *reinterpret_cast<const uint32_t*>(&noise[r3]);
+    auto n4 = *reinterpret_cast<const uint32_t*>(&noise[r4]);
+
+    acc = (n4 ^ CMath::rotl3(n3) ^ CMath::rotl2(n2) ^ CMath::rotl1(n1)) + acc;
+
+    seed.rndacc = acc;
+    seed.rndvls = r4 | (r3 << 8) | ((r2 | (r1 << 8)) << 16);
+
+    return acc;
+}

tempest/random/CRandom.hpp

@@ -0,0 +1,16 @@
+#ifndef TEMPEST_RANDOM_C_RANDOM_HPP
+#define TEMPEST_RANDOM_C_RANDOM_HPP
+
+#include "tempest/random/CRndSeed.hpp"
+#include <cstdint>
+
+extern const uint32_t gnoise32[];
+
+class CRandom {
+    public:
+        // Static functions
+        static uint32_t dice(uint32_t sides, CRndSeed& seed);
+        static uint32_t uint32(CRndSeed& seed);
+};
+
+#endif

tempest/random/CRndSeed.cpp

@@ -0,0 +1,10 @@
+#include "tempest/random/CRndSeed.hpp"
+
+CRndSeed::CRndSeed(uint32_t seed) {
+    this->SetSeed(seed);
+}
+
+void CRndSeed::SetSeed(uint32_t seed) {
+    this->rndacc = seed;
+    this->rndvls = (4 * (seed % 0x3D)) | ((seed % 0x3B) << 10) | ((seed % 0x35) << 18) | ((seed % 0x2F) << 26);
+}

tempest/random/CRndSeed.hpp

@@ -0,0 +1,20 @@
+#ifndef TEMPEST_RANDOM_C_RND_SEED_HPP
+#define TEMPEST_RANDOM_C_RND_SEED_HPP
+
+#include <cstdint>
+
+class CRndSeed {
+    friend class CRandom;
+
+    public:
+        // Public member functions
+        CRndSeed(uint32_t seed);
+        void SetSeed(uint32_t seed);
+
+    protected:
+        // Protected member variables
+        uint32_t rndacc;
+        uint32_t rndvls;
+};
+
+#endif

tempest/vector/C2iVector.hpp

@@ -8,6 +8,17 @@ class C2iVector {
     // Member variables
     int32_t x;
     int32_t y;
+
+    // Member functions
+    C2iVector()
+        : C2iVector(0) {};
+    C2iVector(int32_t a)
+        : C2iVector(a, a) {};
+    C2iVector(const C2iVector& v)
+        : C2iVector(v.x, v.y) {};
+    C2iVector(int32_t x, int32_t y)
+        : x(x)
+        , y(y) {};
 };
 
 #endif

tempest/vector/C3Vector.cpp

@@ -114,7 +114,7 @@ C3Vector& C3Vector::operator/=(const C3Vector& a) {
     return *this;
 }
 
-C3Vector C3Vector::operator-() {
+C3Vector C3Vector::operator-() const {
     return { -this->x, -this->y, -this->z };
 }
 
@@ -234,10 +234,18 @@ C3Vector operator*(float l, const C3Vector& r) {
     return { l * r.x, l * r.y, l * r.z };
 }
 
+C3Vector operator*(const C3Vector& l, const C33Matrix& r) {
+    float x = l.x * r.a0 + l.y * r.b0 + l.z * r.c0;
+    float y = l.x * r.a1 + l.y * r.b1 + l.z * r.c1;
+    float z = l.x * r.a2 + l.y * r.b2 + l.z * r.c2;
+
+    return { x, y, z };
+}
+
 C3Vector operator*(const C3Vector& l, const C44Matrix& r) {
-    float x = r.c0 * l.z + r.b0 * l.y + r.a0 * l.x + r.d0;
+    float x = l.x * r.a0 + l.y * r.b0 + l.z * r.c0 + r.d0;
-    float y = r.c1 * l.z + r.b1 * l.y + r.a1 * l.x + r.d1;
+    float y = l.x * r.a1 + l.y * r.b1 + l.z * r.c1 + r.d1;
-    float z = r.c2 * l.z + r.b2 * l.y + r.a2 * l.x + r.d2;
+    float z = l.x * r.a2 + l.y * r.b2 + l.z * r.c2 + r.d2;
 
     return { x, y, z };
 }

tempest/vector/C3Vector.hpp

@@ -4,6 +4,7 @@
 #include "tempest/vector/C2Vector.hpp"
 #include "tempest/vector/CImVector.hpp"
 
+class C33Matrix;
 class C44Matrix;
 
 class C3Vector {
@@ -63,7 +64,7 @@ class C3Vector {
     C3Vector& operator/=(float a);
     C3Vector& operator/=(const C3Vector& a);
 
-    C3Vector operator-();
+    C3Vector operator-() const;
 
     float& operator[](uint32_t sub);
     const float& operator[](uint32_t sub) const;
@@ -85,6 +86,7 @@ C3Vector operator-(const C3Vector& l, const C3Vector& r);
 
 C3Vector operator*(const C3Vector& l, float r);
 C3Vector operator*(float l, const C3Vector& r);
+C3Vector operator*(const C3Vector& l, const C33Matrix& r);
 C3Vector operator*(const C3Vector& l, const C44Matrix& r);
 
 bool operator==(const C3Vector& l, const C3Vector& r);

tempest/vector/C3iVector.hpp

@@ -0,0 +1,26 @@
+#ifndef TEMPEST_VECTOR_C_3IVECTOR_HPP
+#define TEMPEST_VECTOR_C_3IVECTOR_HPP
+
+#include <cstdint>
+
+class C3iVector {
+    public:
+    // Member variables
+    int32_t x;
+    int32_t y;
+    int32_t z;
+
+    // Member functions
+    C3iVector()
+        : C3iVector(0) {};
+    C3iVector(int32_t a)
+        : C3iVector(a, a, a) {};
+    C3iVector(const C3iVector& v)
+        : C3iVector(v.x, v.y, v.z) {};
+    C3iVector(int32_t x, int32_t y, int32_t z)
+        : x(x)
+        , y(y)
+        , z(z) {};
+};
+
+#endif

test/Random.cpp

@@ -0,0 +1,37 @@
+#include "tempest/Random.hpp"
+#include "test/Test.hpp"
+
+TEST_CASE("CRandom::dice", "[random]") {
+    SECTION("generates numbers within range") {
+        auto seed = CRndSeed(0xC0D1E2F3);
+
+        for (int32_t i = 0; i < 200; i++) {
+            auto result1 = CRandom::dice(0, seed);
+            REQUIRE(result1 == 0);
+
+            auto result2 = CRandom::dice(1, seed);
+            REQUIRE(result2 == 0);
+
+            auto result3 = CRandom::dice(2, seed);
+            REQUIRE(result3 < 2);
+
+            auto result4 = CRandom::dice(10, seed);
+            REQUIRE(result4 < 10);
+
+            auto result5 = CRandom::dice(100, seed);
+            REQUIRE(result5 < 100);
+        }
+    }
+}
+
+TEST_CASE("CRandom::uint32", "[random]") {
+    SECTION("generates expected numbers for given seed") {
+        auto seed = CRndSeed(0x294823);
+
+        auto result1 = CRandom::uint32(seed);
+        REQUIRE(result1 == 0xC0D0F1AE);
+
+        auto result2 = CRandom::uint32(seed);
+        REQUIRE(result2 == 0x5CA56410);
+    }
+}

test/Vector.cpp

@@ -47,6 +47,16 @@ TEST_CASE("C3Vector", "[vector]") {
     }
 }
 
+TEST_CASE("C3Vector::operator-", "[vector]") {
+    SECTION("returns vector with negated values") {
+        auto vector = C3Vector(1.0f, 2.0f, 3.0f);
+        auto negated = -vector;
+        CHECK(negated.x == -1.0f);
+        CHECK(negated.y == -2.0f);
+        CHECK(negated.z == -3.0f);
+    }
+}
+
 TEST_CASE("C3Vector::operator*=", "[vector]") {
     SECTION("multiplies by 2.0f") {
         auto vector = C3Vector(1.0f, 2.0f, 3.0f);