Add CTransform::quickInverse, toCStyleMatrix

- Fixes m128d dot3
- Simplifies negate operations

CMakeLists.txt

@@ -37,9 +37,11 @@ add_library(zeus
     include/zeus/zeus.hpp
     include/zeus/CVector2i.hpp
     include/zeus/CVector2f.hpp
+    include/zeus/CVector2d.hpp
     include/zeus/CVector3f.hpp
     include/zeus/CVector3d.hpp
     include/zeus/CVector4f.hpp
+    include/zeus/CVector4d.hpp
     include/zeus/CRectangle.hpp
     include/zeus/CMatrix4f.hpp
     include/zeus/CFrustum.hpp
@@ -58,10 +60,5 @@ add_library(zeus
     include/zeus/simd/parallelism_v2_simd.hpp)
 
 target_include_directories(zeus PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>)
-if(TARGET athena-core)
-  target_link_libraries(zeus PUBLIC athena-core)
-  target_compile_definitions(zeus PUBLIC ZE_ATHENA_TYPES=1)
-endif()
-
 add_subdirectory(test)
 

include/zeus/CAABox.hpp

@@ -9,10 +9,6 @@
 #include "zeus/CVector3f.hpp"
 #include "zeus/Math.hpp"
 
-#if ZE_ATHENA_TYPES
-#include <athena/IStreamReader.hpp>
-#endif
-
 namespace zeus {
 class CAABox {
 public:
@@ -40,21 +36,6 @@ public:
   constexpr CAABox(float minX, float minY, float minZ, float maxX, float maxY, float maxZ)
   : min(minX, minY, minZ), max(maxX, maxY, maxZ) {}
 
-#if ZE_ATHENA_TYPES
-
-  void readBoundingBoxBig(athena::io::IStreamReader& in) {
-    min.readBig(in);
-    max.readBig(in);
-  }
-
-  [[nodiscard]] static CAABox ReadBoundingBoxBig(athena::io::IStreamReader& in) {
-    CAABox ret;
-    ret.readBoundingBoxBig(in);
-    return ret;
-  }
-
-#endif
-
   [[nodiscard]] bool intersects(const CAABox& other) const {
     const auto mmax = max >= other.min;
     const auto mmin = min <= other.max;

include/zeus/CColor.hpp

@@ -8,20 +8,15 @@
 #include "zeus/Global.hpp"
 #include "zeus/Math.hpp"
 
-#if ZE_ATHENA_TYPES
-#include <athena/FileReader.hpp>
-#include <athena/FileWriter.hpp>
-#endif
-
 #undef min
 #undef max
 
-#if BYTE_ORDER == __ORDER_LITTLE_ENDIAN__
+#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
 #define COLOR(rgba)                                                                                                    \
-  (unsigned)(((rgba)&0x000000FF) << 24 | ((rgba)&0x0000FF00) << 8 | ((rgba)&0x00FF0000) >> 8 |                         \
-             ((rgba)&0xFF000000) >> 24)
+  (zeus::Comp32)(((rgba)&0x000000FF) << 24 | ((rgba)&0x0000FF00) << 8 | ((rgba)&0x00FF0000) >> 8 |                     \
+                 ((rgba)&0xFF000000) >> 24)
 #else
-#define COLOR(rgba) rgba
+#define COLOR(rgba) (zeus::Comp32)(rgba)
 #endif
 
 namespace zeus {
@@ -48,21 +43,12 @@ public:
 
   constexpr CColor(float r, float g, float b, float a = 1.0f) : mSimd(r, g, b, a) {}
 
-#if ZE_ATHENA_TYPES
+  constexpr CColor(Comp32 rgba)
+  : mSimd(((COLOR(rgba) >> 0) & 0xff) * OneOver255, ((COLOR(rgba) >> 8) & 0xff) * OneOver255,
+          ((COLOR(rgba) >> 16) & 0xff) * OneOver255, ((COLOR(rgba) >> 24) & 0xff) * OneOver255) {}
 
-  constexpr CColor(const atVec4f& vec) : mSimd(vec.simd) {}
-
-#endif
-
-  constexpr CColor(Comp32 rgba) : mSimd(((COLOR(rgba) >> 0) & 0xff) * OneOver255,
-                                        ((COLOR(rgba) >> 8) & 0xff) * OneOver255,
-                                        ((COLOR(rgba) >> 16) & 0xff) * OneOver255,
-                                        ((COLOR(rgba) >> 24) & 0xff) * OneOver255) {}
-
-  constexpr CColor(const Comp8* rgba) : mSimd(rgba[0] * OneOver255,
-                                              rgba[1] * OneOver255,
-                                              rgba[2] * OneOver255,
-                                              rgba[3] * OneOver255) {}
+  constexpr CColor(const Comp8* rgba)
+  : mSimd(rgba[0] * OneOver255, rgba[1] * OneOver255, rgba[2] * OneOver255, rgba[3] * OneOver255) {}
 
   constexpr CColor(const CVector4f& other) : mSimd(other.mSimd) {}
 
@@ -259,6 +245,28 @@ public:
     ao = Comp8(a() * 255);
   }
 
+  Comp32 toRGBA() const {
+    RGBA32 ret;
+    ret.r = r() * 255;
+    ret.g = g() * 255;
+    ret.b = b() * 255;
+    ret.a = a() * 255;
+    return COLOR(ret.rgba);
+  }
+
+  [[nodiscard]] unsigned short toRGB5A3() const {
+    Comp8 r;
+    Comp8 g;
+    Comp8 b;
+    Comp8 a;
+    toRGBA8(r, g, b, a);
+    if (a == 255) {
+      return static_cast<unsigned short>((r & 0xf8) << 7 | (g & 0xf8) << 2 | b >> 3 | 0x8000);
+    }
+
+    return static_cast<unsigned short>((r & 0xf0) << 4 | (g & 0xf0) | b >> 4 | (a & 0xe0) << 7);
+  }
+
   /**
    * @brief Assigns rgba from hsv
    * @param h[0-1] The hue percentage of the color.
@@ -351,4 +359,4 @@ struct hash<zeus::CColor> {
     return ret;
   }
 };
-}
+} // namespace std

include/zeus/CEulerAngles.hpp

@@ -9,8 +9,9 @@ class CTransform;
 class CEulerAngles : public CVector3f {
 public:
   constexpr CEulerAngles(float x, float y, float z) : CVector3f(x, y, z) {}
-  CEulerAngles(const CQuaternion& quat);
-  CEulerAngles(const CTransform& xf);
+  explicit CEulerAngles(const CQuaternion& quat);
+  explicit CEulerAngles(const CTransform& xf);
+  explicit CEulerAngles(const CVector3f& vec) : CVector3f(vec) {}
 };
 
 } // namespace zeus

include/zeus/CMatrix3f.hpp

@@ -23,13 +23,13 @@ public:
   constexpr CMatrix3f(float m00, float m01, float m02, float m10, float m11, float m12, float m20, float m21, float m22)
   : m{{{m00, m10, m20}, {m01, m11, m21}, {m02, m12, m22}}} {}
 
-  constexpr CMatrix3f(const CVector3f& scaleVec) {
+  CMatrix3f(const CVector3f& scaleVec) {
     m[0][0] = scaleVec[0];
     m[1][1] = scaleVec[1];
     m[2][2] = scaleVec[2];
   }
 
-  constexpr CMatrix3f(float scale) : CMatrix3f(CVector3f(scale)) {}
+  CMatrix3f(float scale) : CMatrix3f(CVector3f(scale)) {}
 
   constexpr CMatrix3f(const CVector3f& r0, const CVector3f& r1, const CVector3f& r2) : m{{r0, r1, r2}} {}
 

include/zeus/CMatrix4f.hpp

@@ -51,6 +51,14 @@ public:
 
   constexpr CMatrix4f& operator=(const CMatrix4f& other) = default;
 
+  [[nodiscard]] bool operator==(const CMatrix4f& rhs) const {
+    return m[0] == rhs.m[0] && m[1] == rhs.m[1] && m[2] == rhs.m[2] && m[3] == rhs.m[3];
+  }
+
+  [[nodiscard]] bool operator!=(const CMatrix4f& rhs) const {
+    return m[0] != rhs.m[0] || m[1] != rhs.m[1] || m[2] != rhs.m[2] || m[3] != rhs.m[3];
+  }
+
   [[nodiscard]] CVector4f operator*(const CVector4f& other) const {
     return m[0].mSimd * other.mSimd.shuffle<0, 0, 0, 0>() + m[1].mSimd * other.mSimd.shuffle<1, 1, 1, 1>() +
            m[2].mSimd * other.mSimd.shuffle<2, 2, 2, 2>() + m[3].mSimd * other.mSimd.shuffle<3, 3, 3, 3>();

include/zeus/COBBox.hpp

@@ -4,28 +4,9 @@
 #include "zeus/CTransform.hpp"
 #include "zeus/CVector3f.hpp"
 
-#if ZE_ATHENA_TYPES
-#include <athena/IStreamReader.hpp>
-#endif
-
 namespace zeus {
 class COBBox {
 public:
-#if ZE_ATHENA_TYPES
-
-  void readBig(athena::io::IStreamReader& in) {
-    transform.read34RowMajor(in);
-    extents.readBig(in);
-  }
-
-  [[nodiscard]] static COBBox ReadBig(athena::io::IStreamReader& in) {
-    COBBox out;
-    out.readBig(in);
-    return out;
-  }
-
-#endif
-
   CTransform transform;
   CVector3f extents;
 

include/zeus/CPlane.hpp

@@ -38,7 +38,7 @@ public:
     mSimd[3] = nd * mag;
   }
 
-  [[nodiscard]] float pointToPlaneDist(const CVector3f& pos) const { return pos.dot(normal()) - d(); }
+  [[nodiscard]] float pointToPlaneDist(const CVector3f& pos) const { return normal().dot(pos) - d(); }
 
   [[nodiscard]] bool rayPlaneIntersection(const CVector3f& from, const CVector3f& to, CVector3f& point) const;
 

include/zeus/CQuaternion.hpp

@@ -10,10 +10,6 @@
 #include "zeus/CVector4f.hpp"
 #include "zeus/Global.hpp"
 
-#if ZE_ATHENA_TYPES
-#include <athena/IStreamReader.hpp>
-#endif
-
 namespace zeus {
 
 class CNUQuaternion;
@@ -32,25 +28,6 @@ public:
   template <typename T>
   constexpr CQuaternion(const simd<T>& s) : mSimd(s) {}
 
-#if ZE_ATHENA_TYPES
-
-  void readBig(athena::io::IStreamReader& input) {
-    simd_floats f;
-    f[0] = input.readFloatBig();
-    f[1] = input.readFloatBig();
-    f[2] = input.readFloatBig();
-    f[3] = input.readFloatBig();
-    mSimd.copy_from(f);
-  }
-
-  constexpr CQuaternion(const atVec4f& vec) : mSimd(vec.simd) {}
-
-  operator atVec4f&() { return *reinterpret_cast<atVec4f*>(this); }
-
-  operator const atVec4f&() const { return *reinterpret_cast<const atVec4f*>(this); }
-
-#endif
-
   CQuaternion(const CMatrix3f& mat);
 
   CQuaternion(const CVector3f& vec) { fromVector3f(vec); }
@@ -75,9 +52,16 @@ public:
 
   [[nodiscard]] CQuaternion operator-(const CQuaternion& q) const { return mSimd - q.mSimd; }
 
-  [[nodiscard]] CQuaternion operator*(const CQuaternion& q) const;
+  [[nodiscard]] CQuaternion operator*(const CQuaternion& q) const {
+    return CQuaternion(w() * q.w() - CVector3f(x(), y(), z()).dot({q.x(), q.y(), q.z()}),
+                       y() * q.z() - z() * q.y() + w() * q.x() + x() * q.w(),
+                       z() * q.x() - x() * q.z() + w() * q.y() + y() * q.w(),
+                       x() * q.y() - y() * q.x() + w() * q.z() + z() * q.w());
+  }
 
-  [[nodiscard]] CQuaternion operator/(const CQuaternion& q) const;
+  [[nodiscard]] CQuaternion operator/(const CQuaternion& q) const {
+    return *this * q.inverse();
+  }
 
   [[nodiscard]] CQuaternion operator*(float scale) const { return mSimd * simd<float>(scale); }
 
@@ -234,9 +218,9 @@ public:
  */
 class CNUQuaternion {
 public:
-  CNUQuaternion() : mSimd(1.f, 0.f, 0.f, 0.f) {}
+  constexpr CNUQuaternion() : mSimd(1.f, 0.f, 0.f, 0.f) {}
 
-  CNUQuaternion(float wi, float xi, float yi, float zi) : mSimd(wi, xi, yi, zi) {}
+  constexpr CNUQuaternion(float wi, float xi, float yi, float zi) : mSimd(wi, xi, yi, zi) {}
 
   CNUQuaternion(float win, const zeus::CVector3f& vec) : mSimd(vec.mSimd.shuffle<0, 0, 1, 2>()) { w() = win; }
 
@@ -244,7 +228,7 @@ public:
 
   CNUQuaternion(const CMatrix3f& mtx) : CNUQuaternion(CQuaternion(mtx)) {}
 
-  CNUQuaternion(const simd<float>& s) : mSimd(s) {}
+  constexpr CNUQuaternion(const simd<float>& s) : mSimd(s) {}
 
   [[nodiscard]] static CNUQuaternion fromAxisAngle(const CUnitVector3f& axis, const CRelAngle& angle) {
     return CNUQuaternion(CQuaternion::fromAxisAngle(axis, angle));
@@ -264,7 +248,12 @@ public:
     return mSimd * simd<float>(magDiv);
   }
 
-  [[nodiscard]] CNUQuaternion operator*(const CNUQuaternion& q) const;
+  [[nodiscard]] CNUQuaternion operator*(const CNUQuaternion& q) const {
+    return CNUQuaternion(w() * q.w() - CVector3f(x(), y(), z()).dot({q.x(), q.y(), q.z()}),
+                         y() * q.z() - z() * q.y() + w() * q.x() + x() * q.w(),
+                         z() * q.x() - x() * q.z() + w() * q.y() + y() * q.w(),
+                         x() * q.y() - y() * q.x() + w() * q.z() + z() * q.w());
+  }
 
   [[nodiscard]] CNUQuaternion operator*(float f) const { return mSimd * simd<float>(f); }
 

include/zeus/CTransform.hpp

@@ -14,24 +14,7 @@ class CTransform {
 public:
   constexpr CTransform() : basis(false) {}
 
-  constexpr CTransform(const CMatrix3f& basis, const CVector3f& offset = {})
-  : basis(basis), origin(offset) {}
-
-#if ZE_ATHENA_TYPES
-
-  CTransform(const atVec4f* mtx)
-  : basis(mtx[0], mtx[1], mtx[2]), origin(mtx[0].simd[3], mtx[1].simd[3], mtx[2].simd[3]) {}
-
-  void read34RowMajor(athena::io::IStreamReader& r) {
-    atVec4f r0 = r.readVec4fBig();
-    atVec4f r1 = r.readVec4fBig();
-    atVec4f r2 = r.readVec4fBig();
-    basis = CMatrix3f(r0, r1, r2);
-    basis.transpose();
-    origin = CVector3f(r0.simd[3], r1.simd[3], r2.simd[3]);
-  }
-
-#endif
+  constexpr CTransform(const CMatrix3f& basis, const CVector3f& offset = {}) : basis(basis), origin(offset) {}
 
   /* Column constructor */
   constexpr CTransform(const CVector3f& c0, const CVector3f& c1, const CVector3f& c2, const CVector3f& c3)
@@ -52,6 +35,15 @@ public:
     return CTransform(inv, inv * -origin);
   }
 
+  [[nodiscard]] CTransform quickInverse() const {
+    return CTransform{basis.transposed(),
+                      CVector3f{
+                          basis[0][0] * -origin.x() - basis[0][1] * origin.y() - basis[0][2] * origin.z(),
+                          basis[1][0] * -origin.x() - basis[1][1] * origin.y() - basis[1][2] * origin.z(),
+                          basis[2][0] * -origin.x() - basis[2][1] * origin.y() - basis[2][2] * origin.z(),
+                      }};
+  }
+
   [[nodiscard]] static CTransform Translate(const CVector3f& position) { return {CMatrix3f(), position}; }
 
   [[nodiscard]] static CTransform Translate(float x, float y, float z) { return Translate({x, y, z}); }
@@ -194,6 +186,24 @@ public:
     return ret;
   }
 
+  /**
+   * Outputs the matrix to a C-style array (column-major, GX style)
+   */
+  void toCStyleMatrix(float mtx[3][4]) const {
+    mtx[0][0] = basis[0][0];
+    mtx[0][1] = basis[1][0];
+    mtx[0][2] = basis[2][0];
+    mtx[0][3] = origin.x();
+    mtx[1][0] = basis[0][1];
+    mtx[1][1] = basis[1][1];
+    mtx[1][2] = basis[2][1];
+    mtx[1][3] = origin.y();
+    mtx[2][0] = basis[0][2];
+    mtx[2][1] = basis[1][2];
+    mtx[2][2] = basis[2][2];
+    mtx[2][3] = origin.z();
+  }
+
   [[nodiscard]] CVector3f upVector() const { return basis.m[2]; }
 
   [[nodiscard]] CVector3f frontVector() const { return basis.m[1]; }

include/zeus/CVector2d.hpp

@@ -0,0 +1,191 @@
+#pragma once
+
+#include <cassert>
+
+#include "zeus/Global.hpp"
+#include "zeus/Math.hpp"
+
+namespace zeus {
+class CVector2d {
+public:
+  simd<double> mSimd;
+  constexpr CVector2d() : mSimd() {}
+
+  template <typename T>
+  constexpr CVector2d(const simd<T>& s) : mSimd(s) {}
+
+  explicit constexpr CVector2d(double xy) : mSimd(xy) {}
+
+  constexpr void assign(double x, double y) {
+    mSimd.set(x, y);
+  }
+
+  constexpr CVector2d(double x, double y) : mSimd(x, y, 0.f, 0.f) {}
+
+  [[nodiscard]] bool operator==(const CVector2d& rhs) const {
+    return mSimd[0] == rhs.mSimd[0] && mSimd[1] == rhs.mSimd[1];
+  }
+
+  [[nodiscard]] bool operator!=(const CVector2d& rhs) const {
+    return mSimd[0] != rhs.mSimd[0] || mSimd[1] != rhs.mSimd[1];
+  }
+
+  [[nodiscard]] bool operator<(const CVector2d& rhs) const {
+    return mSimd[0] < rhs.mSimd[0] && mSimd[1] < rhs.mSimd[1];
+  }
+
+  [[nodiscard]] bool operator<=(const CVector2d& rhs) const {
+    return mSimd[0] <= rhs.mSimd[0] && mSimd[1] <= rhs.mSimd[1];
+  }
+
+  [[nodiscard]] bool operator>(const CVector2d& rhs) const {
+    return mSimd[0] > rhs.mSimd[0] && mSimd[1] > rhs.mSimd[1];
+  }
+
+  [[nodiscard]] bool operator>=(const CVector2d& rhs) const {
+    return mSimd[0] >= rhs.mSimd[0] && mSimd[1] >= rhs.mSimd[1];
+  }
+
+  [[nodiscard]] CVector2d operator+(const CVector2d& rhs) const { return mSimd + rhs.mSimd; }
+
+  [[nodiscard]] CVector2d operator-(const CVector2d& rhs) const { return mSimd - rhs.mSimd; }
+
+  [[nodiscard]] CVector2d operator-() const { return -mSimd; }
+
+  [[nodiscard]] CVector2d operator*(const CVector2d& rhs) const { return mSimd * rhs.mSimd; }
+
+  [[nodiscard]] CVector2d operator/(const CVector2d& rhs) const { return mSimd / rhs.mSimd; }
+
+  [[nodiscard]] CVector2d operator+(double val) const { return mSimd + simd<double>(val); }
+
+  [[nodiscard]] CVector2d operator-(double val) const { return mSimd - simd<double>(val); }
+
+  [[nodiscard]] CVector2d operator*(double val) const { return mSimd * simd<double>(val); }
+
+  [[nodiscard]] CVector2d operator/(double val) const {
+    double ooval = 1.0 / val;
+    return mSimd * simd<double>(ooval);
+  }
+
+  const CVector2d& operator+=(const CVector2d& rhs) {
+    mSimd += rhs.mSimd;
+    return *this;
+  }
+
+  const CVector2d& operator-=(const CVector2d& rhs) {
+    mSimd -= rhs.mSimd;
+    return *this;
+  }
+
+  const CVector2d& operator*=(const CVector2d& rhs) {
+    mSimd *= rhs.mSimd;
+    return *this;
+  }
+
+  const CVector2d& operator/=(const CVector2d& rhs) {
+    mSimd /= rhs.mSimd;
+    return *this;
+  }
+
+  const CVector2d& operator+=(double rhs) {
+    mSimd += simd<double>(rhs);
+    return *this;
+  }
+
+  const CVector2d& operator-=(double rhs) {
+    mSimd -= simd<double>(rhs);
+    return *this;
+  }
+
+  const CVector2d& operator*=(double rhs) {
+    mSimd *= simd<double>(rhs);
+    return *this;
+  }
+
+  const CVector2d& operator/=(double rhs) {
+    double oorhs = 1.0 / rhs;
+    mSimd *= simd<double>(oorhs);
+    return *this;
+  }
+
+  void normalize() {
+    double mag = magnitude();
+    mag = 1.f / mag;
+    *this *= CVector2d(mag);
+  }
+
+  [[nodiscard]] CVector2d normalized() const {
+    double mag = magnitude();
+    mag = 1.f / mag;
+    return *this * mag;
+  }
+
+  [[nodiscard]] constexpr CVector2d perpendicularVector() const { return {-y(), x()}; }
+
+  [[nodiscard]] constexpr double cross(const CVector2d& rhs) const { return (x() * rhs.y()) - (y() * rhs.x()); }
+
+  [[nodiscard]] constexpr double dot(const CVector2d& rhs) const { return mSimd.dot2(rhs.mSimd); }
+
+  [[nodiscard]] constexpr double magSquared() const { return mSimd.dot2(mSimd); }
+
+  [[nodiscard]] double magnitude() const { return std::sqrt(magSquared()); }
+
+  constexpr void zeroOut() { mSimd = 0.f; }
+
+  constexpr void splat(double xy) { mSimd = xy; }
+
+  [[nodiscard]] static double getAngleDiff(const CVector2d& a, const CVector2d& b);
+
+  [[nodiscard]] static CVector2d lerp(const CVector2d& a, const CVector2d& b, double t) {
+    return zeus::simd<double>(1.f - t) * a.mSimd + b.mSimd * zeus::simd<double>(t);
+  }
+
+  [[nodiscard]] static CVector2d nlerp(const CVector2d& a, const CVector2d& b, double t) {
+    return lerp(a, b, t).normalized();
+  }
+
+  [[nodiscard]] static CVector2d slerp(const CVector2d& a, const CVector2d& b, double t);
+
+  [[nodiscard]] bool isNormalized() const { return std::fabs(1.0 - magSquared()) < 0.01; }
+
+  [[nodiscard]] bool canBeNormalized() const {
+    if (std::isinf(x()) || std::isinf(y()))
+      return false;
+    return std::fabs(x()) >= DBL_EPSILON || std::fabs(y()) >= DBL_EPSILON;
+  }
+
+  [[nodiscard]] bool isZero() const { return mSimd[0] == 0.f && mSimd[1] == 0.f; }
+
+  [[nodiscard]] bool isEqu(const CVector2d& other, double epsilon = FLT_EPSILON) const {
+    const CVector2d diffVec = other - *this;
+    return (diffVec.x() <= epsilon && diffVec.y() <= epsilon);
+  }
+
+  [[nodiscard]] simd<double>::reference operator[](size_t idx) {
+    assert(idx < 2);
+    return mSimd[idx];
+  }
+
+  [[nodiscard]] constexpr double operator[](size_t idx) const {
+    assert(idx < 2);
+    return mSimd[idx];
+  }
+
+  [[nodiscard]] constexpr double x() const { return mSimd[0]; }
+  [[nodiscard]] constexpr double y() const { return mSimd[1]; }
+
+  [[nodiscard]] simd<double>::reference x() { return mSimd[0]; }
+  [[nodiscard]] simd<double>::reference y() { return mSimd[1]; }
+};
+constexpr inline CVector2d skOne2d(1.0);
+constexpr inline CVector2d skNegOne2d(-1.0);
+constexpr inline CVector2d skZero2d(0.0);
+
+[[nodiscard]] inline CVector2d operator+(double lhs, const CVector2d& rhs) { return zeus::simd<double>(lhs) + rhs.mSimd; }
+
+[[nodiscard]] inline CVector2d operator-(double lhs, const CVector2d& rhs) { return zeus::simd<double>(lhs) - rhs.mSimd; }
+
+[[nodiscard]] inline CVector2d operator*(double lhs, const CVector2d& rhs) { return zeus::simd<double>(lhs) * rhs.mSimd; }
+
+[[nodiscard]] inline CVector2d operator/(double lhs, const CVector2d& rhs) { return zeus::simd<double>(lhs) / rhs.mSimd; }
+} // namespace zeus

include/zeus/CVector2f.hpp

@@ -5,10 +5,6 @@
 #include "zeus/Global.hpp"
 #include "zeus/Math.hpp"
 
-#if ZE_ATHENA_TYPES
-#include <athena/IStreamReader.hpp>
-#endif
-
 namespace zeus {
 class CVector2f {
 public:
@@ -18,29 +14,6 @@ public:
   template <typename T>
   constexpr CVector2f(const simd<T>& s) : mSimd(s) {}
 
-#if ZE_ATHENA_TYPES
-
-  constexpr CVector2f(const atVec2f& vec) : mSimd(vec.simd) {}
-
-  operator atVec2f&() { return *reinterpret_cast<atVec2f*>(this); }
-
-  operator const atVec2f&() const { return *reinterpret_cast<const atVec2f*>(this); }
-
-  void readBig(athena::io::IStreamReader& input) {
-    mSimd[0] = input.readFloatBig();
-    mSimd[1] = input.readFloatBig();
-    mSimd[2] = 0.0f;
-    mSimd[3] = 0.0f;
-  }
-
-  [[nodiscard]] static CVector2f ReadBig(athena::io::IStreamReader& input) {
-    CVector2f ret;
-    ret.readBig(input);
-    return ret;
-  }
-
-#endif
-
   explicit constexpr CVector2f(float xy) : mSimd(xy) {}
 
   constexpr void assign(float x, float y) {
@@ -131,7 +104,7 @@ public:
 
   const CVector2f& operator/=(float rhs) {
     float oorhs = 1.f / rhs;
-    mSimd /= simd<float>(oorhs);
+    mSimd *= simd<float>(oorhs);
     return *this;
   }
 
@@ -155,7 +128,7 @@ public:
 
   [[nodiscard]] constexpr float magSquared() const { return mSimd.dot2(mSimd); }
 
-  [[nodiscard]] constexpr float magnitude() const { return std::sqrt(magSquared()); }
+  [[nodiscard]] float magnitude() const { return std::sqrt(magSquared()); }
 
   constexpr void zeroOut() { mSimd = 0.f; }
 
@@ -173,9 +146,9 @@ public:
 
   [[nodiscard]] static CVector2f slerp(const CVector2f& a, const CVector2f& b, float t);
 
-  [[nodiscard]] constexpr bool isNormalized() const { return std::fabs(1.f - magSquared()) < 0.01f; }
+  [[nodiscard]] bool isNormalized() const { return std::fabs(1.f - magSquared()) < 0.01f; }
 
-  [[nodiscard]] constexpr bool canBeNormalized() const {
+  [[nodiscard]] bool canBeNormalized() const {
     if (std::isinf(x()) || std::isinf(y()))
       return false;
     return std::fabs(x()) >= FLT_EPSILON || std::fabs(y()) >= FLT_EPSILON;
@@ -188,7 +161,7 @@ public:
     return (diffVec.x() <= epsilon && diffVec.y() <= epsilon);
   }
 
-  [[nodiscard]] constexpr simd<float>::reference operator[](size_t idx) {
+  [[nodiscard]] simd<float>::reference operator[](size_t idx) {
     assert(idx < 2);
     return mSimd[idx];
   }
@@ -201,8 +174,8 @@ public:
   [[nodiscard]] constexpr float x() const { return mSimd[0]; }
   [[nodiscard]] constexpr float y() const { return mSimd[1]; }
 
-  [[nodiscard]] constexpr simd<float>::reference x() { return mSimd[0]; }
-  [[nodiscard]] constexpr simd<float>::reference y() { return mSimd[1]; }
+  [[nodiscard]] simd<float>::reference x() { return mSimd[0]; }
+  [[nodiscard]] simd<float>::reference y() { return mSimd[1]; }
 };
 constexpr inline CVector2f skOne2f(1.f);
 constexpr inline CVector2f skNegOne2f(-1.f);

include/zeus/CVector3d.hpp

@@ -5,10 +5,6 @@
 #include "zeus/CVector3f.hpp"
 #include "zeus/Global.hpp"
 
-#if ZE_ATHENA_TYPES
-#include <athena/Types.hpp>
-#endif
-
 namespace zeus {
 
 class CVector3d {
@@ -19,10 +15,6 @@ public:
   template <typename T>
   constexpr CVector3d(const simd<T>& s) : mSimd(s) {}
 
-#if ZE_ATHENA_TYPES
-  constexpr CVector3d(const atVec3d& vec) : mSimd(vec.simd) {}
-#endif
-
   explicit constexpr CVector3d(double xyz) : mSimd(xyz) {}
 
   CVector3d(const CVector3f& vec) : mSimd(vec.mSimd) {}

include/zeus/CVector3f.hpp

@@ -8,10 +8,6 @@
 #include "zeus/Global.hpp"
 #include "zeus/Math.hpp"
 
-#if ZE_ATHENA_TYPES
-#include <athena/IStreamReader.hpp>
-#endif
-
 namespace zeus {
 class CVector3d;
 class CRelAngle;
@@ -24,31 +20,6 @@ public:
   template <typename T>
   constexpr CVector3f(const simd<T>& s) : mSimd(s) {}
 
-#if ZE_ATHENA_TYPES
-
-  constexpr CVector3f(const atVec3f& vec) : mSimd(vec.simd) {}
-
-  operator atVec3f&() { return *reinterpret_cast<atVec3f*>(this); }
-
-  operator const atVec3f&() const { return *reinterpret_cast<const atVec3f*>(this); }
-
-  void readBig(athena::io::IStreamReader& input) {
-    simd_floats f;
-    f[0] = input.readFloatBig();
-    f[1] = input.readFloatBig();
-    f[2] = input.readFloatBig();
-    f[3] = 0.0f;
-    mSimd.copy_from(f);
-  }
-
-  [[nodiscard]] static CVector3f ReadBig(athena::io::IStreamReader& input) {
-    CVector3f ret;
-    ret.readBig(input);
-    return ret;
-  }
-
-#endif
-
   inline CVector3f(const CVector3d& vec);
 
   explicit constexpr CVector3f(float xyz) : mSimd(xyz) {}
@@ -139,6 +110,7 @@ public:
   [[nodiscard]] bool isNotInf() const { return !(std::isinf(x()) || std::isinf(y()) || std::isinf(z())); }
 
   [[nodiscard]] bool isMagnitudeSafe() const { return isNotInf() && magSquared() >= 9.9999994e-29; }
+  [[nodiscard]] bool isNaN() const { return std::isnan(x()) || std::isnan(y()) || std::isnan(z()); }
 
   void zeroOut() { mSimd.broadcast(0.f); }
 

include/zeus/CVector4d.hpp

@@ -0,0 +1,189 @@
+#pragma once
+
+#include <cassert>
+#include <cfloat>
+#include <cmath>
+
+#include "zeus/CVector3f.hpp"
+#include "zeus/Global.hpp"
+
+#if ZE_ATHENA_TYPES
+#include <athena/IStreamReader.hpp>
+#endif
+
+namespace zeus {
+class CColor;
+
+class CVector4d {
+public:
+  zeus::simd<double> mSimd;
+
+  constexpr CVector4d() : mSimd() {}
+
+  template <typename T>
+  constexpr CVector4d(const simd<T>& s) : mSimd(s) {}
+
+  explicit constexpr CVector4d(double xyzw) : mSimd(xyzw) {}
+
+  void assign(double x, double y, double z, double w) { mSimd = simd<double>(x, y, z, w); }
+
+  constexpr CVector4d(double x, double y, double z, double w) : mSimd(x, y, z, w) {}
+
+  constexpr CVector4d(const CColor& other);
+
+  CVector4d(const CVector3f& other, double wIn = 1.0) : mSimd(other.mSimd) { mSimd[3] = wIn; }
+
+  [[nodiscard]] static CVector4d ToClip(const zeus::CVector3f& v, double w) { return CVector4d(v * w, w); }
+
+  [[nodiscard]] CVector3f toVec3f() const { return CVector3f(mSimd); }
+
+  constexpr CVector4d& operator=(const CColor& other);
+
+  [[nodiscard]] bool operator==(const CVector4d& rhs) const {
+    auto eq_mask = mSimd == rhs.mSimd;
+    return eq_mask[0] && eq_mask[1] && eq_mask[2] && eq_mask[3];
+  }
+
+  [[nodiscard]] bool operator!=(const CVector4d& rhs) const {
+    auto eq_mask = mSimd != rhs.mSimd;
+    return eq_mask[0] || eq_mask[1] || eq_mask[2] || eq_mask[3];
+  }
+
+  [[nodiscard]] bool operator<(const CVector4d& rhs) const {
+    auto eq_mask = mSimd < rhs.mSimd;
+    return eq_mask[0] && eq_mask[1] && eq_mask[2] && eq_mask[3];
+  }
+
+  [[nodiscard]] bool operator<=(const CVector4d& rhs) const {
+    auto eq_mask = mSimd <= rhs.mSimd;
+    return eq_mask[0] && eq_mask[1] && eq_mask[2] && eq_mask[3];
+  }
+
+  [[nodiscard]] bool operator>(const CVector4d& rhs) const {
+    auto eq_mask = mSimd > rhs.mSimd;
+    return eq_mask[0] && eq_mask[1] && eq_mask[2] && eq_mask[3];
+  }
+
+  [[nodiscard]] bool operator>=(const CVector4d& rhs) const {
+    auto eq_mask = mSimd >= rhs.mSimd;
+    return eq_mask[0] && eq_mask[1] && eq_mask[2] && eq_mask[3];
+  }
+
+  [[nodiscard]] CVector4d operator+(const CVector4d& rhs) const { return mSimd + rhs.mSimd; }
+
+  [[nodiscard]] CVector4d operator-(const CVector4d& rhs) const { return mSimd - rhs.mSimd; }
+
+  [[nodiscard]] CVector4d operator-() const { return -mSimd; }
+
+  [[nodiscard]] CVector4d operator*(const CVector4d& rhs) const { return mSimd * rhs.mSimd; }
+
+  [[nodiscard]] CVector4d operator/(const CVector4d& rhs) const { return mSimd / rhs.mSimd; }
+
+  [[nodiscard]] CVector4d operator+(double val) const { return mSimd + zeus::simd<double>(val); }
+
+  [[nodiscard]] CVector4d operator-(double val) const { return mSimd - zeus::simd<double>(val); }
+
+  [[nodiscard]] CVector4d operator*(double val) const { return mSimd * zeus::simd<double>(val); }
+
+  [[nodiscard]] CVector4d operator/(double val) const {
+    double ooval = 1.0 / val;
+    return mSimd * zeus::simd<double>(ooval);
+  }
+
+  const CVector4d& operator+=(const CVector4d& rhs) {
+    mSimd += rhs.mSimd;
+    return *this;
+  }
+
+  const CVector4d& operator-=(const CVector4d& rhs) {
+    mSimd -= rhs.mSimd;
+    return *this;
+  }
+
+  const CVector4d& operator*=(const CVector4d& rhs) {
+    mSimd *= rhs.mSimd;
+    return *this;
+  }
+
+  const CVector4d& operator/=(const CVector4d& rhs) {
+    mSimd /= rhs.mSimd;
+    return *this;
+  }
+
+  void normalize() {
+    double mag = magnitude();
+    mag = 1.0 / mag;
+    *this *= CVector4d(mag);
+  }
+
+  [[nodiscard]] CVector4d normalized() const {
+    double mag = magnitude();
+    mag = 1.0 / mag;
+    return *this * mag;
+  }
+
+  [[nodiscard]] double dot(const CVector4d& rhs) const { return mSimd.dot4(rhs.mSimd); }
+
+  [[nodiscard]] double magSquared() const { return mSimd.dot4(mSimd); }
+
+  [[nodiscard]] double magnitude() const { return std::sqrt(magSquared()); }
+
+  void zeroOut() { mSimd = zeus::simd<double>(0.0); }
+
+  void splat(double xyzw) { mSimd = zeus::simd<double>(xyzw); }
+
+  [[nodiscard]] static CVector4d lerp(const CVector4d& a, const CVector4d& b, double t) {
+    return zeus::simd<double>(1.0 - t) * a.mSimd + b.mSimd * zeus::simd<double>(t);
+  }
+
+  [[nodiscard]] static CVector4d nlerp(const CVector4d& a, const CVector4d& b, double t) {
+    return lerp(a, b, t).normalized();
+  }
+
+  [[nodiscard]] bool isNormalized() const { return std::fabs(1.0 - magSquared()) < 0.01f; }
+
+  [[nodiscard]] bool canBeNormalized() const {
+    if (std::isinf(x()) || std::isinf(y()) || std::isinf(z()) || std::isinf(w()))
+      return false;
+    return std::fabs(x()) >= FLT_EPSILON || std::fabs(y()) >= FLT_EPSILON || std::fabs(z()) >= FLT_EPSILON ||
+           std::fabs(w()) >= FLT_EPSILON;
+  }
+
+  [[nodiscard]] bool isEqu(const CVector4d& other, double epsilon = FLT_EPSILON) const {
+    const CVector4d diffVec = other - *this;
+    return (diffVec.x() <= epsilon && diffVec.y() <= epsilon && diffVec.z() <= epsilon && diffVec.w() <= epsilon);
+  }
+
+  [[nodiscard]] simd<double>::reference operator[](size_t idx) {
+    assert(idx < 4);
+    return mSimd[idx];
+  }
+
+  [[nodiscard]] double operator[](size_t idx) const {
+    assert(idx < 4);
+    return mSimd[idx];
+  }
+
+  [[nodiscard]] double x() const { return mSimd[0]; }
+  [[nodiscard]] double y() const { return mSimd[1]; }
+  [[nodiscard]] double z() const { return mSimd[2]; }
+  [[nodiscard]] double w() const { return mSimd[3]; }
+
+  [[nodiscard]] simd<double>::reference x() { return mSimd[0]; }
+  [[nodiscard]] simd<double>::reference y() { return mSimd[1]; }
+  [[nodiscard]] simd<double>::reference z() { return mSimd[2]; }
+  [[nodiscard]] simd<double>::reference w() { return mSimd[3]; }
+};
+constexpr CVector4d skOne4d(1.0);
+constexpr CVector4d skNegOne4d(-1.0);
+constexpr CVector4d skZero4d(0.0);
+
+[[nodiscard]] inline CVector4d operator+(double lhs, const CVector4d& rhs) { return zeus::simd<double>(lhs) + rhs.mSimd; }
+
+[[nodiscard]] inline CVector4d operator-(double lhs, const CVector4d& rhs) { return zeus::simd<double>(lhs) - rhs.mSimd; }
+
+[[nodiscard]] inline CVector4d operator*(double lhs, const CVector4d& rhs) { return zeus::simd<double>(lhs) * rhs.mSimd; }
+
+[[nodiscard]] inline CVector4d operator/(double lhs, const CVector4d& rhs) { return zeus::simd<double>(lhs) / rhs.mSimd; }
+
+} // namespace zeus

include/zeus/CVector4f.hpp

@@ -23,25 +23,6 @@ public:
   template <typename T>
   constexpr CVector4f(const simd<T>& s) : mSimd(s) {}
 
-#if ZE_ATHENA_TYPES
-
-  constexpr CVector4f(const atVec4f& vec) : mSimd(vec.simd) {}
-
-  operator atVec4f&() { return *reinterpret_cast<atVec4f*>(this); }
-
-  operator const atVec4f&() const { return *reinterpret_cast<const atVec4f*>(this); }
-
-  void readBig(athena::io::IStreamReader& input) {
-    simd_floats f;
-    f[0] = input.readFloatBig();
-    f[1] = input.readFloatBig();
-    f[2] = input.readFloatBig();
-    f[3] = input.readFloatBig();
-    mSimd.copy_from(f);
-  }
-
-#endif
-
   explicit constexpr CVector4f(float xyzw) : mSimd(xyzw) {}
 
   void assign(float x, float y, float z, float w) { mSimd = simd<float>(x, y, z, w); }

include/zeus/Global.hpp

@@ -1,20 +1,8 @@
 #pragma once
 
-#if ZE_ATHENA_TYPES
-#include "athena/IStreamReader.hpp"
-#include "athena/simd/simd.hpp"
-#else
 #include "simd/simd.hpp"
-#endif
 
 namespace zeus {
-#if ZE_ATHENA_TYPES
-template <typename T>
-using simd = athena::simd<T>;
-using simd_floats = athena::simd_floats;
-using simd_doubles = athena::simd_doubles;
-#endif
-
 template <typename SizeT>
 constexpr void hash_combine_impl(SizeT& seed, SizeT value) noexcept {
   seed ^= value + 0x9e3779b9 + (seed << 6) + (seed >> 2);

include/zeus/Math.hpp

@@ -184,10 +184,10 @@ template <typename E>
 [[nodiscard]] bool close_enough(const CVector2f& a, const CVector2f& b, float epsilon = FLT_EPSILON);
 
 [[nodiscard]] inline bool close_enough(float a, float b, double epsilon = FLT_EPSILON) {
-  return std::fabs(a - b) < epsilon;
+  return std::fabs(a - b) <= epsilon;
 }
 
 [[nodiscard]] inline bool close_enough(double a, double b, double epsilon = FLT_EPSILON) {
-  return std::fabs(a - b) < epsilon;
+  return std::fabs(a - b) <= epsilon;
 }
 } // namespace zeus

include/zeus/simd/parallelism_v2_simd.hpp

@@ -1477,6 +1477,9 @@ private:
   friend class simd_mask;
 
 public:
+  constexpr __simd_storage(_Tp __rv) : __storage_{__rv, __rv, __rv, __rv} {}
+  constexpr __simd_storage(_Tp a, _Tp b, _Tp c, _Tp d) : __storage_{a, b, c, d} {}
+
   constexpr _Tp __get(size_t __index) const noexcept { return __storage_[__index]; };
   constexpr void __set(size_t __index, _Tp __val) noexcept { __storage_[__index] = __val; }
   constexpr std::enable_if_t<__num_element >= 4> __set4(float a, float b, float c, float d) noexcept {

include/zeus/simd/simd.hpp

@@ -18,7 +18,7 @@ using namespace std;
 #elif __ARM_NEON
 #include "simd_neon.hpp"
 #else
-namespace simd_abi {
+namespace zeus::_simd::simd_abi {
 template <typename T>
 struct zeus_native {};
 template <>
@@ -29,7 +29,8 @@ template <>
 struct zeus_native<double> {
   using type = fixed_size<4>;
 };
-} // namespace simd_abi
+} // namespace zeus::_simd::simd_abi
+#include "simd_none.hpp"
 #endif
 #ifdef __GNUC__
 #pragma GCC diagnostic pop

include/zeus/simd/simd_neon.hpp

@@ -18,7 +18,7 @@ class __simd_storage<float, m128_abi> {
 public:
   using storage_type = float32x4_t;
   storage_type __storage_{};
-  [[nodiscard]] constexpr float __get(size_t __index) const noexcept { return __storage_[__index]; }
+  [[nodiscard]] float __get(size_t __index) const noexcept { return __storage_[__index]; }
   inline void __set(size_t __index, float __val) noexcept { __storage_[__index] = __val; }
   constexpr __simd_storage(float a, float b, float c, float d) : __storage_{a, b, c, d} {}
   constexpr void __set4(float a, float b, float c, float d) noexcept { __storage_ = storage_type{a, b, c, d}; }
@@ -35,12 +35,7 @@ public:
   }
   template <int x, int y, int z, int w>
   [[nodiscard]] inline __simd_storage __shuffle() const noexcept {
-    storage_type ret;
-    ret = vmovq_n_f32(vgetq_lane_f32(__storage_, x));
-    ret = vsetq_lane_f32(vgetq_lane_f32(__storage_, y), ret, 1);
-    ret = vsetq_lane_f32(vgetq_lane_f32(__storage_, z), ret, 2);
-    ret = vsetq_lane_f32(vgetq_lane_f32(__storage_, w), ret, 3);
-    return __simd_storage(ret);
+    return __simd_storage{__storage_[x], __storage_[y], __storage_[z], __storage_[w]};
   }
 
   inline void __copy_from(const simd_data<simd<float, m128_abi>>& __buffer) noexcept {
@@ -71,8 +66,7 @@ public:
 
 template <>
 inline simd<float, m128_abi> simd<float, m128_abi>::operator-() const {
-  return vreinterpretq_f32_s32(
-      veorq_s32(vreinterpretq_s32_f32(__s_.__storage_), vreinterpretq_s32_f32(vdupq_n_f32(-0.f))));
+  return vnegq_f32(__s_.__storage_);
 }
 
 inline simd<float, m128_abi> operator+(const simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
@@ -172,7 +166,7 @@ public:
   [[nodiscard]] inline double __dot3(const __simd_storage<double, m128d_abi>& other) const noexcept {
     const vector_type mul1 = vmulq_f64(__storage_.val[0], other.__storage_.val[0]);
     const vector_type mul2 = vmulq_f64(__storage_.val[1], other.__storage_.val[1]);
-    return vaddvq_f64(vcombine_f64(vcreate_f64(vaddvq_f64(mul1)), vget_low_f64(mul2)));
+    return vaddvq_f64(vcombine_f64(vdup_n_f64(vaddvq_f64(mul1)), vget_low_f64(mul2)));
   }
   [[nodiscard]] inline double __dot4(const __simd_storage<double, m128d_abi>& other) const noexcept {
     const vector_type mul1 = vmulq_f64(__storage_.val[0], other.__storage_.val[0]);
@@ -215,8 +209,7 @@ template <>
 inline simd<double, m128d_abi> simd<double, m128d_abi>::operator-() const {
   simd<double, m128d_abi> ret;
   for (int i = 0; i < 2; ++i)
-    ret.__s_.__storage_.val[i] = vreinterpretq_f64_s64(
-        veorq_s64(vreinterpretq_s64_f64(__s_.__storage_.val[i]), vreinterpretq_s64_f64(vdupq_n_f64(-0.0))));
+    ret.__s_.__storage_.val[i] = vnegq_f64(__s_.__storage_.val[i]);
   return ret;
 }
 

include/zeus/simd/simd_none.hpp

@@ -0,0 +1,224 @@
+#pragma once
+#ifndef _ZEUS_SIMD_INCLUDED
+#error simd_none.hpp must not be included directly. Include simd.hpp instead.
+#endif
+namespace zeus::_simd {
+using m128_abi = __simd_abi<_StorageKind::_Array, 4>;
+using m128d_abi = __simd_abi<_StorageKind::_Array, 4>;
+
+// m128 ABI
+template <>
+inline simd<float, m128_abi> simd<float, m128_abi>::operator-() const {
+  return {-__s_.__storage_[0], -__s_.__storage_[1], -__s_.__storage_[2], -__s_.__storage_[3]};
+}
+
+inline simd<float, m128_abi> operator+(const simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  return {a[0] + b[0], a[1] + b[1], a[2] + b[2], a[3] + b[3]};
+}
+
+inline simd<float, m128_abi> operator-(const simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  return {a[0] - b[0], a[1] - b[1], a[2] - b[2], a[3] - b[3]};
+}
+
+inline simd<float, m128_abi> operator*(const simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  return {a[0] * b[0], a[1] * b[1], a[2] * b[2], a[3] * b[3]};
+}
+
+inline simd<float, m128_abi> operator/(const simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  return {a[0] / b[0], a[1] / b[1], a[2] / b[2], a[3] / b[3]};
+}
+
+inline simd<float, m128_abi>& operator+=(simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  a[0] += b[0];
+  a[1] += b[1];
+  a[2] += b[2];
+  a[3] += b[3];
+  return a;
+}
+
+inline simd<float, m128_abi>& operator-=(simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  a[0] -= b[0];
+  a[1] -= b[1];
+  a[2] -= b[2];
+  a[3] -= b[3];
+  return a;
+}
+
+inline simd<float, m128_abi>& operator*=(simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  a[0] *= b[0];
+  a[1] *= b[1];
+  a[2] *= b[2];
+  a[3] *= b[3];
+  return a;
+}
+
+inline simd<float, m128_abi>& operator/=(simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  a[0] /= b[0];
+  a[1] /= b[1];
+  a[2] /= b[2];
+  a[3] /= b[3];
+  return a;
+}
+
+inline simd<float, m128_abi>::mask_type operator==(const simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  simd<float, m128_abi>::mask_type ret;
+  ret[0] = a[0] == b[0];
+  ret[1] = a[1] == b[1];
+  ret[2] = a[2] == b[2];
+  ret[3] = a[3] == b[3];
+  return ret;
+}
+
+inline simd<float, m128_abi>::mask_type operator!=(const simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  simd<float, m128_abi>::mask_type ret;
+  ret[0] = a[0] != b[0];
+  ret[1] = a[1] != b[1];
+  ret[2] = a[2] != b[2];
+  ret[3] = a[3] != b[3];
+  return ret;
+}
+
+inline simd<float, m128_abi>::mask_type operator>=(const simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  simd<float, m128_abi>::mask_type ret;
+  ret[0] = a[0] >= b[0];
+  ret[1] = a[1] >= b[1];
+  ret[2] = a[2] >= b[2];
+  ret[3] = a[3] >= b[3];
+  return ret;
+}
+
+inline simd<float, m128_abi>::mask_type operator<=(const simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  simd<float, m128_abi>::mask_type ret;
+  ret[0] = a[0] <= b[0];
+  ret[1] = a[1] <= b[1];
+  ret[2] = a[2] <= b[2];
+  ret[3] = a[3] <= b[3];
+  return ret;
+}
+
+inline simd<float, m128_abi>::mask_type operator>(const simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  simd<float, m128_abi>::mask_type ret;
+  ret[0] = a[0] > b[0];
+  ret[1] = a[1] > b[1];
+  ret[2] = a[2] > b[2];
+  ret[3] = a[3] > b[3];
+  return ret;
+}
+
+inline simd<float, m128_abi>::mask_type operator<(const simd<float, m128_abi>& a, const simd<float, m128_abi>& b) {
+  simd<float, m128_abi>::mask_type ret;
+  ret[0] = a[0] < b[0];
+  ret[1] = a[1] < b[1];
+  ret[2] = a[2] < b[2];
+  ret[3] = a[3] < b[3];
+  return ret;
+}
+
+// m128d ABI
+template <>
+inline simd<double, m128d_abi> simd<double, m128d_abi>::operator-() const {
+  return {-__s_.__storage_[0], -__s_.__storage_[1], -__s_.__storage_[2], -__s_.__storage_[3]};
+}
+
+inline simd<double, m128d_abi> operator+(const simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  return {a[0] + b[0], a[1] + b[1], a[2] + b[2], a[3] + b[3]};
+}
+
+inline simd<double, m128d_abi> operator-(const simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  return {a[0] - b[0], a[1] - b[1], a[2] - b[2], a[3] - b[3]};
+}
+
+inline simd<double, m128d_abi> operator*(const simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  return {a[0] * b[0], a[1] * b[1], a[2] * b[2], a[3] * b[3]};
+}
+
+inline simd<double, m128d_abi> operator/(const simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  return {a[0] / b[0], a[1] / b[1], a[2] / b[2], a[3] / b[3]};
+}
+
+inline simd<double, m128d_abi>& operator+=(simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  a[0] += b[0];
+  a[1] += b[1];
+  a[2] += b[2];
+  a[3] += b[3];
+  return a;
+}
+
+inline simd<double, m128d_abi>& operator-=(simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  a[0] -= b[0];
+  a[1] -= b[1];
+  a[2] -= b[2];
+  a[3] -= b[3];
+  return a;
+}
+
+inline simd<double, m128d_abi>& operator*=(simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  a[0] *= b[0];
+  a[1] *= b[1];
+  a[2] *= b[2];
+  a[3] *= b[3];
+  return a;
+}
+
+inline simd<double, m128d_abi>& operator/=(simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  a[0] /= b[0];
+  a[1] /= b[1];
+  a[2] /= b[2];
+  a[3] /= b[3];
+  return a;
+}
+
+inline simd<double, m128d_abi>::mask_type operator==(const simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  simd<double, m128d_abi>::mask_type ret;
+  ret[0] = a[0] == b[0];
+  ret[1] = a[1] == b[1];
+  ret[2] = a[2] == b[2];
+  ret[3] = a[3] == b[3];
+  return ret;
+}
+
+inline simd<double, m128d_abi>::mask_type operator!=(const simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  simd<double, m128d_abi>::mask_type ret;
+  ret[0] = a[0] != b[0];
+  ret[1] = a[1] != b[1];
+  ret[2] = a[2] != b[2];
+  ret[3] = a[3] != b[3];
+  return ret;
+}
+
+inline simd<double, m128d_abi>::mask_type operator>=(const simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  simd<double, m128d_abi>::mask_type ret;
+  ret[0] = a[0] >= b[0];
+  ret[1] = a[1] >= b[1];
+  ret[2] = a[2] >= b[2];
+  ret[3] = a[3] >= b[3];
+  return ret;
+}
+
+inline simd<double, m128d_abi>::mask_type operator<=(const simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  simd<double, m128d_abi>::mask_type ret;
+  ret[0] = a[0] <= b[0];
+  ret[1] = a[1] <= b[1];
+  ret[2] = a[2] <= b[2];
+  ret[3] = a[3] <= b[3];
+  return ret;
+}
+
+inline simd<double, m128d_abi>::mask_type operator>(const simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  simd<double, m128d_abi>::mask_type ret;
+  ret[0] = a[0] > b[0];
+  ret[1] = a[1] > b[1];
+  ret[2] = a[2] > b[2];
+  ret[3] = a[3] > b[3];
+  return ret;
+}
+
+inline simd<double, m128d_abi>::mask_type operator<(const simd<double, m128d_abi>& a, const simd<double, m128d_abi>& b) {
+  simd<double, m128d_abi>::mask_type ret;
+  ret[0] = a[0] < b[0];
+  ret[1] = a[1] < b[1];
+  ret[2] = a[2] < b[2];
+  ret[3] = a[3] < b[3];
+  return ret;
+}
+} // namespace zeus::_simd

include/zeus/simd/simd_sse.hpp

@@ -1,4 +1,5 @@
 #pragma once
+#include <cstdint>
 #ifndef _ZEUS_SIMD_INCLUDED
 #error simd_sse.hpp must not be included directly. Include simd.hpp instead.
 #endif
@@ -30,8 +31,25 @@ class __simd_storage<float, m128_abi> {
 public:
   using storage_type = __m128;
   storage_type __storage_{};
-  [[nodiscard]] inline float __get(size_t __index) const noexcept { return __storage_[__index]; }
-  inline void __set(size_t __index, float __val) noexcept { __storage_[__index] = __val; }
+  [[nodiscard]] inline float __get(size_t __index) const noexcept {
+#if _MSC_VER && !defined(__clang__)
+    alignas(16) std::array<float, 4> sse_data;
+    _mm_store_ps(sse_data.data(), __storage_);
+    return sse_data[__index];
+#else
+    return __storage_[__index];
+#endif
+  }
+  inline void __set(size_t __index, float __val) noexcept {
+#if _MSC_VER && !defined(__clang__)
+    alignas(16) std::array<float, 4> sse_data;
+    _mm_store_ps(sse_data.data(), __storage_);
+    sse_data[__index] = __val;
+    __storage_ = _mm_load_ps(sse_data.data());
+#else
+    __storage_[__index] = __val;
+#endif
+  }
   constexpr __simd_storage(float a, float b, float c, float d) : __storage_{a, b, c, d} {}
   constexpr void __set4(float a, float b, float c, float d) noexcept { __storage_ = storage_type{a, b, c, d}; }
   constexpr explicit __simd_storage(float rv) : __storage_{rv, rv, rv, rv} {}
@@ -192,8 +210,25 @@ class __simd_storage<double, m128d_abi> {
 public:
   using storage_type = std::array<__m128d, 2>;
   storage_type __storage_{};
-  [[nodiscard]] inline double __get(size_t __index) const noexcept { return __storage_[__index / 2][__index % 2]; }
-  inline void __set(size_t __index, double __val) noexcept { __storage_[__index / 2][__index % 2] = __val; }
+  [[nodiscard]] inline double __get(size_t __index) const noexcept {
+#if _MSC_VER && !defined(__clang__)
+    alignas(16) std::array<double, 2> sse_data;
+    _mm_store_pd(sse_data.data(), __storage_[__index / 2]);
+    return sse_data[__index % 2];
+#else
+    return __storage_[__index / 2][__index % 2];
+#endif
+  }
+  inline void __set(size_t __index, double __val) noexcept {
+#if _MSC_VER && !defined(__clang__)
+    alignas(16) std::array<double, 2> sse_data;
+    _mm_store_pd(sse_data.data(), __storage_[__index / 2]);
+    sse_data[__index % 2] = __val;
+    __storage_[__index / 2] = _mm_load_pd(sse_data.data());
+#else
+    __storage_[__index / 2][__index % 2] = __val;
+#endif
+  }
   // Make GCC happy
   static constexpr storage_type __make_array(__m128d a, __m128d b) { return {a, b}; }
   constexpr __simd_storage(double a, double b, double c, double d) : __storage_(__make_array(__m128d{a, b}, __m128d{c, d})) {}

include/zeus/zeus.hpp

@@ -18,8 +18,11 @@
 #include "zeus/CTransform.hpp"
 #include "zeus/CUnitVector.hpp"
 #include "zeus/CVector2f.hpp"
-#include "zeus/CVector3d.hpp"
+#include "zeus/CVector2i.hpp"
+#include "zeus/CVector2d.hpp"
 #include "zeus/CVector3f.hpp"
+#include "zeus/CVector3d.hpp"
 #include "zeus/CVector4f.hpp"
+#include "zeus/CVector4d.hpp"
 #include "zeus/Global.hpp"
 #include "zeus/Math.hpp"

src/CMatrix3f.cpp

@@ -4,8 +4,7 @@
 
 namespace zeus {
 
-CMatrix3f::CMatrix3f(const CQuaternion& quat) {
-  CQuaternion nq = quat.normalized();
+CMatrix3f::CMatrix3f(const CQuaternion& nq) {
   float x2 = nq.x() * nq.x();
   float y2 = nq.y() * nq.y();
   float z2 = nq.z() * nq.z();
@@ -47,15 +46,15 @@ void CMatrix3f::transpose() {
   float tmp;
 
   tmp = m[0][1];
-  m[0][1] = m[1][0];
+  m[0][1] = m[1][0].operator float();
   m[1][0] = tmp;
 
   tmp = m[0][2];
-  m[0][2] = m[2][0];
+  m[0][2] = m[2][0].operator float();
   m[2][0] = tmp;
 
   tmp = m[1][2];
-  m[1][2] = m[2][1];
+  m[1][2] = m[2][1].operator float();
   m[2][1] = tmp;
 #endif
 }
@@ -81,15 +80,15 @@ CMatrix3f CMatrix3f::transposed() const {
   float tmp;
 
   tmp = ret.m[0][1];
-  ret.m[0][1] = ret.m[1][0];
+  ret.m[0][1] = ret.m[1][0].operator float();
   ret.m[1][0] = tmp;
 
   tmp = m[0][2];
-  ret.m[0][2] = ret.m[2][0];
+  ret.m[0][2] = ret.m[2][0].operator float();
   ret.m[2][0] = tmp;
 
   tmp = m[1][2];
-  ret.m[1][2] = ret.m[2][1];
+  ret.m[1][2] = ret.m[2][1].operator float();
   ret.m[2][1] = tmp;
 
   return ret;

src/CQuaternion.cpp

@@ -79,25 +79,6 @@ void CQuaternion::fromVector3f(const CVector3f& vec) {
   mSimd.copy_from(f);
 }
 
-CQuaternion CQuaternion::operator*(const CQuaternion& q) const {
-  return CQuaternion(w() * q.w() - CVector3f(x(), y(), z()).dot({q.x(), q.y(), q.z()}),
-                     y() * q.z() - z() * q.y() + w() * q.x() + x() * q.w(),
-                     z() * q.x() - x() * q.z() + w() * q.y() + y() * q.w(),
-                     x() * q.y() - y() * q.x() + w() * q.z() + z() * q.w());
-}
-
-CNUQuaternion CNUQuaternion::operator*(const CNUQuaternion& q) const {
-  return CNUQuaternion(w() * q.w() - CVector3f(x(), y(), z()).dot({q.x(), q.y(), q.z()}),
-                       y() * q.z() - z() * q.y() + w() * q.x() + x() * q.w(),
-                       z() * q.x() - x() * q.z() + w() * q.y() + y() * q.w(),
-                       x() * q.y() - y() * q.x() + w() * q.z() + z() * q.w());
-}
-
-CQuaternion CQuaternion::operator/(const CQuaternion& q) const {
-  CQuaternion p(q);
-  p.invert();
-  return *this * p;
-}
 
 const CQuaternion& CQuaternion::operator*=(const CQuaternion& q) {
   CQuaternion orig = *this;

src/Math.cpp

@@ -21,7 +21,7 @@ static CPUInfo g_cpuFeatures = {};
 static CPUInfo g_missingFeatures = {};
 
 void getCpuInfo(int eax, int regs[4]) {
-#if __x86_64__
+#if defined(__x86_64__) || defined(_M_X64)
   #if _WIN32
   __cpuid(regs, eax);
 #else
@@ -31,7 +31,7 @@ void getCpuInfo(int eax, int regs[4]) {
 }
 
 void getCpuInfoEx(int eax, int ecx, int regs[4]) {
-#if __x86_64__
+#if defined(__x86_64__) || defined(_M_X64)
   #if _WIN32
   __cpuidex(regs, eax, ecx);
 #else
@@ -41,7 +41,7 @@ void getCpuInfoEx(int eax, int ecx, int regs[4]) {
 }
 
 void detectCPU() {
-#if __x86_64__
+#if defined(__x86_64__) || defined(_M_X64)
   if (isCPUInit)
     return;
 
@@ -183,19 +183,14 @@ CTransform lookAt(const CVector3f& pos, const CVector3f& lookPos, const CVector3
 
 CVector3f getBezierPoint(const CVector3f& a, const CVector3f& b, const CVector3f& c, const CVector3f& d, float t) {
   const float omt = 1.f - t;
-  return ((a * omt + b * t) * omt + (b * omt + c * t) * t) * omt +
+  return (((a * omt) + b * t) * omt + (b * omt + c * t) * t) * omt +
          ((b * omt + c * t) * omt + (c * omt + d * t) * t) * t;
 }
 
 int floorPowerOfTwo(int x) {
   if (x == 0)
     return 0;
-  /*
-   * we want to ensure that we always get the previous power,
-   * but if we have values like 256, we'll always get the same value,
-   * x-1 ensures that we always get the previous power.
-   */
-  x = (x - 1) | (x >> 1);
+  x = x | (x >> 1);
   x = x | (x >> 2);
   x = x | (x >> 4);
   x = x | (x >> 8);
@@ -277,11 +272,11 @@ CVector3f baryToWorld(const CVector3f& p0, const CVector3f& p1, const CVector3f&
 }
 
 bool close_enough(const CVector3f& a, const CVector3f& b, float epsilon) {
-  return std::fabs(a.x() - b.x()) < epsilon && std::fabs(a.y() - b.y()) < epsilon && std::fabs(a.z() - b.z()) < epsilon;
+  return std::fabs(a.x() - b.x()) <= epsilon && std::fabs(a.y() - b.y()) <= epsilon && std::fabs(a.z() - b.z()) <= epsilon;
 }
 
 bool close_enough(const CVector2f& a, const CVector2f& b, float epsilon) {
-  return std::fabs(a.x() - b.x()) < epsilon && std::fabs(a.y() - b.y()) < epsilon;
+  return std::fabs(a.x() - b.x()) <= epsilon && std::fabs(a.y() - b.y()) <= epsilon;
 }
 
 template <>