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