1/**********************************************************************************************
   2*
   3*   raymath v2.0 - Math functions to work with Vector2, Vector3, Matrix and Quaternions
   4*
   5*   CONVENTIONS:
   6*     - Matrix structure is defined as row-major (memory layout) but parameters naming AND all
   7*       math operations performed by the library consider the structure as it was column-major
   8*       It is like transposed versions of the matrices are used for all the maths
   9*       It benefits some functions making them cache-friendly and also avoids matrix
  10*       transpositions sometimes required by OpenGL
  11*       Example: In memory order, row0 is [m0 m4 m8 m12] but in semantic math row0 is [m0 m1 m2 m3]
  12*     - Functions are always self-contained, no function use another raymath function inside,
  13*       required code is directly re-implemented inside
  14*     - Functions input parameters are always received by value (2 unavoidable exceptions)
  15*     - Functions use always a "result" variable for return (except C++ operators)
  16*     - Functions are always defined inline
  17*     - Angles are always in radians (DEG2RAD/RAD2DEG macros provided for convenience)
  18*     - No compound literals used to make sure the library is compatible with C++
  19*
  20*   CONFIGURATION:
  21*       #define RAYMATH_IMPLEMENTATION
  22*           Generates the implementation of the library into the included file
  23*           If not defined, the library is in header only mode and can be included in other headers
  24*           or source files without problems. But only ONE file should hold the implementation
  25*
  26*       #define RAYMATH_STATIC_INLINE
  27*           Define static inline functions code, so #include header suffices for use
  28*           This may use up lots of memory
  29*
  30*       #define RAYMATH_DISABLE_CPP_OPERATORS
  31*           Disables C++ operator overloads for raymath types.
  32*
  33*       #define RAYMATH_USE_SIMD_INTRINSICS   1
  34*           Try to enable SIMD intrinsics for MatrixMultiply()
  35*           Note that users enabling it must be aware of the target platform where application will
  36*           run to support the selected SIMD intrinsic, for now, only SSE is supported
  37*
  38*   LICENSE: zlib/libpng
  39*
  40*   Copyright (c) 2015-2026 Ramon Santamaria (@raysan5)
  41*
  42*   This software is provided "as-is", without any express or implied warranty. In no event
  43*   will the authors be held liable for any damages arising from the use of this software.
  44*
  45*   Permission is granted to anyone to use this software for any purpose, including commercial
  46*   applications, and to alter it and redistribute it freely, subject to the following restrictions:
  47*
  48*     1. The origin of this software must not be misrepresented; you must not claim that you
  49*     wrote the original software. If you use this software in a product, an acknowledgment
  50*     in the product documentation would be appreciated but is not required.
  51*
  52*     2. Altered source versions must be plainly marked as such, and must not be misrepresented
  53*     as being the original software.
  54*
  55*     3. This notice may not be removed or altered from any source distribution.
  56*
  57**********************************************************************************************/
  58
  59#ifndef RAYMATH_H
  60#define RAYMATH_H
  61
  62#if defined(RAYMATH_IMPLEMENTATION) && defined(RAYMATH_STATIC_INLINE)
  63    #error "Specifying both RAYMATH_IMPLEMENTATION and RAYMATH_STATIC_INLINE is contradictory"
  64#endif
  65
  66// Function specifiers definition
  67#if defined(RAYMATH_IMPLEMENTATION)
  68    #if defined(_WIN32) && defined(BUILD_LIBTYPE_SHARED)
  69        #define RMAPI __declspec(dllexport) extern inline    // Building raylib as a Win32 shared library (.dll)
  70    #elif defined(BUILD_LIBTYPE_SHARED)
  71        #define RMAPI __attribute__((visibility("default"))) // Building raylib as a Unix shared library (.so/.dylib)
  72    #elif defined(_WIN32) && defined(USE_LIBTYPE_SHARED)
  73        #define RMAPI __declspec(dllimport)                  // Using raylib as a Win32 shared library (.dll)
  74    #else
  75        #define RMAPI extern inline // Provide external definition
  76    #endif
  77#elif defined(RAYMATH_STATIC_INLINE)
  78    #define RMAPI static inline // Functions may be inlined, no external out-of-line definition
  79#else
  80    #if defined(__TINYC__)
  81        #define RMAPI static inline // plain inline not supported by tinycc (See issue #435)
  82    #else
  83        #define RMAPI inline        // Functions may be inlined or external definition used
  84    #endif
  85#endif
  86
  87//----------------------------------------------------------------------------------
  88// Defines and Macros
  89//----------------------------------------------------------------------------------
  90#ifndef PI
  91    #define PI 3.14159265358979323846f
  92#endif
  93
  94#ifndef EPSILON
  95    #define EPSILON 0.000001f
  96#endif
  97
  98#ifndef DEG2RAD
  99    #define DEG2RAD (PI/180.0f)
 100#endif
 101
 102#ifndef RAD2DEG
 103    #define RAD2DEG (180.0f/PI)
 104#endif
 105
 106// Get float vector for Matrix
 107#ifndef MatrixToFloat
 108    #define MatrixToFloat(mat) (MatrixToFloatV(mat).v)
 109#endif
 110
 111// Get float vector for Vector3
 112#ifndef Vector3ToFloat
 113    #define Vector3ToFloat(vec) (Vector3ToFloatV(vec).v)
 114#endif
 115
 116//----------------------------------------------------------------------------------
 117// Types and Structures Definition
 118//----------------------------------------------------------------------------------
 119#if !defined(RL_VECTOR2_TYPE)
 120// Vector2 type
 121typedef struct Vector2 {
 122    float x;
 123    float y;
 124} Vector2;
 125#define RL_VECTOR2_TYPE
 126#endif
 127
 128#if !defined(RL_VECTOR3_TYPE)
 129// Vector3 type
 130typedef struct Vector3 {
 131    float x;
 132    float y;
 133    float z;
 134} Vector3;
 135#define RL_VECTOR3_TYPE
 136#endif
 137
 138#if !defined(RL_VECTOR4_TYPE)
 139// Vector4 type
 140typedef struct Vector4 {
 141    float x;
 142    float y;
 143    float z;
 144    float w;
 145} Vector4;
 146#define RL_VECTOR4_TYPE
 147#endif
 148
 149#if !defined(RL_QUATERNION_TYPE)
 150// Quaternion type
 151typedef Vector4 Quaternion;
 152#define RL_QUATERNION_TYPE
 153#endif
 154
 155#if !defined(RL_MATRIX_TYPE)
 156// Matrix type (OpenGL style 4x4 - right handed, column major)
 157typedef struct Matrix {
 158    float m0, m4, m8, m12;      // Matrix first row (4 components)
 159    float m1, m5, m9, m13;      // Matrix second row (4 components)
 160    float m2, m6, m10, m14;     // Matrix third row (4 components)
 161    float m3, m7, m11, m15;     // Matrix fourth row (4 components)
 162} Matrix;
 163#define RL_MATRIX_TYPE
 164#endif
 165
 166// NOTE: Helper types to be used instead of array return types for *ToFloat functions
 167#if !defined(RL_FLOAT3_TYPE)
 168typedef struct float3 {
 169    float v[3];
 170} float3;
 171#define RL_FLOAT3_TYPE
 172#endif
 173
 174#if !defined(RL_FLOAT16_TYPE)
 175typedef struct float16 {
 176    float v[16];
 177} float16;
 178#define RL_FLOAT16_TYPE
 179#endif
 180
 181#include <math.h>       // Required for: sinf(), cosf(), tan(), atan2f(), sqrtf(), floor(), fminf(), fmaxf(), fabsf()
 182
 183#if RAYMATH_USE_SIMD_INTRINSICS
 184    // SIMD is used on the most costly raymath function MatrixMultiply()
 185    // NOTE: Only SSE intrinsics support implemented
 186    // TODO: Consider support for other SIMD intrinsics:
 187    //  - SSEx, AVX, AVX2, FMA, NEON, RVV
 188    /*
 189    #if defined(__SSE4_2__)
 190        #include <nmmintrin.h>
 191        #define RAYMATH_SSE42_ENABLED
 192    #elif defined(__SSE4_1__)
 193        #include <smmintrin.h>
 194        #define RAYMATH_SSE41_ENABLED
 195    #elif defined(__SSSE3__)
 196        #include <tmmintrin.h>
 197        #define RAYMATH_SSSE3_ENABLED
 198    #elif defined(__SSE3__)
 199        #include <pmmintrin.h>
 200        #define RAYMATH_SSE3_ENABLED
 201    #elif defined(__SSE2__) || (defined(_M_AMD64) || defined(_M_X64)) // SSE2 x64
 202        #include <emmintrin.h>
 203        #define RAYMATH_SSE2_ENABLED
 204    #endif
 205    */
 206    #if defined(__SSE__) || defined(_M_X64) || (defined(_M_IX86_FP) && (_M_IX86_FP >= 1))
 207        #include <xmmintrin.h>
 208        #define RAYMATH_SSE_ENABLED
 209    #endif
 210#endif
 211
 212//----------------------------------------------------------------------------------
 213// Module Functions Definition - Utils math
 214//----------------------------------------------------------------------------------
 215
 216// Clamp float value
 217RMAPI float Clamp(float value, float min, float max)
 218{
 219    float result = (value < min)? min : value;
 220
 221    if (result > max) result = max;
 222
 223    return result;
 224}
 225
 226// Calculate linear interpolation between two floats
 227RMAPI float Lerp(float start, float end, float amount)
 228{
 229    float result = start + amount*(end - start);
 230
 231    return result;
 232}
 233
 234// Normalize input value within input range
 235RMAPI float Normalize(float value, float start, float end)
 236{
 237    float result = (value - start)/(end - start);
 238
 239    return result;
 240}
 241
 242// Remap input value within input range to output range
 243RMAPI float Remap(float value, float inputStart, float inputEnd, float outputStart, float outputEnd)
 244{
 245    float result = (value - inputStart)/(inputEnd - inputStart)*(outputEnd - outputStart) + outputStart;
 246
 247    return result;
 248}
 249
 250// Wrap input value from min to max
 251RMAPI float Wrap(float value, float min, float max)
 252{
 253    float result = value - (max - min)*floorf((value - min)/(max - min));
 254
 255    return result;
 256}
 257
 258// Check whether two given floats are almost equal
 259RMAPI int FloatEquals(float x, float y)
 260{
 261#if !defined(EPSILON)
 262    #define EPSILON 0.000001f
 263#endif
 264
 265    int result = (fabsf(x - y)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(x), fabsf(y))));
 266
 267    return result;
 268}
 269
 270//----------------------------------------------------------------------------------
 271// Module Functions Definition - Vector2 math
 272//----------------------------------------------------------------------------------
 273
 274// Vector with components value 0.0f
 275RMAPI Vector2 Vector2Zero(void)
 276{
 277    Vector2 result = { 0.0f, 0.0f };
 278
 279    return result;
 280}
 281
 282// Vector with components value 1.0f
 283RMAPI Vector2 Vector2One(void)
 284{
 285    Vector2 result = { 1.0f, 1.0f };
 286
 287    return result;
 288}
 289
 290// Add two vectors (v1 + v2)
 291RMAPI Vector2 Vector2Add(Vector2 v1, Vector2 v2)
 292{
 293    Vector2 result = { v1.x + v2.x, v1.y + v2.y };
 294
 295    return result;
 296}
 297
 298// Add vector and float value
 299RMAPI Vector2 Vector2AddValue(Vector2 v, float add)
 300{
 301    Vector2 result = { v.x + add, v.y + add };
 302
 303    return result;
 304}
 305
 306// Subtract two vectors (v1 - v2)
 307RMAPI Vector2 Vector2Subtract(Vector2 v1, Vector2 v2)
 308{
 309    Vector2 result = { v1.x - v2.x, v1.y - v2.y };
 310
 311    return result;
 312}
 313
 314// Subtract vector by float value
 315RMAPI Vector2 Vector2SubtractValue(Vector2 v, float sub)
 316{
 317    Vector2 result = { v.x - sub, v.y - sub };
 318
 319    return result;
 320}
 321
 322// Calculate vector length
 323RMAPI float Vector2Length(Vector2 v)
 324{
 325    float result = sqrtf((v.x*v.x) + (v.y*v.y));
 326
 327    return result;
 328}
 329
 330// Calculate vector square length
 331RMAPI float Vector2LengthSqr(Vector2 v)
 332{
 333    float result = (v.x*v.x) + (v.y*v.y);
 334
 335    return result;
 336}
 337
 338// Calculate two vectors dot product
 339RMAPI float Vector2DotProduct(Vector2 v1, Vector2 v2)
 340{
 341    float result = (v1.x*v2.x + v1.y*v2.y);
 342
 343    return result;
 344}
 345
 346// Calculate two vectors cross product
 347RMAPI float Vector2CrossProduct(Vector2 v1, Vector2 v2)
 348{
 349    float result = (v1.x*v2.y - v1.y*v2.x);
 350
 351    return result;
 352}
 353
 354// Calculate distance between two vectors
 355RMAPI float Vector2Distance(Vector2 v1, Vector2 v2)
 356{
 357    float result = sqrtf((v1.x - v2.x)*(v1.x - v2.x) + (v1.y - v2.y)*(v1.y - v2.y));
 358
 359    return result;
 360}
 361
 362// Calculate square distance between two vectors
 363RMAPI float Vector2DistanceSqr(Vector2 v1, Vector2 v2)
 364{
 365    float result = ((v1.x - v2.x)*(v1.x - v2.x) + (v1.y - v2.y)*(v1.y - v2.y));
 366
 367    return result;
 368}
 369
 370// Calculate the signed angle from v1 to v2, relative to the origin (0, 0)
 371// NOTE: Coordinate system convention: positive X right, positive Y down
 372// positive angles appear clockwise, and negative angles appear counterclockwise
 373RMAPI float Vector2Angle(Vector2 v1, Vector2 v2)
 374{
 375    float result = 0.0f;
 376
 377    float dot = v1.x*v2.x + v1.y*v2.y;
 378    float det = v1.x*v2.y - v1.y*v2.x;
 379
 380    result = atan2f(det, dot);
 381
 382    return result;
 383}
 384
 385// Calculate angle defined by a two vectors line
 386// NOTE: Parameters need to be normalized
 387// Current implementation should be aligned with glm::angle
 388RMAPI float Vector2LineAngle(Vector2 start, Vector2 end)
 389{
 390    float result = 0.0f;
 391
 392    // TODO(10/9/2023): Currently angles move clockwise, determine if this is wanted behavior
 393    result = -atan2f(end.y - start.y, end.x - start.x);
 394
 395    return result;
 396}
 397
 398// Scale vector (multiply by value)
 399RMAPI Vector2 Vector2Scale(Vector2 v, float scale)
 400{
 401    Vector2 result = { v.x*scale, v.y*scale };
 402
 403    return result;
 404}
 405
 406// Multiply vector by vector
 407RMAPI Vector2 Vector2Multiply(Vector2 v1, Vector2 v2)
 408{
 409    Vector2 result = { v1.x*v2.x, v1.y*v2.y };
 410
 411    return result;
 412}
 413
 414// Negate vector
 415RMAPI Vector2 Vector2Negate(Vector2 v)
 416{
 417    Vector2 result = { -v.x, -v.y };
 418
 419    return result;
 420}
 421
 422// Divide vector by vector
 423RMAPI Vector2 Vector2Divide(Vector2 v1, Vector2 v2)
 424{
 425    Vector2 result = { v1.x/v2.x, v1.y/v2.y };
 426
 427    return result;
 428}
 429
 430// Normalize provided vector
 431RMAPI Vector2 Vector2Normalize(Vector2 v)
 432{
 433    Vector2 result = { 0 };
 434    float length = sqrtf((v.x*v.x) + (v.y*v.y));
 435
 436    if (length > 0)
 437    {
 438        float ilength = 1.0f/length;
 439        result.x = v.x*ilength;
 440        result.y = v.y*ilength;
 441    }
 442
 443    return result;
 444}
 445
 446// Transforms a Vector2 by a given Matrix
 447RMAPI Vector2 Vector2Transform(Vector2 v, Matrix mat)
 448{
 449    Vector2 result = { 0 };
 450
 451    float x = v.x;
 452    float y = v.y;
 453    float z = 0;
 454
 455    result.x = mat.m0*x + mat.m4*y + mat.m8*z + mat.m12;
 456    result.y = mat.m1*x + mat.m5*y + mat.m9*z + mat.m13;
 457
 458    return result;
 459}
 460
 461// Calculate linear interpolation between two vectors
 462RMAPI Vector2 Vector2Lerp(Vector2 v1, Vector2 v2, float amount)
 463{
 464    Vector2 result = { 0 };
 465
 466    result.x = v1.x + amount*(v2.x - v1.x);
 467    result.y = v1.y + amount*(v2.y - v1.y);
 468
 469    return result;
 470}
 471
 472// Calculate reflected vector to normal
 473RMAPI Vector2 Vector2Reflect(Vector2 v, Vector2 normal)
 474{
 475    Vector2 result = { 0 };
 476
 477    float dotProduct = (v.x*normal.x + v.y*normal.y); // Dot product
 478
 479    result.x = v.x - (2.0f*normal.x)*dotProduct;
 480    result.y = v.y - (2.0f*normal.y)*dotProduct;
 481
 482    return result;
 483}
 484
 485// Get min value for each pair of components
 486RMAPI Vector2 Vector2Min(Vector2 v1, Vector2 v2)
 487{
 488    Vector2 result = { 0 };
 489
 490    result.x = fminf(v1.x, v2.x);
 491    result.y = fminf(v1.y, v2.y);
 492
 493    return result;
 494}
 495
 496// Get max value for each pair of components
 497RMAPI Vector2 Vector2Max(Vector2 v1, Vector2 v2)
 498{
 499    Vector2 result = { 0 };
 500
 501    result.x = fmaxf(v1.x, v2.x);
 502    result.y = fmaxf(v1.y, v2.y);
 503
 504    return result;
 505}
 506
 507// Rotate vector by angle
 508RMAPI Vector2 Vector2Rotate(Vector2 v, float angle)
 509{
 510    Vector2 result = { 0 };
 511
 512    float cosres = cosf(angle);
 513    float sinres = sinf(angle);
 514
 515    result.x = v.x*cosres - v.y*sinres;
 516    result.y = v.x*sinres + v.y*cosres;
 517
 518    return result;
 519}
 520
 521// Move Vector towards target
 522RMAPI Vector2 Vector2MoveTowards(Vector2 v, Vector2 target, float maxDistance)
 523{
 524    Vector2 result = { 0 };
 525
 526    float dx = target.x - v.x;
 527    float dy = target.y - v.y;
 528    float value = (dx*dx) + (dy*dy);
 529
 530    if ((value == 0) || ((maxDistance >= 0) && (value <= maxDistance*maxDistance))) return target;
 531
 532    float dist = sqrtf(value);
 533
 534    result.x = v.x + dx/dist*maxDistance;
 535    result.y = v.y + dy/dist*maxDistance;
 536
 537    return result;
 538}
 539
 540// Invert the given vector
 541RMAPI Vector2 Vector2Invert(Vector2 v)
 542{
 543    Vector2 result = { 1.0f/v.x, 1.0f/v.y };
 544
 545    return result;
 546}
 547
 548// Clamp the components of the vector between
 549// min and max values specified by the given vectors
 550RMAPI Vector2 Vector2Clamp(Vector2 v, Vector2 min, Vector2 max)
 551{
 552    Vector2 result = { 0 };
 553
 554    result.x = fminf(max.x, fmaxf(min.x, v.x));
 555    result.y = fminf(max.y, fmaxf(min.y, v.y));
 556
 557    return result;
 558}
 559
 560// Clamp the magnitude of the vector between two min and max values
 561RMAPI Vector2 Vector2ClampValue(Vector2 v, float min, float max)
 562{
 563    Vector2 result = v;
 564
 565    float length = (v.x*v.x) + (v.y*v.y);
 566    if (length > 0.0f)
 567    {
 568        length = sqrtf(length);
 569
 570        float scale = 1; // By default, 1 as the neutral element
 571        if (length < min) scale = min/length;
 572        else if (length > max) scale = max/length;
 573
 574        result.x = v.x*scale;
 575        result.y = v.y*scale;
 576    }
 577
 578    return result;
 579}
 580
 581// Check whether two given vectors are almost equal
 582RMAPI int Vector2Equals(Vector2 p, Vector2 q)
 583{
 584#if !defined(EPSILON)
 585    #define EPSILON 0.000001f
 586#endif
 587
 588    int result = ((fabsf(p.x - q.x)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.x), fabsf(q.x))))) &&
 589                  ((fabsf(p.y - q.y)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.y), fabsf(q.y)))));
 590
 591    return result;
 592}
 593
 594// Compute the direction of a refracted ray
 595// v: normalized direction of the incoming ray
 596// n: normalized normal vector of the interface of two optical media
 597// r: ratio of the refractive index of the medium from where the ray comes
 598// to the refractive index of the medium on the other side of the surface
 599RMAPI Vector2 Vector2Refract(Vector2 v, Vector2 n, float r)
 600{
 601    Vector2 result = { 0 };
 602
 603    float dot = v.x*n.x + v.y*n.y;
 604    float d = 1.0f - r*r*(1.0f - dot*dot);
 605
 606    if (d >= 0.0f)
 607    {
 608        d = sqrtf(d);
 609        v.x = r*v.x - (r*dot + d)*n.x;
 610        v.y = r*v.y - (r*dot + d)*n.y;
 611
 612        result = v;
 613    }
 614
 615    return result;
 616}
 617
 618
 619//----------------------------------------------------------------------------------
 620// Module Functions Definition - Vector3 math
 621//----------------------------------------------------------------------------------
 622
 623// Vector with components value 0.0f
 624RMAPI Vector3 Vector3Zero(void)
 625{
 626    Vector3 result = { 0.0f, 0.0f, 0.0f };
 627
 628    return result;
 629}
 630
 631// Vector with components value 1.0f
 632RMAPI Vector3 Vector3One(void)
 633{
 634    Vector3 result = { 1.0f, 1.0f, 1.0f };
 635
 636    return result;
 637}
 638
 639// Add two vectors
 640RMAPI Vector3 Vector3Add(Vector3 v1, Vector3 v2)
 641{
 642    Vector3 result = { v1.x + v2.x, v1.y + v2.y, v1.z + v2.z };
 643
 644    return result;
 645}
 646
 647// Add vector and float value
 648RMAPI Vector3 Vector3AddValue(Vector3 v, float add)
 649{
 650    Vector3 result = { v.x + add, v.y + add, v.z + add };
 651
 652    return result;
 653}
 654
 655// Subtract two vectors
 656RMAPI Vector3 Vector3Subtract(Vector3 v1, Vector3 v2)
 657{
 658    Vector3 result = { v1.x - v2.x, v1.y - v2.y, v1.z - v2.z };
 659
 660    return result;
 661}
 662
 663// Subtract vector by float value
 664RMAPI Vector3 Vector3SubtractValue(Vector3 v, float sub)
 665{
 666    Vector3 result = { v.x - sub, v.y - sub, v.z - sub };
 667
 668    return result;
 669}
 670
 671// Multiply vector by scalar
 672RMAPI Vector3 Vector3Scale(Vector3 v, float scalar)
 673{
 674    Vector3 result = { v.x*scalar, v.y*scalar, v.z*scalar };
 675
 676    return result;
 677}
 678
 679// Multiply vector by vector
 680RMAPI Vector3 Vector3Multiply(Vector3 v1, Vector3 v2)
 681{
 682    Vector3 result = { v1.x*v2.x, v1.y*v2.y, v1.z*v2.z };
 683
 684    return result;
 685}
 686
 687// Calculate two vectors cross product
 688RMAPI Vector3 Vector3CrossProduct(Vector3 v1, Vector3 v2)
 689{
 690    Vector3 result = { v1.y*v2.z - v1.z*v2.y, v1.z*v2.x - v1.x*v2.z, v1.x*v2.y - v1.y*v2.x };
 691
 692    return result;
 693}
 694
 695// Calculate one vector perpendicular vector
 696RMAPI Vector3 Vector3Perpendicular(Vector3 v)
 697{
 698    Vector3 result = { 0 };
 699
 700    float min = fabsf(v.x);
 701    Vector3 cardinalAxis = {1.0f, 0.0f, 0.0f};
 702
 703    if (fabsf(v.y) < min)
 704    {
 705        min = fabsf(v.y);
 706        Vector3 tmp = {0.0f, 1.0f, 0.0f};
 707        cardinalAxis = tmp;
 708    }
 709
 710    if (fabsf(v.z) < min)
 711    {
 712        Vector3 tmp = {0.0f, 0.0f, 1.0f};
 713        cardinalAxis = tmp;
 714    }
 715
 716    // Cross product between vectors
 717    result.x = v.y*cardinalAxis.z - v.z*cardinalAxis.y;
 718    result.y = v.z*cardinalAxis.x - v.x*cardinalAxis.z;
 719    result.z = v.x*cardinalAxis.y - v.y*cardinalAxis.x;
 720
 721    return result;
 722}
 723
 724// Calculate vector length
 725RMAPI float Vector3Length(const Vector3 v)
 726{
 727    float result = sqrtf(v.x*v.x + v.y*v.y + v.z*v.z);
 728
 729    return result;
 730}
 731
 732// Calculate vector square length
 733RMAPI float Vector3LengthSqr(const Vector3 v)
 734{
 735    float result = v.x*v.x + v.y*v.y + v.z*v.z;
 736
 737    return result;
 738}
 739
 740// Calculate two vectors dot product
 741RMAPI float Vector3DotProduct(Vector3 v1, Vector3 v2)
 742{
 743    float result = (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z);
 744
 745    return result;
 746}
 747
 748// Calculate distance between two vectors
 749RMAPI float Vector3Distance(Vector3 v1, Vector3 v2)
 750{
 751    float result = 0.0f;
 752
 753    float dx = v2.x - v1.x;
 754    float dy = v2.y - v1.y;
 755    float dz = v2.z - v1.z;
 756    result = sqrtf(dx*dx + dy*dy + dz*dz);
 757
 758    return result;
 759}
 760
 761// Calculate square distance between two vectors
 762RMAPI float Vector3DistanceSqr(Vector3 v1, Vector3 v2)
 763{
 764    float result = 0.0f;
 765
 766    float dx = v2.x - v1.x;
 767    float dy = v2.y - v1.y;
 768    float dz = v2.z - v1.z;
 769    result = dx*dx + dy*dy + dz*dz;
 770
 771    return result;
 772}
 773
 774// Calculate angle between two vectors
 775RMAPI float Vector3Angle(Vector3 v1, Vector3 v2)
 776{
 777    float result = 0.0f;
 778
 779    Vector3 cross = { v1.y*v2.z - v1.z*v2.y, v1.z*v2.x - v1.x*v2.z, v1.x*v2.y - v1.y*v2.x };
 780    float len = sqrtf(cross.x*cross.x + cross.y*cross.y + cross.z*cross.z);
 781    float dot = (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z);
 782    result = atan2f(len, dot);
 783
 784    return result;
 785}
 786
 787// Negate provided vector (invert direction)
 788RMAPI Vector3 Vector3Negate(Vector3 v)
 789{
 790    Vector3 result = { -v.x, -v.y, -v.z };
 791
 792    return result;
 793}
 794
 795// Divide vector by vector
 796RMAPI Vector3 Vector3Divide(Vector3 v1, Vector3 v2)
 797{
 798    Vector3 result = { v1.x/v2.x, v1.y/v2.y, v1.z/v2.z };
 799
 800    return result;
 801}
 802
 803// Normalize provided vector
 804RMAPI Vector3 Vector3Normalize(Vector3 v)
 805{
 806    Vector3 result = v;
 807
 808    float length = sqrtf(v.x*v.x + v.y*v.y + v.z*v.z);
 809    if (length != 0.0f)
 810    {
 811        float ilength = 1.0f/length;
 812
 813        result.x *= ilength;
 814        result.y *= ilength;
 815        result.z *= ilength;
 816    }
 817
 818    return result;
 819}
 820
 821//Calculate the projection of the vector v1 on to v2
 822RMAPI Vector3 Vector3Project(Vector3 v1, Vector3 v2)
 823{
 824    Vector3 result = { 0 };
 825
 826    float v1dv2 = (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z);
 827    float v2dv2 = (v2.x*v2.x + v2.y*v2.y + v2.z*v2.z);
 828
 829    float mag = v1dv2/v2dv2;
 830
 831    result.x = v2.x*mag;
 832    result.y = v2.y*mag;
 833    result.z = v2.z*mag;
 834
 835    return result;
 836}
 837
 838//Calculate the rejection of the vector v1 on to v2
 839RMAPI Vector3 Vector3Reject(Vector3 v1, Vector3 v2)
 840{
 841    Vector3 result = { 0 };
 842
 843    float v1dv2 = (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z);
 844    float v2dv2 = (v2.x*v2.x + v2.y*v2.y + v2.z*v2.z);
 845
 846    float mag = v1dv2/v2dv2;
 847
 848    result.x = v1.x - (v2.x*mag);
 849    result.y = v1.y - (v2.y*mag);
 850    result.z = v1.z - (v2.z*mag);
 851
 852    return result;
 853}
 854
 855// Orthonormalize provided vectors
 856// Makes vectors normalized and orthogonal to each other
 857// Gram-Schmidt function implementation
 858RMAPI void Vector3OrthoNormalize(Vector3 *v1, Vector3 *v2)
 859{
 860    float length = 0.0f;
 861    float ilength = 0.0f;
 862
 863    // Vector3Normalize(*v1);
 864    Vector3 v = *v1;
 865    length = sqrtf(v.x*v.x + v.y*v.y + v.z*v.z);
 866    if (length == 0.0f) length = 1.0f;
 867    ilength = 1.0f/length;
 868    v1->x *= ilength;
 869    v1->y *= ilength;
 870    v1->z *= ilength;
 871
 872    // Vector3CrossProduct(*v1, *v2)
 873    Vector3 vn1 = { v1->y*v2->z - v1->z*v2->y, v1->z*v2->x - v1->x*v2->z, v1->x*v2->y - v1->y*v2->x };
 874
 875    // Vector3Normalize(vn1);
 876    v = vn1;
 877    length = sqrtf(v.x*v.x + v.y*v.y + v.z*v.z);
 878    if (length == 0.0f) length = 1.0f;
 879    ilength = 1.0f/length;
 880    vn1.x *= ilength;
 881    vn1.y *= ilength;
 882    vn1.z *= ilength;
 883
 884    // Vector3CrossProduct(vn1, *v1)
 885    Vector3 vn2 = { vn1.y*v1->z - vn1.z*v1->y, vn1.z*v1->x - vn1.x*v1->z, vn1.x*v1->y - vn1.y*v1->x };
 886
 887    *v2 = vn2;
 888}
 889
 890// Transforms a Vector3 by a given Matrix
 891RMAPI Vector3 Vector3Transform(Vector3 v, Matrix mat)
 892{
 893    Vector3 result = { 0 };
 894
 895    float x = v.x;
 896    float y = v.y;
 897    float z = v.z;
 898
 899    result.x = mat.m0*x + mat.m4*y + mat.m8*z + mat.m12;
 900    result.y = mat.m1*x + mat.m5*y + mat.m9*z + mat.m13;
 901    result.z = mat.m2*x + mat.m6*y + mat.m10*z + mat.m14;
 902
 903    return result;
 904}
 905
 906// Transform a vector by quaternion rotation
 907RMAPI Vector3 Vector3RotateByQuaternion(Vector3 v, Quaternion q)
 908{
 909    Vector3 result = { 0 };
 910
 911    result.x = v.x*(q.x*q.x + q.w*q.w - q.y*q.y - q.z*q.z) + v.y*(2*q.x*q.y - 2*q.w*q.z) + v.z*(2*q.x*q.z + 2*q.w*q.y);
 912    result.y = v.x*(2*q.w*q.z + 2*q.x*q.y) + v.y*(q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z) + v.z*(-2*q.w*q.x + 2*q.y*q.z);
 913    result.z = v.x*(-2*q.w*q.y + 2*q.x*q.z) + v.y*(2*q.w*q.x + 2*q.y*q.z)+ v.z*(q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z);
 914
 915    return result;
 916}
 917
 918// Rotates a vector around an axis
 919RMAPI Vector3 Vector3RotateByAxisAngle(Vector3 v, Vector3 axis, float angle)
 920{
 921    // Using Euler-Rodrigues Formula
 922    // Ref.: https://en.wikipedia.org/w/index.php?title=Euler%E2%80%93Rodrigues_formula
 923
 924    Vector3 result = v;
 925
 926    // Vector3Normalize(axis);
 927    float length = sqrtf(axis.x*axis.x + axis.y*axis.y + axis.z*axis.z);
 928    if (length == 0.0f) length = 1.0f;
 929    float ilength = 1.0f/length;
 930    axis.x *= ilength;
 931    axis.y *= ilength;
 932    axis.z *= ilength;
 933
 934    angle /= 2.0f;
 935    float a = sinf(angle);
 936    float b = axis.x*a;
 937    float c = axis.y*a;
 938    float d = axis.z*a;
 939    a = cosf(angle);
 940    Vector3 w = { b, c, d };
 941
 942    // Vector3CrossProduct(w, v)
 943    Vector3 wv = { w.y*v.z - w.z*v.y, w.z*v.x - w.x*v.z, w.x*v.y - w.y*v.x };
 944
 945    // Vector3CrossProduct(w, wv)
 946    Vector3 wwv = { w.y*wv.z - w.z*wv.y, w.z*wv.x - w.x*wv.z, w.x*wv.y - w.y*wv.x };
 947
 948    // Vector3Scale(wv, 2*a)
 949    a *= 2;
 950    wv.x *= a;
 951    wv.y *= a;
 952    wv.z *= a;
 953
 954    // Vector3Scale(wwv, 2)
 955    wwv.x *= 2;
 956    wwv.y *= 2;
 957    wwv.z *= 2;
 958
 959    result.x += wv.x;
 960    result.y += wv.y;
 961    result.z += wv.z;
 962
 963    result.x += wwv.x;
 964    result.y += wwv.y;
 965    result.z += wwv.z;
 966
 967    return result;
 968}
 969
 970// Move Vector towards target
 971RMAPI Vector3 Vector3MoveTowards(Vector3 v, Vector3 target, float maxDistance)
 972{
 973    Vector3 result = { 0 };
 974
 975    float dx = target.x - v.x;
 976    float dy = target.y - v.y;
 977    float dz = target.z - v.z;
 978    float value = (dx*dx) + (dy*dy) + (dz*dz);
 979
 980    if ((value == 0) || ((maxDistance >= 0) && (value <= maxDistance*maxDistance))) return target;
 981
 982    float dist = sqrtf(value);
 983
 984    result.x = v.x + dx/dist*maxDistance;
 985    result.y = v.y + dy/dist*maxDistance;
 986    result.z = v.z + dz/dist*maxDistance;
 987
 988    return result;
 989}
 990
 991// Calculate linear interpolation between two vectors
 992RMAPI Vector3 Vector3Lerp(Vector3 v1, Vector3 v2, float amount)
 993{
 994    Vector3 result = { 0 };
 995
 996    result.x = v1.x + amount*(v2.x - v1.x);
 997    result.y = v1.y + amount*(v2.y - v1.y);
 998    result.z = v1.z + amount*(v2.z - v1.z);
 999
1000    return result;
1001}
1002
1003// Calculate cubic hermite interpolation between two vectors and their tangents
1004// as described in the GLTF 2.0 specification: https://registry.khronos.org/glTF/specs/2.0/glTF-2.0.html#interpolation-cubic
1005RMAPI Vector3 Vector3CubicHermite(Vector3 v1, Vector3 tangent1, Vector3 v2, Vector3 tangent2, float amount)
1006{
1007    Vector3 result = { 0 };
1008
1009    float amountPow2 = amount*amount;
1010    float amountPow3 = amount*amount*amount;
1011
1012    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;
1013    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;
1014    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;
1015
1016    return result;
1017}
1018
1019// Calculate reflected vector to normal
1020RMAPI Vector3 Vector3Reflect(Vector3 v, Vector3 normal)
1021{
1022    Vector3 result = { 0 };
1023
1024    // I is the original vector
1025    // N is the normal of the incident plane
1026    // R = I - (2*N*(DotProduct[I, N]))
1027
1028    float dotProduct = (v.x*normal.x + v.y*normal.y + v.z*normal.z);
1029
1030    result.x = v.x - (2.0f*normal.x)*dotProduct;
1031    result.y = v.y - (2.0f*normal.y)*dotProduct;
1032    result.z = v.z - (2.0f*normal.z)*dotProduct;
1033
1034    return result;
1035}
1036
1037// Get min value for each pair of components
1038RMAPI Vector3 Vector3Min(Vector3 v1, Vector3 v2)
1039{
1040    Vector3 result = { 0 };
1041
1042    result.x = fminf(v1.x, v2.x);
1043    result.y = fminf(v1.y, v2.y);
1044    result.z = fminf(v1.z, v2.z);
1045
1046    return result;
1047}
1048
1049// Get max value for each pair of components
1050RMAPI Vector3 Vector3Max(Vector3 v1, Vector3 v2)
1051{
1052    Vector3 result = { 0 };
1053
1054    result.x = fmaxf(v1.x, v2.x);
1055    result.y = fmaxf(v1.y, v2.y);
1056    result.z = fmaxf(v1.z, v2.z);
1057
1058    return result;
1059}
1060
1061// Compute barycenter coordinates (u, v, w) for point p with respect to triangle (a, b, c)
1062// NOTE: Assumes P is on the plane of the triangle
1063RMAPI Vector3 Vector3Barycenter(Vector3 p, Vector3 a, Vector3 b, Vector3 c)
1064{
1065    Vector3 result = { 0 };
1066
1067    Vector3 v0 = { b.x - a.x, b.y - a.y, b.z - a.z };   // Vector3Subtract(b, a)
1068    Vector3 v1 = { c.x - a.x, c.y - a.y, c.z - a.z };   // Vector3Subtract(c, a)
1069    Vector3 v2 = { p.x - a.x, p.y - a.y, p.z - a.z };   // Vector3Subtract(p, a)
1070    float d00 = (v0.x*v0.x + v0.y*v0.y + v0.z*v0.z);    // Vector3DotProduct(v0, v0)
1071    float d01 = (v0.x*v1.x + v0.y*v1.y + v0.z*v1.z);    // Vector3DotProduct(v0, v1)
1072    float d11 = (v1.x*v1.x + v1.y*v1.y + v1.z*v1.z);    // Vector3DotProduct(v1, v1)
1073    float d20 = (v2.x*v0.x + v2.y*v0.y + v2.z*v0.z);    // Vector3DotProduct(v2, v0)
1074    float d21 = (v2.x*v1.x + v2.y*v1.y + v2.z*v1.z);    // Vector3DotProduct(v2, v1)
1075
1076    float denom = d00*d11 - d01*d01;
1077
1078    result.y = (d11*d20 - d01*d21)/denom;
1079    result.z = (d00*d21 - d01*d20)/denom;
1080    result.x = 1.0f - (result.z + result.y);
1081
1082    return result;
1083}
1084
1085// Projects a Vector3 from screen space into object space
1086// NOTE: Self-contained function, no other raymath functions are called
1087RMAPI Vector3 Vector3Unproject(Vector3 source, Matrix projection, Matrix view)
1088{
1089    Vector3 result = { 0 };
1090
1091    // Calculate unprojected matrix (multiply view matrix by projection matrix) and invert it
1092    Matrix matViewProj = {      // MatrixMultiply(view, projection);
1093        view.m0*projection.m0 + view.m1*projection.m4 + view.m2*projection.m8 + view.m3*projection.m12,
1094        view.m0*projection.m1 + view.m1*projection.m5 + view.m2*projection.m9 + view.m3*projection.m13,
1095        view.m0*projection.m2 + view.m1*projection.m6 + view.m2*projection.m10 + view.m3*projection.m14,
1096        view.m0*projection.m3 + view.m1*projection.m7 + view.m2*projection.m11 + view.m3*projection.m15,
1097        view.m4*projection.m0 + view.m5*projection.m4 + view.m6*projection.m8 + view.m7*projection.m12,
1098        view.m4*projection.m1 + view.m5*projection.m5 + view.m6*projection.m9 + view.m7*projection.m13,
1099        view.m4*projection.m2 + view.m5*projection.m6 + view.m6*projection.m10 + view.m7*projection.m14,
1100        view.m4*projection.m3 + view.m5*projection.m7 + view.m6*projection.m11 + view.m7*projection.m15,
1101        view.m8*projection.m0 + view.m9*projection.m4 + view.m10*projection.m8 + view.m11*projection.m12,
1102        view.m8*projection.m1 + view.m9*projection.m5 + view.m10*projection.m9 + view.m11*projection.m13,
1103        view.m8*projection.m2 + view.m9*projection.m6 + view.m10*projection.m10 + view.m11*projection.m14,
1104        view.m8*projection.m3 + view.m9*projection.m7 + view.m10*projection.m11 + view.m11*projection.m15,
1105        view.m12*projection.m0 + view.m13*projection.m4 + view.m14*projection.m8 + view.m15*projection.m12,
1106        view.m12*projection.m1 + view.m13*projection.m5 + view.m14*projection.m9 + view.m15*projection.m13,
1107        view.m12*projection.m2 + view.m13*projection.m6 + view.m14*projection.m10 + view.m15*projection.m14,
1108        view.m12*projection.m3 + view.m13*projection.m7 + view.m14*projection.m11 + view.m15*projection.m15 };
1109
1110    // Calculate inverted matrix -> MatrixInvert(matViewProj);
1111    // Cache the matrix values (speed optimization)
1112    float a00 = matViewProj.m0, a01 = matViewProj.m1, a02 = matViewProj.m2, a03 = matViewProj.m3;
1113    float a10 = matViewProj.m4, a11 = matViewProj.m5, a12 = matViewProj.m6, a13 = matViewProj.m7;
1114    float a20 = matViewProj.m8, a21 = matViewProj.m9, a22 = matViewProj.m10, a23 = matViewProj.m11;
1115    float a30 = matViewProj.m12, a31 = matViewProj.m13, a32 = matViewProj.m14, a33 = matViewProj.m15;
1116
1117    float b00 = a00*a11 - a01*a10;
1118    float b01 = a00*a12 - a02*a10;
1119    float b02 = a00*a13 - a03*a10;
1120    float b03 = a01*a12 - a02*a11;
1121    float b04 = a01*a13 - a03*a11;
1122    float b05 = a02*a13 - a03*a12;
1123    float b06 = a20*a31 - a21*a30;
1124    float b07 = a20*a32 - a22*a30;
1125    float b08 = a20*a33 - a23*a30;
1126    float b09 = a21*a32 - a22*a31;
1127    float b10 = a21*a33 - a23*a31;
1128    float b11 = a22*a33 - a23*a32;
1129
1130    // Calculate the invert determinant (inlined to avoid double-caching)
1131    float invDet = 1.0f/(b00*b11 - b01*b10 + b02*b09 + b03*b08 - b04*b07 + b05*b06);
1132
1133    Matrix matViewProjInv = {
1134        (a11*b11 - a12*b10 + a13*b09)*invDet,
1135        (-a01*b11 + a02*b10 - a03*b09)*invDet,
1136        (a31*b05 - a32*b04 + a33*b03)*invDet,
1137        (-a21*b05 + a22*b04 - a23*b03)*invDet,
1138        (-a10*b11 + a12*b08 - a13*b07)*invDet,
1139        (a00*b11 - a02*b08 + a03*b07)*invDet,
1140        (-a30*b05 + a32*b02 - a33*b01)*invDet,
1141        (a20*b05 - a22*b02 + a23*b01)*invDet,
1142        (a10*b10 - a11*b08 + a13*b06)*invDet,
1143        (-a00*b10 + a01*b08 - a03*b06)*invDet,
1144        (a30*b04 - a31*b02 + a33*b00)*invDet,
1145        (-a20*b04 + a21*b02 - a23*b00)*invDet,
1146        (-a10*b09 + a11*b07 - a12*b06)*invDet,
1147        (a00*b09 - a01*b07 + a02*b06)*invDet,
1148        (-a30*b03 + a31*b01 - a32*b00)*invDet,
1149        (a20*b03 - a21*b01 + a22*b00)*invDet };
1150
1151    // Create quaternion from source point
1152    Quaternion quat = { source.x, source.y, source.z, 1.0f };
1153
1154    // Multiply quat point by unprojected matrix
1155    Quaternion qtransformed = {     // QuaternionTransform(quat, matViewProjInv)
1156        matViewProjInv.m0*quat.x + matViewProjInv.m4*quat.y + matViewProjInv.m8*quat.z + matViewProjInv.m12*quat.w,
1157        matViewProjInv.m1*quat.x + matViewProjInv.m5*quat.y + matViewProjInv.m9*quat.z + matViewProjInv.m13*quat.w,
1158        matViewProjInv.m2*quat.x + matViewProjInv.m6*quat.y + matViewProjInv.m10*quat.z + matViewProjInv.m14*quat.w,
1159        matViewProjInv.m3*quat.x + matViewProjInv.m7*quat.y + matViewProjInv.m11*quat.z + matViewProjInv.m15*quat.w };
1160
1161    // Normalized world points in vectors
1162    result.x = qtransformed.x/qtransformed.w;
1163    result.y = qtransformed.y/qtransformed.w;
1164    result.z = qtransformed.z/qtransformed.w;
1165
1166    return result;
1167}
1168
1169// Get Vector3 as float array
1170RMAPI float3 Vector3ToFloatV(Vector3 v)
1171{
1172    float3 buffer = { 0 };
1173
1174    buffer.v[0] = v.x;
1175    buffer.v[1] = v.y;
1176    buffer.v[2] = v.z;
1177
1178    return buffer;
1179}
1180
1181// Invert the given vector
1182RMAPI Vector3 Vector3Invert(Vector3 v)
1183{
1184    Vector3 result = { 1.0f/v.x, 1.0f/v.y, 1.0f/v.z };
1185
1186    return result;
1187}
1188
1189// Clamp the components of the vector between
1190// min and max values specified by the given vectors
1191RMAPI Vector3 Vector3Clamp(Vector3 v, Vector3 min, Vector3 max)
1192{
1193    Vector3 result = { 0 };
1194
1195    result.x = fminf(max.x, fmaxf(min.x, v.x));
1196    result.y = fminf(max.y, fmaxf(min.y, v.y));
1197    result.z = fminf(max.z, fmaxf(min.z, v.z));
1198
1199    return result;
1200}
1201
1202// Clamp the magnitude of the vector between two values
1203RMAPI Vector3 Vector3ClampValue(Vector3 v, float min, float max)
1204{
1205    Vector3 result = v;
1206
1207    float length = (v.x*v.x) + (v.y*v.y) + (v.z*v.z);
1208    if (length > 0.0f)
1209    {
1210        length = sqrtf(length);
1211
1212        float scale = 1; // By default, 1 as the neutral element
1213        if (length < min) scale = min/length;
1214        else if (length > max) scale = max/length;
1215
1216        result.x = v.x*scale;
1217        result.y = v.y*scale;
1218        result.z = v.z*scale;
1219    }
1220
1221    return result;
1222}
1223
1224// Check whether two given vectors are almost equal
1225RMAPI int Vector3Equals(Vector3 p, Vector3 q)
1226{
1227#if !defined(EPSILON)
1228    #define EPSILON 0.000001f
1229#endif
1230
1231    int result = ((fabsf(p.x - q.x)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.x), fabsf(q.x))))) &&
1232                 ((fabsf(p.y - q.y)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.y), fabsf(q.y))))) &&
1233                 ((fabsf(p.z - q.z)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.z), fabsf(q.z)))));
1234
1235    return result;
1236}
1237
1238// Compute the direction of a refracted ray
1239// v: normalized direction of the incoming ray
1240// n: normalized normal vector of the interface of two optical media
1241// r: ratio of the refractive index of the medium from where the ray comes
1242// to the refractive index of the medium on the other side of the surface
1243RMAPI Vector3 Vector3Refract(Vector3 v, Vector3 n, float r)
1244{
1245    Vector3 result = { 0 };
1246
1247    float dot = v.x*n.x + v.y*n.y + v.z*n.z;
1248    float d = 1.0f - r*r*(1.0f - dot*dot);
1249
1250    if (d >= 0.0f)
1251    {
1252        d = sqrtf(d);
1253        v.x = r*v.x - (r*dot + d)*n.x;
1254        v.y = r*v.y - (r*dot + d)*n.y;
1255        v.z = r*v.z - (r*dot + d)*n.z;
1256
1257        result = v;
1258    }
1259
1260    return result;
1261}
1262
1263
1264//----------------------------------------------------------------------------------
1265// Module Functions Definition - Vector4 math
1266//----------------------------------------------------------------------------------
1267// Get  vector zero
1268RMAPI Vector4 Vector4Zero(void)
1269{
1270    Vector4 result = { 0.0f, 0.0f, 0.0f, 0.0f };
1271    return result;
1272}
1273
1274// Get vector one
1275RMAPI Vector4 Vector4One(void)
1276{
1277    Vector4 result = { 1.0f, 1.0f, 1.0f, 1.0f };
1278    return result;
1279}
1280
1281// Add two vectors
1282RMAPI Vector4 Vector4Add(Vector4 v1, Vector4 v2)
1283{
1284    Vector4 result = {
1285        v1.x + v2.x,
1286        v1.y + v2.y,
1287        v1.z + v2.z,
1288        v1.w + v2.w
1289    };
1290    return result;
1291}
1292
1293// Add value to vector components
1294RMAPI Vector4 Vector4AddValue(Vector4 v, float add)
1295{
1296    Vector4 result = {
1297        v.x + add,
1298        v.y + add,
1299        v.z + add,
1300        v.w + add
1301    };
1302    return result;
1303}
1304
1305// Substract vectors
1306RMAPI Vector4 Vector4Subtract(Vector4 v1, Vector4 v2)
1307{
1308    Vector4 result = {
1309        v1.x - v2.x,
1310        v1.y - v2.y,
1311        v1.z - v2.z,
1312        v1.w - v2.w
1313    };
1314    return result;
1315}
1316
1317// Substract value from vector components
1318RMAPI Vector4 Vector4SubtractValue(Vector4 v, float add)
1319{
1320    Vector4 result = {
1321        v.x - add,
1322        v.y - add,
1323        v.z - add,
1324        v.w - add
1325    };
1326    return result;
1327}
1328
1329// Vector length
1330RMAPI float Vector4Length(Vector4 v)
1331{
1332    float result = sqrtf((v.x*v.x) + (v.y*v.y) + (v.z*v.z) + (v.w*v.w));
1333    return result;
1334}
1335
1336// Vector square length
1337RMAPI float Vector4LengthSqr(Vector4 v)
1338{
1339    float result = (v.x*v.x) + (v.y*v.y) + (v.z*v.z) + (v.w*v.w);
1340    return result;
1341}
1342
1343// Vectors dot product
1344RMAPI float Vector4DotProduct(Vector4 v1, Vector4 v2)
1345{
1346    float result = (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z + v1.w*v2.w);
1347    return result;
1348}
1349
1350// Calculate distance between two vectors
1351RMAPI float Vector4Distance(Vector4 v1, Vector4 v2)
1352{
1353    float result = sqrtf(
1354        (v1.x - v2.x)*(v1.x - v2.x) + (v1.y - v2.y)*(v1.y - v2.y) +
1355        (v1.z - v2.z)*(v1.z - v2.z) + (v1.w - v2.w)*(v1.w - v2.w));
1356    return result;
1357}
1358
1359// Calculate square distance between two vectors
1360RMAPI float Vector4DistanceSqr(Vector4 v1, Vector4 v2)
1361{
1362    float result =
1363        (v1.x - v2.x)*(v1.x - v2.x) + (v1.y - v2.y)*(v1.y - v2.y) +
1364        (v1.z - v2.z)*(v1.z - v2.z) + (v1.w - v2.w)*(v1.w - v2.w);
1365
1366    return result;
1367}
1368
1369// Scale vector components by value (multiply)
1370RMAPI Vector4 Vector4Scale(Vector4 v, float scale)
1371{
1372    Vector4 result = { v.x*scale, v.y*scale, v.z*scale, v.w*scale };
1373    return result;
1374}
1375
1376// Multiply vector by vector
1377RMAPI Vector4 Vector4Multiply(Vector4 v1, Vector4 v2)
1378{
1379    Vector4 result = { v1.x*v2.x, v1.y*v2.y, v1.z*v2.z, v1.w*v2.w };
1380    return result;
1381}
1382
1383// Negate vector
1384RMAPI Vector4 Vector4Negate(Vector4 v)
1385{
1386    Vector4 result = { -v.x, -v.y, -v.z, -v.w };
1387    return result;
1388}
1389
1390// Divide vector by vector
1391RMAPI Vector4 Vector4Divide(Vector4 v1, Vector4 v2)
1392{
1393    Vector4 result = { v1.x/v2.x, v1.y/v2.y, v1.z/v2.z, v1.w/v2.w };
1394    return result;
1395}
1396
1397// Normalize provided vector
1398RMAPI Vector4 Vector4Normalize(Vector4 v)
1399{
1400    Vector4 result = { 0 };
1401    float length = sqrtf((v.x*v.x) + (v.y*v.y) + (v.z*v.z) + (v.w*v.w));
1402
1403    if (length > 0)
1404    {
1405        float ilength = 1.0f/length;
1406        result.x = v.x*ilength;
1407        result.y = v.y*ilength;
1408        result.z = v.z*ilength;
1409        result.w = v.w*ilength;
1410    }
1411
1412    return result;
1413}
1414
1415// Get min value for each pair of components
1416RMAPI Vector4 Vector4Min(Vector4 v1, Vector4 v2)
1417{
1418    Vector4 result = { 0 };
1419
1420    result.x = fminf(v1.x, v2.x);
1421    result.y = fminf(v1.y, v2.y);
1422    result.z = fminf(v1.z, v2.z);
1423    result.w = fminf(v1.w, v2.w);
1424
1425    return result;
1426}
1427
1428// Get max value for each pair of components
1429RMAPI Vector4 Vector4Max(Vector4 v1, Vector4 v2)
1430{
1431    Vector4 result = { 0 };
1432
1433    result.x = fmaxf(v1.x, v2.x);
1434    result.y = fmaxf(v1.y, v2.y);
1435    result.z = fmaxf(v1.z, v2.z);
1436    result.w = fmaxf(v1.w, v2.w);
1437
1438    return result;
1439}
1440
1441// Calculate linear interpolation between two vectors
1442RMAPI Vector4 Vector4Lerp(Vector4 v1, Vector4 v2, float amount)
1443{
1444    Vector4 result = { 0 };
1445
1446    result.x = v1.x + amount*(v2.x - v1.x);
1447    result.y = v1.y + amount*(v2.y - v1.y);
1448    result.z = v1.z + amount*(v2.z - v1.z);
1449    result.w = v1.w + amount*(v2.w - v1.w);
1450
1451    return result;
1452}
1453
1454// Move Vector towards target
1455RMAPI Vector4 Vector4MoveTowards(Vector4 v, Vector4 target, float maxDistance)
1456{
1457    Vector4 result = { 0 };
1458
1459    float dx = target.x - v.x;
1460    float dy = target.y - v.y;
1461    float dz = target.z - v.z;
1462    float dw = target.w - v.w;
1463    float value = (dx*dx) + (dy*dy) + (dz*dz) + (dw*dw);
1464
1465    if ((value == 0) || ((maxDistance >= 0) && (value <= maxDistance*maxDistance))) return target;
1466
1467    float dist = sqrtf(value);
1468
1469    result.x = v.x + dx/dist*maxDistance;
1470    result.y = v.y + dy/dist*maxDistance;
1471    result.z = v.z + dz/dist*maxDistance;
1472    result.w = v.w + dw/dist*maxDistance;
1473
1474    return result;
1475}
1476
1477// Invert the given vector
1478RMAPI Vector4 Vector4Invert(Vector4 v)
1479{
1480    Vector4 result = { 1.0f/v.x, 1.0f/v.y, 1.0f/v.z, 1.0f/v.w };
1481    return result;
1482}
1483
1484// Check whether two given vectors are almost equal
1485RMAPI int Vector4Equals(Vector4 p, Vector4 q)
1486{
1487#if !defined(EPSILON)
1488    #define EPSILON 0.000001f
1489#endif
1490
1491    int result = ((fabsf(p.x - q.x)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.x), fabsf(q.x))))) &&
1492                  ((fabsf(p.y - q.y)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.y), fabsf(q.y))))) &&
1493                  ((fabsf(p.z - q.z)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.z), fabsf(q.z))))) &&
1494                  ((fabsf(p.w - q.w)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.w), fabsf(q.w)))));
1495    return result;
1496}
1497
1498
1499//----------------------------------------------------------------------------------
1500// Module Functions Definition - Matrix math
1501//----------------------------------------------------------------------------------
1502
1503// Compute matrix determinant
1504RMAPI float MatrixDeterminant(Matrix mat)
1505{
1506    float result = 0.0f;
1507/*
1508    // Cache the matrix values (speed optimization)
1509    float a00 = mat.m0, a01 = mat.m1, a02 = mat.m2, a03 = mat.m3;
1510    float a10 = mat.m4, a11 = mat.m5, a12 = mat.m6, a13 = mat.m7;
1511    float a20 = mat.m8, a21 = mat.m9, a22 = mat.m10, a23 = mat.m11;
1512    float a30 = mat.m12, a31 = mat.m13, a32 = mat.m14, a33 = mat.m15;
1513
1514    // NOTE: It takes 72 multiplication to calculate 4x4 matrix determinant
1515    result = a30*a21*a12*a03 - a20*a31*a12*a03 - a30*a11*a22*a03 + a10*a31*a22*a03 +
1516             a20*a11*a32*a03 - a10*a21*a32*a03 - a30*a21*a02*a13 + a20*a31*a02*a13 +
1517             a30*a01*a22*a13 - a00*a31*a22*a13 - a20*a01*a32*a13 + a00*a21*a32*a13 +
1518             a30*a11*a02*a23 - a10*a31*a02*a23 - a30*a01*a12*a23 + a00*a31*a12*a23 +
1519             a10*a01*a32*a23 - a00*a11*a32*a23 - a20*a11*a02*a33 + a10*a21*a02*a33 +
1520             a20*a01*a12*a33 - a00*a21*a12*a33 - a10*a01*a22*a33 + a00*a11*a22*a33;
1521*/
1522    // Using Laplace expansion (https://en.wikipedia.org/wiki/Laplace_expansion),
1523    // previous operation can be simplified to 40 multiplications, decreasing matrix
1524    // size from 4x4 to 2x2 using minors
1525
1526    // Cache the matrix values (speed optimization)
1527    float m0 = mat.m0, m1 = mat.m1, m2 = mat.m2, m3 = mat.m3;
1528    float m4 = mat.m4, m5 = mat.m5, m6 = mat.m6, m7 = mat.m7;
1529    float m8 = mat.m8, m9 = mat.m9, m10 = mat.m10, m11 = mat.m11;
1530    float m12 = mat.m12, m13 = mat.m13, m14 = mat.m14, m15 = mat.m15;
1531
1532    result = (m0*((m5*(m10*m15 - m11*m14) - m9*(m6*m15 - m7*m14) + m13*(m6*m11 - m7*m10))) -
1533        m4*((m1*(m10*m15 - m11*m14) - m9*(m2*m15 - m3*m14) + m13*(m2*m11 - m3*m10))) +
1534        m8*((m1*(m6*m15 - m7*m14) - m5*(m2*m15 - m3*m14) + m13*(m2*m7 - m3*m6))) -
1535        m12*((m1*(m6*m11 - m7*m10) - m5*(m2*m11 - m3*m10) + m9*(m2*m7 - m3*m6))));
1536
1537    return result;
1538}
1539
1540// Get the trace of the matrix (sum of the values along the diagonal)
1541RMAPI float MatrixTrace(Matrix mat)
1542{
1543    float result = (mat.m0 + mat.m5 + mat.m10 + mat.m15);
1544
1545    return result;
1546}
1547
1548// Transposes provided matrix
1549RMAPI Matrix MatrixTranspose(Matrix mat)
1550{
1551    Matrix result = { 0 };
1552
1553    result.m0 = mat.m0;
1554    result.m1 = mat.m4;
1555    result.m2 = mat.m8;
1556    result.m3 = mat.m12;
1557    result.m4 = mat.m1;
1558    result.m5 = mat.m5;
1559    result.m6 = mat.m9;
1560    result.m7 = mat.m13;
1561    result.m8 = mat.m2;
1562    result.m9 = mat.m6;
1563    result.m10 = mat.m10;
1564    result.m11 = mat.m14;
1565    result.m12 = mat.m3;
1566    result.m13 = mat.m7;
1567    result.m14 = mat.m11;
1568    result.m15 = mat.m15;
1569
1570    return result;
1571}
1572
1573// Invert provided matrix
1574RMAPI Matrix MatrixInvert(Matrix mat)
1575{
1576    Matrix result = { 0 };
1577
1578    // Cache the matrix values (speed optimization)
1579    float a00 = mat.m0, a01 = mat.m1, a02 = mat.m2, a03 = mat.m3;
1580    float a10 = mat.m4, a11 = mat.m5, a12 = mat.m6, a13 = mat.m7;
1581    float a20 = mat.m8, a21 = mat.m9, a22 = mat.m10, a23 = mat.m11;
1582    float a30 = mat.m12, a31 = mat.m13, a32 = mat.m14, a33 = mat.m15;
1583
1584    float b00 = a00*a11 - a01*a10;
1585    float b01 = a00*a12 - a02*a10;
1586    float b02 = a00*a13 - a03*a10;
1587    float b03 = a01*a12 - a02*a11;
1588    float b04 = a01*a13 - a03*a11;
1589    float b05 = a02*a13 - a03*a12;
1590    float b06 = a20*a31 - a21*a30;
1591    float b07 = a20*a32 - a22*a30;
1592    float b08 = a20*a33 - a23*a30;
1593    float b09 = a21*a32 - a22*a31;
1594    float b10 = a21*a33 - a23*a31;
1595    float b11 = a22*a33 - a23*a32;
1596
1597    // Calculate the invert determinant (inlined to avoid double-caching)
1598    float invDet = 1.0f/(b00*b11 - b01*b10 + b02*b09 + b03*b08 - b04*b07 + b05*b06);
1599
1600    result.m0 = (a11*b11 - a12*b10 + a13*b09)*invDet;
1601    result.m1 = (-a01*b11 + a02*b10 - a03*b09)*invDet;
1602    result.m2 = (a31*b05 - a32*b04 + a33*b03)*invDet;
1603    result.m3 = (-a21*b05 + a22*b04 - a23*b03)*invDet;
1604    result.m4 = (-a10*b11 + a12*b08 - a13*b07)*invDet;
1605    result.m5 = (a00*b11 - a02*b08 + a03*b07)*invDet;
1606    result.m6 = (-a30*b05 + a32*b02 - a33*b01)*invDet;
1607    result.m7 = (a20*b05 - a22*b02 + a23*b01)*invDet;
1608    result.m8 = (a10*b10 - a11*b08 + a13*b06)*invDet;
1609    result.m9 = (-a00*b10 + a01*b08 - a03*b06)*invDet;
1610    result.m10 = (a30*b04 - a31*b02 + a33*b00)*invDet;
1611    result.m11 = (-a20*b04 + a21*b02 - a23*b00)*invDet;
1612    result.m12 = (-a10*b09 + a11*b07 - a12*b06)*invDet;
1613    result.m13 = (a00*b09 - a01*b07 + a02*b06)*invDet;
1614    result.m14 = (-a30*b03 + a31*b01 - a32*b00)*invDet;
1615    result.m15 = (a20*b03 - a21*b01 + a22*b00)*invDet;
1616
1617    return result;
1618}
1619
1620// Get identity matrix
1621RMAPI Matrix MatrixIdentity(void)
1622{
1623    Matrix result = { 1.0f, 0.0f, 0.0f, 0.0f,
1624                      0.0f, 1.0f, 0.0f, 0.0f,
1625                      0.0f, 0.0f, 1.0f, 0.0f,
1626                      0.0f, 0.0f, 0.0f, 1.0f };
1627
1628    return result;
1629}
1630
1631// Add two matrices
1632RMAPI Matrix MatrixAdd(Matrix left, Matrix right)
1633{
1634    Matrix result = { 0 };
1635
1636    result.m0 = left.m0 + right.m0;
1637    result.m1 = left.m1 + right.m1;
1638    result.m2 = left.m2 + right.m2;
1639    result.m3 = left.m3 + right.m3;
1640    result.m4 = left.m4 + right.m4;
1641    result.m5 = left.m5 + right.m5;
1642    result.m6 = left.m6 + right.m6;
1643    result.m7 = left.m7 + right.m7;
1644    result.m8 = left.m8 + right.m8;
1645    result.m9 = left.m9 + right.m9;
1646    result.m10 = left.m10 + right.m10;
1647    result.m11 = left.m11 + right.m11;
1648    result.m12 = left.m12 + right.m12;
1649    result.m13 = left.m13 + right.m13;
1650    result.m14 = left.m14 + right.m14;
1651    result.m15 = left.m15 + right.m15;
1652
1653    return result;
1654}
1655
1656// Subtract two matrices (left - right)
1657RMAPI Matrix MatrixSubtract(Matrix left, Matrix right)
1658{
1659    Matrix result = { 0 };
1660
1661    result.m0 = left.m0 - right.m0;
1662    result.m1 = left.m1 - right.m1;
1663    result.m2 = left.m2 - right.m2;
1664    result.m3 = left.m3 - right.m3;
1665    result.m4 = left.m4 - right.m4;
1666    result.m5 = left.m5 - right.m5;
1667    result.m6 = left.m6 - right.m6;
1668    result.m7 = left.m7 - right.m7;
1669    result.m8 = left.m8 - right.m8;
1670    result.m9 = left.m9 - right.m9;
1671    result.m10 = left.m10 - right.m10;
1672    result.m11 = left.m11 - right.m11;
1673    result.m12 = left.m12 - right.m12;
1674    result.m13 = left.m13 - right.m13;
1675    result.m14 = left.m14 - right.m14;
1676    result.m15 = left.m15 - right.m15;
1677
1678    return result;
1679}
1680
1681// Get two matrix multiplication
1682// NOTE: When multiplying matrices... the order matters!
1683RMAPI Matrix MatrixMultiply(Matrix left, Matrix right)
1684{
1685    Matrix result = { 0 };
1686
1687#if defined(RAYMATH_SSE_ENABLED)
1688    // Load left side and right side
1689    __m128 c0 = _mm_set_ps(right.m12, right.m8,  right.m4,  right.m0);
1690    __m128 c1 = _mm_set_ps(right.m13, right.m9,  right.m5,  right.m1);
1691    __m128 c2 = _mm_set_ps(right.m14, right.m10, right.m6,  right.m2);
1692    __m128 c3 = _mm_set_ps(right.m15, right.m11, right.m7,  right.m3);
1693
1694    // Transpose so c0..c3 become *rows* of the right matrix in semantic order
1695    _MM_TRANSPOSE4_PS(c0, c1, c2, c3);
1696
1697    float tmp[4] = { 0 };
1698    __m128 row;
1699
1700    // Row 0 of result: [m0, m1, m2, m3]
1701    row  = _mm_mul_ps(_mm_set1_ps(left.m0),  c0);
1702    row  = _mm_add_ps(row, _mm_mul_ps(_mm_set1_ps(left.m1),  c1));
1703    row  = _mm_add_ps(row, _mm_mul_ps(_mm_set1_ps(left.m2),  c2));
1704    row  = _mm_add_ps(row, _mm_mul_ps(_mm_set1_ps(left.m3),  c3));
1705    _mm_storeu_ps(tmp, row);
1706    result.m0 = tmp[0];
1707    result.m1 = tmp[1];
1708    result.m2 = tmp[2];
1709    result.m3 = tmp[3];
1710
1711    // Row 1 of result: [m4, m5, m6, m7]
1712    row  = _mm_mul_ps(_mm_set1_ps(left.m4),  c0);
1713    row  = _mm_add_ps(row, _mm_mul_ps(_mm_set1_ps(left.m5),  c1));
1714    row  = _mm_add_ps(row, _mm_mul_ps(_mm_set1_ps(left.m6),  c2));
1715    row  = _mm_add_ps(row, _mm_mul_ps(_mm_set1_ps(left.m7),  c3));
1716    _mm_storeu_ps(tmp, row);
1717    result.m4 = tmp[0];
1718    result.m5 = tmp[1];
1719    result.m6 = tmp[2];
1720    result.m7 = tmp[3];
1721
1722    // Row 2 of result: [m8, m9, m10, m11]
1723    row  = _mm_mul_ps(_mm_set1_ps(left.m8),  c0);
1724    row  = _mm_add_ps(row, _mm_mul_ps(_mm_set1_ps(left.m9),  c1));
1725    row  = _mm_add_ps(row, _mm_mul_ps(_mm_set1_ps(left.m10), c2));
1726    row  = _mm_add_ps(row, _mm_mul_ps(_mm_set1_ps(left.m11), c3));
1727    _mm_storeu_ps(tmp, row);
1728    result.m8  = tmp[0];
1729    result.m9  = tmp[1];
1730    result.m10 = tmp[2];
1731    result.m11 = tmp[3];
1732
1733    // Row 3 of result: [m12, m13, m14, m15]
1734    row  = _mm_mul_ps(_mm_set1_ps(left.m12), c0);
1735    row  = _mm_add_ps(row, _mm_mul_ps(_mm_set1_ps(left.m13), c1));
1736    row  = _mm_add_ps(row, _mm_mul_ps(_mm_set1_ps(left.m14), c2));
1737    row  = _mm_add_ps(row, _mm_mul_ps(_mm_set1_ps(left.m15), c3));
1738    _mm_storeu_ps(tmp, row);
1739    result.m12 = tmp[0];
1740    result.m13 = tmp[1];
1741    result.m14 = tmp[2];
1742    result.m15 = tmp[3];
1743#else
1744    result.m0 = left.m0*right.m0 + left.m1*right.m4 + left.m2*right.m8 + left.m3*right.m12;
1745    result.m1 = left.m0*right.m1 + left.m1*right.m5 + left.m2*right.m9 + left.m3*right.m13;
1746    result.m2 = left.m0*right.m2 + left.m1*right.m6 + left.m2*right.m10 + left.m3*right.m14;
1747    result.m3 = left.m0*right.m3 + left.m1*right.m7 + left.m2*right.m11 + left.m3*right.m15;
1748    result.m4 = left.m4*right.m0 + left.m5*right.m4 + left.m6*right.m8 + left.m7*right.m12;
1749    result.m5 = left.m4*right.m1 + left.m5*right.m5 + left.m6*right.m9 + left.m7*right.m13;
1750    result.m6 = left.m4*right.m2 + left.m5*right.m6 + left.m6*right.m10 + left.m7*right.m14;
1751    result.m7 = left.m4*right.m3 + left.m5*right.m7 + left.m6*right.m11 + left.m7*right.m15;
1752    result.m8 = left.m8*right.m0 + left.m9*right.m4 + left.m10*right.m8 + left.m11*right.m12;
1753    result.m9 = left.m8*right.m1 + left.m9*right.m5 + left.m10*right.m9 + left.m11*right.m13;
1754    result.m10 = left.m8*right.m2 + left.m9*right.m6 + left.m10*right.m10 + left.m11*right.m14;
1755    result.m11 = left.m8*right.m3 + left.m9*right.m7 + left.m10*right.m11 + left.m11*right.m15;
1756    result.m12 = left.m12*right.m0 + left.m13*right.m4 + left.m14*right.m8 + left.m15*right.m12;
1757    result.m13 = left.m12*right.m1 + left.m13*right.m5 + left.m14*right.m9 + left.m15*right.m13;
1758    result.m14 = left.m12*right.m2 + left.m13*right.m6 + left.m14*right.m10 + left.m15*right.m14;
1759    result.m15 = left.m12*right.m3 + left.m13*right.m7 + left.m14*right.m11 + left.m15*right.m15;
1760#endif
1761
1762    return result;
1763}
1764
1765// Multiply matrix components by value
1766RMAPI Matrix MatrixMultiplyValue(Matrix left, float value)
1767{
1768    Matrix result = {
1769        left.m0*value, left.m4*value, left.m8*value, left.m12*value,
1770        left.m1*value, left.m5*value, left.m9*value, left.m13*value,
1771        left.m2*value, left.m6*value, left.m10*value, left.m14*value,
1772        left.m3*value, left.m7*value, left.m11*value, left.m15*value
1773    };
1774
1775    return result;
1776}
1777
1778// Get translation matrix
1779RMAPI Matrix MatrixTranslate(float x, float y, float z)
1780{
1781    Matrix result = { 1.0f, 0.0f, 0.0f, x,
1782                      0.0f, 1.0f, 0.0f, y,
1783                      0.0f, 0.0f, 1.0f, z,
1784                      0.0f, 0.0f, 0.0f, 1.0f };
1785
1786    return result;
1787}
1788
1789// Create rotation matrix from axis and angle
1790// NOTE: Angle should be provided in radians
1791RMAPI Matrix MatrixRotate(Vector3 axis, float angle)
1792{
1793    Matrix result = { 0 };
1794
1795    float x = axis.x, y = axis.y, z = axis.z;
1796
1797    float lengthSquared = x*x + y*y + z*z;
1798
1799    if ((lengthSquared != 1.0f) && (lengthSquared != 0.0f))
1800    {
1801        float ilength = 1.0f/sqrtf(lengthSquared);
1802        x *= ilength;
1803        y *= ilength;
1804        z *= ilength;
1805    }
1806
1807    float sinres = sinf(angle);
1808    float cosres = cosf(angle);
1809    float t = 1.0f - cosres;
1810
1811    result.m0 = x*x*t + cosres;
1812    result.m1 = y*x*t + z*sinres;
1813    result.m2 = z*x*t - y*sinres;
1814    result.m3 = 0.0f;
1815
1816    result.m4 = x*y*t - z*sinres;
1817    result.m5 = y*y*t + cosres;
1818    result.m6 = z*y*t + x*sinres;
1819    result.m7 = 0.0f;
1820
1821    result.m8 = x*z*t + y*sinres;
1822    result.m9 = y*z*t - x*sinres;
1823    result.m10 = z*z*t + cosres;
1824    result.m11 = 0.0f;
1825
1826    result.m12 = 0.0f;
1827    result.m13 = 0.0f;
1828    result.m14 = 0.0f;
1829    result.m15 = 1.0f;
1830
1831    return result;
1832}
1833
1834// Get x-rotation matrix
1835// NOTE: Angle must be provided in radians
1836RMAPI Matrix MatrixRotateX(float angle)
1837{
1838    Matrix result = { 1.0f, 0.0f, 0.0f, 0.0f,
1839                      0.0f, 1.0f, 0.0f, 0.0f,
1840                      0.0f, 0.0f, 1.0f, 0.0f,
1841                      0.0f, 0.0f, 0.0f, 1.0f }; // MatrixIdentity()
1842
1843    float cosres = cosf(angle);
1844    float sinres = sinf(angle);
1845
1846    result.m5 = cosres;
1847    result.m6 = sinres;
1848    result.m9 = -sinres;
1849    result.m10 = cosres;
1850
1851    return result;
1852}
1853
1854// Get y-rotation matrix
1855// NOTE: Angle must be provided in radians
1856RMAPI Matrix MatrixRotateY(float angle)
1857{
1858    Matrix result = { 1.0f, 0.0f, 0.0f, 0.0f,
1859                      0.0f, 1.0f, 0.0f, 0.0f,
1860                      0.0f, 0.0f, 1.0f, 0.0f,
1861                      0.0f, 0.0f, 0.0f, 1.0f }; // MatrixIdentity()
1862
1863    float cosres = cosf(angle);
1864    float sinres = sinf(angle);
1865
1866    result.m0 = cosres;
1867    result.m2 = -sinres;
1868    result.m8 = sinres;
1869    result.m10 = cosres;
1870
1871    return result;
1872}
1873
1874// Get z-rotation matrix
1875// NOTE: Angle must be provided in radians
1876RMAPI Matrix MatrixRotateZ(float angle)
1877{
1878    Matrix result = { 1.0f, 0.0f, 0.0f, 0.0f,
1879                      0.0f, 1.0f, 0.0f, 0.0f,
1880                      0.0f, 0.0f, 1.0f, 0.0f,
1881                      0.0f, 0.0f, 0.0f, 1.0f }; // MatrixIdentity()
1882
1883    float cosres = cosf(angle);
1884    float sinres = sinf(angle);
1885
1886    result.m0 = cosres;
1887    result.m1 = sinres;
1888    result.m4 = -sinres;
1889    result.m5 = cosres;
1890
1891    return result;
1892}
1893
1894
1895// Get xyz-rotation matrix
1896// NOTE: Angle must be provided in radians
1897RMAPI Matrix MatrixRotateXYZ(Vector3 angle)
1898{
1899    Matrix result = { 1.0f, 0.0f, 0.0f, 0.0f,
1900                      0.0f, 1.0f, 0.0f, 0.0f,
1901                      0.0f, 0.0f, 1.0f, 0.0f,
1902                      0.0f, 0.0f, 0.0f, 1.0f }; // MatrixIdentity()
1903
1904    float cosz = cosf(-angle.z);
1905    float sinz = sinf(-angle.z);
1906    float cosy = cosf(-angle.y);
1907    float siny = sinf(-angle.y);
1908    float cosx = cosf(-angle.x);
1909    float sinx = sinf(-angle.x);
1910
1911    result.m0 = cosz*cosy;
1912    result.m1 = (cosz*siny*sinx) - (sinz*cosx);
1913    result.m2 = (cosz*siny*cosx) + (sinz*sinx);
1914
1915    result.m4 = sinz*cosy;
1916    result.m5 = (sinz*siny*sinx) + (cosz*cosx);
1917    result.m6 = (sinz*siny*cosx) - (cosz*sinx);
1918
1919    result.m8 = -siny;
1920    result.m9 = cosy*sinx;
1921    result.m10= cosy*cosx;
1922
1923    return result;
1924}
1925
1926// Get zyx-rotation matrix
1927// NOTE: Angle must be provided in radians
1928RMAPI Matrix MatrixRotateZYX(Vector3 angle)
1929{
1930    Matrix result = { 0 };
1931
1932    float cz = cosf(angle.z);
1933    float sz = sinf(angle.z);
1934    float cy = cosf(angle.y);
1935    float sy = sinf(angle.y);
1936    float cx = cosf(angle.x);
1937    float sx = sinf(angle.x);
1938
1939    result.m0 = cz*cy;
1940    result.m4 = cz*sy*sx - cx*sz;
1941    result.m8 = sz*sx + cz*cx*sy;
1942    result.m12 = 0;
1943
1944    result.m1 = cy*sz;
1945    result.m5 = cz*cx + sz*sy*sx;
1946    result.m9 = cx*sz*sy - cz*sx;
1947    result.m13 = 0;
1948
1949    result.m2 = -sy;
1950    result.m6 = cy*sx;
1951    result.m10 = cy*cx;
1952    result.m14 = 0;
1953
1954    result.m3 = 0;
1955    result.m7 = 0;
1956    result.m11 = 0;
1957    result.m15 = 1;
1958
1959    return result;
1960}
1961
1962// Get scaling matrix
1963RMAPI Matrix MatrixScale(float x, float y, float z)
1964{
1965    Matrix result = { x, 0.0f, 0.0f, 0.0f,
1966                      0.0f, y, 0.0f, 0.0f,
1967                      0.0f, 0.0f, z, 0.0f,
1968                      0.0f, 0.0f, 0.0f, 1.0f };
1969
1970    return result;
1971}
1972
1973// Get perspective projection matrix
1974RMAPI Matrix MatrixFrustum(double left, double right, double bottom, double top, double nearPlane, double farPlane)
1975{
1976    Matrix result = { 0 };
1977
1978    float rl = (float)(right - left);
1979    float tb = (float)(top - bottom);
1980    float fn = (float)(farPlane - nearPlane);
1981
1982    result.m0 = ((float)nearPlane*2.0f)/rl;
1983    result.m1 = 0.0f;
1984    result.m2 = 0.0f;
1985    result.m3 = 0.0f;
1986
1987    result.m4 = 0.0f;
1988    result.m5 = ((float)nearPlane*2.0f)/tb;
1989    result.m6 = 0.0f;
1990    result.m7 = 0.0f;
1991
1992    result.m8 = ((float)right + (float)left)/rl;
1993    result.m9 = ((float)top + (float)bottom)/tb;
1994    result.m10 = -((float)farPlane + (float)nearPlane)/fn;
1995    result.m11 = -1.0f;
1996
1997    result.m12 = 0.0f;
1998    result.m13 = 0.0f;
1999    result.m14 = -((float)farPlane*(float)nearPlane*2.0f)/fn;
2000    result.m15 = 0.0f;
2001
2002    return result;
2003}
2004
2005// Get perspective projection matrix
2006// NOTE: Fovy angle must be provided in radians
2007RMAPI Matrix MatrixPerspective(double fovY, double aspect, double nearPlane, double farPlane)
2008{
2009    Matrix result = { 0 };
2010
2011    double top = nearPlane*tan(fovY*0.5);
2012    double bottom = -top;
2013    double right = top*aspect;
2014    double left = -right;
2015
2016    // MatrixFrustum(-right, right, -top, top, near, far);
2017    float rl = (float)(right - left);
2018    float tb = (float)(top - bottom);
2019    float fn = (float)(farPlane - nearPlane);
2020
2021    result.m0 = ((float)nearPlane*2.0f)/rl;
2022    result.m5 = ((float)nearPlane*2.0f)/tb;
2023    result.m8 = ((float)right + (float)left)/rl;
2024    result.m9 = ((float)top + (float)bottom)/tb;
2025    result.m10 = -((float)farPlane + (float)nearPlane)/fn;
2026    result.m11 = -1.0f;
2027    result.m14 = -((float)farPlane*(float)nearPlane*2.0f)/fn;
2028
2029    return result;
2030}
2031
2032// Get orthographic projection matrix
2033RMAPI Matrix MatrixOrtho(double left, double right, double bottom, double top, double nearPlane, double farPlane)
2034{
2035    Matrix result = { 0 };
2036
2037    float rl = (float)(right - left);
2038    float tb = (float)(top - bottom);
2039    float fn = (float)(farPlane - nearPlane);
2040
2041    result.m0 = 2.0f/rl;
2042    result.m1 = 0.0f;
2043    result.m2 = 0.0f;
2044    result.m3 = 0.0f;
2045    result.m4 = 0.0f;
2046    result.m5 = 2.0f/tb;
2047    result.m6 = 0.0f;
2048    result.m7 = 0.0f;
2049    result.m8 = 0.0f;
2050    result.m9 = 0.0f;
2051    result.m10 = -2.0f/fn;
2052    result.m11 = 0.0f;
2053    result.m12 = -((float)left + (float)right)/rl;
2054    result.m13 = -((float)top + (float)bottom)/tb;
2055    result.m14 = -((float)farPlane + (float)nearPlane)/fn;
2056    result.m15 = 1.0f;
2057
2058    return result;
2059}
2060
2061// Get camera look-at matrix (view matrix)
2062RMAPI Matrix MatrixLookAt(Vector3 eye, Vector3 target, Vector3 up)
2063{
2064    Matrix result = { 0 };
2065
2066    float length = 0.0f;
2067    float ilength = 0.0f;
2068
2069    // Vector3Subtract(eye, target)
2070    Vector3 vz = { eye.x - target.x, eye.y - target.y, eye.z - target.z };
2071
2072    // Vector3Normalize(vz)
2073    Vector3 v = vz;
2074    length = sqrtf(v.x*v.x + v.y*v.y + v.z*v.z);
2075    if (length == 0.0f) length = 1.0f;
2076    ilength = 1.0f/length;
2077    vz.x *= ilength;
2078    vz.y *= ilength;
2079    vz.z *= ilength;
2080
2081    // Vector3CrossProduct(up, vz)
2082    Vector3 vx = { up.y*vz.z - up.z*vz.y, up.z*vz.x - up.x*vz.z, up.x*vz.y - up.y*vz.x };
2083
2084    // Vector3Normalize(x)
2085    v = vx;
2086    length = sqrtf(v.x*v.x + v.y*v.y + v.z*v.z);
2087    if (length == 0.0f) length = 1.0f;
2088    ilength = 1.0f/length;
2089    vx.x *= ilength;
2090    vx.y *= ilength;
2091    vx.z *= ilength;
2092
2093    // Vector3CrossProduct(vz, vx)
2094    Vector3 vy = { vz.y*vx.z - vz.z*vx.y, vz.z*vx.x - vz.x*vx.z, vz.x*vx.y - vz.y*vx.x };
2095
2096    result.m0 = vx.x;
2097    result.m1 = vy.x;
2098    result.m2 = vz.x;
2099    result.m3 = 0.0f;
2100    result.m4 = vx.y;
2101    result.m5 = vy.y;
2102    result.m6 = vz.y;
2103    result.m7 = 0.0f;
2104    result.m8 = vx.z;
2105    result.m9 = vy.z;
2106    result.m10 = vz.z;
2107    result.m11 = 0.0f;
2108    result.m12 = -(vx.x*eye.x + vx.y*eye.y + vx.z*eye.z);   // Vector3DotProduct(vx, eye)
2109    result.m13 = -(vy.x*eye.x + vy.y*eye.y + vy.z*eye.z);   // Vector3DotProduct(vy, eye)
2110    result.m14 = -(vz.x*eye.x + vz.y*eye.y + vz.z*eye.z);   // Vector3DotProduct(vz, eye)
2111    result.m15 = 1.0f;
2112
2113    return result;
2114}
2115
2116// Get float array of matrix data
2117RMAPI float16 MatrixToFloatV(Matrix mat)
2118{
2119    float16 result = { 0 };
2120
2121    result.v[0] = mat.m0;
2122    result.v[1] = mat.m1;
2123    result.v[2] = mat.m2;
2124    result.v[3] = mat.m3;
2125    result.v[4] = mat.m4;
2126    result.v[5] = mat.m5;
2127    result.v[6] = mat.m6;
2128    result.v[7] = mat.m7;
2129    result.v[8] = mat.m8;
2130    result.v[9] = mat.m9;
2131    result.v[10] = mat.m10;
2132    result.v[11] = mat.m11;
2133    result.v[12] = mat.m12;
2134    result.v[13] = mat.m13;
2135    result.v[14] = mat.m14;
2136    result.v[15] = mat.m15;
2137
2138    return result;
2139}
2140
2141//----------------------------------------------------------------------------------
2142// Module Functions Definition - Quaternion math
2143//----------------------------------------------------------------------------------
2144
2145// Add two quaternions
2146RMAPI Quaternion QuaternionAdd(Quaternion q1, Quaternion q2)
2147{
2148    Quaternion result = {q1.x + q2.x, q1.y + q2.y, q1.z + q2.z, q1.w + q2.w};
2149
2150    return result;
2151}
2152
2153// Add quaternion and float value
2154RMAPI Quaternion QuaternionAddValue(Quaternion q, float add)
2155{
2156    Quaternion result = {q.x + add, q.y + add, q.z + add, q.w + add};
2157
2158    return result;
2159}
2160
2161// Subtract two quaternions
2162RMAPI Quaternion QuaternionSubtract(Quaternion q1, Quaternion q2)
2163{
2164    Quaternion result = {q1.x - q2.x, q1.y - q2.y, q1.z - q2.z, q1.w - q2.w};
2165
2166    return result;
2167}
2168
2169// Subtract quaternion and float value
2170RMAPI Quaternion QuaternionSubtractValue(Quaternion q, float sub)
2171{
2172    Quaternion result = {q.x - sub, q.y - sub, q.z - sub, q.w - sub};
2173
2174    return result;
2175}
2176
2177// Get identity quaternion
2178RMAPI Quaternion QuaternionIdentity(void)
2179{
2180    Quaternion result = { 0.0f, 0.0f, 0.0f, 1.0f };
2181
2182    return result;
2183}
2184
2185// Computes the length of a quaternion
2186RMAPI float QuaternionLength(Quaternion q)
2187{
2188    float result = sqrtf(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w);
2189
2190    return result;
2191}
2192
2193// Normalize provided quaternion
2194RMAPI Quaternion QuaternionNormalize(Quaternion q)
2195{
2196    Quaternion result = { 0 };
2197
2198    float length = sqrtf(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w);
2199    if (length == 0.0f) length = 1.0f;
2200    float ilength = 1.0f/length;
2201
2202    result.x = q.x*ilength;
2203    result.y = q.y*ilength;
2204    result.z = q.z*ilength;
2205    result.w = q.w*ilength;
2206
2207    return result;
2208}
2209
2210// Invert provided quaternion
2211RMAPI Quaternion QuaternionInvert(Quaternion q)
2212{
2213    Quaternion result = q;
2214
2215    float lengthSq = q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w;
2216
2217    if (lengthSq != 0.0f)
2218    {
2219        float invLength = 1.0f/lengthSq;
2220
2221        result.x *= -invLength;
2222        result.y *= -invLength;
2223        result.z *= -invLength;
2224        result.w *= invLength;
2225    }
2226
2227    return result;
2228}
2229
2230// Calculate two quaternion multiplication
2231RMAPI Quaternion QuaternionMultiply(Quaternion q1, Quaternion q2)
2232{
2233    Quaternion result = { 0 };
2234
2235    float qax = q1.x, qay = q1.y, qaz = q1.z, qaw = q1.w;
2236    float qbx = q2.x, qby = q2.y, qbz = q2.z, qbw = q2.w;
2237
2238    result.x = qax*qbw + qaw*qbx + qay*qbz - qaz*qby;
2239    result.y = qay*qbw + qaw*qby + qaz*qbx - qax*qbz;
2240    result.z = qaz*qbw + qaw*qbz + qax*qby - qay*qbx;
2241    result.w = qaw*qbw - qax*qbx - qay*qby - qaz*qbz;
2242
2243    return result;
2244}
2245
2246// Scale quaternion by float value
2247RMAPI Quaternion QuaternionScale(Quaternion q, float mul)
2248{
2249    Quaternion result = { 0 };
2250
2251    result.x = q.x*mul;
2252    result.y = q.y*mul;
2253    result.z = q.z*mul;
2254    result.w = q.w*mul;
2255
2256    return result;
2257}
2258
2259// Divide two quaternions
2260RMAPI Quaternion QuaternionDivide(Quaternion q1, Quaternion q2)
2261{
2262    Quaternion result = { q1.x/q2.x, q1.y/q2.y, q1.z/q2.z, q1.w/q2.w };
2263
2264    return result;
2265}
2266
2267// Calculate linear interpolation between two quaternions
2268RMAPI Quaternion QuaternionLerp(Quaternion q1, Quaternion q2, float amount)
2269{
2270    Quaternion result = { 0 };
2271
2272    result.x = q1.x + amount*(q2.x - q1.x);
2273    result.y = q1.y + amount*(q2.y - q1.y);
2274    result.z = q1.z + amount*(q2.z - q1.z);
2275    result.w = q1.w + amount*(q2.w - q1.w);
2276
2277    return result;
2278}
2279
2280// Calculate slerp-optimized interpolation between two quaternions
2281RMAPI Quaternion QuaternionNlerp(Quaternion q1, Quaternion q2, float amount)
2282{
2283    Quaternion result = { 0 };
2284
2285    // QuaternionLerp(q1, q2, amount)
2286    result.x = q1.x + amount*(q2.x - q1.x);
2287    result.y = q1.y + amount*(q2.y - q1.y);
2288    result.z = q1.z + amount*(q2.z - q1.z);
2289    result.w = q1.w + amount*(q2.w - q1.w);
2290
2291    // QuaternionNormalize(q);
2292    Quaternion q = result;
2293    float length = sqrtf(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w);
2294    if (length == 0.0f) length = 1.0f;
2295    float ilength = 1.0f/length;
2296
2297    result.x = q.x*ilength;
2298    result.y = q.y*ilength;
2299    result.z = q.z*ilength;
2300    result.w = q.w*ilength;
2301
2302    return result;
2303}
2304
2305// Calculates spherical linear interpolation between two quaternions
2306RMAPI Quaternion QuaternionSlerp(Quaternion q1, Quaternion q2, float amount)
2307{
2308    Quaternion result = { 0 };
2309
2310#if !defined(EPSILON)
2311    #define EPSILON 0.000001f
2312#endif
2313
2314    float cosHalfTheta = q1.x*q2.x + q1.y*q2.y + q1.z*q2.z + q1.w*q2.w;
2315
2316    if (cosHalfTheta < 0)
2317    {
2318        q2.x = -q2.x; q2.y = -q2.y; q2.z = -q2.z; q2.w = -q2.w;
2319        cosHalfTheta = -cosHalfTheta;
2320    }
2321
2322    if (fabsf(cosHalfTheta) >= 1.0f) result = q1;
2323    else if (cosHalfTheta > 0.95f) result = QuaternionNlerp(q1, q2, amount);
2324    else
2325    {
2326        float halfTheta = acosf(cosHalfTheta);
2327        float sinHalfTheta = sqrtf(1.0f - cosHalfTheta*cosHalfTheta);
2328
2329        if (fabsf(sinHalfTheta) < EPSILON)
2330        {
2331            result.x = (q1.x*0.5f + q2.x*0.5f);
2332            result.y = (q1.y*0.5f + q2.y*0.5f);
2333            result.z = (q1.z*0.5f + q2.z*0.5f);
2334            result.w = (q1.w*0.5f + q2.w*0.5f);
2335        }
2336        else
2337        {
2338            float ratioA = sinf((1 - amount)*halfTheta)/sinHalfTheta;
2339            float ratioB = sinf(amount*halfTheta)/sinHalfTheta;
2340
2341            result.x = (q1.x*ratioA + q2.x*ratioB);
2342            result.y = (q1.y*ratioA + q2.y*ratioB);
2343            result.z = (q1.z*ratioA + q2.z*ratioB);
2344            result.w = (q1.w*ratioA + q2.w*ratioB);
2345        }
2346    }
2347
2348    return result;
2349}
2350
2351// Calculate quaternion cubic spline interpolation using Cubic Hermite Spline algorithm
2352// as described in the GLTF 2.0 specification: https://registry.khronos.org/glTF/specs/2.0/glTF-2.0.html#interpolation-cubic
2353RMAPI Quaternion QuaternionCubicHermiteSpline(Quaternion q1, Quaternion outTangent1, Quaternion q2, Quaternion inTangent2, float t)
2354{
2355    float t2 = t*t;
2356    float t3 = t2*t;
2357    float h00 = 2*t3 - 3*t2 + 1;
2358    float h10 = t3 - 2*t2 + t;
2359    float h01 = -2*t3 + 3*t2;
2360    float h11 = t3 - t2;
2361
2362    Quaternion p0 = QuaternionScale(q1, h00);
2363    Quaternion m0 = QuaternionScale(outTangent1, h10);
2364    Quaternion p1 = QuaternionScale(q2, h01);
2365    Quaternion m1 = QuaternionScale(inTangent2, h11);
2366
2367    Quaternion result = { 0 };
2368
2369    result = QuaternionAdd(p0, m0);
2370    result = QuaternionAdd(result, p1);
2371    result = QuaternionAdd(result, m1);
2372    result = QuaternionNormalize(result);
2373
2374    return result;
2375}
2376
2377// Calculate quaternion based on the rotation from one vector to another
2378RMAPI Quaternion QuaternionFromVector3ToVector3(Vector3 from, Vector3 to)
2379{
2380    Quaternion result = { 0 };
2381
2382    float cos2Theta = (from.x*to.x + from.y*to.y + from.z*to.z); // Vector3DotProduct(from, to)
2383    Vector3 cross = { from.y*to.z - from.z*to.y, from.z*to.x - from.x*to.z, from.x*to.y - from.y*to.x }; // Vector3CrossProduct(from, to)
2384
2385    result.x = cross.x;
2386    result.y = cross.y;
2387    result.z = cross.z;
2388    result.w = sqrtf(cross.x*cross.x + cross.y*cross.y + cross.z*cross.z + cos2Theta*cos2Theta) + cos2Theta;
2389
2390    // QuaternionNormalize(q);
2391    // NOTE: Normalize to essentially nlerp the original and identity to 0.5
2392    Quaternion q = result;
2393    float length = sqrtf(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w);
2394    if (length == 0.0f) length = 1.0f;
2395    float ilength = 1.0f/length;
2396
2397    result.x = q.x*ilength;
2398    result.y = q.y*ilength;
2399    result.z = q.z*ilength;
2400    result.w = q.w*ilength;
2401
2402    return result;
2403}
2404
2405// Get a quaternion for a given rotation matrix
2406RMAPI Quaternion QuaternionFromMatrix(Matrix mat)
2407{
2408    Quaternion result = { 0 };
2409
2410    float fourWSquaredMinus1 = mat.m0  + mat.m5 + mat.m10;
2411    float fourXSquaredMinus1 = mat.m0  - mat.m5 - mat.m10;
2412    float fourYSquaredMinus1 = mat.m5  - mat.m0 - mat.m10;
2413    float fourZSquaredMinus1 = mat.m10 - mat.m0 - mat.m5;
2414
2415    int biggestIndex = 0;
2416    float fourBiggestSquaredMinus1 = fourWSquaredMinus1;
2417    if (fourXSquaredMinus1 > fourBiggestSquaredMinus1)
2418    {
2419        fourBiggestSquaredMinus1 = fourXSquaredMinus1;
2420        biggestIndex = 1;
2421    }
2422
2423    if (fourYSquaredMinus1 > fourBiggestSquaredMinus1)
2424    {
2425        fourBiggestSquaredMinus1 = fourYSquaredMinus1;
2426        biggestIndex = 2;
2427    }
2428
2429    if (fourZSquaredMinus1 > fourBiggestSquaredMinus1)
2430    {
2431        fourBiggestSquaredMinus1 = fourZSquaredMinus1;
2432        biggestIndex = 3;
2433    }
2434
2435    float biggestVal = sqrtf(fourBiggestSquaredMinus1 + 1.0f)*0.5f;
2436    float mult = 0.25f/biggestVal;
2437
2438    switch (biggestIndex)
2439    {
2440        case 0:
2441            result.w = biggestVal;
2442            result.x = (mat.m6 - mat.m9)*mult;
2443            result.y = (mat.m8 - mat.m2)*mult;
2444            result.z = (mat.m1 - mat.m4)*mult;
2445            break;
2446        case 1:
2447            result.x = biggestVal;
2448            result.w = (mat.m6 - mat.m9)*mult;
2449            result.y = (mat.m1 + mat.m4)*mult;
2450            result.z = (mat.m8 + mat.m2)*mult;
2451            break;
2452        case 2:
2453            result.y = biggestVal;
2454            result.w = (mat.m8 - mat.m2)*mult;
2455            result.x = (mat.m1 + mat.m4)*mult;
2456            result.z = (mat.m6 + mat.m9)*mult;
2457            break;
2458        case 3:
2459            result.z = biggestVal;
2460            result.w = (mat.m1 - mat.m4)*mult;
2461            result.x = (mat.m8 + mat.m2)*mult;
2462            result.y = (mat.m6 + mat.m9)*mult;
2463            break;
2464    }
2465
2466    return result;
2467}
2468
2469// Get a matrix for a given quaternion
2470RMAPI Matrix QuaternionToMatrix(Quaternion q)
2471{
2472    Matrix result = { 1.0f, 0.0f, 0.0f, 0.0f,
2473                      0.0f, 1.0f, 0.0f, 0.0f,
2474                      0.0f, 0.0f, 1.0f, 0.0f,
2475                      0.0f, 0.0f, 0.0f, 1.0f }; // MatrixIdentity()
2476
2477    float a2 = q.x*q.x;
2478    float b2 = q.y*q.y;
2479    float c2 = q.z*q.z;
2480    float ac = q.x*q.z;
2481    float ab = q.x*q.y;
2482    float bc = q.y*q.z;
2483    float ad = q.w*q.x;
2484    float bd = q.w*q.y;
2485    float cd = q.w*q.z;
2486
2487    result.m0 = 1 - 2*(b2 + c2);
2488    result.m1 = 2*(ab + cd);
2489    result.m2 = 2*(ac - bd);
2490
2491    result.m4 = 2*(ab - cd);
2492    result.m5 = 1 - 2*(a2 + c2);
2493    result.m6 = 2*(bc + ad);
2494
2495    result.m8 = 2*(ac + bd);
2496    result.m9 = 2*(bc - ad);
2497    result.m10 = 1 - 2*(a2 + b2);
2498
2499    return result;
2500}
2501
2502// Get rotation quaternion for an angle and axis
2503// NOTE: Angle must be provided in radians
2504RMAPI Quaternion QuaternionFromAxisAngle(Vector3 axis, float angle)
2505{
2506    Quaternion result = { 0.0f, 0.0f, 0.0f, 1.0f };
2507
2508    float length = sqrtf(axis.x*axis.x + axis.y*axis.y + axis.z*axis.z);
2509
2510    if (length != 0.0f)
2511    {
2512        angle *= 0.5f;
2513
2514        // Vector3Normalize(axis)
2515        float ilength = 1.0f/length;
2516        axis.x *= ilength;
2517        axis.y *= ilength;
2518        axis.z *= ilength;
2519
2520        float sinres = sinf(angle);
2521        float cosres = cosf(angle);
2522
2523        result.x = axis.x*sinres;
2524        result.y = axis.y*sinres;
2525        result.z = axis.z*sinres;
2526        result.w = cosres;
2527
2528        // QuaternionNormalize(q);
2529        Quaternion q = result;
2530        length = sqrtf(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w);
2531        if (length == 0.0f) length = 1.0f;
2532        ilength = 1.0f/length;
2533        result.x = q.x*ilength;
2534        result.y = q.y*ilength;
2535        result.z = q.z*ilength;
2536        result.w = q.w*ilength;
2537    }
2538
2539    return result;
2540}
2541
2542// Get the rotation angle and axis for a given quaternion
2543RMAPI void QuaternionToAxisAngle(Quaternion q, Vector3 *outAxis, float *outAngle)
2544{
2545    if (fabsf(q.w) > 1.0f)
2546    {
2547        // QuaternionNormalize(q);
2548        float length = sqrtf(q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w);
2549        if (length == 0.0f) length = 1.0f;
2550        float ilength = 1.0f/length;
2551
2552        q.x = q.x*ilength;
2553        q.y = q.y*ilength;
2554        q.z = q.z*ilength;
2555        q.w = q.w*ilength;
2556    }
2557
2558    Vector3 resAxis = { 0.0f, 0.0f, 0.0f };
2559    float resAngle = 2.0f*acosf(q.w);
2560    float den = sqrtf(1.0f - q.w*q.w);
2561
2562    if (den > EPSILON)
2563    {
2564        resAxis.x = q.x/den;
2565        resAxis.y = q.y/den;
2566        resAxis.z = q.z/den;
2567    }
2568    else
2569    {
2570        // This occurs when the angle is zero
2571        // Not a problem, set an arbitrary normalized axis
2572        resAxis.x = 1.0f;
2573    }
2574
2575    *outAxis = resAxis;
2576    *outAngle = resAngle;
2577}
2578
2579// Get the quaternion equivalent to Euler angles
2580// NOTE: Rotation order is ZYX
2581RMAPI Quaternion QuaternionFromEuler(float pitch, float yaw, float roll)
2582{
2583    Quaternion result = { 0 };
2584
2585    float x0 = cosf(pitch*0.5f);
2586    float x1 = sinf(pitch*0.5f);
2587    float y0 = cosf(yaw*0.5f);
2588    float y1 = sinf(yaw*0.5f);
2589    float z0 = cosf(roll*0.5f);
2590    float z1 = sinf(roll*0.5f);
2591
2592    result.x = x1*y0*z0 - x0*y1*z1;
2593    result.y = x0*y1*z0 + x1*y0*z1;
2594    result.z = x0*y0*z1 - x1*y1*z0;
2595    result.w = x0*y0*z0 + x1*y1*z1;
2596
2597    return result;
2598}
2599
2600// Get the Euler angles equivalent to quaternion (roll, pitch, yaw)
2601// NOTE: Angles are returned in a Vector3 struct in radians
2602RMAPI Vector3 QuaternionToEuler(Quaternion q)
2603{
2604    Vector3 result = { 0 };
2605
2606    // Roll (x-axis rotation)
2607    float x0 = 2.0f*(q.w*q.x + q.y*q.z);
2608    float x1 = 1.0f - 2.0f*(q.x*q.x + q.y*q.y);
2609    result.x = atan2f(x0, x1);
2610
2611    // Pitch (y-axis rotation)
2612    float y0 = 2.0f*(q.w*q.y - q.z*q.x);
2613    y0 = y0 > 1.0f ? 1.0f : y0;
2614    y0 = y0 < -1.0f ? -1.0f : y0;
2615    result.y = asinf(y0);
2616
2617    // Yaw (z-axis rotation)
2618    float z0 = 2.0f*(q.w*q.z + q.x*q.y);
2619    float z1 = 1.0f - 2.0f*(q.y*q.y + q.z*q.z);
2620    result.z = atan2f(z0, z1);
2621
2622    return result;
2623}
2624
2625// Transform a quaternion given a transformation matrix
2626RMAPI Quaternion QuaternionTransform(Quaternion q, Matrix mat)
2627{
2628    Quaternion result = { 0 };
2629
2630    result.x = mat.m0*q.x + mat.m4*q.y + mat.m8*q.z + mat.m12*q.w;
2631    result.y = mat.m1*q.x + mat.m5*q.y + mat.m9*q.z + mat.m13*q.w;
2632    result.z = mat.m2*q.x + mat.m6*q.y + mat.m10*q.z + mat.m14*q.w;
2633    result.w = mat.m3*q.x + mat.m7*q.y + mat.m11*q.z + mat.m15*q.w;
2634
2635    return result;
2636}
2637
2638// Check whether two given quaternions are almost equal
2639RMAPI int QuaternionEquals(Quaternion p, Quaternion q)
2640{
2641#if !defined(EPSILON)
2642    #define EPSILON 0.000001f
2643#endif
2644
2645    int result = (((fabsf(p.x - q.x)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.x), fabsf(q.x))))) &&
2646                  ((fabsf(p.y - q.y)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.y), fabsf(q.y))))) &&
2647                  ((fabsf(p.z - q.z)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.z), fabsf(q.z))))) &&
2648                  ((fabsf(p.w - q.w)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.w), fabsf(q.w)))))) ||
2649                 (((fabsf(p.x + q.x)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.x), fabsf(q.x))))) &&
2650                  ((fabsf(p.y + q.y)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.y), fabsf(q.y))))) &&
2651                  ((fabsf(p.z + q.z)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.z), fabsf(q.z))))) &&
2652                  ((fabsf(p.w + q.w)) <= (EPSILON*fmaxf(1.0f, fmaxf(fabsf(p.w), fabsf(q.w))))));
2653
2654    return result;
2655}
2656
2657// Compose a transformation matrix from rotational, translational and scaling components
2658// TODO: This function is not following raymath conventions defined in header: NOT self-contained
2659RMAPI Matrix MatrixCompose(Vector3 translation, Quaternion rotation, Vector3 scale)
2660{
2661    // Initialize vectors
2662    Vector3 right = { 1.0f, 0.0f, 0.0f };
2663    Vector3 up = { 0.0f, 1.0f, 0.0f };
2664    Vector3 forward = { 0.0f, 0.0f, 1.0f };
2665
2666    // Scale vectors
2667    right = Vector3Scale(right, scale.x);
2668    up = Vector3Scale(up, scale.y);
2669    forward = Vector3Scale(forward , scale.z);
2670
2671    // Rotate vectors
2672    right = Vector3RotateByQuaternion(right, rotation);
2673    up = Vector3RotateByQuaternion(up, rotation);
2674    forward = Vector3RotateByQuaternion(forward, rotation);
2675
2676    // Set result matrix output
2677    Matrix result = {
2678        right.x, up.x, forward.x, translation.x,
2679        right.y, up.y, forward.y, translation.y,
2680        right.z, up.z, forward.z, translation.z,
2681        0.0f, 0.0f, 0.0f, 1.0f
2682    };
2683
2684    return result;
2685}
2686
2687// Decompose a transformation matrix into its rotational, translational and scaling components and remove shear
2688// TODO: This function is not following raymath conventions defined in header: NOT self-contained
2689RMAPI void MatrixDecompose(Matrix mat, Vector3 *translation, Quaternion *rotation, Vector3 *scale)
2690{
2691    float eps = (float)1e-9;
2692
2693    // Extract Translation
2694    translation->x = mat.m12;
2695    translation->y = mat.m13;
2696    translation->z = mat.m14;
2697
2698    // Matrix Columns - Rotation will be extracted into here
2699    Vector3 matColumns[3] = {{ mat.m0, mat.m4, mat.m8 },
2700                             { mat.m1, mat.m5, mat.m9 },
2701                             { mat.m2, mat.m6, mat.m10 }};
2702
2703    // Shear Parameters XY, XZ, and YZ (extract and ignored)
2704    float shear[3] = { 0 };
2705
2706    // Normalized Scale Parameters
2707    Vector3 scl = { 0 };
2708
2709    // Max-Normalizing helps numerical stability
2710    float stabilizer = eps;
2711    for (int i = 0; i < 3; i++)
2712    {
2713        stabilizer = fmaxf(stabilizer, fabsf(matColumns[i].x));
2714        stabilizer = fmaxf(stabilizer, fabsf(matColumns[i].y));
2715        stabilizer = fmaxf(stabilizer, fabsf(matColumns[i].z));
2716    }
2717    matColumns[0] = Vector3Scale(matColumns[0], 1.0f / stabilizer);
2718    matColumns[1] = Vector3Scale(matColumns[1], 1.0f / stabilizer);
2719    matColumns[2] = Vector3Scale(matColumns[2], 1.0f / stabilizer);
2720
2721    // X Scale
2722    scl.x = Vector3Length(matColumns[0]);
2723    if (scl.x > eps) matColumns[0] = Vector3Scale(matColumns[0], 1.0f / scl.x);
2724
2725    // Compute XY shear and make col2 orthogonal
2726    shear[0] = Vector3DotProduct(matColumns[0], matColumns[1]);
2727    matColumns[1] = Vector3Subtract(matColumns[1], Vector3Scale(matColumns[0], shear[0]));
2728
2729    // Y Scale
2730    scl.y = Vector3Length(matColumns[1]);
2731    if (scl.y > eps)
2732    {
2733        matColumns[1] = Vector3Scale(matColumns[1], 1.0f / scl.y);
2734        shear[0] /= scl.y; // Correct XY shear
2735    }
2736
2737    // Compute XZ and YZ shears and make col3 orthogonal
2738    shear[1] = Vector3DotProduct(matColumns[0], matColumns[2]);
2739    matColumns[2] = Vector3Subtract(matColumns[2], Vector3Scale(matColumns[0], shear[1]));
2740    shear[2] = Vector3DotProduct(matColumns[1], matColumns[2]);
2741    matColumns[2] = Vector3Subtract(matColumns[2], Vector3Scale(matColumns[1], shear[2]));
2742
2743    // Z Scale
2744    scl.z = Vector3Length(matColumns[2]);
2745    if (scl.z > eps)
2746    {
2747        matColumns[2] = Vector3Scale(matColumns[2], 1.0f / scl.z);
2748        shear[1] /= scl.z; // Correct XZ shear
2749        shear[2] /= scl.z; // Correct YZ shear
2750    }
2751
2752    // matColumns are now orthonormal in O(3). Now ensure its in SO(3) by enforcing det = 1
2753    if (Vector3DotProduct(matColumns[0], Vector3CrossProduct(matColumns[1], matColumns[2])) < 0)
2754    {
2755        scl = Vector3Negate(scl);
2756        matColumns[0] = Vector3Negate(matColumns[0]);
2757        matColumns[1] = Vector3Negate(matColumns[1]);
2758        matColumns[2] = Vector3Negate(matColumns[2]);
2759    }
2760
2761    // Set Scale
2762    *scale = Vector3Scale(scl, stabilizer);
2763
2764    // Extract Rotation
2765    Matrix rotationMatrix = { matColumns[0].x, matColumns[0].y, matColumns[0].z, 0,
2766                             matColumns[1].x, matColumns[1].y, matColumns[1].z, 0,
2767                             matColumns[2].x, matColumns[2].y, matColumns[2].z, 0,
2768                             0, 0, 0, 1 };
2769    *rotation = QuaternionFromMatrix(rotationMatrix);
2770}
2771
2772#if defined(__cplusplus) && !defined(RAYMATH_DISABLE_CPP_OPERATORS)
2773
2774// Optional C++ math operators
2775//-------------------------------------------------------------------------------
2776
2777// Vector2 operators
2778static constexpr Vector2 Vector2Zeros = { 0, 0 };
2779static constexpr Vector2 Vector2Ones = { 1, 1 };
2780static constexpr Vector2 Vector2UnitX = { 1, 0 };
2781static constexpr Vector2 Vector2UnitY = { 0, 1 };
2782
2783inline Vector2 operator + (const Vector2& lhs, const Vector2& rhs)
2784{
2785    return Vector2Add(lhs, rhs);
2786}
2787
2788inline const Vector2& operator += (Vector2& lhs, const Vector2& rhs)
2789{
2790    lhs = Vector2Add(lhs, rhs);
2791    return lhs;
2792}
2793
2794inline Vector2 operator - (const Vector2& lhs, const Vector2& rhs)
2795{
2796    return Vector2Subtract(lhs, rhs);
2797}
2798
2799inline const Vector2& operator -= (Vector2& lhs, const Vector2& rhs)
2800{
2801    lhs = Vector2Subtract(lhs, rhs);
2802    return lhs;
2803}
2804
2805inline Vector2 operator * (const Vector2& lhs, const float& rhs)
2806{
2807    return Vector2Scale(lhs, rhs);
2808}
2809
2810inline const Vector2& operator *= (Vector2& lhs, const float& rhs)
2811{
2812    lhs = Vector2Scale(lhs, rhs);
2813    return lhs;
2814}
2815
2816inline Vector2 operator * (const Vector2& lhs, const Vector2& rhs)
2817{
2818    return Vector2Multiply(lhs, rhs);
2819}
2820
2821inline const Vector2& operator *= (Vector2& lhs, const Vector2& rhs)
2822{
2823    lhs = Vector2Multiply(lhs, rhs);
2824    return lhs;
2825}
2826
2827inline Vector2 operator * (const Vector2& lhs, const Matrix& rhs)
2828{
2829    return Vector2Transform(lhs, rhs);
2830}
2831
2832inline const Vector2& operator *= (Vector2& lhs, const Matrix& rhs)
2833{
2834    lhs = Vector2Transform(lhs, rhs);
2835    return lhs;
2836}
2837
2838inline Vector2 operator / (const Vector2& lhs, const float& rhs)
2839{
2840    return Vector2Scale(lhs, 1.0f/rhs);
2841}
2842
2843inline const Vector2& operator /= (Vector2& lhs, const float& rhs)
2844{
2845    lhs = Vector2Scale(lhs, 1.0f/rhs);
2846    return lhs;
2847}
2848
2849inline Vector2 operator / (const Vector2& lhs, const Vector2& rhs)
2850{
2851    return Vector2Divide(lhs, rhs);
2852}
2853
2854inline const Vector2& operator /= (Vector2& lhs, const Vector2& rhs)
2855{
2856    lhs = Vector2Divide(lhs, rhs);
2857    return lhs;
2858}
2859
2860inline bool operator == (const Vector2& lhs, const Vector2& rhs)
2861{
2862    return FloatEquals(lhs.x, rhs.x) && FloatEquals(lhs.y, rhs.y);
2863}
2864
2865inline bool operator != (const Vector2& lhs, const Vector2& rhs)
2866{
2867    return !FloatEquals(lhs.x, rhs.x) || !FloatEquals(lhs.y, rhs.y);
2868}
2869
2870// Vector3 operators
2871static constexpr Vector3 Vector3Zeros = { 0, 0, 0 };
2872static constexpr Vector3 Vector3Ones = { 1, 1, 1 };
2873static constexpr Vector3 Vector3UnitX = { 1, 0, 0 };
2874static constexpr Vector3 Vector3UnitY = { 0, 1, 0 };
2875static constexpr Vector3 Vector3UnitZ = { 0, 0, 1 };
2876
2877inline Vector3 operator + (const Vector3& lhs, const Vector3& rhs)
2878{
2879    return Vector3Add(lhs, rhs);
2880}
2881
2882inline const Vector3& operator += (Vector3& lhs, const Vector3& rhs)
2883{
2884    lhs = Vector3Add(lhs, rhs);
2885    return lhs;
2886}
2887
2888inline Vector3 operator - (const Vector3& lhs, const Vector3& rhs)
2889{
2890    return Vector3Subtract(lhs, rhs);
2891}
2892
2893inline const Vector3& operator -= (Vector3& lhs, const Vector3& rhs)
2894{
2895    lhs = Vector3Subtract(lhs, rhs);
2896    return lhs;
2897}
2898
2899inline Vector3 operator * (const Vector3& lhs, const float& rhs)
2900{
2901    return Vector3Scale(lhs, rhs);
2902}
2903
2904inline const Vector3& operator *= (Vector3& lhs, const float& rhs)
2905{
2906    lhs = Vector3Scale(lhs, rhs);
2907    return lhs;
2908}
2909
2910inline Vector3 operator * (const Vector3& lhs, const Vector3& rhs)
2911{
2912    return Vector3Multiply(lhs, rhs);
2913}
2914
2915inline const Vector3& operator *= (Vector3& lhs, const Vector3& rhs)
2916{
2917    lhs = Vector3Multiply(lhs, rhs);
2918    return lhs;
2919}
2920
2921inline Vector3 operator * (const Vector3& lhs, const Matrix& rhs)
2922{
2923    return Vector3Transform(lhs, rhs);
2924}
2925
2926inline const Vector3& operator *= (Vector3& lhs, const Matrix& rhs)
2927{
2928    lhs = Vector3Transform(lhs, rhs);
2929    return lhs;
2930}
2931
2932inline Vector3 operator / (const Vector3& lhs, const float& rhs)
2933{
2934    return Vector3Scale(lhs, 1.0f/rhs);
2935}
2936
2937inline const Vector3& operator /= (Vector3& lhs, const float& rhs)
2938{
2939    lhs = Vector3Scale(lhs, 1.0f/rhs);
2940    return lhs;
2941}
2942
2943inline Vector3 operator / (const Vector3& lhs, const Vector3& rhs)
2944{
2945    return Vector3Divide(lhs, rhs);
2946}
2947
2948inline const Vector3& operator /= (Vector3& lhs, const Vector3& rhs)
2949{
2950    lhs = Vector3Divide(lhs, rhs);
2951    return lhs;
2952}
2953
2954inline bool operator == (const Vector3& lhs, const Vector3& rhs)
2955{
2956    return FloatEquals(lhs.x, rhs.x) && FloatEquals(lhs.y, rhs.y) && FloatEquals(lhs.z, rhs.z);
2957}
2958
2959inline bool operator != (const Vector3& lhs, const Vector3& rhs)
2960{
2961    return !FloatEquals(lhs.x, rhs.x) || !FloatEquals(lhs.y, rhs.y) || !FloatEquals(lhs.z, rhs.z);
2962}
2963
2964// Vector4 operators
2965static constexpr Vector4 Vector4Zeros = { 0, 0, 0, 0 };
2966static constexpr Vector4 Vector4Ones = { 1, 1, 1, 1 };
2967static constexpr Vector4 Vector4UnitX = { 1, 0, 0, 0 };
2968static constexpr Vector4 Vector4UnitY = { 0, 1, 0, 0 };
2969static constexpr Vector4 Vector4UnitZ = { 0, 0, 1, 0 };
2970static constexpr Vector4 Vector4UnitW = { 0, 0, 0, 1 };
2971
2972inline Vector4 operator + (const Vector4& lhs, const Vector4& rhs)
2973{
2974    return Vector4Add(lhs, rhs);
2975}
2976
2977inline const Vector4& operator += (Vector4& lhs, const Vector4& rhs)
2978{
2979    lhs = Vector4Add(lhs, rhs);
2980    return lhs;
2981}
2982
2983inline Vector4 operator - (const Vector4& lhs, const Vector4& rhs)
2984{
2985    return Vector4Subtract(lhs, rhs);
2986}
2987
2988inline const Vector4& operator -= (Vector4& lhs, const Vector4& rhs)
2989{
2990    lhs = Vector4Subtract(lhs, rhs);
2991    return lhs;
2992}
2993
2994inline Vector4 operator * (const Vector4& lhs, const float& rhs)
2995{
2996    return Vector4Scale(lhs, rhs);
2997}
2998
2999inline const Vector4& operator *= (Vector4& lhs, const float& rhs)
3000{
3001    lhs = Vector4Scale(lhs, rhs);
3002    return lhs;
3003}
3004
3005inline Vector4 operator * (const Vector4& lhs, const Vector4& rhs)
3006{
3007    return Vector4Multiply(lhs, rhs);
3008}
3009
3010inline const Vector4& operator *= (Vector4& lhs, const Vector4& rhs)
3011{
3012    lhs = Vector4Multiply(lhs, rhs);
3013    return lhs;
3014}
3015
3016inline Vector4 operator / (const Vector4& lhs, const float& rhs)
3017{
3018    return Vector4Scale(lhs, 1.0f/rhs);
3019}
3020
3021inline const Vector4& operator /= (Vector4& lhs, const float& rhs)
3022{
3023    lhs = Vector4Scale(lhs, 1.0f/rhs);
3024    return lhs;
3025}
3026
3027inline Vector4 operator / (const Vector4& lhs, const Vector4& rhs)
3028{
3029    return Vector4Divide(lhs, rhs);
3030}
3031
3032inline const Vector4& operator /= (Vector4& lhs, const Vector4& rhs)
3033{
3034    lhs = Vector4Divide(lhs, rhs);
3035    return lhs;
3036}
3037
3038inline bool operator == (const Vector4& lhs, const Vector4& rhs)
3039{
3040    return FloatEquals(lhs.x, rhs.x) && FloatEquals(lhs.y, rhs.y) && FloatEquals(lhs.z, rhs.z) && FloatEquals(lhs.w, rhs.w);
3041}
3042
3043inline bool operator != (const Vector4& lhs, const Vector4& rhs)
3044{
3045    return !FloatEquals(lhs.x, rhs.x) || !FloatEquals(lhs.y, rhs.y) || !FloatEquals(lhs.z, rhs.z) || !FloatEquals(lhs.w, rhs.w);
3046}
3047
3048// Quaternion operators
3049static constexpr Quaternion QuaternionZeros = { 0, 0, 0, 0 };
3050static constexpr Quaternion QuaternionOnes = { 1, 1, 1, 1 };
3051static constexpr Quaternion QuaternionUnitX = { 0, 0, 0, 1 };
3052
3053inline Quaternion operator + (const Quaternion& lhs, const float& rhs)
3054{
3055    return QuaternionAddValue(lhs, rhs);
3056}
3057
3058inline const Quaternion& operator += (Quaternion& lhs, const float& rhs)
3059{
3060    lhs = QuaternionAddValue(lhs, rhs);
3061    return lhs;
3062}
3063
3064inline Quaternion operator - (const Quaternion& lhs, const float& rhs)
3065{
3066    return QuaternionSubtractValue(lhs, rhs);
3067}
3068
3069inline const Quaternion& operator -= (Quaternion& lhs, const float& rhs)
3070{
3071    lhs = QuaternionSubtractValue(lhs, rhs);
3072    return lhs;
3073}
3074
3075inline Quaternion operator * (const Quaternion& lhs, const Matrix& rhs)
3076{
3077    return QuaternionTransform(lhs, rhs);
3078}
3079
3080inline const Quaternion& operator *= (Quaternion& lhs, const Matrix& rhs)
3081{
3082    lhs = QuaternionTransform(lhs, rhs);
3083    return lhs;
3084}
3085
3086// Matrix operators
3087static constexpr Matrix MatrixUnit = { 1, 0, 0, 0,
3088                                       0, 1, 0, 0,
3089                                       0, 0, 1, 0,
3090                                       0, 0, 0, 1 };
3091
3092inline Matrix operator + (const Matrix& lhs, const Matrix& rhs)
3093{
3094    return MatrixAdd(lhs, rhs);
3095}
3096
3097inline const Matrix& operator += (Matrix& lhs, const Matrix& rhs)
3098{
3099    lhs = MatrixAdd(lhs, rhs);
3100    return lhs;
3101}
3102
3103inline Matrix operator - (const Matrix& lhs, const Matrix& rhs)
3104{
3105    return MatrixSubtract(lhs, rhs);
3106}
3107
3108inline const Matrix& operator -= (Matrix& lhs, const Matrix& rhs)
3109{
3110    lhs = MatrixSubtract(lhs, rhs);
3111    return lhs;
3112}
3113
3114inline Matrix operator * (const Matrix& lhs, const Matrix& rhs)
3115{
3116    return MatrixMultiply(lhs, rhs);
3117}
3118
3119inline const Matrix& operator *= (Matrix& lhs, const Matrix& rhs)
3120{
3121    lhs = MatrixMultiply(lhs, rhs);
3122    return lhs;
3123}
3124
3125inline Matrix operator * (const Matrix& lhs, const float value)
3126{
3127    return MatrixMultiplyValue(lhs, value);
3128}
3129
3130inline const Matrix& operator *= (Matrix& lhs, const float value)
3131{
3132    lhs = MatrixMultiplyValue(lhs, value);
3133    return lhs;
3134}
3135
3136//-------------------------------------------------------------------------------
3137#endif // C++ operators
3138
3139#endif // RAYMATH_H