Initial switch to raylib 5.5.

1 files changed, 784 insertions, 33 deletions

diff --git a/include/raymath.h b/include/raymath.h
index ff60170..e522113 100644
--- a/include/raymath.h
+++ b/include/raymath.h
@@ -1,6 +1,6 @@
 /**********************************************************************************************
 *
-*   raymath v1.5 - Math functions to work with Vector2, Vector3, Matrix and Quaternions
+*   raymath v2.0 - Math functions to work with Vector2, Vector3, Matrix and Quaternions
 *
 *   CONVENTIONS:
 *     - Matrix structure is defined as row-major (memory layout) but parameters naming AND all
@@ -12,7 +12,7 @@
 *     - Functions are always self-contained, no function use another raymath function inside,
 *       required code is directly re-implemented inside
 *     - Functions input parameters are always received by value (2 unavoidable exceptions)
-*     - Functions use always a "result" variable for return
+*     - Functions use always a "result" variable for return (except C++ operators)
 *     - Functions are always defined inline
 *     - Angles are always in radians (DEG2RAD/RAD2DEG macros provided for convenience)
 *     - No compound literals used to make sure libray is compatible with C++
@@ -27,10 +27,12 @@
 *           Define static inline functions code, so #include header suffices for use.
 *           This may use up lots of memory.
 *
+*       #define RAYMATH_DISABLE_CPP_OPERATORS
+*           Disables C++ operator overloads for raymath types.
 *
 *   LICENSE: zlib/libpng
 *
-*   Copyright (c) 2015-2023 Ramon Santamaria (@raysan5)
+*   Copyright (c) 2015-2024 Ramon Santamaria (@raysan5)
 *
 *   This software is provided "as-is", without any express or implied warranty. In no event
 *   will the authors be held liable for any damages arising from the use of this software.
@@ -59,7 +61,9 @@
 // Function specifiers definition
 #if defined(RAYMATH_IMPLEMENTATION)
     #if defined(_WIN32) && defined(BUILD_LIBTYPE_SHARED)
-        #define RMAPI __declspec(dllexport) extern inline // We are building raylib as a Win32 shared library (.dll).
+        #define RMAPI __declspec(dllexport) extern inline // We are building raylib as a Win32 shared library (.dll)
+    #elif defined(BUILD_LIBTYPE_SHARED)
+        #define RMAPI __attribute__((visibility("default"))) // We are building raylib as a Unix shared library (.so/.dylib)
     #elif defined(_WIN32) && defined(USE_LIBTYPE_SHARED)
         #define RMAPI __declspec(dllimport)         // We are using raylib as a Win32 shared library (.dll)
     #else
@@ -75,6 +79,7 @@
     #endif
 #endif
 
+
 //----------------------------------------------------------------------------------
 // Defines and Macros
 //----------------------------------------------------------------------------------
@@ -163,7 +168,7 @@ typedef struct float16 {
     float v[16];
 } float16;
 
-#include <math.h>       // Required for: sinf(), cosf(), tan(), atan2f(), sqrtf(), floor(), fminf(), fmaxf(), fabs()
+#include <math.h>       // Required for: sinf(), cosf(), tan(), atan2f(), sqrtf(), floor(), fminf(), fmaxf(), fabsf()
 
 //----------------------------------------------------------------------------------
 // Module Functions Definition - Utils math
@@ -429,6 +434,28 @@ RMAPI Vector2 Vector2Reflect(Vector2 v, Vector2 normal)
     return result;
 }
 
+// Get min value for each pair of components
+RMAPI Vector2 Vector2Min(Vector2 v1, Vector2 v2)
+{
+    Vector2 result = { 0 };
+
+    result.x = fminf(v1.x, v2.x);
+    result.y = fminf(v1.y, v2.y);
+
+    return result;
+}
+
+// Get max value for each pair of components
+RMAPI Vector2 Vector2Max(Vector2 v1, Vector2 v2)
+{
+    Vector2 result = { 0 };
+
+    result.x = fmaxf(v1.x, v2.x);
+    result.y = fmaxf(v1.y, v2.y);
+
+    return result;
+}
+
 // Rotate vector by angle
 RMAPI Vector2 Vector2Rotate(Vector2 v, float angle)
 {
@@ -492,18 +519,18 @@ RMAPI Vector2 Vector2ClampValue(Vector2 v, float min, float max)
     {
         length = sqrtf(length);
 
+        float scale = 1;    // By default, 1 as the neutral element.
         if (length < min)
         {
-            float scale = min/length;
-            result.x = v.x*scale;
-            result.y = v.y*scale;
+            scale = min/length;
         }
         else if (length > max)
         {
-            float scale = max/length;
-            result.x = v.x*scale;
-            result.y = v.y*scale;
+            scale = max/length;
         }
+
+        result.x = v.x*scale;
+        result.y = v.y*scale;
     }
 
     return result;
@@ -522,6 +549,31 @@ RMAPI int Vector2Equals(Vector2 p, Vector2 q)
     return result;
 }
 
+// Compute the direction of a refracted ray
+// v: normalized direction of the incoming ray
+// n: normalized normal vector of the interface of two optical media
+// r: ratio of the refractive index of the medium from where the ray comes
+//    to the refractive index of the medium on the other side of the surface
+RMAPI Vector2 Vector2Refract(Vector2 v, Vector2 n, float r)
+{
+    Vector2 result = { 0 };
+
+    float dot = v.x*n.x + v.y*n.y;
+    float d = 1.0f - r*r*(1.0f - dot*dot);
+
+    if (d >= 0.0f)
+    {
+        d = sqrtf(d);
+        v.x = r*v.x - (r*dot + d)*n.x;
+        v.y = r*v.y - (r*dot + d)*n.y;
+
+        result = v;
+    }
+
+    return result;
+}
+
+
 //----------------------------------------------------------------------------------
 // Module Functions Definition - Vector3 math
 //----------------------------------------------------------------------------------
@@ -603,12 +655,12 @@ RMAPI Vector3 Vector3Perpendicular(Vector3 v)
 {
     Vector3 result = { 0 };
 
-    float min = (float) fabs(v.x);
+    float min = fabsf(v.x);
     Vector3 cardinalAxis = {1.0f, 0.0f, 0.0f};
 
     if (fabsf(v.y) < min)
     {
-        min = (float) fabs(v.y);
+        min = fabsf(v.y);
         Vector3 tmp = {0.0f, 1.0f, 0.0f};
         cardinalAxis = tmp;
     }
@@ -728,7 +780,7 @@ RMAPI Vector3 Vector3Normalize(Vector3 v)
 RMAPI Vector3 Vector3Project(Vector3 v1, Vector3 v2)
 {
     Vector3 result = { 0 };
-    
+
     float v1dv2 = (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z);
     float v2dv2 = (v2.x*v2.x + v2.y*v2.y + v2.z*v2.z);
 
@@ -745,7 +797,7 @@ RMAPI Vector3 Vector3Project(Vector3 v1, Vector3 v2)
 RMAPI Vector3 Vector3Reject(Vector3 v1, Vector3 v2)
 {
     Vector3 result = { 0 };
-    
+
     float v1dv2 = (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z);
     float v2dv2 = (v2.x*v2.x + v2.y*v2.y + v2.z*v2.z);
 
@@ -832,7 +884,7 @@ RMAPI Vector3 Vector3RotateByAxisAngle(Vector3 v, Vector3 axis, float angle)
     // Vector3Normalize(axis);
     float length = sqrtf(axis.x*axis.x + axis.y*axis.y + axis.z*axis.z);
     if (length == 0.0f) length = 1.0f;
-    float ilength = 1.0f / length;
+    float ilength = 1.0f/length;
     axis.x *= ilength;
     axis.y *= ilength;
     axis.z *= ilength;
@@ -873,6 +925,27 @@ RMAPI Vector3 Vector3RotateByAxisAngle(Vector3 v, Vector3 axis, float angle)
     return result;
 }
 
+// Move Vector towards target
+RMAPI Vector3 Vector3MoveTowards(Vector3 v, Vector3 target, float maxDistance)
+{
+    Vector3 result = { 0 };
+
+    float dx = target.x - v.x;
+    float dy = target.y - v.y;
+    float dz = target.z - v.z;
+    float value = (dx*dx) + (dy*dy) + (dz*dz);
+
+    if ((value == 0) || ((maxDistance >= 0) && (value <= maxDistance*maxDistance))) return target;
+
+    float dist = sqrtf(value);
+
+    result.x = v.x + dx/dist*maxDistance;
+    result.y = v.y + dy/dist*maxDistance;
+    result.z = v.z + dz/dist*maxDistance;
+
+    return result;
+}
+
 // Calculate linear interpolation between two vectors
 RMAPI Vector3 Vector3Lerp(Vector3 v1, Vector3 v2, float amount)
 {
@@ -885,6 +958,22 @@ RMAPI Vector3 Vector3Lerp(Vector3 v1, Vector3 v2, float amount)
     return result;
 }
 
+// Calculate cubic hermite interpolation between two vectors and their tangents
+// as described in the GLTF 2.0 specification: https://registry.khronos.org/glTF/specs/2.0/glTF-2.0.html#interpolation-cubic
+RMAPI Vector3 Vector3CubicHermite(Vector3 v1, Vector3 tangent1, Vector3 v2, Vector3 tangent2, float amount)
+{
+    Vector3 result = { 0 };
+
+    float amountPow2 = amount*amount;
+    float amountPow3 = amount*amount*amount;
+
+    result.x = (2*amountPow3 - 3*amountPow2 + 1)*v1.x + (amountPow3 - 2*amountPow2 + amount)*tangent1.x + (-2*amountPow3 + 3*amountPow2)*v2.x + (amountPow3 - amountPow2)*tangent2.x;
+    result.y = (2*amountPow3 - 3*amountPow2 + 1)*v1.y + (amountPow3 - 2*amountPow2 + amount)*tangent1.y + (-2*amountPow3 + 3*amountPow2)*v2.y + (amountPow3 - amountPow2)*tangent2.y;
+    result.z = (2*amountPow3 - 3*amountPow2 + 1)*v1.z + (amountPow3 - 2*amountPow2 + amount)*tangent1.z + (-2*amountPow3 + 3*amountPow2)*v2.z + (amountPow3 - amountPow2)*tangent2.z;
+
+    return result;
+}
+
 // Calculate reflected vector to normal
 RMAPI Vector3 Vector3Reflect(Vector3 v, Vector3 normal)
 {
@@ -1078,20 +1167,19 @@ RMAPI Vector3 Vector3ClampValue(Vector3 v, float min, float max)
     {
         length = sqrtf(length);
 
+        float scale = 1;    // By default, 1 as the neutral element.
         if (length < min)
         {
-            float scale = min/length;
-            result.x = v.x*scale;
-            result.y = v.y*scale;
-            result.z = v.z*scale;
+            scale = min/length;
         }
         else if (length > max)
         {
-            float scale = max/length;
-            result.x = v.x*scale;
-            result.y = v.y*scale;
-            result.z = v.z*scale;
+            scale = max/length;
         }
+
+        result.x = v.x*scale;
+        result.y = v.y*scale;
+        result.z = v.z*scale;
     }
 
     return result;
@@ -1136,6 +1224,233 @@ RMAPI Vector3 Vector3Refract(Vector3 v, Vector3 n, float r)
     return result;
 }
 
+
+//----------------------------------------------------------------------------------
+// Module Functions Definition - Vector4 math
+//----------------------------------------------------------------------------------
+
+RMAPI Vector4 Vector4Zero(void)
+{
+    Vector4 result = { 0.0f, 0.0f, 0.0f, 0.0f };
+    return result;
+}
+
+RMAPI Vector4 Vector4One(void)
+{
+    Vector4 result = { 1.0f, 1.0f, 1.0f, 1.0f };
+    return result;
+}
+
+RMAPI Vector4 Vector4Add(Vector4 v1, Vector4 v2)
+{
+    Vector4 result = {
+        v1.x + v2.x,
+        v1.y + v2.y,
+        v1.z + v2.z,
+        v1.w + v2.w
+    };
+    return result;
+}
+
+RMAPI Vector4 Vector4AddValue(Vector4 v, float add)
+{
+    Vector4 result = {
+        v.x + add,
+        v.y + add,
+        v.z + add,
+        v.w + add
+    };
+    return result;
+}
+
+RMAPI Vector4 Vector4Subtract(Vector4 v1, Vector4 v2)
+{
+    Vector4 result = {
+        v1.x - v2.x,
+        v1.y - v2.y,
+        v1.z - v2.z,
+        v1.w - v2.w
+    };
+    return result;
+}
+
+RMAPI Vector4 Vector4SubtractValue(Vector4 v, float add)
+{
+    Vector4 result = {
+        v.x - add,
+        v.y - add,
+        v.z - add,
+        v.w - add
+    };
+    return result;
+}
+
+RMAPI float Vector4Length(Vector4 v)
+{
+    float result = sqrtf((v.x*v.x) + (v.y*v.y) + (v.z*v.z) + (v.w*v.w));
+    return result;
+}
+
+RMAPI float Vector4LengthSqr(Vector4 v)
+{
+    float result = (v.x*v.x) + (v.y*v.y) + (v.z*v.z) + (v.w*v.w);
+    return result;
+}
+
+RMAPI float Vector4DotProduct(Vector4 v1, Vector4 v2)
+{
+    float result = (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z + v1.w*v2.w);
+    return result;
+}
+
+// Calculate distance between two vectors
+RMAPI float Vector4Distance(Vector4 v1, Vector4 v2)
+{
+    float result = sqrtf(
+        (v1.x - v2.x)*(v1.x - v2.x) + (v1.y - v2.y)*(v1.y - v2.y) +
+        (v1.z - v2.z)*(v1.z - v2.z) + (v1.w - v2.w)*(v1.w - v2.w));
+    return result;
+}
+
+// Calculate square distance between two vectors
+RMAPI float Vector4DistanceSqr(Vector4 v1, Vector4 v2)
+{
+    float result =
+        (v1.x - v2.x)*(v1.x - v2.x) + (v1.y - v2.y)*(v1.y - v2.y) +
+        (v1.z - v2.z)*(v1.z - v2.z) + (v1.w - v2.w)*(v1.w - v2.w);
+
+    return result;
+}
+
+RMAPI Vector4 Vector4Scale(Vector4 v, float scale)
+{
+    Vector4 result = { v.x*scale, v.y*scale, v.z*scale, v.w*scale };
+    return result;
+}
+
+// Multiply vector by vector
+RMAPI Vector4 Vector4Multiply(Vector4 v1, Vector4 v2)
+{
+    Vector4 result = { v1.x*v2.x, v1.y*v2.y, v1.z*v2.z, v1.w*v2.w };
+    return result;
+}
+
+// Negate vector
+RMAPI Vector4 Vector4Negate(Vector4 v)
+{
+    Vector4 result = { -v.x, -v.y, -v.z, -v.w };
+    return result;
+}
+
+// Divide vector by vector
+RMAPI Vector4 Vector4Divide(Vector4 v1, Vector4 v2)
+{
+    Vector4 result = { v1.x/v2.x, v1.y/v2.y, v1.z/v2.z, v1.w/v2.w };
+    return result;
+}
+
+// Normalize provided vector
+RMAPI Vector4 Vector4Normalize(Vector4 v)
+{
+    Vector4 result = { 0 };
+    float length = sqrtf((v.x*v.x) + (v.y*v.y) + (v.z*v.z) + (v.w*v.w));
+
+    if (length > 0)
+    {
+        float ilength = 1.0f/length;
+        result.x = v.x*ilength;
+        result.y = v.y*ilength;
+        result.z = v.z*ilength;
+        result.w = v.w*ilength;
+    }
+
+    return result;
+}
+
+// Get min value for each pair of components
+RMAPI Vector4 Vector4Min(Vector4 v1, Vector4 v2)
+{
+    Vector4 result = { 0 };
+
+    result.x = fminf(v1.x, v2.x);
+    result.y = fminf(v1.y, v2.y);
+    result.z = fminf(v1.z, v2.z);
+    result.w = fminf(v1.w, v2.w);
+
+    return result;
+}
+
+// Get max value for each pair of components
+RMAPI Vector4 Vector4Max(Vector4 v1, Vector4 v2)
+{
+    Vector4 result = { 0 };
+
+    result.x = fmaxf(v1.x, v2.x);
+    result.y = fmaxf(v1.y, v2.y);
+    result.z = fmaxf(v1.z, v2.z);
+    result.w = fmaxf(v1.w, v2.w);
+
+    return result;
+}
+
+// Calculate linear interpolation between two vectors
+RMAPI Vector4 Vector4Lerp(Vector4 v1, Vector4 v2, float amount)
+{
+    Vector4 result = { 0 };
+
+    result.x = v1.x + amount*(v2.x - v1.x);
+    result.y = v1.y + amount*(v2.y - v1.y);
+    result.z = v1.z + amount*(v2.z - v1.z);
+    result.w = v1.w + amount*(v2.w - v1.w);
+
+    return result;
+}
+
+// Move Vector towards target
+RMAPI Vector4 Vector4MoveTowards(Vector4 v, Vector4 target, float maxDistance)
+{
+    Vector4 result = { 0 };
+
+    float dx = target.x - v.x;
+    float dy = target.y - v.y;
+    float dz = target.z - v.z;
+    float dw = target.w - v.w;
+    float value = (dx*dx) + (dy*dy) + (dz*dz) + (dw*dw);
+
+    if ((value == 0) || ((maxDistance >= 0) && (value <= maxDistance*maxDistance))) return target;
+
+    float dist = sqrtf(value);
+
+    result.x = v.x + dx/dist*maxDistance;
+    result.y = v.y + dy/dist*maxDistance;
+    result.z = v.z + dz/dist*maxDistance;
+    result.w = v.w + dw/dist*maxDistance;
+
+    return result;
+}
+
+// Invert the given vector
+RMAPI Vector4 Vector4Invert(Vector4 v)
+{
+    Vector4 result = { 1.0f/v.x, 1.0f/v.y, 1.0f/v.z, 1.0f/v.w };
+    return result;
+}
+
+// Check whether two given vectors are almost equal
+RMAPI int Vector4Equals(Vector4 p, Vector4 q)
+{
+#if !defined(EPSILON)
+    #define EPSILON 0.000001f
+#endif
+
+    int result = ((fabsf(p.x - q.x)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.x), fabsf(q.x))))) &&
+                  ((fabsf(p.y - q.y)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.y), fabsf(q.y))))) &&
+                  ((fabsf(p.z - q.z)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.z), fabsf(q.z))))) &&
+                  ((fabsf(p.w - q.w)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.w), fabsf(q.w)))));
+    return result;
+}
+
+
 //----------------------------------------------------------------------------------
 // Module Functions Definition - Matrix math
 //----------------------------------------------------------------------------------
@@ -1524,32 +1839,32 @@ RMAPI Matrix MatrixScale(float x, float y, float z)
 }
 
 // Get perspective projection matrix
-RMAPI Matrix MatrixFrustum(double left, double right, double bottom, double top, double near, double far)
+RMAPI Matrix MatrixFrustum(double left, double right, double bottom, double top, double nearPlane, double farPlane)
 {
     Matrix result = { 0 };
 
     float rl = (float)(right - left);
     float tb = (float)(top - bottom);
-    float fn = (float)(far - near);
+    float fn = (float)(farPlane - nearPlane);
 
-    result.m0 = ((float)near*2.0f)/rl;
+    result.m0 = ((float)nearPlane*2.0f)/rl;
     result.m1 = 0.0f;
     result.m2 = 0.0f;
     result.m3 = 0.0f;
 
     result.m4 = 0.0f;
-    result.m5 = ((float)near*2.0f)/tb;
+    result.m5 = ((float)nearPlane*2.0f)/tb;
     result.m6 = 0.0f;
     result.m7 = 0.0f;
 
     result.m8 = ((float)right + (float)left)/rl;
     result.m9 = ((float)top + (float)bottom)/tb;
-    result.m10 = -((float)far + (float)near)/fn;
+    result.m10 = -((float)farPlane + (float)nearPlane)/fn;
     result.m11 = -1.0f;
 
     result.m12 = 0.0f;
     result.m13 = 0.0f;
-    result.m14 = -((float)far*(float)near*2.0f)/fn;
+    result.m14 = -((float)farPlane*(float)nearPlane*2.0f)/fn;
     result.m15 = 0.0f;
 
     return result;
@@ -1901,6 +2216,32 @@ RMAPI Quaternion QuaternionSlerp(Quaternion q1, Quaternion q2, float amount)
     return result;
 }
 
+// Calculate quaternion cubic spline interpolation using Cubic Hermite Spline algorithm
+// as described in the GLTF 2.0 specification: https://registry.khronos.org/glTF/specs/2.0/glTF-2.0.html#interpolation-cubic
+RMAPI Quaternion QuaternionCubicHermiteSpline(Quaternion q1, Quaternion outTangent1, Quaternion q2, Quaternion inTangent2, float t)
+{
+    float t2 = t*t;
+    float t3 = t2*t;
+    float h00 = 2*t3 - 3*t2 + 1;
+    float h10 = t3 - 2*t2 + t;
+    float h01 = -2*t3 + 3*t2;
+    float h11 = t3 - t2;
+
+    Quaternion p0 = QuaternionScale(q1, h00);
+    Quaternion m0 = QuaternionScale(outTangent1, h10);
+    Quaternion p1 = QuaternionScale(q2, h01);
+    Quaternion m1 = QuaternionScale(inTangent2, h11);
+
+    Quaternion result = { 0 };
+
+    result = QuaternionAdd(p0, m0);
+    result = QuaternionAdd(result, p1);
+    result = QuaternionAdd(result, m1);
+    result = QuaternionNormalize(result);
+
+    return result;
+}
+
 // Calculate quaternion based on the rotation from one vector to another
 RMAPI Quaternion QuaternionFromVector3ToVector3(Vector3 from, Vector3 to)
 {
@@ -1960,7 +2301,7 @@ RMAPI Quaternion QuaternionFromMatrix(Matrix mat)
     }
 
     float biggestVal = sqrtf(fourBiggestSquaredMinus1 + 1.0f)*0.5f;
-    float mult = 0.25f / biggestVal;
+    float mult = 0.25f/biggestVal;
 
     switch (biggestIndex)
     {
@@ -2042,8 +2383,7 @@ RMAPI Quaternion QuaternionFromAxisAngle(Vector3 axis, float angle)
         float ilength = 0.0f;
 
         // Vector3Normalize(axis)
-        Vector3 v = axis;
-        length = sqrtf(v.x*v.x + v.y*v.y + v.z*v.z);
+        length = axisLength;
         if (length == 0.0f) length = 1.0f;
         ilength = 1.0f/length;
         axis.x *= ilength;
@@ -2187,4 +2527,415 @@ RMAPI int QuaternionEquals(Quaternion p, Quaternion q)
     return result;
 }
 
+// Decompose a transformation matrix into its rotational, translational and scaling components
+RMAPI void MatrixDecompose(Matrix mat, Vector3 *translation, Quaternion *rotation, Vector3 *scale)
+{
+    // Extract translation.
+    translation->x = mat.m12;
+    translation->y = mat.m13;
+    translation->z = mat.m14;
+
+    // Extract upper-left for determinant computation
+    const float a = mat.m0;
+    const float b = mat.m4;
+    const float c = mat.m8;
+    const float d = mat.m1;
+    const float e = mat.m5;
+    const float f = mat.m9;
+    const float g = mat.m2;
+    const float h = mat.m6;
+    const float i = mat.m10;
+    const float A = e*i - f*h;
+    const float B = f*g - d*i;
+    const float C = d*h - e*g;
+
+    // Extract scale
+    const float det = a*A + b*B + c*C;
+    Vector3 abc = { a, b, c };
+    Vector3 def = { d, e, f };
+    Vector3 ghi = { g, h, i };
+
+    float scalex = Vector3Length(abc);
+    float scaley = Vector3Length(def);
+    float scalez = Vector3Length(ghi);
+    Vector3 s = { scalex, scaley, scalez };
+
+    if (det < 0) s = Vector3Negate(s);
+
+    *scale = s;
+
+    // Remove scale from the matrix if it is not close to zero
+    Matrix clone = mat;
+    if (!FloatEquals(det, 0))
+    {
+        clone.m0 /= s.x;
+        clone.m4 /= s.x;
+        clone.m8 /= s.x;
+        clone.m1 /= s.y;
+        clone.m5 /= s.y;
+        clone.m9 /= s.y;
+        clone.m2 /= s.z;
+        clone.m6 /= s.z;
+        clone.m10 /= s.z;
+
+        // Extract rotation
+        *rotation = QuaternionFromMatrix(clone);
+    }
+    else
+    {
+        // Set to identity if close to zero
+        *rotation = QuaternionIdentity();
+    }
+}
+
+#if defined(__cplusplus) && !defined(RAYMATH_DISABLE_CPP_OPERATORS)
+
+// Optional C++ math operators
+//-------------------------------------------------------------------------------
+
+// Vector2 operators
+static constexpr Vector2 Vector2Zeros = { 0, 0 };
+static constexpr Vector2 Vector2Ones = { 1, 1 };
+static constexpr Vector2 Vector2UnitX = { 1, 0 };
+static constexpr Vector2 Vector2UnitY = { 0, 1 };
+
+inline Vector2 operator + (const Vector2& lhs, const Vector2& rhs)
+{
+    return Vector2Add(lhs, rhs);
+}
+
+inline const Vector2& operator += (Vector2& lhs, const Vector2& rhs)
+{
+    lhs = Vector2Add(lhs, rhs);
+    return lhs;
+}
+
+inline Vector2 operator - (const Vector2& lhs, const Vector2& rhs)
+{
+    return Vector2Subtract(lhs, rhs);
+}
+
+inline const Vector2& operator -= (Vector2& lhs, const Vector2& rhs)
+{
+    lhs = Vector2Subtract(lhs, rhs);
+    return lhs;
+}
+
+inline Vector2 operator * (const Vector2& lhs, const float& rhs)
+{
+    return Vector2Scale(lhs, rhs);
+}
+
+inline const Vector2& operator *= (Vector2& lhs, const float& rhs)
+{
+    lhs = Vector2Scale(lhs, rhs);
+    return lhs;
+}
+
+inline Vector2 operator * (const Vector2& lhs, const Vector2& rhs)
+{
+    return Vector2Multiply(lhs, rhs);
+}
+
+inline const Vector2& operator *= (Vector2& lhs, const Vector2& rhs)
+{
+    lhs = Vector2Multiply(lhs, rhs);
+    return lhs;
+}
+
+inline Vector2 operator * (const Vector2& lhs, const Matrix& rhs)
+{
+    return Vector2Transform(lhs, rhs);
+}
+
+inline const Vector2& operator -= (Vector2& lhs, const Matrix& rhs)
+{
+    lhs = Vector2Transform(lhs, rhs);
+    return lhs;
+}
+
+inline Vector2 operator / (const Vector2& lhs, const float& rhs)
+{
+    return Vector2Scale(lhs, 1.0f / rhs);
+}
+
+inline const Vector2& operator /= (Vector2& lhs, const float& rhs)
+{
+    lhs = Vector2Scale(lhs, rhs);
+    return lhs;
+}
+
+inline Vector2 operator / (const Vector2& lhs, const Vector2& rhs)
+{
+    return Vector2Divide(lhs, rhs);
+}
+
+inline const Vector2& operator /= (Vector2& lhs, const Vector2& rhs)
+{
+    lhs = Vector2Divide(lhs, rhs);
+    return lhs;
+}
+
+inline bool operator == (const Vector2& lhs, const Vector2& rhs)
+{
+    return FloatEquals(lhs.x, rhs.x) && FloatEquals(lhs.y, rhs.y);
+}
+
+inline bool operator != (const Vector2& lhs, const Vector2& rhs)
+{
+    return !FloatEquals(lhs.x, rhs.x) || !FloatEquals(lhs.y, rhs.y);
+}
+
+// Vector3 operators
+static constexpr Vector3 Vector3Zeros = { 0, 0, 0 };
+static constexpr Vector3 Vector3Ones = { 1, 1, 1 };
+static constexpr Vector3 Vector3UnitX = { 1, 0, 0 };
+static constexpr Vector3 Vector3UnitY = { 0, 1, 0 };
+static constexpr Vector3 Vector3UnitZ = { 0, 0, 1 };
+
+inline Vector3 operator + (const Vector3& lhs, const Vector3& rhs)
+{
+    return Vector3Add(lhs, rhs);
+}
+
+inline const Vector3& operator += (Vector3& lhs, const Vector3& rhs)
+{
+    lhs = Vector3Add(lhs, rhs);
+    return lhs;
+}
+
+inline Vector3 operator - (const Vector3& lhs, const Vector3& rhs)
+{
+    return Vector3Subtract(lhs, rhs);
+}
+
+inline const Vector3& operator -= (Vector3& lhs, const Vector3& rhs)
+{
+    lhs = Vector3Subtract(lhs, rhs);
+    return lhs;
+}
+
+inline Vector3 operator * (const Vector3& lhs, const float& rhs)
+{
+    return Vector3Scale(lhs, rhs);
+}
+
+inline const Vector3& operator *= (Vector3& lhs, const float& rhs)
+{
+    lhs = Vector3Scale(lhs, rhs);
+    return lhs;
+}
+
+inline Vector3 operator * (const Vector3& lhs, const Vector3& rhs)
+{
+    return Vector3Multiply(lhs, rhs);
+}
+
+inline const Vector3& operator *= (Vector3& lhs, const Vector3& rhs)
+{
+    lhs = Vector3Multiply(lhs, rhs);
+    return lhs;
+}
+
+inline Vector3 operator * (const Vector3& lhs, const Matrix& rhs)
+{
+    return Vector3Transform(lhs, rhs);
+}
+
+inline const Vector3& operator -= (Vector3& lhs, const Matrix& rhs)
+{
+    lhs = Vector3Transform(lhs, rhs);
+    return lhs;
+}
+
+inline Vector3 operator / (const Vector3& lhs, const float& rhs)
+{
+    return Vector3Scale(lhs, 1.0f / rhs);
+}
+
+inline const Vector3& operator /= (Vector3& lhs, const float& rhs)
+{
+    lhs = Vector3Scale(lhs, rhs);
+    return lhs;
+}
+
+inline Vector3 operator / (const Vector3& lhs, const Vector3& rhs)
+{
+    return Vector3Divide(lhs, rhs);
+}
+
+inline const Vector3& operator /= (Vector3& lhs, const Vector3& rhs)
+{
+    lhs = Vector3Divide(lhs, rhs);
+    return lhs;
+}
+
+inline bool operator == (const Vector3& lhs, const Vector3& rhs)
+{
+    return FloatEquals(lhs.x, rhs.x) && FloatEquals(lhs.y, rhs.y) && FloatEquals(lhs.z, rhs.z);
+}
+
+inline bool operator != (const Vector3& lhs, const Vector3& rhs)
+{
+    return !FloatEquals(lhs.x, rhs.x) || !FloatEquals(lhs.y, rhs.y) || !FloatEquals(lhs.z, rhs.z);
+}
+
+// Vector4 operators
+static constexpr Vector4 Vector4Zeros = { 0, 0, 0, 0 };
+static constexpr Vector4 Vector4Ones = { 1, 1, 1, 1 };
+static constexpr Vector4 Vector4UnitX = { 1, 0, 0, 0 };
+static constexpr Vector4 Vector4UnitY = { 0, 1, 0, 0 };
+static constexpr Vector4 Vector4UnitZ = { 0, 0, 1, 0 };
+static constexpr Vector4 Vector4UnitW = { 0, 0, 0, 1 };
+
+inline Vector4 operator + (const Vector4& lhs, const Vector4& rhs)
+{
+    return Vector4Add(lhs, rhs);
+}
+
+inline const Vector4& operator += (Vector4& lhs, const Vector4& rhs)
+{
+    lhs = Vector4Add(lhs, rhs);
+    return lhs;
+}
+
+inline Vector4 operator - (const Vector4& lhs, const Vector4& rhs)
+{
+    return Vector4Subtract(lhs, rhs);
+}
+
+inline const Vector4& operator -= (Vector4& lhs, const Vector4& rhs)
+{
+    lhs = Vector4Subtract(lhs, rhs);
+    return lhs;
+}
+
+inline Vector4 operator * (const Vector4& lhs, const float& rhs)
+{
+    return Vector4Scale(lhs, rhs);
+}
+
+inline const Vector4& operator *= (Vector4& lhs, const float& rhs)
+{
+    lhs = Vector4Scale(lhs, rhs);
+    return lhs;
+}
+
+inline Vector4 operator * (const Vector4& lhs, const Vector4& rhs)
+{
+    return Vector4Multiply(lhs, rhs);
+}
+
+inline const Vector4& operator *= (Vector4& lhs, const Vector4& rhs)
+{
+    lhs = Vector4Multiply(lhs, rhs);
+    return lhs;
+}
+
+inline Vector4 operator / (const Vector4& lhs, const float& rhs)
+{
+    return Vector4Scale(lhs, 1.0f / rhs);
+}
+
+inline const Vector4& operator /= (Vector4& lhs, const float& rhs)
+{
+    lhs = Vector4Scale(lhs, rhs);
+    return lhs;
+}
+
+inline Vector4 operator / (const Vector4& lhs, const Vector4& rhs)
+{
+    return Vector4Divide(lhs, rhs);
+}
+
+inline const Vector4& operator /= (Vector4& lhs, const Vector4& rhs)
+{
+    lhs = Vector4Divide(lhs, rhs);
+    return lhs;
+}
+
+inline bool operator == (const Vector4& lhs, const Vector4& rhs)
+{
+    return FloatEquals(lhs.x, rhs.x) && FloatEquals(lhs.y, rhs.y) && FloatEquals(lhs.z, rhs.z) && FloatEquals(lhs.w, rhs.w);
+}
+
+inline bool operator != (const Vector4& lhs, const Vector4& rhs)
+{
+    return !FloatEquals(lhs.x, rhs.x) || !FloatEquals(lhs.y, rhs.y) || !FloatEquals(lhs.z, rhs.z) || !FloatEquals(lhs.w, rhs.w);
+}
+
+// Quaternion operators
+static constexpr Quaternion QuaternionZeros = { 0, 0, 0, 0 };
+static constexpr Quaternion QuaternionOnes = { 1, 1, 1, 1 };
+static constexpr Quaternion QuaternionUnitX = { 0, 0, 0, 1 };
+
+inline Quaternion operator + (const Quaternion& lhs, const float& rhs)
+{
+    return QuaternionAddValue(lhs, rhs);
+}
+
+inline const Quaternion& operator += (Quaternion& lhs, const float& rhs)
+{
+    lhs = QuaternionAddValue(lhs, rhs);
+    return lhs;
+}
+
+inline Quaternion operator - (const Quaternion& lhs, const float& rhs)
+{
+    return QuaternionSubtractValue(lhs, rhs);
+}
+
+inline const Quaternion& operator -= (Quaternion& lhs, const float& rhs)
+{
+    lhs = QuaternionSubtractValue(lhs, rhs);
+    return lhs;
+}
+
+inline Quaternion operator * (const Quaternion& lhs, const Matrix& rhs)
+{
+    return QuaternionTransform(lhs, rhs);
+}
+
+inline const Quaternion& operator *= (Quaternion& lhs, const Matrix& rhs)
+{
+    lhs = QuaternionTransform(lhs, rhs);
+    return lhs;
+}
+
+// Matrix operators
+inline Matrix operator + (const Matrix& lhs, const Matrix& rhs)
+{
+    return MatrixAdd(lhs, rhs);
+}
+
+inline const Matrix& operator += (Matrix& lhs, const Matrix& rhs)
+{
+    lhs = MatrixAdd(lhs, rhs);
+    return lhs;
+}
+
+inline Matrix operator - (const Matrix& lhs, const Matrix& rhs)
+{
+    return MatrixSubtract(lhs, rhs);
+}
+
+inline const Matrix& operator -= (Matrix& lhs, const Matrix& rhs)
+{
+    lhs = MatrixSubtract(lhs, rhs);
+    return lhs;
+}
+
+inline Matrix operator * (const Matrix& lhs, const Matrix& rhs)
+{
+    return MatrixMultiply(lhs, rhs);
+}
+
+inline const Matrix& operator *= (Matrix& lhs, const Matrix& rhs)
+{
+    lhs = MatrixMultiply(lhs, rhs);
+    return lhs;
+}
+//-------------------------------------------------------------------------------
+#endif  // C++ operators
+
 #endif  // RAYMATH_H