text: // Copyright (c) 2017 The Khronos Group Inc.
  // Copyright (c) 2016 Oculus VR, LLC.
  //
  // SPDX-License-Identifier: Apache-2.0
  //
  // Licensed under the Apache License, Version 2.0 (the "License");
  // you may not use this file except in compliance with the License.
  // You may obtain a copy of the License at
  //
  //     http://www.apache.org/licenses/LICENSE-2.0
  //
  // Unless required by applicable law or agreed to in writing, software
  // distributed under the License is distributed on an "AS IS" BASIS,
  // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  // See the License for the specific language governing permissions and
  // limitations under the License.
  //
  // Author: J.M.P. van Waveren
  //
  
  #ifndef XR_LINEAR_H_
  #define XR_LINEAR_H_
  
  #if defined(OS_LINUX_XCB) || defined(OS_LINUX_XCB_GLX) || defined(OS_LINUX_WAYLAND)
  #pragma GCC diagnostic ignored "-Wunused-function"
  #pragma clang diagnostic ignored "-Wunused-function"
  #endif
  
  #include "openxr.h"
  
  /*
  ================================================================================================
  
  Description     :       Vector, matrix and quaternion math.
  Author          :       J.M.P. van Waveren
  Date            :       12/10/2016
  Language        :       C99
  Format          :       Indent 4 spaces - no tabs.
  Copyright       :       Copyright (c) 2016 Oculus VR, LLC. All Rights reserved.
  
  
  DESCRIPTION
  ===========
  
  All matrices are column-major.
  
  INTERFACE
  =========
  
  XrVector2f
  XrVector3f
  XrVector4f
  XrQuaternionf
  XrMatrix4x4f
  
  inline static void XrVector3f_Set(XrVector3f* v, const float value);
  inline static void XrVector3f_Add(XrVector3f* result, const XrVector3f* a, const XrVector3f* b);
  inline static void XrVector3f_Sub(XrVector3f* result, const XrVector3f* a, const XrVector3f* b);
  inline static void XrVector3f_Min(XrVector3f* result, const XrVector3f* a, const XrVector3f* b);
  inline static void XrVector3f_Max(XrVector3f* result, const XrVector3f* a, const XrVector3f* b);
  inline static void XrVector3f_Decay(XrVector3f* result, const XrVector3f* a, const float value);
  inline static void XrVector3f_Lerp(XrVector3f* result, const XrVector3f* a, const XrVector3f* b, const float fraction);
  inline static void XrVector3f_Scale(XrVector3f* result, const XrVector3f* a, const float scaleFactor);
  inline static void XrVector3f_Normalize(XrVector3f* v);
  inline static float XrVector3f_Length(const XrVector3f* v);
  
  inline static void XrQuaternionf_Lerp(XrQuaternionf* result, const XrQuaternionf* a, const XrQuaternionf* b, const float fraction);
  inline static void XrQuaternionf_Multiply(XrQuaternionf* result, const XrQuaternionf* a, const XrQuaternionf* b;
  
  inline static void XrMatrix4x4f_CreateIdentity(XrMatrix4x4f* result);
  inline static void XrMatrix4x4f_CreateTranslation(XrMatrix4x4f* result, const float x, const float y, const float z);
  inline static void XrMatrix4x4f_CreateRotation(XrMatrix4x4f* result, const float degreesX, const float degreesY,
                                                 const float degreesZ);
  inline static void XrMatrix4x4f_CreateScale(XrMatrix4x4f* result, const float x, const float y, const float z);
  inline static void XrMatrix4x4f_CreateTranslationRotationScale(XrMatrix4x4f* result, const XrVector3f* translation,
                                                                 const XrQuaternionf* rotation, const XrVector3f* scale);
  inline static void XrMatrix4x4f_CreateProjection(XrMatrix4x4f* result, const float tanAngleLeft, const float tanAngleRight,
                                                   const float tanAngleUp, float const tanAngleDown, const float nearZ,
                                                   const float farZ);
  inline static void XrMatrix4x4f_CreateProjectionFov(XrMatrix4x4f* result, const float fovDegreesLeft, const float fovDegreesRight,
                                                      const float fovDegreeUp, const float fovDegreesDown, const float nearZ,
                                                      const float farZ);
  inline static void XrMatrix4x4f_CreateFromQuaternion(XrMatrix4x4f* result, const XrQuaternionf* src);
  inline static void XrMatrix4x4f_CreateOffsetScaleForBounds(XrMatrix4x4f* result, const XrMatrix4x4f* matrix, const XrVector3f* mins,
                                                             const XrVector3f* maxs);
  
  inline static bool XrMatrix4x4f_IsAffine(const XrMatrix4x4f* matrix, const float epsilon);
  inline static bool XrMatrix4x4f_IsOrthogonal(const XrMatrix4x4f* matrix, const float epsilon);
  inline static bool XrMatrix4x4f_IsOrthonormal(const XrMatrix4x4f* matrix, const float epsilon);
  inline static bool XrMatrix4x4f_IsRigidBody(const XrMatrix4x4f* matrix, const float epsilon);
  
  inline static void XrMatrix4x4f_GetTranslation(XrVector3f* result, const XrMatrix4x4f* src);
  inline static void XrMatrix4x4f_GetRotation(XrQuaternionf* result, const XrMatrix4x4f* src);
  inline static void XrMatrix4x4f_GetScale(XrVector3f* result, const XrMatrix4x4f* src);
  
  inline static void XrMatrix4x4f_Multiply(XrMatrix4x4f* result, const XrMatrix4x4f* a, const XrMatrix4x4f* b);
  inline static void XrMatrix4x4f_Transpose(XrMatrix4x4f* result, const XrMatrix4x4f* src);
  inline static void XrMatrix4x4f_Invert(XrMatrix4x4f* result, const XrMatrix4x4f* src);
  inline static void XrMatrix4x4f_InvertRigidBody(XrMatrix4x4f* result, const XrMatrix4x4f* src);
  
  inline static void XrMatrix4x4f_TransformVector3f(XrVector3f* result, const XrMatrix4x4f* m, const XrVector3f* v);
  inline static void XrMatrix4x4f_TransformVector4f(XrVector4f* result, const XrMatrix4x4f* m, const XrVector4f* v);
  
  inline static void XrMatrix4x4f_TransformBounds(XrVector3f* resultMins, XrVector3f* resultMaxs, const XrMatrix4x4f* matrix,
                                                  const XrVector3f* mins, const XrVector3f* maxs);
  inline static bool XrMatrix4x4f_CullBounds(const XrMatrix4x4f* mvp, const XrVector3f* mins, const XrVector3f* maxs);
  
  ================================================================================================
  */
  
  #include <assert.h>
  #include <math.h>
  #include <stdbool.h>
  
  #define MATH_PI 3.14159265358979323846f
  
  #define DEFAULT_NEAR_Z 0.015625f  // exact floating point representation
  #define INFINITE_FAR_Z 0.0f
  
  static const XrColor4f XrColorRed = {1.0f, 0.0f, 0.0f, 1.0f};
  static const XrColor4f XrColorGreen = {0.0f, 1.0f, 0.0f, 1.0f};
  static const XrColor4f XrColorBlue = {0.0f, 0.0f, 1.0f, 1.0f};
  static const XrColor4f XrColorYellow = {1.0f, 1.0f, 0.0f, 1.0f};
  static const XrColor4f XrColorPurple = {1.0f, 0.0f, 1.0f, 1.0f};
  static const XrColor4f XrColorCyan = {0.0f, 1.0f, 1.0f, 1.0f};
  static const XrColor4f XrColorLightGrey = {0.7f, 0.7f, 0.7f, 1.0f};
  static const XrColor4f XrColorDarkGrey = {0.3f, 0.3f, 0.3f, 1.0f};
  
  enum GraphicsAPI { GRAPHICS_VULKAN, GRAPHICS_OPENGL, GRAPHICS_OPENGL_ES, GRAPHICS_D3D };
  
  // Column-major, pre-multiplied. This type does not exist in the OpenXR API and is provided for convenience.
  struct XrMatrix4x4f {
      float m[16];
  };
  
  inline static float XrRcpSqrt(const float x) {
      const float SMALLEST_NON_DENORMAL = 1.1754943508222875e-038f;  // ( 1U << 23 )
      const float rcp = (x >= SMALLEST_NON_DENORMAL) ? 1.0f / sqrtf(x) : 1.0f;
      return rcp;
  }
  
  inline static void XrVector3f_Set(XrVector3f* v, const float value) {
      v->x = value;
      v->y = value;
      v->z = value;
  }
  
  inline static void XrVector3f_Add(XrVector3f* result, const XrVector3f* a, const XrVector3f* b) {
      result->x = a->x + b->x;
      result->y = a->y + b->y;
      result->z = a->z + b->z;
  }
  
  inline static void XrVector3f_Sub(XrVector3f* result, const XrVector3f* a, const XrVector3f* b) {
      result->x = a->x - b->x;
      result->y = a->y - b->y;
      result->z = a->z - b->z;
  }
  
  inline static void XrVector3f_Min(XrVector3f* result, const XrVector3f* a, const XrVector3f* b) {
      result->x = (a->x < b->x) ? a->x : b->x;
      result->y = (a->y < b->y) ? a->y : b->y;
      result->z = (a->z < b->z) ? a->z : b->z;
  }
  
  inline static void XrVector3f_Max(XrVector3f* result, const XrVector3f* a, const XrVector3f* b) {
      result->x = (a->x > b->x) ? a->x : b->x;
      result->y = (a->y > b->y) ? a->y : b->y;
      result->z = (a->z > b->z) ? a->z : b->z;
  }
  
  inline static void XrVector3f_Decay(XrVector3f* result, const XrVector3f* a, const float value) {
      result->x = (fabsf(a->x) > value) ? ((a->x > 0.0f) ? (a->x - value) : (a->x + value)) : 0.0f;
      result->y = (fabsf(a->y) > value) ? ((a->y > 0.0f) ? (a->y - value) : (a->y + value)) : 0.0f;
      result->z = (fabsf(a->z) > value) ? ((a->z > 0.0f) ? (a->z - value) : (a->z + value)) : 0.0f;
  }
  
  inline static void XrVector3f_Lerp(XrVector3f* result, const XrVector3f* a, const XrVector3f* b, const float fraction) {
      result->x = a->x + fraction * (b->x - a->x);
      result->y = a->y + fraction * (b->y - a->y);
      result->z = a->z + fraction * (b->z - a->z);
  }
  
  inline static void XrVector3f_Scale(XrVector3f* result, const XrVector3f* a, const float scaleFactor) {
      result->x = a->x * scaleFactor;
      result->y = a->y * scaleFactor;
      result->z = a->z * scaleFactor;
  }
  
  inline static float XrVector3f_Dot(const XrVector3f* a, const XrVector3f* b) { return a->x * b->x + a->y * b->y + a->z * b->z; }
  
  // Compute cross product, which generates a normal vector.
  // Direction vector can be determined by right-hand rule: Pointing index finder in
  // direction a and middle finger in direction b, thumb will point in Cross(a, b).
  inline static void XrVector3f_Cross(XrVector3f* result, const XrVector3f* a, const XrVector3f* b) {
      result->x = a->y * b->z - a->z * b->y;
      result->y = a->z * b->x - a->x * b->z;
      result->x = a->x * b->y - a->y * b->x;
  }
  
  inline static void XrVector3f_Normalize(XrVector3f* v) {
      const float lengthRcp = XrRcpSqrt(v->x * v->x + v->y * v->y + v->z * v->z);
      v->x *= lengthRcp;
      v->y *= lengthRcp;
      v->z *= lengthRcp;
  }
  
  inline static float XrVector3f_Length(const XrVector3f* v) { return sqrtf(v->x * v->x + v->y * v->y + v->z * v->z); }
  
  inline static void XrQuaternionf_CreateFromAxisAngle(XrQuaternionf* result, const XrVector3f* axis, const float angleInRadians) {
      float s = sinf(angleInRadians / 2.0f);
      float lengthRcp = XrRcpSqrt(axis->x * axis->x + axis->y * axis->y + axis->z * axis->z);
      result->x = s * axis->x * lengthRcp;
      result->y = s * axis->y * lengthRcp;
      result->z = s * axis->z * lengthRcp;
      result->w = cosf(angleInRadians / 2.0f);
  }
  
  inline static void XrQuaternionf_Lerp(XrQuaternionf* result, const XrQuaternionf* a, const XrQuaternionf* b, const float fraction) {
      const float s = a->x * b->x + a->y * b->y + a->z * b->z + a->w * b->w;
      const float fa = 1.0f - fraction;
      const float fb = (s < 0.0f) ? -fraction : fraction;
      const float x = a->x * fa + b->x * fb;
      const float y = a->y * fa + b->y * fb;
      const float z = a->z * fa + b->z * fb;
      const float w = a->w * fa + b->w * fb;
      const float lengthRcp = XrRcpSqrt(x * x + y * y + z * z + w * w);
      result->x = x * lengthRcp;
      result->y = y * lengthRcp;
      result->z = z * lengthRcp;
      result->w = w * lengthRcp;
  }
  
  inline static void XrQuaternionf_Multiply(XrQuaternionf* result, const XrQuaternionf* a, const XrQuaternionf* b) {
      result->x = (b->w * a->x) + (b->x * a->w) + (b->y * a->z) - (b->z * a->y);
      result->y = (b->w * a->y) - (b->x * a->z) + (b->y * a->w) + (b->z * a->x);
      result->z = (b->w * a->z) + (b->x * a->y) - (b->y * a->x) + (b->z * a->w);
      result->w = (b->w * a->w) - (b->x * a->x) - (b->y * a->y) - (b->z * a->z);
  }
  
  // Use left-multiplication to accumulate transformations.
  inline static void XrMatrix4x4f_Multiply(XrMatrix4x4f* result, const XrMatrix4x4f* a, const XrMatrix4x4f* b) {
      result->m[0] = a->m[0] * b->m[0] + a->m[4] * b->m[1] + a->m[8] * b->m[2] + a->m[12] * b->m[3];
      result->m[1] = a->m[1] * b->m[0] + a->m[5] * b->m[1] + a->m[9] * b->m[2] + a->m[13] * b->m[3];
      result->m[2] = a->m[2] * b->m[0] + a->m[6] * b->m[1] + a->m[10] * b->m[2] + a->m[14] * b->m[3];
      result->m[3] = a->m[3] * b->m[0] + a->m[7] * b->m[1] + a->m[11] * b->m[2] + a->m[15] * b->m[3];
  
      result->m[4] = a->m[0] * b->m[4] + a->m[4] * b->m[5] + a->m[8] * b->m[6] + a->m[12] * b->m[7];
      result->m[5] = a->m[1] * b->m[4] + a->m[5] * b->m[5] + a->m[9] * b->m[6] + a->m[13] * b->m[7];
      result->m[6] = a->m[2] * b->m[4] + a->m[6] * b->m[5] + a->m[10] * b->m[6] + a->m[14] * b->m[7];
      result->m[7] = a->m[3] * b->m[4] + a->m[7] * b->m[5] + a->m[11] * b->m[6] + a->m[15] * b->m[7];
  
      result->m[8] = a->m[0] * b->m[8] + a->m[4] * b->m[9] + a->m[8] * b->m[10] + a->m[12] * b->m[11];
      result->m[9] = a->m[1] * b->m[8] + a->m[5] * b->m[9] + a->m[9] * b->m[10] + a->m[13] * b->m[11];
      result->m[10] = a->m[2] * b->m[8] + a->m[6] * b->m[9] + a->m[10] * b->m[10] + a->m[14] * b->m[11];
      result->m[11] = a->m[3] * b->m[8] + a->m[7] * b->m[9] + a->m[11] * b->m[10] + a->m[15] * b->m[11];
  
      result->m[12] = a->m[0] * b->m[12] + a->m[4] * b->m[13] + a->m[8] * b->m[14] + a->m[12] * b->m[15];
      result->m[13] = a->m[1] * b->m[12] + a->m[5] * b->m[13] + a->m[9] * b->m[14] + a->m[13] * b->m[15];
      result->m[14] = a->m[2] * b->m[12] + a->m[6] * b->m[13] + a->m[10] * b->m[14] + a->m[14] * b->m[15];
      result->m[15] = a->m[3] * b->m[12] + a->m[7] * b->m[13] + a->m[11] * b->m[14] + a->m[15] * b->m[15];
  }
  
  // Creates the transpose of the given matrix.
  inline static void XrMatrix4x4f_Transpose(XrMatrix4x4f* result, const XrMatrix4x4f* src) {
      result->m[0] = src->m[0];
      result->m[1] = src->m[4];
      result->m[2] = src->m[8];
      result->m[3] = src->m[12];
  
      result->m[4] = src->m[1];
      result->m[5] = src->m[5];
      result->m[6] = src->m[9];
      result->m[7] = src->m[13];
  
      result->m[8] = src->m[2];
      result->m[9] = src->m[6];
      result->m[10] = src->m[10];
      result->m[11] = src->m[14];
  
      result->m[12] = src->m[3];
      result->m[13] = src->m[7];
      result->m[14] = src->m[11];
      result->m[15] = src->m[15];
  }
  
  // Returns a 3x3 minor of a 4x4 matrix.
  inline static float XrMatrix4x4f_Minor(const XrMatrix4x4f* matrix, int r0, int r1, int r2, int c0, int c1, int c2) {
      return matrix->m[4 * r0 + c0] *
                 (matrix->m[4 * r1 + c1] * matrix->m[4 * r2 + c2] - matrix->m[4 * r2 + c1] * matrix->m[4 * r1 + c2]) -
             matrix->m[4 * r0 + c1] *
                 (matrix->m[4 * r1 + c0] * matrix->m[4 * r2 + c2] - matrix->m[4 * r2 + c0] * matrix->m[4 * r1 + c2]) +
             matrix->m[4 * r0 + c2] *
                 (matrix->m[4 * r1 + c0] * matrix->m[4 * r2 + c1] - matrix->m[4 * r2 + c0] * matrix->m[4 * r1 + c1]);
  }
  
  // Calculates the inverse of a 4x4 matrix.
  inline static void XrMatrix4x4f_Invert(XrMatrix4x4f* result, const XrMatrix4x4f* src) {
      const float rcpDet =
          1.0f / (src->m[0] * XrMatrix4x4f_Minor(src, 1, 2, 3, 1, 2, 3) - src->m[1] * XrMatrix4x4f_Minor(src, 1, 2, 3, 0, 2, 3) +
                  src->m[2] * XrMatrix4x4f_Minor(src, 1, 2, 3, 0, 1, 3) - src->m[3] * XrMatrix4x4f_Minor(src, 1, 2, 3, 0, 1, 2));
  
      result->m[0] = XrMatrix4x4f_Minor(src, 1, 2, 3, 1, 2, 3) * rcpDet;
      result->m[1] = -XrMatrix4x4f_Minor(src, 0, 2, 3, 1, 2, 3) * rcpDet;
      result->m[2] = XrMatrix4x4f_Minor(src, 0, 1, 3, 1, 2, 3) * rcpDet;
      result->m[3] = -XrMatrix4x4f_Minor(src, 0, 1, 2, 1, 2, 3) * rcpDet;
      result->m[4] = -XrMatrix4x4f_Minor(src, 1, 2, 3, 0, 2, 3) * rcpDet;
      result->m[5] = XrMatrix4x4f_Minor(src, 0, 2, 3, 0, 2, 3) * rcpDet;
      result->m[6] = -XrMatrix4x4f_Minor(src, 0, 1, 3, 0, 2, 3) * rcpDet;
      result->m[7] = XrMatrix4x4f_Minor(src, 0, 1, 2, 0, 2, 3) * rcpDet;
      result->m[8] = XrMatrix4x4f_Minor(src, 1, 2, 3, 0, 1, 3) * rcpDet;
      result->m[9] = -XrMatrix4x4f_Minor(src, 0, 2, 3, 0, 1, 3) * rcpDet;
      result->m[10] = XrMatrix4x4f_Minor(src, 0, 1, 3, 0, 1, 3) * rcpDet;
      result->m[11] = -XrMatrix4x4f_Minor(src, 0, 1, 2, 0, 1, 3) * rcpDet;
      result->m[12] = -XrMatrix4x4f_Minor(src, 1, 2, 3, 0, 1, 2) * rcpDet;
      result->m[13] = XrMatrix4x4f_Minor(src, 0, 2, 3, 0, 1, 2) * rcpDet;
      result->m[14] = -XrMatrix4x4f_Minor(src, 0, 1, 3, 0, 1, 2) * rcpDet;
      result->m[15] = XrMatrix4x4f_Minor(src, 0, 1, 2, 0, 1, 2) * rcpDet;
  }
  
  // Calculates the inverse of a rigid body transform.
  inline static void XrMatrix4x4f_InvertRigidBody(XrMatrix4x4f* result, const XrMatrix4x4f* src) {
      result->m[0] = src->m[0];
      result->m[1] = src->m[4];
      result->m[2] = src->m[8];
      result->m[3] = 0.0f;
      result->m[4] = src->m[1];
      result->m[5] = src->m[5];
      result->m[6] = src->m[9];
      result->m[7] = 0.0f;
      result->m[8] = src->m[2];
      result->m[9] = src->m[6];
      result->m[10] = src->m[10];
      result->m[11] = 0.0f;
      result->m[12] = -(src->m[0] * src->m[12] + src->m[1] * src->m[13] + src->m[2] * src->m[14]);
      result->m[13] = -(src->m[4] * src->m[12] + src->m[5] * src->m[13] + src->m[6] * src->m[14]);
      result->m[14] = -(src->m[8] * src->m[12] + src->m[9] * src->m[13] + src->m[10] * src->m[14]);
      result->m[15] = 1.0f;
  }
  
  // Creates an identity matrix.
  inline static void XrMatrix4x4f_CreateIdentity(XrMatrix4x4f* result) {
      result->m[0] = 1.0f;
      result->m[1] = 0.0f;
      result->m[2] = 0.0f;
      result->m[3] = 0.0f;
      result->m[4] = 0.0f;
      result->m[5] = 1.0f;
      result->m[6] = 0.0f;
      result->m[7] = 0.0f;
      result->m[8] = 0.0f;
      result->m[9] = 0.0f;
      result->m[10] = 1.0f;
      result->m[11] = 0.0f;
      result->m[12] = 0.0f;
      result->m[13] = 0.0f;
      result->m[14] = 0.0f;
      result->m[15] = 1.0f;
  }
  
  // Creates a translation matrix.
  inline static void XrMatrix4x4f_CreateTranslation(XrMatrix4x4f* result, const float x, const float y, const float z) {
      result->m[0] = 1.0f;
      result->m[1] = 0.0f;
      result->m[2] = 0.0f;
      result->m[3] = 0.0f;
      result->m[4] = 0.0f;
      result->m[5] = 1.0f;
      result->m[6] = 0.0f;
      result->m[7] = 0.0f;
      result->m[8] = 0.0f;
      result->m[9] = 0.0f;
      result->m[10] = 1.0f;
      result->m[11] = 0.0f;
      result->m[12] = x;
      result->m[13] = y;
      result->m[14] = z;
      result->m[15] = 1.0f;
  }
  
  // Creates a rotation matrix.
  // If -Z=forward, +Y=up, +X=right, then degreesX=pitch, degreesY=yaw, degreesZ=roll.
  inline static void XrMatrix4x4f_CreateRotation(XrMatrix4x4f* result, const float degreesX, const float degreesY,
                                                 const float degreesZ) {
      const float sinX = sinf(degreesX * (MATH_PI / 180.0f));
      const float cosX = cosf(degreesX * (MATH_PI / 180.0f));
      const XrMatrix4x4f rotationX = {{1, 0, 0, 0, 0, cosX, sinX, 0, 0, -sinX, cosX, 0, 0, 0, 0, 1}};
      const float sinY = sinf(degreesY * (MATH_PI / 180.0f));
      const float cosY = cosf(degreesY * (MATH_PI / 180.0f));
      const XrMatrix4x4f rotationY = {{cosY, 0, -sinY, 0, 0, 1, 0, 0, sinY, 0, cosY, 0, 0, 0, 0, 1}};
      const float sinZ = sinf(degreesZ * (MATH_PI / 180.0f));
      const float cosZ = cosf(degreesZ * (MATH_PI / 180.0f));
      const XrMatrix4x4f rotationZ = {{cosZ, sinZ, 0, 0, -sinZ, cosZ, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1}};
      XrMatrix4x4f rotationXY;
      XrMatrix4x4f_Multiply(&rotationXY, &rotationY, &rotationX);
      XrMatrix4x4f_Multiply(result, &rotationZ, &rotationXY);
  }
  
  // Creates a scale matrix.
  inline static void XrMatrix4x4f_CreateScale(XrMatrix4x4f* result, const float x, const float y, const float z) {
      result->m[0] = x;
      result->m[1] = 0.0f;
      result->m[2] = 0.0f;
      result->m[3] = 0.0f;
      result->m[4] = 0.0f;
      result->m[5] = y;
      result->m[6] = 0.0f;
      result->m[7] = 0.0f;
      result->m[8] = 0.0f;
      result->m[9] = 0.0f;
      result->m[10] = z;
      result->m[11] = 0.0f;
      result->m[12] = 0.0f;
      result->m[13] = 0.0f;
      result->m[14] = 0.0f;
      result->m[15] = 1.0f;
  }
  
  // Creates a matrix from a quaternion.
  inline static void XrMatrix4x4f_CreateFromQuaternion(XrMatrix4x4f* result, const XrQuaternionf* quat) {
      const float x2 = quat->x + quat->x;
      const float y2 = quat->y + quat->y;
      const float z2 = quat->z + quat->z;
  
      const float xx2 = quat->x * x2;
      const float yy2 = quat->y * y2;
      const float zz2 = quat->z * z2;
  
      const float yz2 = quat->y * z2;
      const float wx2 = quat->w * x2;
      const float xy2 = quat->x * y2;
      const float wz2 = quat->w * z2;
      const float xz2 = quat->x * z2;
      const float wy2 = quat->w * y2;
  
      result->m[0] = 1.0f - yy2 - zz2;
      result->m[1] = xy2 + wz2;
      result->m[2] = xz2 - wy2;
      result->m[3] = 0.0f;
  
      result->m[4] = xy2 - wz2;
      result->m[5] = 1.0f - xx2 - zz2;
      result->m[6] = yz2 + wx2;
      result->m[7] = 0.0f;
  
      result->m[8] = xz2 + wy2;
      result->m[9] = yz2 - wx2;
      result->m[10] = 1.0f - xx2 - yy2;
      result->m[11] = 0.0f;
  
      result->m[12] = 0.0f;
      result->m[13] = 0.0f;
      result->m[14] = 0.0f;
      result->m[15] = 1.0f;
  }
  
  // Creates a combined translation(rotation(scale(object))) matrix.
  inline static void XrMatrix4x4f_CreateTranslationRotationScale(XrMatrix4x4f* result, const XrVector3f* translation,
                                                                 const XrQuaternionf* rotation, const XrVector3f* scale) {
      XrMatrix4x4f scaleMatrix;
      XrMatrix4x4f_CreateScale(&scaleMatrix, scale->x, scale->y, scale->z);
  
      XrMatrix4x4f rotationMatrix;
      XrMatrix4x4f_CreateFromQuaternion(&rotationMatrix, rotation);
  
      XrMatrix4x4f translationMatrix;
      XrMatrix4x4f_CreateTranslation(&translationMatrix, translation->x, translation->y, translation->z);
  
      XrMatrix4x4f combinedMatrix;
      XrMatrix4x4f_Multiply(&combinedMatrix, &rotationMatrix, &scaleMatrix);
      XrMatrix4x4f_Multiply(result, &translationMatrix, &combinedMatrix);
  }
  
  // Creates a projection matrix based on the specified dimensions.
  // The projection matrix transforms -Z=forward, +Y=up, +X=right to the appropriate clip space for the graphics API.
  // The far plane is placed at infinity if farZ <= nearZ.
  // An infinite projection matrix is preferred for rasterization because, except for
  // things *right* up against the near plane, it always provides better precision:
  //              "Tightening the Precision of Perspective Rendering"
  //              Paul Upchurch, Mathieu Desbrun
  //              Journal of Graphics Tools, Volume 16, Issue 1, 2012
  inline static void XrMatrix4x4f_CreateProjection(XrMatrix4x4f* result, GraphicsAPI graphicsApi, const float tanAngleLeft,
                                                   const float tanAngleRight, const float tanAngleUp, float const tanAngleDown,
                                                   const float nearZ, const float farZ) {
      const float tanAngleWidth = tanAngleRight - tanAngleLeft;
  
      // Set to tanAngleDown - tanAngleUp for a clip space with positive Y down (Vulkan).
      // Set to tanAngleUp - tanAngleDown for a clip space with positive Y up (OpenGL / D3D / Metal).
      const float tanAngleHeight = graphicsApi == GRAPHICS_VULKAN ? (tanAngleDown - tanAngleUp) : (tanAngleUp - tanAngleDown);
  
      // Set to nearZ for a [-1,1] Z clip space (OpenGL / OpenGL ES).
      // Set to zero for a [0,1] Z clip space (Vulkan / D3D / Metal).
      const float offsetZ = (graphicsApi == GRAPHICS_OPENGL || graphicsApi == GRAPHICS_OPENGL_ES) ? nearZ : 0;
  
      if (farZ <= nearZ) {
          // place the far plane at infinity
          result->m[0] = 2.0f / tanAngleWidth;
          result->m[4] = 0.0f;
          result->m[8] = (tanAngleRight + tanAngleLeft) / tanAngleWidth;
          result->m[12] = 0.0f;
  
          result->m[1] = 0.0f;
          result->m[5] = 2.0f / tanAngleHeight;
          result->m[9] = (tanAngleUp + tanAngleDown) / tanAngleHeight;
          result->m[13] = 0.0f;
  
          result->m[2] = 0.0f;
          result->m[6] = 0.0f;
          result->m[10] = -1.0f;
          result->m[14] = -(nearZ + offsetZ);
  
          result->m[3] = 0.0f;
          result->m[7] = 0.0f;
          result->m[11] = -1.0f;
          result->m[15] = 0.0f;
      } else {
          // normal projection
          result->m[0] = 2.0f / tanAngleWidth;
          result->m[4] = 0.0f;
          result->m[8] = (tanAngleRight + tanAngleLeft) / tanAngleWidth;
          result->m[12] = 0.0f;
  
          result->m[1] = 0.0f;
          result->m[5] = 2.0f / tanAngleHeight;
          result->m[9] = (tanAngleUp + tanAngleDown) / tanAngleHeight;
          result->m[13] = 0.0f;
  
          result->m[2] = 0.0f;
          result->m[6] = 0.0f;
          result->m[10] = -(farZ + offsetZ) / (farZ - nearZ);
          result->m[14] = -(farZ * (nearZ + offsetZ)) / (farZ - nearZ);
  
          result->m[3] = 0.0f;
          result->m[7] = 0.0f;
          result->m[11] = -1.0f;
          result->m[15] = 0.0f;
      }
  }
  
  // Creates a projection matrix based on the specified FOV.
  inline static void XrMatrix4x4f_CreateProjectionFov(XrMatrix4x4f* result, GraphicsAPI graphicsApi, const XrFovf fov,
                                                      const float nearZ, const float farZ) {
      const float tanLeft = tanf(fov.angleLeft);
      const float tanRight = tanf(fov.angleRight);
  
      const float tanDown = tanf(fov.angleDown);
      const float tanUp = tanf(fov.angleUp);
  
      XrMatrix4x4f_CreateProjection(result, graphicsApi, tanLeft, tanRight, tanUp, tanDown, nearZ, farZ);
  }
  
  // Creates a matrix that transforms the -1 to 1 cube to cover the given 'mins' and 'maxs' transformed with the given 'matrix'.
  inline static void XrMatrix4x4f_CreateOffsetScaleForBounds(XrMatrix4x4f* result, const XrMatrix4x4f* matrix, const XrVector3f* mins,
                                                             const XrVector3f* maxs) {
      const XrVector3f offset = {(maxs->x + mins->x) * 0.5f, (maxs->y + mins->y) * 0.5f, (maxs->z + mins->z) * 0.5f};
      const XrVector3f scale = {(maxs->x - mins->x) * 0.5f, (maxs->y - mins->y) * 0.5f, (maxs->z - mins->z) * 0.5f};
  
      result->m[0] = matrix->m[0] * scale.x;
      result->m[1] = matrix->m[1] * scale.x;
      result->m[2] = matrix->m[2] * scale.x;
      result->m[3] = matrix->m[3] * scale.x;
  
      result->m[4] = matrix->m[4] * scale.y;
      result->m[5] = matrix->m[5] * scale.y;
      result->m[6] = matrix->m[6] * scale.y;
      result->m[7] = matrix->m[7] * scale.y;
  
      result->m[8] = matrix->m[8] * scale.z;
      result->m[9] = matrix->m[9] * scale.z;
      result->m[10] = matrix->m[10] * scale.z;
      result->m[11] = matrix->m[11] * scale.z;
  
      result->m[12] = matrix->m[12] + matrix->m[0] * offset.x + matrix->m[4] * offset.y + matrix->m[8] * offset.z;
      result->m[13] = matrix->m[13] + matrix->m[1] * offset.x + matrix->m[5] * offset.y + matrix->m[9] * offset.z;
      result->m[14] = matrix->m[14] + matrix->m[2] * offset.x + matrix->m[6] * offset.y + matrix->m[10] * offset.z;
      result->m[15] = matrix->m[15] + matrix->m[3] * offset.x + matrix->m[7] * offset.y + matrix->m[11] * offset.z;
  }
  
  // Returns true if the given matrix is affine.
  inline static bool XrMatrix4x4f_IsAffine(const XrMatrix4x4f* matrix, const float epsilon) {
      return fabsf(matrix->m[3]) <= epsilon && fabsf(matrix->m[7]) <= epsilon && fabsf(matrix->m[11]) <= epsilon &&
             fabsf(matrix->m[15] - 1.0f) <= epsilon;
  }
  
  // Returns true if the given matrix is orthogonal.
  inline static bool XrMatrix4x4f_IsOrthogonal(const XrMatrix4x4f* matrix, const float epsilon) {
      for (int i = 0; i < 3; i++) {
          for (int j = 0; j < 3; j++) {
              if (i != j) {
                  if (fabsf(matrix->m[4 * i + 0] * matrix->m[4 * j + 0] + matrix->m[4 * i + 1] * matrix->m[4 * j + 1] +
                            matrix->m[4 * i + 2] * matrix->m[4 * j + 2]) > epsilon) {
                      return false;
                  }
                  if (fabsf(matrix->m[4 * 0 + i] * matrix->m[4 * 0 + j] + matrix->m[4 * 1 + i] * matrix->m[4 * 1 + j] +
                            matrix->m[4 * 2 + i] * matrix->m[4 * 2 + j]) > epsilon) {
                      return false;
                  }
              }
          }
      }
      return true;
  }
  
  // Returns true if the given matrix is orthonormal.
  inline static bool XrMatrix4x4f_IsOrthonormal(const XrMatrix4x4f* matrix, const float epsilon) {
      for (int i = 0; i < 3; i++) {
          for (int j = 0; j < 3; j++) {
              const float kd = (i == j) ? 1.0f : 0.0f;  // Kronecker delta
              if (fabsf(kd - (matrix->m[4 * i + 0] * matrix->m[4 * j + 0] + matrix->m[4 * i + 1] * matrix->m[4 * j + 1] +
                              matrix->m[4 * i + 2] * matrix->m[4 * j + 2])) > epsilon) {
                  return false;
              }
              if (fabsf(kd - (matrix->m[4 * 0 + i] * matrix->m[4 * 0 + j] + matrix->m[4 * 1 + i] * matrix->m[4 * 1 + j] +
                              matrix->m[4 * 2 + i] * matrix->m[4 * 2 + j])) > epsilon) {
                  return false;
              }
          }
      }
      return true;
  }
  
  // Returns true if the given matrix is a rigid body transform.
  inline static bool XrMatrix4x4f_IsRigidBody(const XrMatrix4x4f* matrix, const float epsilon) {
      return XrMatrix4x4f_IsAffine(matrix, epsilon) && XrMatrix4x4f_IsOrthonormal(matrix, epsilon);
  }
  
  // Get the translation from a combined translation(rotation(scale(object))) matrix.
  inline static void XrMatrix4x4f_GetTranslation(XrVector3f* result, const XrMatrix4x4f* src) {
      assert(XrMatrix4x4f_IsAffine(src, 1e-4f));
      assert(XrMatrix4x4f_IsOrthogonal(src, 1e-4f));
  
      result->x = src->m[12];
      result->y = src->m[13];
      result->z = src->m[14];
  }
  
  // Get the rotation from a combined translation(rotation(scale(object))) matrix.
  inline static void XrMatrix4x4f_GetRotation(XrQuaternionf* result, const XrMatrix4x4f* src) {
      assert(XrMatrix4x4f_IsAffine(src, 1e-4f));
      assert(XrMatrix4x4f_IsOrthogonal(src, 1e-4f));
  
      const float rcpScaleX = XrRcpSqrt(src->m[0] * src->m[0] + src->m[1] * src->m[1] + src->m[2] * src->m[2]);
      const float rcpScaleY = XrRcpSqrt(src->m[4] * src->m[4] + src->m[5] * src->m[5] + src->m[6] * src->m[6]);
      const float rcpScaleZ = XrRcpSqrt(src->m[8] * src->m[8] + src->m[9] * src->m[9] + src->m[10] * src->m[10]);
      const float m[9] = {src->m[0] * rcpScaleX, src->m[1] * rcpScaleX, src->m[2] * rcpScaleX,
                          src->m[4] * rcpScaleY, src->m[5] * rcpScaleY, src->m[6] * rcpScaleY,
                          src->m[8] * rcpScaleZ, src->m[9] * rcpScaleZ, src->m[10] * rcpScaleZ};
      if (m[0 * 3 + 0] + m[1 * 3 + 1] + m[2 * 3 + 2] > 0.0f) {
          float t = +m[0 * 3 + 0] + m[1 * 3 + 1] + m[2 * 3 + 2] + 1.0f;
          float s = XrRcpSqrt(t) * 0.5f;
          result->w = s * t;
          result->z = (m[0 * 3 + 1] - m[1 * 3 + 0]) * s;
          result->y = (m[2 * 3 + 0] - m[0 * 3 + 2]) * s;
          result->x = (m[1 * 3 + 2] - m[2 * 3 + 1]) * s;
      } else if (m[0 * 3 + 0] > m[1 * 3 + 1] && m[0 * 3 + 0] > m[2 * 3 + 2]) {
          float t = +m[0 * 3 + 0] - m[1 * 3 + 1] - m[2 * 3 + 2] + 1.0f;
          float s = XrRcpSqrt(t) * 0.5f;
          result->x = s * t;
          result->y = (m[0 * 3 + 1] + m[1 * 3 + 0]) * s;
          result->z = (m[2 * 3 + 0] + m[0 * 3 + 2]) * s;
          result->w = (m[1 * 3 + 2] - m[2 * 3 + 1]) * s;
      } else if (m[1 * 3 + 1] > m[2 * 3 + 2]) {
          float t = -m[0 * 3 + 0] + m[1 * 3 + 1] - m[2 * 3 + 2] + 1.0f;
          float s = XrRcpSqrt(t) * 0.5f;
          result->y = s * t;
          result->x = (m[0 * 3 + 1] + m[1 * 3 + 0]) * s;
          result->w = (m[2 * 3 + 0] - m[0 * 3 + 2]) * s;
          result->z = (m[1 * 3 + 2] + m[2 * 3 + 1]) * s;
      } else {
          float t = -m[0 * 3 + 0] - m[1 * 3 + 1] + m[2 * 3 + 2] + 1.0f;
          float s = XrRcpSqrt(t) * 0.5f;
          result->z = s * t;
          result->w = (m[0 * 3 + 1] - m[1 * 3 + 0]) * s;
          result->x = (m[2 * 3 + 0] + m[0 * 3 + 2]) * s;
          result->y = (m[1 * 3 + 2] + m[2 * 3 + 1]) * s;
      }
  }
  
  // Get the scale from a combined translation(rotation(scale(object))) matrix.
  inline static void XrMatrix4x4f_GetScale(XrVector3f* result, const XrMatrix4x4f* src) {
      assert(XrMatrix4x4f_IsAffine(src, 1e-4f));
      assert(XrMatrix4x4f_IsOrthogonal(src, 1e-4f));
  
      result->x = sqrtf(src->m[0] * src->m[0] + src->m[1] * src->m[1] + src->m[2] * src->m[2]);
      result->y = sqrtf(src->m[4] * src->m[4] + src->m[5] * src->m[5] + src->m[6] * src->m[6]);
      result->z = sqrtf(src->m[8] * src->m[8] + src->m[9] * src->m[9] + src->m[10] * src->m[10]);
  }
  
  // Transforms a 3D vector.
  inline static void XrMatrix4x4f_TransformVector3f(XrVector3f* result, const XrMatrix4x4f* m, const XrVector3f* v) {
      const float w = m->m[3] * v->x + m->m[7] * v->y + m->m[11] * v->z + m->m[15];
      const float rcpW = 1.0f / w;
      result->x = (m->m[0] * v->x + m->m[4] * v->y + m->m[8] * v->z + m->m[12]) * rcpW;
      result->y = (m->m[1] * v->x + m->m[5] * v->y + m->m[9] * v->z + m->m[13]) * rcpW;
      result->z = (m->m[2] * v->x + m->m[6] * v->y + m->m[10] * v->z + m->m[14]) * rcpW;
  }
  
  // Transforms a 4D vector.
  inline static void XrMatrix4x4f_TransformVector4f(XrVector4f* result, const XrMatrix4x4f* m, const XrVector4f* v) {
      result->x = m->m[0] * v->x + m->m[4] * v->y + m->m[8] * v->z + m->m[12] * v->w;
      result->y = m->m[1] * v->x + m->m[5] * v->y + m->m[9] * v->z + m->m[13] * v->w;
      result->z = m->m[2] * v->x + m->m[6] * v->y + m->m[10] * v->z + m->m[14] * v->w;
      result->w = m->m[3] * v->x + m->m[7] * v->y + m->m[11] * v->z + m->m[15] * v->w;
  }
  
  // Transforms the 'mins' and 'maxs' bounds with the given 'matrix'.
  inline static void XrMatrix4x4f_TransformBounds(XrVector3f* resultMins, XrVector3f* resultMaxs, const XrMatrix4x4f* matrix,
                                                  const XrVector3f* mins, const XrVector3f* maxs) {
      assert(XrMatrix4x4f_IsAffine(matrix, 1e-4f));
  
      const XrVector3f center = {(mins->x + maxs->x) * 0.5f, (mins->y + maxs->y) * 0.5f, (mins->z + maxs->z) * 0.5f};
      const XrVector3f extents = {maxs->x - center.x, maxs->y - center.y, maxs->z - center.z};
      const XrVector3f newCenter = {matrix->m[0] * center.x + matrix->m[4] * center.y + matrix->m[8] * center.z + matrix->m[12],
                                    matrix->m[1] * center.x + matrix->m[5] * center.y + matrix->m[9] * center.z + matrix->m[13],
                                    matrix->m[2] * center.x + matrix->m[6] * center.y + matrix->m[10] * center.z + matrix->m[14]};
      const XrVector3f newExtents = {
          fabsf(extents.x * matrix->m[0]) + fabsf(extents.y * matrix->m[4]) + fabsf(extents.z * matrix->m[8]),
          fabsf(extents.x * matrix->m[1]) + fabsf(extents.y * matrix->m[5]) + fabsf(extents.z * matrix->m[9]),
          fabsf(extents.x * matrix->m[2]) + fabsf(extents.y * matrix->m[6]) + fabsf(extents.z * matrix->m[10])};
      XrVector3f_Sub(resultMins, &newCenter, &newExtents);
      XrVector3f_Add(resultMaxs, &newCenter, &newExtents);
  }
  
  // Returns true if the 'mins' and 'maxs' bounds is completely off to one side of the projection matrix.
  inline static bool XrMatrix4x4f_CullBounds(const XrMatrix4x4f* mvp, const XrVector3f* mins, const XrVector3f* maxs) {
      if (maxs->x <= mins->x && maxs->y <= mins->y && maxs->z <= mins->z) {
          return false;
      }
  
      XrVector4f c[8];
      for (int i = 0; i < 8; i++) {
          const XrVector4f corner = {(i & 1) != 0 ? maxs->x : mins->x, (i & 2) != 0 ? maxs->y : mins->y,
                                     (i & 4) != 0 ? maxs->z : mins->z, 1.0f};
          XrMatrix4x4f_TransformVector4f(&c[i], mvp, &corner);
      }
  
      int i;
      for (i = 0; i < 8; i++) {
          if (c[i].x > -c[i].w) {
              break;
          }
      }
      if (i == 8) {
          return true;
      }
      for (i = 0; i < 8; i++) {
          if (c[i].x < c[i].w) {
              break;
          }
      }
      if (i == 8) {
          return true;
      }
  
      for (i = 0; i < 8; i++) {
          if (c[i].y > -c[i].w) {
              break;
          }
      }
      if (i == 8) {
          return true;
      }
      for (i = 0; i < 8; i++) {
          if (c[i].y < c[i].w) {
              break;
          }
      }
      if (i == 8) {
          return true;
      }
      for (i = 0; i < 8; i++) {
          if (c[i].z > -c[i].w) {
              break;
          }
      }
      if (i == 8) {
          return true;
      }
      for (i = 0; i < 8; i++) {
          if (c[i].z < c[i].w) {
              break;
          }
      }
      return i == 8;
  }
  
  #endif  // XR_LINEAR_H_
  