/*
 * This file is part of Afterwarp Framework (http://afterwarp.io).
 * Copyright (c) 2015 - 2025 Dr. Yuriy Kotsarenko. All rights reserved.
 *
 * 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.
 */

// Afterwarp.Vectors.h

#ifndef AFTERWARP_VECTORS_IMPLEMENTATION
  #define AFTERWARP_VECTORS_INTERFACE
#endif

#ifdef AFTERWARP_VECTORS_INTERFACE

#ifndef AFTERWARP_VECTORS_H
#define AFTERWARP_VECTORS_H

#include "Afterwarp.Types.h"
#include "Afterwarp.Structs.h"

// Global macros.

#ifndef MAX
  #define MAX(x, y) (((x) > (y)) ? (x) : (y))
#endif
#ifndef MIN
  #define MIN(x, y) (((x) < (y)) ? (x) : (y))
#endif

/// A special value that determines absolute precision limit when comparing vectors and coordinates.
/// This assumes that the corresponding vectors have their values normalized to [0, 1] range.
extern float const vectorEpsilon;

/// Default font effect parameters provided for calling convenience.
extern FontEffect const fontEffectDefault;

// Color functions.

/// Predefined constant for opaque Black color.
extern Color const colorBlack;

/// Predefined constant for opaque White color.
extern Color const colorWhite;

/// Predefined constant for translucent Black color.
extern Color const colorTranslucentBlack;

/// Predefined constant for translucent White color.
extern Color const colorTranslucentWhite;

// ColorPair constants and functions.

/// Predefined constant for a pair of opaque Black colors.
extern ColorPair const colorPairBlack;

/// Predefined constant for a pair of opaque White colors.
extern ColorPair const colorPairWhite;

/// Predefined constant for a pair of translucent Black colors.
extern ColorPair const colorPairTraslucentBlack;

/// Predefined constant for a pair of translucent White colors.
extern ColorPair const colorPairTraslucentWhite;

/// Creates color pair with the given parameters.
extern ColorPair makeColorPair(Color first, Color second);

/// Creates color pair with both colors equal to the given one.
extern ColorPair colorPair(Color color);

// ColorRect constants.

/// Predefined constant for four opaque Black colors.
extern ColorRect const colorRectBlack;

/// Predefined constant for four opaque White colors.
extern ColorRect const colorRectWhite;

/// Predefined constant for four translucent Black colors.
extern ColorRect const colorRectTraslucentBlack;

/// Predefined constant for four translucent White colors.
extern ColorRect const colorRectTraslucentWhite;

// FloatColorRGB constants and functions.

/// Predefined constant for opaque Black color.
extern FloatColorRGB const floatColorBlackRGB;

/// Predefined constant for opaque White color.
extern FloatColorRGB const floatColorWhiteRGB;

/// Creates floating-point color from the given parameters.
extern FloatColorRGB makeFloatColorRGB(float red, float green, float blue);

/// Creates floating-point color from 32-bit RGBA color value.
extern FloatColorRGB floatColorRGB(Color color);

/// Creates floating-point color with all components equal to the given value.
extern FloatColorRGB floatColorFromValueRGB(float value);

/// Returns color with all components negated.
extern FloatColorRGB floatColorNegateRGB(FloatColorRGB color);

/// Returns color with all components inverted (one minus value).
extern FloatColorRGB floatColorInverseRGB(FloatColorRGB color);

/// Takes floating-point color and clamps each component within [0, 1] range.
extern FloatColorRGB floatColorSaturateRGB(FloatColorRGB color);

/// Add components of the given two colors.
extern FloatColorRGB floatColorAddRGB(FloatColorRGB color1, FloatColorRGB color2);

/// Subtracts components of the second color from the first one.
extern FloatColorRGB floatColorSubtractRGB(FloatColorRGB color1, FloatColorRGB color2);

/// Multiplies components of the given two colors.
extern FloatColorRGB floatColorMultiplyRGB(FloatColorRGB color1, FloatColorRGB color2);

/// Divides components of the first color by the components of the second one.
extern FloatColorRGB floatColorDivideRGB(FloatColorRGB color1, FloatColorRGB color2);

/// Computes average between the two given floating-point color values.
extern FloatColorRGB floatColorAverageRGB(FloatColorRGB color1, FloatColorRGB color2);

/// Interpolates between current and the specified color; theta should be in range of [0, 1].
extern FloatColorRGB floatColorLerpRGB(FloatColorRGB color1, FloatColorRGB color2, float theta);

/// Tests whether the components of two given colors are indistinguishable.
extern bool floatColorEqualsRGB(FloatColorRGB color1, FloatColorRGB color2);

/// Tests whether the components of given color are nearly zero.
extern bool floatColorEmptyRGB(FloatColorRGB color);

/// Returns grayscale value in range of [0, 1] from the given floating-point color value. The resulting
/// value can be considered the color's Luma. Alpha-channel is ignored.
extern float floatColorGrayRGB(FloatColorRGB color);

/// Converts current floating-point color to 32-bit RGBA integer representation.
extern Color floatToColorRGB(FloatColorRGB color);

// FloatColor constants and functions.

/// Predefined constant for opaque Black color.
extern FloatColor const floatColorBlack;

/// Predefined constant for opaque White color.
extern FloatColor const floatColorWhite;

/// Predefined constant for translucent Black color.
extern FloatColor const floatColorTraslucentBlack;

/// Predefined constant for translucent White color.
extern FloatColor const floatColorTraslucentWhite;

/// Creates floating-point color from the given parameters.
extern FloatColor makeFloatColor(float red, float green, float blue, float alpha);

/// Creates floating-point color from 32-bit RGBA color value.
extern FloatColor floatColor(Color color);

/// Creates floating-point color with all components equal to the given value.
extern FloatColor floatColorFromValue(float value);

/// Returns color with all components negated.
extern FloatColor floatColorNegate(FloatColor color);

/// Returns color with all components inverted (one minus value).
extern FloatColor floatColorInverse(FloatColor color);

/// Takes floating-point color with unpremultiplied alpha and multiplies each of red, green, and blue
/// components by its alpha-channel, resulting in premultiplied alpha color.
extern FloatColor floatColorPremultiply(FloatColor color);

/// Takes floating-point color with premultiplied alpha-channel and divides each of its red, green, and
/// blue components by alpha, resulting in unpremultiplied alpha color.
extern FloatColor floatColorUnpremultiply(FloatColor color);

/// Takes floating-point color and clamps each component within [0, 1] range.
extern FloatColor floatColorSaturate(FloatColor color);

/// Add components of the given two colors.
extern FloatColor floatColorAdd(FloatColor color1, FloatColor color2);

/// Subtracts components of the second color from the first one.
extern FloatColor floatColorSubtract(FloatColor color1, FloatColor color2);

/// Multiplies components of the given two colors.
extern FloatColor floatColorMultiply(FloatColor color1, FloatColor color2);

/// Divides components of the first color by the components of the second one.
extern FloatColor floatColorDivide(FloatColor color1, FloatColor color2);

/// Computes average between the two given floating-point color values.
extern FloatColor floatColorAverage(FloatColor color1, FloatColor color2);

/// Interpolates between current and the specified color; theta should be in range of [0, 1].
extern FloatColor floatColorLerp(FloatColor color1, FloatColor color2, float theta);

/// Tests whether the components of two given colors are indistinguishable.
extern bool floatColorEquals(FloatColor color1, FloatColor color2);

/// Tests whether the components of given color are nearly zero.
extern bool floatColorEmpty(FloatColor color);

/// Returns grayscale value in range of [0, 1] from the given floating-point color value. The resulting
/// value can be considered the color's Luma. Alpha-channel is ignored.
extern float floatColorGray(FloatColor color);

/// Converts current floating-point color to 32-bit RGBA integer representation.
extern Color floatToColor(FloatColor color);

/// Returns three components from an existing 4-component floating-point color.
extern FloatColorRGB floatColorToRGB(FloatColor color);

// Point constants and functions.

/// Predefined constant, where X and Y are zero.
extern Point const pointZero;

/// Predefined constant, where X and Y are one.
extern Point const pointUnity;

/// Predefined constant, where X = 1 and Y = 0.
extern Point const pointAxisX;

/// Predefined constant, where X = 0 and Y = 1.
extern Point const pointAxisY;

/// Creates 2D integer point from the given parameters.
extern Point makePoint(int32_t x, int32_t y);

/// Creates 2D integer point with both X and Y set to the given value.
extern Point pointValue(int32_t value);

/// Returns 2D integer point with both components negated.
extern Point pointNegate(Point point);

/// Adds components of two 2D integer points and returns the result.
extern Point pointAdd(Point point1, Point point2);

/// Subtracts components of the second 2D integer point from the first one and returns the result.
extern Point pointSubtract(Point point1, Point point2);

/// Multiplies components of two 2D integer points and returns the result.
extern Point pointMultiply(Point point1, Point point2);

/// Divides components of the first 2D integer point by the components of the second one.
extern Point pointDivide(Point point1, Point point2);

/// Multiplies components of 2D integer point by the given coefficient and returns the result.
extern Point pointRescale(Point point, int32_t theta);

/// Returns length of the given 2D integer point.
extern float pointLength(Point point);

/// Returns distance between two 2D integer points.
extern float pointDistance(Point point1, Point point2);

/// Returns angle (in radians) at which the given 2D integer vector is pointing at, in range of [-Pi, Pi].
extern float pointAngle(Point point);

/// Calculates dot-product between two given 2D integer vectors. The dot-product is an indirect measure
/// of the angle between two vectors.
extern int32_t pointDot(Point point1, Point point2);

/// Calculates a so-called "cross-product" between two given 2D integer vectors, or analog of thereof.
extern int32_t pointCross(Point point1, Point point2);

/// Interpolates between first and second specified 2D integer vectors; theta should be in range of [0, 1].
extern Point pointLerp(Point point1, Point point2, float theta);

/// Tests whether two given 2D integer points match.
extern bool pointEquals(Point point1, Point point2);

/// Tests whether the given 2D integer point has both components equal to zero.
extern bool pointEmpty(Point point);

// PointF constants and functions.

/// Predefined constant, where X and Y are zero.
extern PointF const pointZeroF;

/// Predefined constant, where X and Y are one.
extern PointF const pointUnityF;

/// Predefined constant, where X = 1 and Y = 0.
extern PointF const pointAxisXF;

/// Predefined constant, where X = 0 and Y = 1.
extern PointF const pointAxisYF;

/// Creates floating-point 2D point from the given parameters.
extern PointF makePointF(float x, float y);

/// Creates floating-point 2D point from integer point.
extern PointF pointF(Point point);

/// Creates 2D floating-point vector with both X and Y set to the given value.
extern PointF pointValueF(float value);

/// Returns 2D floating-point vector with both components negated.
extern PointF pointNegateF(PointF point);

/// Adds components of two 2D floating-point vectors and returns the result.
extern PointF pointAddF(PointF point1, PointF point2);

/// Subtracts components of the second 2D floating-point vector from the first one and returns the result.
extern PointF pointSubtractF(PointF point1, PointF point2);

/// Multiplies components of two 2D floating-point vectors and returns the result.
extern PointF pointMultiplyF(PointF point1, PointF point2);

/// Divides components of the first 2D floating-point vector by the components of the second one.
extern PointF pointDivideF(PointF point1, PointF point2);

/// Multiplies components of 2D floating-point vector by the given coefficient and returns the result.
extern PointF pointRescaleF(PointF point, float theta);

/// Returns length of the given 2D floating-point vector.
extern float pointLengthF(PointF point);

/// Returns distance between two 2D floating-point vectors.
extern float pointDistanceF(PointF point1, PointF point2);

/// Returns angle (in radians) at which the given 2D floating-point vector is pointing at, in range
/// of [-Pi, Pi].
float pointAngleF(PointF point);

/// Calculates dot-product between two given 2D floating-point vectors. The dot-product is an indirect
/// measure of the angle between two vectors.
float pointDotF(PointF point1, PointF point2);

/// Calculates a so-called "cross-product" between two given 2D floating-point vectors, or analog of thereof.
float pointCrossF(PointF point1, PointF point2);

/// Normalizes given 2D floating-point vector to unity length. If the vector has zero length, it will be
/// returned unchanged.
PointF pointNormalizeF(PointF point);

/// Interpolates between first and second specified 2D floating-point vectors; theta should be in range
/// of [0, 1].
PointF pointLerpF(PointF point1, PointF point2, float theta);

/// Tests whether two given 2D floating-point vectors are indistinguishable.
bool pointEqualsF(PointF point1, PointF point2);

/// Tests whether the given 2D floating-point vector has both components nearly zero.
bool pointEmptyF(PointF point);

/// Multiplies given vector treated as 3D vector with Z = 1, by 3x2 matrix, treated as 3x3 with right side
/// matching identity matrix.
PointF pointTransformF(PointF point, Matrix3x2 matrix);

// Vector constants and functions.

/// Predefined constant, where X, Y and Z are zero.
extern Vector const vectorZero;

/// Predefined constant, where X, Y and Z are one.
extern Vector const vectorUnity;

/// Predefined constant, where X = 1, Y = 0 and Z = 0.
extern Vector const vectorAxisX;

/// Predefined constant, where X = 0, Y = 1 and Z = 0.
extern Vector const vectorAxisY;

/// Predefined constant, where X = 0, Y = 0 and Z = 1.
extern Vector const vectorAxisZ;

/// Creates 3D vector from the given parameters.
Vector makeVector(float x, float y, float z);

/// Creates 3D vector with X, Y and Z set to the given value.
Vector vectorValue(float value);

/// Returns 3D vector with all components negated.
Vector vectorNegate(Vector vector);

/// Adds components of two 3D vectors and returns the result.
Vector vectorAdd(Vector vector1, Vector vector2);

/// Subtracts components of the second 3D vector from the first one and returns the result.
Vector vectorSubtract(Vector vector1, Vector vector2);

/// Multiplies components of two 3D vectors and returns the result.
Vector vectorMultiply(Vector vector1, Vector vector2);

/// Divides components of the first 3D vector by the components of the second one.
Vector vectorDivide(Vector vector1, Vector vector2);

/// Multiplies components of 3D vector by the given coefficient and returns the result.
Vector vectorRescale(Vector vector, float theta);

/// Calculates dot-product between two given 3D vectors. The dot-product is an indirect measure of the angle
/// between two vectors.
float vectorDot(Vector vector1, Vector vector2);

/// Calculates cross product between two given 3D vectors. The resulting vector is perpendicular to both
/// vectors and normal to the plane containing them.
Vector vectorCross(Vector vector1, Vector vector2);

/// Returns length of the given 3D vector.
float vectorLength(Vector vector);

/// Returns distance between two 3D vectors.
float vectorDistance(Vector vector1, Vector vector2);

/// Calculates angle between two 3D vectors. The returned value has range of [0, Pi].
float vectorAngle(Vector vector1, Vector vector2);

/// Normalizes given 3D vector to unity length. If the vector has zero length, it will be returned unchanged.
Vector vectorNormalize(Vector vector);

/// Calculates portion of given 3D vector that is parallel to the direction vector.
Vector vectorParallel(Vector vector, Vector direction);

/// Calculates portion of given 3D vector that is perpendicular to the direction vector.
Vector vectorPerpendicular(Vector vector, Vector direction);

/// Calculates 3D vector that is a reflection of given vector from surface given by the specified normal.
Vector vectorReflect(Vector vector, Vector normal);

/// Interpolates between first and second specified 3D vectors; theta should be in range of [0, 1].
Vector vectorLerp(Vector vector1, Vector vector2, float theta);

/// Tests whether two given 3D vectors are indistinguishable.
bool vectorEquals(Vector vector1, Vector vector2);

/// Tests whether the given 3D vector has all components nearly zero.
bool vectorEmpty(Vector vector);

/// Multiplies given 3D vector as if it was 4D (3D + W) vector by 4x4 matrix and returns the result,
/// discarding the fourth element.
Vector vectorTransform(Vector vector, float w, Matrix matrix);

/// Multiplies given 3D vector as if it was 4D (3D + W) vector by 4x4 matrix, then divides each of the
/// resulting components by the resulting value of W (which is then discarded), returning X, Y and Z.
Vector vectorProject(Vector vector, float w, Matrix matrix);

/// Multiplies given 3D vector as if it was 4D (3D + 1) vector by 4x4 matrix, then divides each of the
/// resulting components by the resulting value of W (which is then discarded), and converts absolute
/// coordinates to screen pixels by using the given target size.
PointF vectorProjectTarget(Vector vector, PointF targetSize, Matrix matrix);

// Vector4 constants and functions.

/// Predefined constant, where X, Y, Z and W are zero.
extern Vector4 const vectorZero3D;

/// Predefined constant, where X, Y, Z and W are one.
extern Vector4 const vectorUnity3D;

/// Predefined constant, where X = 1, Y = 0, Z = 0 and W = 0.
extern Vector4 const vectorAxisX3D;

/// Predefined constant, where X = 0, Y = 1, Z = 0 and W = 0.
extern Vector4 const vectorAxisY3D;

/// Predefined constant, where X = 0, Y = 0, Z = 1 and W = 0.
extern Vector4 const vectorAxisZ3D;

/// Predefined constant, where X = 0, Y = 0, Z = 0 and W = 1.
extern Vector4 const vectorAxisW3D;

/// Creates 4-component vector from the given parameters.
Vector4 makeVector3D(float x, float y, float z, float w);

/// Creates 4-component vector with X, Y, Z and W set to the given value.
Vector4 vectorFromValue3D(float value);

/// Returns 4-component vector with all components negated.
Vector4 vectorNegate3D(Vector4 vector);

/// Adds components of two 4-component vectors and returns the result.
Vector4 vectorAdd3D(Vector4 vector1, Vector4 vector2);

/// Subtracts components of the second 4-component vector from the first one and returns the result.
Vector4 vectorSubtract3D(Vector4 vector1, Vector4 vector2);

/// Multiplies components of two 4-component vectors and returns the result.
Vector4 vectorMultiply3D(Vector4 vector1, Vector4 vector2);

/// Divides components of the first 4-component vector by the components of the second one.
Vector4 vectorDivide3D(Vector4 vector1, Vector4 vector2);

/// Multiplies components of 4-component vector by the given coefficient and returns the result.
Vector4 vectorRescale3D(Vector4 vector, float theta);

/// Calculates dot-product between two given 4-component vectors. The dot-product is an indirect measure of
/// the angle between two vectors.
float vectorDot3D(Vector4 vector1, Vector4 vector2);

/// Returns length of the given 4-component vector.
float vectorLength3D(Vector4 vector);

/// Returns distance between two 4-component vectors.
float vectorDistance3D(Vector4 vector1, Vector4 vector2);

/// Normalizes given 4-component vector to unity length. If the vector has zero length, it will be returned
/// unchanged.
Vector4 vectorNormalize3D(Vector4 vector);

/// Interpolates between first and second specified 4-component vectors; theta should be in range of [0, 1].
Vector4 vectorLerp3D(Vector4 vector1, Vector4 vector2, float theta);

/// Tests whether two given 4-component vectors are indistinguishable.
bool vectorEquals3D(Vector4 vector1, Vector4 vector2);

/// Tests whether the given 4-component vector has all components nearly zero.
bool vectorEmpty3D(Vector4 vector);

/// Multiplies given 4-component vector by the given 4x4 matrix.
Vector4 vectorTransform3D(Vector4 vector, Matrix matrix);

// Matrix3x2 constants and functions.

/// Predefined constant, where all matrix values are zero.
extern Matrix3x2 const matrixZero3x2;

/// Predefined constant with values corresponding to an identity matrix.
extern Matrix3x2 const matrixIdentity3x2;

/// Multiplies two 3x2 matrices and returns the result.
Matrix3x2 matrixMultiply3x2(Matrix3x2 matrix1, Matrix3x2 matrix2);

/// Multiplies 3x2 matrix by a coefficient and returns the result.
Matrix3x2 matrixRescale3x2(Matrix3x2 matrix, float theta);

/// Calculates determinant of current matrix (as if it was 3x2 with right column being identity).
float matrixDeterminant3x2(Matrix3x2 matrix);

/// Calculates inverse matrix of the current matrix (as if it was 3x2 with right column being identity).
Matrix3x2 matrixInverse3x2(Matrix3x2 matrix);

/// Calculates inverse and transpose matrix of the current matrix as if it was 2x2.
Matrix3x2 matrixInverseTranspose3x2(Matrix3x2 matrix);

/// Creates 2D scale 3x2 matrix with the specified coefficients.
Matrix3x2 matrixScale3x2(PointF scale);

/// Creates 2D scale 3x2 matrix around a specific center point.
Matrix3x2 matrixScaleBy3x2(PointF center, PointF scale);

/// Creates 2D translation 3x2 matrix with the specified offset.
Matrix3x2 matrixTranslate3x2(PointF offset);

/// Creates 2D rotation 3x2 matrix with specified angle (in radians).
Matrix3x2 matrixRotate3x2(float angle);

/// Creates 2D rotation 3x2 matrix with specified angle (in radians) around a specific center point.
Matrix3x2 matrixRotateBy3x2(PointF center, float angle);

/// Tests whether both 3x2 matrices are indistinguishable.
bool matrixEquals3x2(Matrix3x2 matrix1, Matrix3x2 matrix2);

/// Tests whether given 3x2 matrix has all components nearly zero.
bool matrixEmpty3x2(Matrix3x2 matrix);

// Matrix constants and functions.

/// Predefined constant, where all matrix values are zero.
extern Matrix const matrixZero;

/// Predefined constant with values corresponding to an identity matrix.
extern Matrix const matrixIdentity;

/// Retrieves 3x2 portion from the given matrix.
Matrix3x2 matrixGet3x2(Matrix matrix);

/// Multiplies two 4x4 matrices and returns the result.
Matrix matrixMultiply(Matrix matrix1, Matrix matrix2);

/// Multiplies 4x4 matrix by a coefficient and returns the result.
Matrix matrixRescale(Matrix matrix, float theta);

/// Calculates determinant of the given 4x4 matrix.
float matrixDeterminant(Matrix matrix);

/// Calculates transposed matrix from the given 4x4 matrix.
Matrix matrixTranspose(Matrix matrix);

/// Calculates adjoint matrix from the given 4x4 matrix.
Matrix matrixAdjoint(Matrix matrix);

/// Calculates inverse matrix from the given 4x4 matrix.
Matrix matrixInverse(Matrix matrix);

/// Copies 3x3 portion of source matrix while leaving remaining portion equal to identity 4x4 matrix.
Matrix matrixSub3x3(Matrix matrix);

/// Assuming that the given 4x4 matrix is a view matrix, calculates the 3D position where the camera
/// (or "eye") is supposedly located.
Vector matrixEyePosition(Matrix matrix);

/// Assuming that the given 4x4 matrix is a world matrix, calculates the 3D position where the object
/// (or "world") is supposedly located.
Vector matrixWorldPosition(Matrix matrix);

/// Creates 3D translation 4x4 matrix with specified offset.
Matrix matrixTranslate(Vector offset);

/// Creates 3D scale 4x4 matrix with specified coefficients.
Matrix matrixScale(Vector scale);

/// Creates 3D rotation 4x4 matrix around X axis with the specified angle.
Matrix matrixRotateX(float angle);

/// Creates 3D rotation 4x4 matrix around Y axis with the specified angle.
Matrix matrixRotateY(float angle);

/// Creates 3D rotation 4x4 matrix around Z axis with the specified angle.
Matrix matrixRotateZ(float angle);

/// Creates 3D rotation 4x4 matrix around specified axis and angle (in radians).
Matrix matrixRotate(Vector axis, float angle);

/// Creates 3D rotation 4x4 matrix based on parameters similar to flight dynamics, specifically Heading,
/// Pitch and Bank. The components are taken from the specified vector with Y corresponding to Heading,
/// X to Pitch and Z to Bank.
Matrix matrixHeadingPitchBank(Vector angles);

/// Creates 3D rotation 4x4 matrix based on parameters similar to flight dynamics, specifically Yaw, Pitch
/// and Roll. The components are taken from the specified vector with Y corresponding to Yaw, X to Pitch and
/// Z to Roll.
Matrix matrixYawPitchRoll(Vector angles);

/// Creates a reflection matrix specified by the given vector defining orientation of that reflection.
Matrix matrixReflect(Vector axis);

/// Creates a view matrix that is defined by the camera's position, its target and vertical axis
/// or "ceiling".
Matrix matrixLookAt(Vector origin, Vector target, Vector ceiling);

/// Creates perspective projection matrix defined by a field of view on Y axis. This is a common way for
/// typical 3D applications. In 3D shooters special care is to be taken because on wide-screen monitors
/// the visible area will be bigger when using this method. The parameters that define the viewed range
/// are important for defining the precision of the depth transformation or a depth-buffer.
///   FieldOfView  -  The camera's field of view in radians. For example Pi/4.
///   AspectRatio  -  The screen's aspect ratio. Can be calculated as y/x.
///   MinRange     -  The closest range at which the scene will be viewed.
///   MaxRange     -  The farthest range at which the scene will be viewed.
Matrix matrixPerspectiveFOVY(float fieldOfView, float aspectRatio, float minRange, float maxRange,
  bool negativeDepthRange);

/// Creates perspective projection matrix defined by a field of view on X axis. In 3D shooters the field
/// of view needs to be adjusted to allow more visible area on wide-screen monitors. The parameters that
/// define the viewed range are important for defining the precision of the depth transformation or a
/// depth-buffer.
///   FieldOfView  -  The camera's field of view in radians. For example Pi/4.
///   AspectRatio  -  The screen's aspect ratio. Can be calculated as x/y.
///   MinRange     -  The closest range at which the scene will be viewed.
///   MaxRange     -  The farthest range at which the scene will be viewed.
Matrix matrixPerspectiveFOVX(float fieldOfView, float aspectRatio, float minRange, float maxRange,
  bool negativeDepthRange);

/// Creates perspective projection matrix defined by the viewing volume in 3D space.
Matrix matrixPerspectiveVOL(float width, float height, float minRange, float maxRange);

/// Creates perspective projection matrix defined by the individual axis's boundaries.
Matrix matrixPerspectiveBDS(float left, float top, float right, float bottom, float minRange,
  float maxRange);

/// Creates orthogonal projection matrix defined by the viewing volume in 3D space.
Matrix matrixOrthogonalVOL(float width, float height, float minRange, float maxRange);

/// Creates orthogonal projection matrix defined by the individual axis's boundaries.
Matrix matrixOrthogonalBDS(float left, float top, float right, float bottom, float minRange,
  float maxRange);

/// Tests whether both 4x4 matrices are indistinguishable.
bool matrixEquals(Matrix matrix1, Matrix matrix2);

/// Tests whether given 4x4 matrix has all components nearly zero.
bool matrixEmpty(Matrix matrix);

// Quaternion constants and functions.

/// Predefined constant with values corresponding to an identity quaternion.
extern Quaternion const quaternionIdentity;

/// Multiplies quaternion by another quaternion, effectively combining their transformations.
Quaternion quaternionMultiply(Quaternion quat1, Quaternion quat2);

/// Computes dot product between two given quaternions.
float quaternionDot(Quaternion quat1, Quaternion quat2);

/// Returns the magnitude of the given quaternion.
float quaternionLength(Quaternion quat);

/// Returns rotational angle that is present in the given quaternion.
float quaternionAngle(Quaternion quat);

/// Returns rotational axis that is present in the given quaternion.
Vector quaternionAxis(Quaternion quat);

/// Normalizes the given quaternion. Note that normally quaternions are always normalized (of course, within
/// limits of numerical precision). This function is provided mainly to combat floating point "error creep",
/// which occurs after many successive quaternion operations.
Quaternion quaternionNormalize(Quaternion quat);

/// Computes conjugate of the given quaternion. The resulting quaternion has opposite rotation.
Quaternion quaternionConjugate(Quaternion quat);

/// Computes exponentiation of the given quaternion.
Quaternion quaternionExponentiate(Quaternion quat, float exponent);

/// Applies spherical linear interpolation between two given quaternions.
Quaternion quaternionSlerp(Quaternion quat1, Quaternion quat2, float theta);

/// Converts the given quaternion into 4x4 matrix representation.
Matrix quaternionMatrix(Quaternion quat);

/// Creates quaternion with transformation converted from the given matrix.
Quaternion matrixQuaternion(Matrix matrix);

/// Creates 3D quaternion containing rotation around X axis with given angle (in radians).
Quaternion quaternionRotateX(float angle);

/// Creates 3D quaternion containing rotation around Y axis with given angle (in radians).
Quaternion quaternionRotateY(float angle);

/// Creates 3D quaternion containing rotation around Z axis with given angle (in radians).
Quaternion quaternionRotateZ(float angle);

/// Creates 3D quaternion containing rotation around an arbitrary axis with given angle (in radians).
Quaternion quaternionRotate(Vector axis, float angle);

/// Creates 3D quaternion setup to perform Object-To-Inertial rotation using the angles specified in
/// Euler format. The components are taken from the specified vector with Y corresponding to Heading,
/// X to Pitch and Z to Bank.
Quaternion quaternionRotateObjectToInertial(Vector angles);

/// Creates 3D quaternion setup to perform Inertial-To-Object rotation using the angles specified in
/// Euler format. The components are taken from the specified vector with Y corresponding to Heading,
/// X to Pitch and Z to Bank.
Quaternion quaternionRotateInertialToObject(Vector angles);

// Rect constants and functions.

/// Empty rectangle with all members set to zero.
extern Rect const rectZero;

/// Creates rectangle with the given parameters.
Rect makeRect(int32_t left, int32_t top, int32_t width, int32_t height);

/// Creates rectangle from the given boundaries.
Rect makeBounds(int32_t left, int32_t top, int32_t right, int32_t bottom);

/// Returns the size of given rectangle as 2D integer vector.
Point rectGetSize(Rect rect);

/// Specifies new size for given rectangle from values of 2D integer vector.
Rect rectSetSize(Rect rect, Point size);

/// Returns right edge of the given rectangle.
int32_t rectGetRight(Rect rect);

/// Specifies new right edge for the given rectangle.
Rect rectSetRight(Rect rect, int32_t right);

/// Returns bottom edge of the given rectangle.
int32_t rectGetBottom(Rect rect);

/// Specifies new bottom edge for the given rectangle.
Rect rectSetBottom(Rect rect, int32_t bottom);

/// Returns top and left of the given rectangle as 2D integer vector.
Point rectGetTopLeft(Rect rect);

/// Specifies Top and Left for the given rectangle from values of 2D integer vector.
Rect rectSetTopLeft(Rect rect, Point topLeft);

/// Returns Bottom and Right of the given rectangle as 2D integer vector.
Point rectGetBottomRight(Rect rect);

/// Specifies Bottom and Right for the given rectangle from values of 2D integer vector.
Rect rectSetBottomRight(Rect rect, Point bottomRight);

/// Returns top and right of the given rectangle as 2D integer vector.
Point rectGetTopRight(Rect rect);

/// Returns bottom and left of the given rectangle as 2D integer vector.
Point rectGetBottomLeft(Rect rect);

/// Tests whether both given rectangles are exactly the same.
bool rectEquals(Rect rect1, Rect rect2);

/// Tests whether the given rectangle is empty, that is, having width or height of zero or less.
bool rectEmpty(Rect rect);

/// Tests whether the specified point is contained within the given rectangle.
bool pointInRect(Point point, Rect rect);

/// Tests whether the specified rectangle is contained within the given rectangle.
bool rectInRect(Rect rect, Rect largerRect);

/// Tests whether two given rectangles overlap.
bool rectOverlap(Rect rect1, Rect rect2);

/// Calculates rectangle that results from intersection between two given rectangles.
Rect rectIntersect(Rect rect1, Rect rect2);

/// Calculates rectangle that results from union between two given rectangles.
Rect rectJoin(Rect rect1, Rect rect2);

/// Calculates displaced rectangle by certain offset.
Rect rectOffset(Rect rect, Point delta);

/// Returns rectangle with left and top decremented, while right and bottom incremented by given offset.
Rect rectInflate(Rect rect, Point delta);

// RectF constants and functions.

/// Empty rectangle with all members set to zero.
extern RectF const rectZeroF;

/// Creates rectangle with the given parameters.
RectF makeRectF(float left, float top, float width, float height);

/// Creates rectangle from the given boundaries.
RectF makeBoundsF(float left, float top, float right, float bottom);

/// Creates rectangle with floating-point coordinates from rectangle with integer coordinates.
RectF rectF(Rect rect);

/// Creates rectangle with integer coordinates from rectangle with floating-point coordinates.
Rect rectToIntRectF(RectF rect);

/// Returns the size of given rectangle as 2D integer vector.
PointF rectGetSizeF(RectF rect);

/// Specifies new size for given rectangle from values of 2D integer vector.
RectF rectSetSizeF(RectF rect, PointF size);

/// Returns right edge of the given rectangle.
float rectGetRightF(RectF rect);

/// Specifies new right edge for the given rectangle.
RectF rectSetRightF(RectF rect, float right);

/// Returns bottom edge of the given rectangle.
float rectGetBottomF(RectF rect);

/// Specifies new bottom edge for the given rectangle.
RectF rectSetBottomF(RectF rect, float bottom);

/// Returns top and left of the given rectangle as 2D integer vector.
PointF rectGetTopLeftF(RectF rect);

/// Specifies Top and Left for the given rectangle from values of 2D integer vector.
RectF rectSetTopLeftF(RectF rect, PointF topLeft);

/// Returns Bottom and Right of the given rectangle as 2D integer vector.
PointF rectGetBottomRightF(RectF rect);

/// Specifies Bottom and Right for the given rectangle from values of 2D integer vector.
RectF rectSetBottomRightF(RectF rect, PointF bottomRight);

/// Returns top and right of the given rectangle as 2D integer vector.
PointF rectGetTopRightF(RectF rect);

/// Returns bottom and left of the given rectangle as 2D integer vector.
PointF rectGetBottomLeftF(RectF rect);

/// Tests whether both given rectangles are nearly same..
bool rectEqualsF(RectF rect1, RectF rect2);

/// Tests whether the given rectangle is empty, that is, having width or height of nearly zero or less.
bool rectEmptyF(RectF rect);

/// Tests whether the specified point is contained within the given rectangle.
bool pointInRectF(PointF point, RectF rect);

/// Tests whether the specified rectangle is contained within the given rectangle.
bool rectInRectF(RectF rect, RectF largerRect);

/// Tests whether two given rectangles overlap.
bool rectOverlapF(RectF rect1, RectF rect2);

/// Calculates rectangle that results from intersection between two given rectangles.
RectF rectIntersectF(RectF rect1, RectF rect2);

/// Calculates rectangle that results from union between two given rectangles.
RectF rectJoinF(RectF rect1, RectF rect2);

/// Calculates displaced rectangle by certain offset.
RectF rectOffsetF(RectF rect, PointF delta);

/// Returns rectangle with left and top decremented, while right and bottom incremented by given offset.
RectF rectInflateF(RectF rect, PointF delta);

// Quad constants and functions.

/// Quadrilateral with all vertices set to zero.
extern Quad const quadZero;

/// Quadrilateral with vertices set at four corners set between (0, 0) and (1, 1).
extern Quad const quadUnity;

/// Creates quadrilateral with the given parameters.
Quad makeQuadWith(PointF topLeft, PointF topRight, PointF bottomRight, PointF bottomLeft);

/// Creates quadrilateral as rectangle with the given parameters.
Quad makeQuad(float left, float top, float width, float height);

/// Creates quadrilateral with the specified top left corner and the given dimensions, which are scaled
/// by the provided coefficient.
Quad makeQuadScaled(float left, float top, float width, float height, float scale, bool centered);

/// Creates quadrilateral specified by its dimensions. The rectangle is then rotated and scaled around
/// the specified middle point (assumed to be inside rectangle's dimensions) and placed in center of the
/// specified origin.
Quad makeQuadRotatedAt(PointF origin, PointF size, PointF center, float angle, float scale);

/// Creates quadrilateral specified by its dimensions. The rectangle is then rotated and scaled around its
/// center and placed at the specified origin.
Quad makeQuadRotated(PointF origin, PointF size, float angle, float scale);

/// Creates quadrilateral specified by top-left corner and size. The rectangle is then rotated and scaled
/// around the specified middle point (assumed to be inside rectangle's dimensions) and placed in the
/// center of the specified origin. The difference between this method and @link(makeQuadRotatedAt) is that
/// the rotation does not preserve centering of the rectangle in case where middle point is not actually
/// located in the middle.
Quad makeQuadRotatedTL(PointF topLeft, PointF size, PointF center, float angle, float scale);

/// Creates a horizontally skewed quadrilateral, oriented around its origin and with the given dimensions.
/// A skew angle defines quadrilateral's top/left corner.
Quad makeQuadSkewedHoriz(PointF origin, PointF size, float angle);

/// Creates a vertically skewed quadrilateral, oriented around its origin and with the given dimensions.
/// A skew angle defines quadrilateral's top/left corner.
Quad makeQuadSkewedVert(PointF origin, PointF size, float angle);

/// Creates quadrilateral from rectangle with integer coordinates.
Quad quadFromRect(Rect rect);

/// Creates quadrilateral from rectangle with floating-point coordinates.
Quad quadFromRectF(RectF rect);

/// Rescales vertices of the given quadrilateral by the provided coefficient, optionally centering them
/// around zero origin.
Quad quadScale(Quad quad, float scale, bool centered);

/// Creates quadrilateral from the given one but having left vertices exchanged with the right ones,
/// effectively mirroring it horizontally.
Quad quadMirror(Quad quad);

/// Creates quadrilateral from the given one but having top vertices exchanged with the bottom ones,
/// effectively flipping it vertically.
Quad quadFlip(Quad quad);

/// Transforms (multiplies) vertices of the given quadrilateral by the specified matrix.
Quad quadTransform(Quad quad, Matrix3x2 matrix);

/// Displaces vertices of the given quadrilateral by the specified offset.
Quad quadOffset(Quad quad, PointF delta);

#endif // AFTERWARP_VECTORS_H
#endif // AFTERWARP_VECTORS_INTERFACE
#ifdef AFTERWARP_VECTORS_IMPLEMENTATION

#include <math.h>
#include <string.h>
#include <stdint.h>

// Forward declarations.

static bool nearlyEqualFloats(float left, float right);
static bool nearlyZeroFloat(float value);
static float saturateFloat(float value, float min, float max);
static float lerpFloat(float value1, float value2, float theta);

// Global constants.

float const vectorEpsilon = 1E-6f;

// Font constants.

FontEffect const fontEffectDefault = {1.0f, 1.0f, FONT_BORDER_NONE, 1.0f, 0.25f, 0.75f, 3.0f,
  {2.0f, 2.0f}, 0.15f, 0.0f, 0.0f};

// Color constants and functions.

Color const colorBlack = 0xFF000000;
Color const colorWhite = 0xFFFFFFFF;
Color const colorTranslucentBlack = 0x00000000;
Color const colorTranslucentWhite = 0x00FFFFFF;

// ColorPair constants and functions.

ColorPair const colorPairBlack = {.first = 0xFF000000, .second = 0xFF000000};
ColorPair const colorPairWhite = {.first = 0xFFFFFFFF, .second = 0xFFFFFFFF};
ColorPair const colorPairTraslucentBlack = {.first = 0x00000000, .second = 0x00000000};
ColorPair const colorPairTraslucentWhite = {.first = 0x00FFFFFF, .second = 0x00FFFFFF};

ColorPair makeColorPair(Color first, Color second)
{
  return (ColorPair){.first = first, .second = second};
}

ColorPair colorPair(Color color)
{
  return (ColorPair){.first = color, .second = color};
}

// ColorRect constants.

ColorRect const colorRectBlack = {0xFF000000, 0xFF000000, 0xFF000000, 0xFF000000};
ColorRect const colorRectWhite = {0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF};
ColorRect const colorRectTraslucentBlack = {0x00000000, 0x00000000, 0x00000000, 0x00000000};
ColorRect const colorRectTraslucentWhite = {0x00FFFFFF, 0x00FFFFFF, 0x00FFFFFF, 0x00FFFFFF};

// FloatColorRGB constants.

FloatColorRGB const floatColorBlackRGB = {.red = 0.0f, .green = 0.0f, .blue = 0.0f};

FloatColorRGB const floatColorWhiteRGB = {.red = 1.0f, .green = 1.0f, .blue = 1.0f};

// FloatColorRGB functions.

FloatColorRGB makeFloatColorRGB(float red, float green, float blue)
{
  return (FloatColorRGB){.red = red, .green = green, .blue = blue};
}

FloatColorRGB floatColorRGB(Color color)
{
  return (FloatColorRGB){
    .red = (color & 0xFFu) / 255.0f,
    .green = ((color >> 8) & 0xFFu) / 255.0f,
    .blue = ((color >> 16) & 0xFFu) / 255.0f};
}

FloatColorRGB floatColorFromValueRGB(float value)
{
  return (FloatColorRGB){
    .red = value,
    .green = value,
    .blue = value};
}

FloatColorRGB floatColorNegateRGB(FloatColorRGB color)
{
  return (FloatColorRGB){
    .red = -color.red,
    .green = -color.green,
    .blue = -color.blue};
}

FloatColorRGB floatColorInverseRGB(FloatColorRGB color)
{
  return (FloatColorRGB){
    .red = 1.0f - color.red,
    .green = 1.0f - color.green,
    .blue = 1.0f - color.blue};
}

FloatColorRGB floatColorSaturateRGB(FloatColorRGB color)
{
  return (FloatColorRGB){
    .red = saturateFloat(color.red, 0.0f, 1.0f),
    .green = saturateFloat(color.green, 0.0f, 1.0f),
    .blue = saturateFloat(color.blue, 0.0f, 1.0f)};
}

FloatColorRGB floatColorAddRGB(FloatColorRGB color1, FloatColorRGB color2)
{
  return (FloatColorRGB){
    .red = color1.red + color2.red,
    .green = color1.green + color2.red,
    .blue = color1.blue + color2.blue};
}

FloatColorRGB floatColorSubtractRGB(FloatColorRGB color1, FloatColorRGB color2)
{
  return (FloatColorRGB){
    .red = color1.red - color2.red,
    .green = color1.green - color2.red,
    .blue = color1.blue - color2.blue};
}

FloatColorRGB floatColorMultiplyRGB(FloatColorRGB color1, FloatColorRGB color2)
{
  return (FloatColorRGB){
    .red = color1.red * color2.red,
    .green = color1.green * color2.red,
    .blue = color1.blue * color2.blue};
}

FloatColorRGB floatColorDivideRGB(FloatColorRGB color1, FloatColorRGB color2)
{
  return (FloatColorRGB){
    .red = color1.red / color2.red,
    .green = color1.green / color2.red,
    .blue = color1.blue / color2.blue};
}

FloatColorRGB floatColorAverageRGB(FloatColorRGB color1, FloatColorRGB color2)
{
  return (FloatColorRGB){
    .red = (color1.red + color2.red) * 0.5f,
    .green = (color1.green + color2.red) * 0.5f,
    .blue = (color1.blue + color2.blue) * 0.5f};
}

FloatColorRGB floatColorLerpRGB(FloatColorRGB color1, FloatColorRGB color2, float theta)
{
  return (FloatColorRGB){
    .red = lerpFloat(color1.red, color2.red, theta),
    .green = lerpFloat(color1.green, color2.red, theta),
    .blue = lerpFloat(color1.blue, color2.blue, theta)};
}

bool floatColorEqualsRGB(FloatColorRGB color1, FloatColorRGB color2)
{
  return
    nearlyEqualFloats(color1.red, color2.red) &&
    nearlyEqualFloats(color1.green, color2.green) &&
    nearlyEqualFloats(color1.blue, color2.blue);
}

bool floatColorEmptyRGB(FloatColorRGB color)
{
  return
    nearlyZeroFloat(color.red) &&
    nearlyZeroFloat(color.green) &&
    nearlyZeroFloat(color.blue);
}

float floatColorGrayRGB(FloatColorRGB color)
{
  return color.red * 0.29889531f + color.green * 0.58662247f + color.blue * 0.11448222f;
}

Color floatToColorRGB(FloatColorRGB color)
{
  return (Color)roundf(color.red * 255.0f) | ((Color)roundf(color.green * 255.0f) << 8) |
    ((Color)roundf(color.blue * 255.0f) << 16) | 0xFF000000u;
}

// FloatColor constants.

FloatColor const floatColorBlack = {.red = 0.0f, .green = 0.0f, .blue = 0.0f, .alpha = 1.0f};

FloatColor const floatColorWhite = {.red = 1.0f, .green = 1.0f, .blue = 1.0f, .alpha = 1.0f};

FloatColor const floatColorTraslucentBlack = {.red = 0.0f, .green = 0.0f, .blue = 0.0f, .alpha = 0.0f};

FloatColor const floatColorTraslucentWhite = {.red = 1.0f, .green = 1.0f, .blue = 1.0f, .alpha = 0.0f};

// FloatColor functions.

FloatColor makeFloatColor(float red, float green, float blue, float alpha)
{
  return (FloatColor){.red = red, .green = green, .blue = blue, .alpha = alpha};
}

FloatColor floatColor(Color color)
{
  return (FloatColor){
    .red = (color & 0xFFu) / 255.0f,
    .green = ((color >> 8) & 0xFFu) / 255.0f,
    .blue = ((color >> 16) & 0xFFu) / 255.0f,
    .alpha = ((color >> 24) & 0xFFu) / 255.0f};
}

FloatColor floatColorFromValue(float value)
{
  return (FloatColor){
    .red = value,
    .green = value,
    .blue = value,
    .alpha = value};
}

FloatColor floatColorNegate(FloatColor color)
{
  return (FloatColor){
    .red = -color.red,
    .green = -color.green,
    .blue = -color.blue,
    .alpha = -color.alpha};
}

FloatColor floatColorInverse(FloatColor color)
{
  return (FloatColor){
    .red = 1.0f - color.red,
    .green = 1.0f - color.green,
    .blue = 1.0f - color.blue,
    .alpha = 1.0f - color.alpha};
}

FloatColor floatColorPremultiply(FloatColor color)
{
  return (FloatColor){
    .red = color.red * color.alpha,
    .green = color.green * color.alpha,
    .blue = color.blue * color.alpha,
    .alpha = color.alpha};
}

FloatColor floatColorUnpremultiply(FloatColor color)
{
  if (fabsf(color.alpha) >= vectorEpsilon)
  {
    float const invAlpha = 1.0f / color.alpha;

    return (FloatColor){
      .red = color.red * invAlpha,
      .green = color.green * invAlpha,
      .blue = color.blue * invAlpha,
      .alpha = color.alpha * invAlpha};
  }
  else
    return color;
}

FloatColor floatColorSaturate(FloatColor color)
{
  return (FloatColor){
    .red = saturateFloat(color.red, 0.0f, 1.0f),
    .green = saturateFloat(color.green, 0.0f, 1.0f),
    .blue = saturateFloat(color.blue, 0.0f, 1.0f),
    .alpha = saturateFloat(color.alpha, 0.0f, 1.0f)};
}

FloatColor floatColorAdd(FloatColor color1, FloatColor color2)
{
  return (FloatColor){
    .red = color1.red + color2.red,
    .green = color1.green + color2.green,
    .blue = color1.blue + color2.blue,
    .alpha = color1.alpha + color2.alpha};
}

FloatColor floatColorSubtract(FloatColor color1, FloatColor color2)
{
  return (FloatColor){
    .red = color1.red - color2.red,
    .green = color1.green - color2.green,
    .blue = color1.blue - color2.blue,
    .alpha = color1.alpha - color2.alpha};
}

FloatColor floatColorMultiply(FloatColor color1, FloatColor color2)
{
  return (FloatColor){
    .red = color1.red * color2.red,
    .green = color1.green * color2.green,
    .blue = color1.blue * color2.blue,
    .alpha = color1.alpha * color2.alpha};
}

FloatColor floatColorDivide(FloatColor color1, FloatColor color2)
{
  return (FloatColor){
    .red = color1.red / color2.red,
    .green = color1.green / color2.green,
    .blue = color1.blue / color2.blue,
    .alpha = color1.alpha / color2.alpha};
}

FloatColor floatColorAverage(FloatColor color1, FloatColor color2)
{
  return (FloatColor){
    .red = (color1.red + color2.red) * 0.5f,
    .green = (color1.green + color2.green) * 0.5f,
    .blue = (color1.blue + color2.blue) * 0.5f,
    .alpha = (color1.alpha + color2.alpha) * 0.5f};
}

FloatColor floatColorLerp(FloatColor color1, FloatColor color2, float theta)
{
  return (FloatColor){
    .red = lerpFloat(color1.red, color2.red, theta),
    .green = lerpFloat(color1.green, color2.green, theta),
    .blue = lerpFloat(color1.blue, color2.blue, theta),
    .alpha = lerpFloat(color1.alpha, color2.alpha, theta)};
}

bool floatColorEquals(FloatColor color1, FloatColor color2)
{
  return
    nearlyEqualFloats(color1.red, color2.red) &&
    nearlyEqualFloats(color1.green, color2.green) &&
    nearlyEqualFloats(color1.blue, color2.blue) &&
    nearlyEqualFloats(color1.alpha, color2.alpha);
}

bool floatColorEmpty(FloatColor color)
{
  return
    nearlyZeroFloat(color.red) &&
    nearlyZeroFloat(color.green) &&
    nearlyZeroFloat(color.blue) &&
    nearlyZeroFloat(color.alpha);
}

float floatColorGray(FloatColor color)
{
  return color.red * 0.29889531f + color.green * 0.58662247f + color.blue * 0.11448222f;
}

Color floatToColor(FloatColor color)
{
  return
    (Color)roundf(color.red * 255.0f) |
    ((Color)roundf(color.green * 255.0f) << 8) |
    ((Color)roundf(color.blue * 255.0f) << 16) |
    ((Color)roundf(color.alpha * 255.0f) << 24);
}

FloatColorRGB floatColorToRGB(FloatColor color)
{
  return (FloatColorRGB){.red = color.red, .green = color.green, .blue = color.blue};
}

// Point constants.

Point const pointZero = {.x = 0, .y = 0};
Point const pointUnity = {.x = 1, .y = 1};
Point const pointAxisX = {.x = 1, .y = 0};
Point const pointAxisY = {.x = 0, .y = 1};

// Point functions.

Point makePoint(int32_t x, int32_t y)
{
  return (Point){.x = x, .y = y};
}

Point pointValue(int32_t value)
{
  return (Point){.x = value, .y = value};
}

Point pointNegate(Point point)
{
  return (Point){.x = -point.x, .y = -point.y};
}

Point pointAdd(Point point1, Point point2)
{
  return (Point){.x = point1.x + point2.x, .y = point1.y + point2.y};
}

Point pointSubtract(Point point1, Point point2)
{
  return (Point){.x = point1.x - point2.x, .y = point1.y - point2.y};
}

Point pointMultiply(Point point1, Point point2)
{
  return (Point){.x = point1.x * point2.x, .y = point1.y * point2.y};
}

Point pointDivide(Point point1, Point point2)
{
  return (Point){.x = point1.x / point2.x, .y = point1.y / point2.y};
}

Point pointRescale(Point point, int32_t theta)
{
  return (Point){.x = point.x * theta, .y = point.y * theta};
}

float pointLength(Point point)
{
  return hypotf((float)point.x, (float)point.y);
}

float pointDistance(Point point1, Point point2)
{
  return pointLength(pointSubtract(point1, point2));
}

float pointAngle(Point point)
{
  return atan2f(-(float)point.y, (float)point.x);
}

int32_t pointDot(Point point1, Point point2)
{
  return (point1.x * point2.x) + (point1.y * point2.y);
}

int32_t pointCross(Point point1, Point point2)
{
  return (point1.x * point2.y) - (point1.y * point2.x);
}

Point pointLerp(Point point1, Point point2, float theta)
{
  return (Point){
    .x = (int32_t)roundf(lerpFloat((float)point1.x, (float)point2.x, theta)),
    .y = (int32_t)roundf(lerpFloat((float)point1.y, (float)point2.y, theta))};
}

bool pointEquals(Point point1, Point point2)
{
  return point1.x == point2.x && point1.y == point2.y;
}

bool pointEmpty(Point point)
{
  return point.x == 0 && point.y == 0;
}

// PointF constants.

PointF const pointZeroF = {.x = 0.0f, .y = 0.0f};
PointF const pointUnityF = {.x = 1.0f, .y = 1.0f};
PointF const pointAxisXF = {.x = 1.0f, .y = 0.0f};
PointF const pointAxisYF = {.x = 0.0f, .y = 1.0f};

// PointF functions.

PointF makePointF(float x, float y)
{
  return (PointF){.x = x, .y = y};
}

PointF pointF(Point point)
{
  return (PointF){.x = (float)point.x, .y = (float)point.y};
}

PointF pointValueF(float value)
{
  return (PointF){.x = value, .y = value};
}

PointF pointNegateF(PointF point)
{
  return (PointF){.x = -point.x, .y = -point.y};
}

PointF pointAddF(PointF point1, PointF point2)
{
  return (PointF){.x = point1.x + point2.x, .y = point1.y + point2.y};
}

PointF pointSubtractF(PointF point1, PointF point2)
{
  return (PointF){.x = point1.x - point2.x, .y = point1.y - point2.y};
}

PointF pointMultiplyF(PointF point1, PointF point2)
{
  return (PointF){.x = point1.x * point2.x, .y = point1.y * point2.y};
}

PointF pointDivideF(PointF point1, PointF point2)
{
  return (PointF){.x = point1.x / point2.x, .y = point1.y / point2.y};
}

PointF pointRescaleF(PointF point, float theta)
{
  return (PointF){.x = point.x * theta, .y = point.y * theta};
}

float pointLengthF(PointF point)
{
  return hypotf(point.x, point.y);
}

float pointDistanceF(PointF point1, PointF point2)
{
  return pointLengthF(pointSubtractF(point1, point2));
}

float pointAngleF(PointF point)
{
  return atan2f(-point.y, point.x);
}

float pointDotF(PointF point1, PointF point2)
{
  return (point1.x * point2.x) + (point1.y * point2.y);
}

float pointCrossF(PointF point1, PointF point2)
{
  return (point1.x * point2.y) - (point1.y * point2.x);
}

PointF pointNormalizeF(PointF point)
{
  return pointRescaleF(point, 1.0f / pointLengthF(point));
}

PointF pointLerpF(PointF point1, PointF point2, float theta)
{
  return (PointF){
    .x = lerpFloat(point1.x, point2.x, theta),
    .y = lerpFloat(point1.y, point2.y, theta)};
}

bool pointEqualsF(PointF point1, PointF point2)
{
  return nearlyEqualFloats(point1.x, point2.x) && nearlyEqualFloats(point1.y, point2.y);
}

bool pointEmptyF(PointF point)
{
  return nearlyZeroFloat(point.x) && nearlyZeroFloat(point.y);
}

PointF pointTransformF(PointF point, Matrix3x2 matrix)
{
  return (PointF){
    .x = (point.x * matrix.data[0][0]) + (point.y * matrix.data[1][0]) + matrix.data[2][0],
    .y = (point.x * matrix.data[0][1]) + (point.y * matrix.data[1][1]) + matrix.data[2][1]};
}

// Vector constants.

Vector const vectorZero = {.x = 0.0f, .y = 0.0f, .z = 0.0f};
Vector const vectorUnity = {.x = 1.0f, .y = 1.0f, .z = 1.0f};
Vector const vectorAxisX = {.x = 1.0f, .y = 0.0f, .z = 0.0f};
Vector const vectorAxisY = {.x = 0.0f, .y = 1.0f, .z = 0.0f};
Vector const vectorAxisZ = {.x = 0.0f, .y = 0.0f, .z = 1.0f};

// Vector functions.

Vector makeVector(float x, float y, float z)
{
  return (Vector){.x = x, .y = y, .z = z};
}

Vector vectorValue(float value)
{
  return (Vector){.x = value, .y = value, .z = value};
}

Vector vectorNegate(Vector vector)
{
  return (Vector){.x = -vector.x, .y = -vector.y, .z = -vector.z};
}

Vector vectorAdd(Vector vector1, Vector vector2)
{
  return (Vector){.x = vector1.x + vector2.x, .y = vector1.y + vector2.y, .z = vector1.z + vector2.z};
}

Vector vectorSubtract(Vector vector1, Vector vector2)
{
  return (Vector){.x = vector1.x - vector2.x, .y = vector1.y - vector2.y, .z = vector1.z - vector2.z};
}

Vector vectorMultiply(Vector vector1, Vector vector2)
{
  return (Vector){.x = vector1.x * vector2.x, .y = vector1.y * vector2.y, .z = vector1.z * vector2.z};
}

Vector vectorDivide(Vector vector1, Vector vector2)
{
  return (Vector){.x = vector1.x / vector2.x, .y = vector1.y / vector2.y, .z = vector1.z / vector2.z};
}

Vector vectorRescale(Vector vector, float theta)
{
  return (Vector){.x = vector.x * theta, .y = vector.y * theta, .z = vector.z * theta};
}

float vectorDot(Vector vector1, Vector vector2)
{
  return (vector1.x * vector2.x) + (vector1.y * vector2.y) + (vector1.z * vector2.z);
}

Vector vectorCross(Vector vector1, Vector vector2)
{
  return (Vector){
    .x = (vector1.y * vector2.z) - (vector1.z * vector2.y),
    .y = (vector1.z * vector2.x) - (vector1.x * vector2.z),
    .z = (vector1.x * vector2.y) - (vector1.y * vector2.x)};
}

float vectorLength(Vector vector)
{
  return sqrtf(vectorDot(vector, vector));
}

float vectorDistance(Vector vector1, Vector vector2)
{
  return vectorLength(vectorSubtract(vector1, vector2));
}

float vectorAngle(Vector vector1, Vector vector2)
{
  float const cosValue = vectorDot(vector1, vector2) / (vectorLength(vector1) * vectorLength(vector2));
  return acosf(saturateFloat(cosValue, -1.0f, 1.0f));
}

Vector vectorNormalize(Vector vector)
{
  float const magnitude = vectorLength(vector);
  return magnitude > vectorEpsilon ? vectorRescale(vector, 1.0f / magnitude) : vector;
}

Vector vectorParallel(Vector vector, Vector direction)
{
  float const dirAmp = vectorLength(direction);
  return vectorRescale(direction, (vectorDot(vector, direction) / (dirAmp * dirAmp)));
}

Vector vectorPerpendicular(Vector vector, Vector direction)
{
  return vectorSubtract(vector, vectorParallel(vector, direction));
}

Vector vectorReflect(Vector vector, Vector normal)
{
  return vectorSubtract(vector, vectorRescale(normal, vectorDot(vector, normal) * 2.0f));
}

Vector vectorLerp(Vector vector1, Vector vector2, float theta)
{
  return (Vector){
    .x = lerpFloat(vector1.x, vector2.x, theta),
    .y = lerpFloat(vector1.y, vector2.y, theta),
    .z = lerpFloat(vector1.z, vector2.z, theta)};
}

bool vectorEquals(Vector vector1, Vector vector2)
{
  return
    nearlyEqualFloats(vector1.x, vector2.x) &&
    nearlyEqualFloats(vector1.y, vector2.y) &&
    nearlyEqualFloats(vector1.z, vector2.z);
}

bool vectorEmpty(Vector vector)
{
  return nearlyZeroFloat(vector.x) && nearlyZeroFloat(vector.y) && nearlyZeroFloat(vector.z);
}

Vector vectorTransform(Vector vector, float w, Matrix matrix)
{
  return (Vector){
    (vector.x * matrix.data[0][0]) + (vector.y * matrix.data[1][0]) + (vector.z * matrix.data[2][0]) +
      w * matrix.data[3][0],
    (vector.x * matrix.data[0][1]) + (vector.y * matrix.data[1][1]) + (vector.z * matrix.data[2][1]) +
      w * matrix.data[3][1],
    (vector.x * matrix.data[0][2]) + (vector.y * matrix.data[1][2]) + (vector.z * matrix.data[2][2]) +
      w * matrix.data[3][2]};
}

Vector vectorProject(Vector vector, float w, Matrix matrix)
{
  float const wf = vector.x * matrix.data[0][3] + vector.y * matrix.data[1][3] +
    vector.z * matrix.data[2][3] + w * matrix.data[3][3];

  return vectorRescale(vectorTransform(vector, w, matrix), 1.0f / wf);
}

PointF vectorProjectTarget(Vector vector, PointF targetSize, Matrix matrix)
{
  float const wf = vector.x * matrix.data[0][3] + vector.y * matrix.data[1][3] +
    vector.z * matrix.data[2][3] + matrix.data[3][3];

  PointF temp = (PointF){
    .x = vector.x * matrix.data[0][0] + vector.y * matrix.data[1][0] + vector.z * matrix.data[2][0] +
      matrix.data[3][0],
    .y = vector.x * matrix.data[0][1] + vector.y * matrix.data[1][1] + vector.z * matrix.data[2][1] +
      matrix.data[3][1]};

  temp = pointRescaleF(temp, 1.0f / wf);

  return (PointF){
    .x = (temp.x * 0.5f + 0.5f) * targetSize.x,
    .y = (0.5f - temp.y * 0.5f) * targetSize.y};
}

// Vector4 constants.

Vector4 const vectorZero3D = {.x = 0.0f, .y = 0.0f, .z = 0.0f, .w = 0.0f};
Vector4 const vectorUnity3D = {.x = 1.0f, .y = 1.0f, .z = 1.0f, .w = 1.0f};
Vector4 const vectorAxisX3D = {.x = 1.0f, .y = 0.0f, .z = 0.0f, .w = 0.0f};
Vector4 const vectorAxisY3D = {.x = 0.0f, .y = 1.0f, .z = 0.0f, .w = 0.0f};
Vector4 const vectorAxisZ3D = {.x = 0.0f, .y = 0.0f, .z = 1.0f, .w = 0.0f};
Vector4 const vectorAxisW3D = {.x = 0.0f, .y = 0.0f, .z = 0.0f, .w = 1.0f};

// Vector4 functions.

Vector4 makeVector3D(float x, float y, float z, float w)
{
  return (Vector4){.x = x, .y = y, .z = z, .w = w};
}

Vector4 vectorFromValue3D(float value)
{
  return (Vector4){.x = value, .y = value, .z = value, .w = value};
}

Vector4 vectorNegate3D(Vector4 vector)
{
  return (Vector4){.x = -vector.x, .y = -vector.y, .z = -vector.z, .w = -vector.w};
}

Vector4 vectorAdd3D(Vector4 vector1, Vector4 vector2)
{
  return (Vector4){.x = vector1.x + vector2.x, .y = vector1.y + vector2.y, .z = vector1.z + vector2.z,
    .w = vector1.w + vector2.w};
}

Vector4 vectorSubtract3D(Vector4 vector1, Vector4 vector2)
{
  return (Vector4){.x = vector1.x - vector2.x, .y = vector1.y - vector2.y, .z = vector1.z - vector2.z,
    .w = vector1.w - vector2.w};
}

Vector4 vectorMultiply3D(Vector4 vector1, Vector4 vector2)
{
  return (Vector4){.x = vector1.x * vector2.x, .y = vector1.y * vector2.y, .z = vector1.z * vector2.z,
    .w = vector1.w * vector2.w};
}

Vector4 vectorDivide3D(Vector4 vector1, Vector4 vector2)
{
  return (Vector4){.x = vector1.x / vector2.x, .y = vector1.y / vector2.y, .z = vector1.z / vector2.z,
    .w = vector1.w / vector2.w};
}

Vector4 vectorRescale3D(Vector4 vector, float theta)
{
  return (Vector4){.x = vector.x * theta, .y = vector.y * theta, .z = vector.z * theta,
    .w = vector.w * theta};
}

float vectorDot3D(Vector4 vector1, Vector4 vector2)
{
  return (vector1.x * vector2.x) + (vector1.y * vector2.y) + (vector1.z * vector2.z) +
    (vector1.w * vector2.w);
}

float vectorLength3D(Vector4 vector)
{
  return sqrtf(vectorDot3D(vector, vector));
}

float vectorDistance3D(Vector4 vector1, Vector4 vector2)
{
  return vectorLength3D(vectorSubtract3D(vector1, vector2));
}

Vector4 vectorNormalize3D(Vector4 vector)
{
  return vectorMultiply3D(vector, vectorFromValue3D(1.0f / vectorLength3D(vector)));
}

Vector4 vectorLerp3D(Vector4 vector1, Vector4 vector2, float theta)
{
  return (Vector4){
    .x = lerpFloat(vector1.x, vector2.x, theta),
    .y = lerpFloat(vector1.y, vector2.y, theta),
    .z = lerpFloat(vector1.z, vector2.z, theta),
    .w = lerpFloat(vector1.w, vector2.w, theta)};
}

bool vectorEquals3D(Vector4 vector1, Vector4 vector2)
{
  return
    nearlyEqualFloats(vector1.x, vector2.x) &&
    nearlyEqualFloats(vector1.y, vector2.y) &&
    nearlyEqualFloats(vector1.z, vector2.z) &&
    nearlyEqualFloats(vector1.w, vector2.w);
}

bool vectorEmpty3D(Vector4 vector)
{
  return nearlyZeroFloat(vector.x) && nearlyZeroFloat(vector.y) && nearlyZeroFloat(vector.z) &&
    nearlyZeroFloat(vector.w);
}

Vector4 vectorTransform3D(Vector4 vector, Matrix matrix)
{
  return (Vector4){
    .x = vector.x * matrix.data[0][0] + vector.y * matrix.data[1][0] +
      vector.z * matrix.data[2][0] + vector.w * matrix.data[3][0],
    .y = vector.x * matrix.data[0][1] + vector.y * matrix.data[1][1] +
      vector.z * matrix.data[2][1] + vector.w * matrix.data[3][1],
    .z = vector.x * matrix.data[0][2] + vector.y * matrix.data[1][2] +
      vector.z * matrix.data[2][2] + vector.w * matrix.data[3][2],
    .w = vector.x * matrix.data[0][3] + vector.y * matrix.data[1][3] +
      vector.z * matrix.data[2][3] + vector.w * matrix.data[3][3]};
}

// Matrix3 constants.

Matrix3x2 const matrixZero3x2 = {{{0.0f, 0.0f}, {0.0f, 0.0f}, {0.0f, 0.0f}}};

Matrix3x2 const matrixIdentity3x2 = {{{1.0f, 0.0f}, {0.0f, 1.0f}, {0.0f, 0.0f}}};

// Matrix3 functions.

Matrix3x2 matrixMultiply3x2(Matrix3x2 matrix1, Matrix3x2 matrix2)
{
  Matrix3x2 res;
  res.data[0][0] = matrix1.data[0][0] * matrix2.data[0][0] + matrix1.data[0][1] * matrix2.data[1][0];
  res.data[0][1] = matrix1.data[0][0] * matrix2.data[0][1] + matrix1.data[0][1] * matrix2.data[1][1];
  res.data[1][0] = matrix1.data[1][0] * matrix2.data[0][0] + matrix1.data[1][1] * matrix2.data[1][0];
  res.data[1][1] = matrix1.data[1][0] * matrix2.data[0][1] + matrix1.data[1][1] * matrix2.data[1][1];
  res.data[2][0] = matrix1.data[2][0] * matrix2.data[0][0] + matrix1.data[2][1] * matrix2.data[1][0] +
    matrix2.data[2][0];
  res.data[2][1] = matrix1.data[2][0] * matrix2.data[0][1] + matrix1.data[2][1] * matrix2.data[1][1] +
    matrix2.data[2][1];
  return res;
}

Matrix3x2 matrixRescale3x2(Matrix3x2 matrix, float theta)
{
  Matrix3x2 res;

  for (int_fast8_t j = 0; j < 3; ++j)
    for (int_fast8_t i = 0; i < 2; ++i)
      res.data[j][i] = matrix.data[j][i] * theta;

  return res;
}

float matrixDeterminant3x2(Matrix3x2 matrix)
{
  return (matrix.data[0][0] * matrix.data[1][1]) - (matrix.data[1][0] * matrix.data[0][1]);
}

Matrix3x2 matrixInverse3x2(Matrix3x2 matrix)
{
  float invDet = 1.0f / matrixDeterminant3x2(matrix);
  Matrix3x2 res;

  res.data[0][0] = matrix.data[1][1] * invDet;
  res.data[0][1] = -matrix.data[0][1] * invDet;
  res.data[1][0] = -matrix.data[1][0] * invDet;
  res.data[1][1] = matrix.data[0][0] * invDet;
  res.data[2][0] = (matrix.data[1][0] * matrix.data[2][1] - matrix.data[2][0] * matrix.data[1][1]) * invDet;
  res.data[2][1] = (matrix.data[2][0] * matrix.data[0][1] - matrix.data[0][0] * matrix.data[2][1]) * invDet;

  return res;
}

Matrix3x2 matrixInverseTranspose3x2(Matrix3x2 matrix)
{
  float invDet = 1.0f / matrixDeterminant3x2(matrix);
  Matrix3x2 res = {0};

  res.data[0][0] = matrix.data[1][1] * invDet;
  res.data[1][0] = -matrix.data[0][1] * invDet;
  res.data[0][1] = -matrix.data[1][0] * invDet;
  res.data[1][1] = matrix.data[0][0] * invDet;

  return res;
}

Matrix3x2 matrixScale3x2(PointF scale)
{
  Matrix3x2 res = matrixIdentity3x2;

  res.data[0][0] = scale.x;
  res.data[1][1] = scale.y;

  return res;
}

Matrix3x2 matrixScaleBy3x2(PointF center, PointF scale)
{
  Matrix3x2 res = matrixIdentity3x2;

  res.data[0][0] = scale.x;
  res.data[1][1] = scale.y;
  res.data[2][0] = center.x * (1.0f - scale.x);
  res.data[2][1] = center.y * (1.0f - scale.y);

  return res;
}

Matrix3x2 matrixTranslate3x2(PointF offset)
{
  Matrix3x2 res = matrixIdentity3x2;

  res.data[2][0] = offset.x;
  res.data[2][1] = offset.y;

  return res;
}

Matrix3x2 matrixRotate3x2(float angle)
{
  Matrix3x2 res = matrixIdentity3x2;

  float const angleSin = sinf(angle);
  float const angleCos = cosf(angle);

  res.data[0][0] = angleCos;
  res.data[0][1] = angleSin;
  res.data[1][0] = -angleSin;
  res.data[1][1] = angleCos;

  return res;
}

Matrix3x2 matrixRotateBy3x2(PointF center, float angle)
{
  float const angleSin = sinf(angle);
  float const angleCos = cosf(angle);

  Matrix3x2 res;

  res.data[0][0] = angleCos;
  res.data[0][1] = angleSin;
  res.data[1][0] = -angleSin;
  res.data[1][1] = angleCos;
  res.data[2][0] = center.x * (1.0f - angleCos) + center.y * angleSin;
  res.data[2][1] = center.y * (1.0f - angleCos) - center.x * angleSin;

  return res;
}

bool matrixEquals3x2(Matrix3x2 matrix1, Matrix3x2 matrix2)
{
  for (int_fast8_t j = 0; j < 3; j++)
    for (int_fast8_t i = 0; i < 3; i++)
      if (!nearlyEqualFloats(matrix1.data[j][i], matrix2.data[j][i]))
        return false;

  return true;
}

bool matrixEmpty3x2(Matrix3x2 matrix)
{
  for (int_fast8_t j = 0; j < 3; j++)
    for (int_fast8_t i = 0; i < 3; i++)
      if (!nearlyZeroFloat(matrix.data[j][i]))
        return false;

  return true;
}

// Matrix constants.

Matrix const matrixZero = {{
  {0.0f, 0.0f, 0.0f, 0.0f},
  {0.0f, 0.0f, 0.0f, 0.0f},
  {0.0f, 0.0f, 0.0f, 0.0f},
  {0.0f, 0.0f, 0.0f, 0.0f}}};

Matrix const matrixIdentity = {{
  {1.0f, 0.0f, 0.0f, 0.0f},
  {0.0f, 1.0f, 0.0f, 0.0f},
  {0.0f, 0.0f, 1.0f, 0.0f},
  {0.0f, 0.0f, 0.0f, 1.0f}}};

// Matrix forward declarations.

float detSub(float a1, float a2, float a3, float b1, float b2, float b3, float c1, float c2, float c3);

// Matrix functions.

Matrix3x2 matrixGet3x2(Matrix matrix)
{
  Matrix3x2 temp;

  for (int_fast8_t j = 0; j < 3; ++j)
    for (int_fast8_t i = 0; i < 2; ++i)
      temp.data[j][i] = matrix.data[j][i];

  return temp;
}

Matrix matrixMultiply(Matrix matrix1, Matrix matrix2)
{
  Matrix temp;

  for (int_fast8_t j = 0; j < 4; ++j)
    for (int_fast8_t i = 0; i < 4; ++i)
      temp.data[j][i] =
        matrix1.data[j][0] * matrix2.data[0][i] +
        matrix1.data[j][1] * matrix2.data[1][i] +
        matrix1.data[j][2] * matrix2.data[2][i] +
        matrix1.data[j][3] * matrix2.data[3][i];

  return temp;
}

Matrix matrixRescale(Matrix matrix, float theta)
{
  Matrix temp;

  for (int_fast8_t j = 0; j < 4; ++j)
    for (int_fast8_t i = 0; i < 4; ++i)
      temp.data[j][i] = matrix.data[j][i] * theta;

  return temp;
}

float matrixDeterminant(Matrix matrix)
{
  return
    matrix.data[0][0] * detSub(
      matrix.data[1][1], matrix.data[2][1], matrix.data[3][1],
      matrix.data[1][2], matrix.data[2][2], matrix.data[3][2],
      matrix.data[1][3], matrix.data[2][3], matrix.data[3][3]) -
    matrix.data[0][1] * detSub(
      matrix.data[1][0], matrix.data[2][0], matrix.data[3][0],
      matrix.data[1][2], matrix.data[2][2], matrix.data[3][2],
      matrix.data[1][3], matrix.data[2][3], matrix.data[3][3]) +
    matrix.data[0][2] * detSub(
      matrix.data[1][0], matrix.data[2][0], matrix.data[3][0],
      matrix.data[1][1], matrix.data[2][1], matrix.data[3][1],
      matrix.data[1][3], matrix.data[2][3], matrix.data[3][3]) -
    matrix.data[0][3] * detSub(
      matrix.data[1][0], matrix.data[2][0], matrix.data[3][0],
      matrix.data[1][1], matrix.data[2][1], matrix.data[3][1],
      matrix.data[1][2], matrix.data[2][2], matrix.data[3][2]);
}

Matrix matrixTranspose(Matrix matrix)
{
  Matrix temp;

  for (int_fast8_t j = 0; j < 4; ++j)
    for (int_fast8_t i = 0; i < 4; ++i)
      temp.data[i][j] = matrix.data[j][i];

  return temp;
}

Matrix matrixAdjoint(Matrix matrix)
{
  Matrix temp;

  temp.data[0][0] = detSub(
    matrix.data[1][1], matrix.data[2][1], matrix.data[3][1],
    matrix.data[1][2], matrix.data[2][2], matrix.data[3][2],
    matrix.data[1][3], matrix.data[2][3], matrix.data[3][3]);

  temp.data[1][0] = -detSub(
    matrix.data[1][0], matrix.data[2][0], matrix.data[3][0],
    matrix.data[1][2], matrix.data[2][2], matrix.data[3][2],
    matrix.data[1][3], matrix.data[2][3], matrix.data[3][3]);

  temp.data[2][0] = detSub(
    matrix.data[1][0], matrix.data[2][0], matrix.data[3][0],
    matrix.data[1][1], matrix.data[2][1], matrix.data[3][1],
    matrix.data[1][3], matrix.data[2][3], matrix.data[3][3]);

  temp.data[3][0] = -detSub(
    matrix.data[1][0], matrix.data[2][0], matrix.data[3][0],
    matrix.data[1][1], matrix.data[2][1], matrix.data[3][1],
    matrix.data[1][2], matrix.data[2][2], matrix.data[3][2]);

  temp.data[0][1] = -detSub(
    matrix.data[0][1], matrix.data[2][1], matrix.data[3][1],
    matrix.data[0][2], matrix.data[2][2], matrix.data[3][2],
    matrix.data[0][3], matrix.data[2][3], matrix.data[3][3]);

  temp.data[1][1] = detSub(
    matrix.data[0][0], matrix.data[2][0], matrix.data[3][0],
    matrix.data[0][2], matrix.data[2][2], matrix.data[3][2],
    matrix.data[0][3], matrix.data[2][3], matrix.data[3][3]);

  temp.data[2][1] = -detSub(
    matrix.data[0][0], matrix.data[2][0], matrix.data[3][0],
    matrix.data[0][1], matrix.data[2][1], matrix.data[3][1],
    matrix.data[0][3], matrix.data[2][3], matrix.data[3][3]);

  temp.data[3][1] = detSub(
    matrix.data[0][0], matrix.data[2][0], matrix.data[3][0],
    matrix.data[0][1], matrix.data[2][1], matrix.data[3][1],
    matrix.data[0][2], matrix.data[2][2], matrix.data[3][2]);

  temp.data[0][2] = detSub(
    matrix.data[0][1], matrix.data[1][1], matrix.data[3][1],
    matrix.data[0][2], matrix.data[1][2], matrix.data[3][2],
    matrix.data[0][3], matrix.data[1][3], matrix.data[3][3]);

  temp.data[1][2] = -detSub(
    matrix.data[0][0], matrix.data[1][0], matrix.data[3][0],
    matrix.data[0][2], matrix.data[1][2], matrix.data[3][2],
    matrix.data[0][3], matrix.data[1][3], matrix.data[3][3]);

  temp.data[2][2] = detSub(
    matrix.data[0][0], matrix.data[1][0], matrix.data[3][0],
    matrix.data[0][1], matrix.data[1][1], matrix.data[3][1],
    matrix.data[0][3], matrix.data[1][3], matrix.data[3][3]);

  temp.data[3][2] = -detSub(
    matrix.data[0][0], matrix.data[1][0], matrix.data[3][0],
    matrix.data[0][1], matrix.data[1][1], matrix.data[3][1],
    matrix.data[0][2], matrix.data[1][2], matrix.data[3][2]);

  temp.data[0][3] = -detSub(
    matrix.data[0][1], matrix.data[1][1], matrix.data[2][1],
    matrix.data[0][2], matrix.data[1][2], matrix.data[2][2],
    matrix.data[0][3], matrix.data[1][3], matrix.data[2][3]);

  temp.data[1][3] = detSub(
    matrix.data[0][0], matrix.data[1][0], matrix.data[2][0],
    matrix.data[0][2], matrix.data[1][2], matrix.data[2][2],
    matrix.data[0][3], matrix.data[1][3], matrix.data[2][3]);

  temp.data[2][3] = -detSub(
    matrix.data[0][0], matrix.data[1][0], matrix.data[2][0],
    matrix.data[0][1], matrix.data[1][1], matrix.data[2][1],
    matrix.data[0][3], matrix.data[1][3], matrix.data[2][3]);

  temp.data[3][3] = detSub(
    matrix.data[0][0], matrix.data[1][0], matrix.data[2][0],
    matrix.data[0][1], matrix.data[1][1], matrix.data[2][1],
    matrix.data[0][2], matrix.data[1][2], matrix.data[2][2]);

  return temp;
}

Matrix matrixInverse(Matrix matrix)
{
  return matrixRescale(matrixAdjoint(matrix), 1.0f / matrixDeterminant(matrix));
}

Matrix matrixSub3x3(Matrix matrix)
{
  Matrix temp = matrixIdentity;

  for (int_fast32_t j = 0; j < 3; ++j)
    for (int_fast32_t i = 0; i < 3; ++i)
      temp.data[j][i] = matrix.data[j][i];

  return temp;
}

Vector matrixEyePosition(Matrix matrix)
{
  return (Vector){
    .x = -matrix.data[0][0] * matrix.data[3][0] - matrix.data[0][1] * matrix.data[3][1] -
      matrix.data[0][2] * matrix.data[3][2],
    .y = -matrix.data[1][0] * matrix.data[3][0] - matrix.data[1][1] * matrix.data[3][1] -
      matrix.data[1][2] * matrix.data[3][2],
    .z = -matrix.data[2][0] * matrix.data[3][0] - matrix.data[2][1] * matrix.data[3][1] -
      matrix.data[2][2] * matrix.data[3][2]};
}

Vector matrixWorldPosition(Matrix matrix)
{
  return (Vector){.x = matrix.data[3][0], .y = matrix.data[3][1], .z = matrix.data[3][2]};
}

Matrix matrixTranslate(Vector offset)
{
  Matrix temp = matrixIdentity;

  temp.data[3][0] = offset.x;
  temp.data[3][1] = offset.y;
  temp.data[3][2] = offset.z;

  return temp;
}

Matrix matrixScale(Vector scale)
{
  Matrix temp = matrixIdentity;

  temp.data[0][0] = scale.x;
  temp.data[1][1] = scale.y;
  temp.data[2][2] = scale.z;

  return temp;
}

Matrix matrixRotateX(float angle)
{
  Matrix temp = matrixIdentity;

  float const sinAngle = sinf(angle);
  float const cosAngle = cosf(angle);

  temp.data[1][1] = cosAngle;
  temp.data[1][2] = sinAngle;
  temp.data[2][1] = -sinAngle;
  temp.data[2][2] = cosAngle;

  return temp;
}

Matrix matrixRotateY(float angle)
{
  Matrix temp = matrixIdentity;

  float const sinAngle = sinf(angle);
  float const cosAngle = cosf(angle);

  temp.data[0][0] = cosAngle;
  temp.data[0][2] = -sinAngle;
  temp.data[2][0] = sinAngle;
  temp.data[2][2] = cosAngle;

  return temp;
}

Matrix matrixRotateZ(float angle)
{
  Matrix temp = matrixIdentity;

  float const sinAngle = sinf(angle);
  float const cosAngle = cosf(angle);

  temp.data[0][0] = cosAngle;
  temp.data[0][1] = sinAngle;
  temp.data[1][0] = -sinAngle;
  temp.data[1][1] = cosAngle;

  return temp;
}

Matrix matrixRotate(Vector axis, float angle)
{
  float const sinAngle = sinf(angle);
  float const cosAngle = cosf(angle);
  float const invCosAngle = 1.0f - cosAngle;

  float const xyMul = axis.x * axis.y * invCosAngle;
  float const xzMul = axis.x * axis.z * invCosAngle;
  float const yzMul = axis.y * axis.z * invCosAngle;

  float const xSin = axis.x * sinAngle;
  float const ySin = axis.y * sinAngle;
  float const zSin = axis.z * sinAngle;

  Matrix temp = matrixIdentity;

  temp.data[0][0] = axis.x * axis.x * invCosAngle + cosAngle;
  temp.data[0][1] = xyMul + zSin;
  temp.data[0][2] = xzMul - ySin;
  temp.data[1][0] = xyMul - zSin;
  temp.data[1][1] = axis.y * axis.y * invCosAngle + cosAngle;
  temp.data[1][2] = yzMul + xSin;
  temp.data[2][0] = xzMul + ySin;
  temp.data[2][1] = yzMul - xSin;
  temp.data[2][2] = axis.z * axis.z * invCosAngle + cosAngle;

  return temp;
}

Matrix matrixHeadingPitchBank(Vector angles)
{
  float const sinHeading = sinf(angles.y);
  float const cosHeading = cosf(angles.y);
  float const sinPitch = sinf(angles.x);
  float const cosPitch = cosf(angles.x);
  float const sinBank = sinf(angles.z);
  float const cosBank = cosf(angles.z);

  Matrix temp = matrixIdentity;

  temp.data[0][0] = cosHeading * cosBank + sinHeading * sinPitch * sinBank;
  temp.data[0][1] = sinHeading * sinPitch * cosBank - cosHeading * sinBank;
  temp.data[0][2] = sinHeading * cosPitch;
  temp.data[1][0] = sinBank * cosPitch;
  temp.data[1][1] = cosBank * cosPitch;
  temp.data[1][2] = -sinPitch;
  temp.data[2][0] = cosHeading * sinPitch * sinBank - sinHeading * cosBank;
  temp.data[2][1] = sinBank * sinHeading + cosHeading * sinPitch * cosBank;
  temp.data[2][2] = cosHeading * cosPitch;

  return temp;
}

Matrix matrixYawPitchRoll(Vector angles)
{
  float const sinYaw = sinf(angles.y);
  float const cosYaw = cosf(angles.y);
  float const sinPitch = sinf(angles.x);
  float const cosPitch = cosf(angles.x);
  float const sinRoll = sinf(angles.z);
  float const cosRoll = cosf(angles.z);

  Matrix temp = matrixIdentity;

  temp.data[0][0] = cosRoll * cosYaw + sinPitch * sinRoll * sinYaw;
  temp.data[0][1] = cosYaw * sinPitch * sinRoll - cosRoll * sinYaw;
  temp.data[0][2] = -cosPitch * sinRoll;
  temp.data[1][0] = cosPitch * sinYaw;
  temp.data[1][1] = cosPitch * cosYaw;
  temp.data[1][2] = sinPitch;
  temp.data[2][0] = cosYaw * sinRoll - cosRoll * sinPitch * sinYaw;
  temp.data[2][1] = -cosRoll * cosYaw * sinPitch - sinRoll * sinYaw;
  temp.data[2][2] = cosPitch * cosRoll;

  return temp;
}

Matrix matrixReflect(Vector axis)
{
  float const xyMul = -2.0f * axis.x * axis.y;
  float const xzMul = -2.0f * axis.x * axis.z;
  float const yzMul = -2.0f * axis.y * axis.z;

  Matrix temp = matrixIdentity;

  temp.data[0][0] = 1.0f - (2.0f * axis.x * axis.x);
  temp.data[0][1] = xyMul;
  temp.data[0][2] = xzMul;
  temp.data[1][0] = xyMul;
  temp.data[1][1] = 1.0f - (2.0f * axis.y * axis.y);
  temp.data[1][2] = yzMul;
  temp.data[2][0] = xzMul;
  temp.data[2][1] = yzMul;
  temp.data[2][2] = 1.0f - (2.0f * axis.z * axis.z);

  return temp;
}

Matrix matrixLookAt(Vector origin, Vector target, Vector ceiling)
{
  Vector const zAxis = vectorNormalize(vectorSubtract(target, origin));
  Vector const xAxis = vectorNormalize(vectorCross(ceiling, zAxis));
  Vector const yAxis = vectorNormalize(vectorCross(zAxis, xAxis));

  Matrix temp;

  temp.data[0][0] = xAxis.x;
  temp.data[0][1] = yAxis.x;
  temp.data[0][2] = zAxis.x;
  temp.data[0][3] = 0.0f;
  temp.data[1][0] = xAxis.y;
  temp.data[1][1] = yAxis.y;
  temp.data[1][2] = zAxis.y;
  temp.data[1][3] = 0.0f;
  temp.data[2][0] = xAxis.z;
  temp.data[2][1] = yAxis.z;
  temp.data[2][2] = zAxis.z;
  temp.data[2][3] = 0.0f;
  temp.data[3][0] = -vectorDot(xAxis, origin);
  temp.data[3][1] = -vectorDot(yAxis, origin);
  temp.data[3][2] = -vectorDot(zAxis, origin);
  temp.data[3][3] = 1.0f;

  return temp;
}

Matrix matrixPerspectiveFOVY(float fieldOfView, float aspectRatio, float minRange, float maxRange,
  bool negativeDepthRange)
{
  float const yScale = 1.0f / tanf(fieldOfView * 0.5f);
  float const xScale = yScale / aspectRatio;

  Matrix temp = matrixZero;

  temp.data[0][0] = xScale;
  temp.data[1][1] = yScale;
  temp.data[2][3] = 1.0f;

  if (negativeDepthRange)
  {
    temp.data[2][2] = (minRange + maxRange) / (maxRange - minRange);
    temp.data[3][2] = -2.0f * minRange * maxRange / (maxRange - minRange);
  }
  else
  {
    temp.data[2][2] = minRange / (minRange - maxRange);
    temp.data[3][2] = -maxRange * temp.data[2][2];
  }
  return temp;
}

Matrix matrixPerspectiveFOVX(float fieldOfView, float aspectRatio, float minRange, float maxRange,
  bool negativeDepthRange)
{
  float const xScale = 1.0f / tanf(fieldOfView * 0.5f);
  float const yScale = xScale / aspectRatio;

  Matrix temp = matrixZero;

  temp.data[0][0] = xScale;
  temp.data[1][1] = yScale;
  temp.data[2][3] = 1.0f;

  if (negativeDepthRange)
  {
    temp.data[2][2] = (minRange + maxRange) / (maxRange - minRange);
    temp.data[3][2] = -2.0f * minRange * maxRange / (maxRange - minRange);
  }
  else
  {
    temp.data[2][2] = minRange / (minRange - maxRange);
    temp.data[3][2] = -maxRange * temp.data[2][2];
  }
  return temp;
}

Matrix matrixPerspectiveVOL(float width, float height, float minRange, float maxRange)
{
  float const rangeDiff = maxRange - minRange;
  Matrix temp = matrixZero;

  temp.data[0][0] = (2.0f * minRange) / width;
  temp.data[1][1] = (2.0f * minRange) / height;
  temp.data[2][2] = maxRange / rangeDiff;
  temp.data[2][3] = 1.0f;
  temp.data[3][2] = -minRange * rangeDiff;

  return temp;
}

Matrix matrixPerspectiveBDS(float left, float top, float right, float bottom, float minRange,
  float maxRange)
{
  Matrix temp = matrixZero;

  temp.data[0][0] = (2.0f * minRange) / (right - left);
  temp.data[1][1] = (2.0f * minRange) / (top - bottom);
  temp.data[2][0] = (left + right) / (left - right);
  temp.data[2][1] = (top + bottom) / (bottom - top);
  temp.data[2][2] = maxRange / (maxRange - minRange);
  temp.data[2][3] = 1.0f;
  temp.data[3][2] = minRange * maxRange / (minRange - maxRange);

  return temp;
}

Matrix matrixOrthogonalVOL(float width, float height, float minRange, float maxRange)
{
  Matrix temp = matrixZero;

  temp.data[0][0] = 2.0f / width;
  temp.data[1][1] = 2.0f / height;
  temp.data[2][2] = 1.0f / (minRange - maxRange);
  temp.data[3][2] = -maxRange / (minRange - maxRange);
  temp.data[3][3] = 1.0f;

  return temp;
}

Matrix matrixOrthogonalBDS(float left, float top, float right, float bottom, float minRange, float maxRange)
{
  Matrix temp = matrixZero;

  temp.data[0][0] = 2.0f / (right - left);
  temp.data[1][1] = 2.0f / (top - bottom);
  temp.data[2][2] = 1.0f / (minRange - maxRange);
  temp.data[3][0] = (left + right) / (left - right);
  temp.data[3][1] = (top + bottom) / (bottom - top);
  temp.data[3][2] = -maxRange / (minRange - maxRange);
  temp.data[3][3] = 1.0f;

  return temp;
}

bool matrixEquals(Matrix matrix1, Matrix matrix2)
{
  for (int_fast8_t j = 0; j < 4; j++)
    for (int_fast8_t i = 0; i < 4; i++)
      if (!nearlyEqualFloats(matrix1.data[j][i], matrix2.data[j][i]))
        return false;

  return true;
}

bool matrixEmpty(Matrix matrix)
{
  for (int_fast8_t j = 0; j < 4; j++)
    for (int_fast8_t i = 0; i < 4; i++)
      if (!nearlyZeroFloat(matrix.data[j][i]))
        return false;

  return true;
}

float detSub(float a1, float a2, float a3, float b1, float b2, float b3, float c1, float c2, float c3)
{
  return a1 * (b2 * c3 - b3 * c2) - b1 * (a2 * c3 - a3 * c2) + c1 * (a2 * b3 - a3 * b2);
}

// Quaternion constants and functions.

Quaternion const quaternionIdentity = {.x = 0.0f, .y = 0.0f, .z = 0.0f, .w = 1.0f};

Quaternion quaternionMultiply(Quaternion quat1, Quaternion quat2)
{
  return (Quaternion){
    .x = quat2.w * quat1.x + quat2.x * quat1.w + quat2.z * quat1.y - quat2.y * quat1.z,
    .y = quat2.w * quat1.y + quat2.y * quat1.w + quat2.x * quat1.z - quat2.z * quat1.x,
    .z = quat2.w * quat1.z + quat2.z * quat1.w + quat2.y * quat1.x - quat2.x * quat1.y,
    .w = quat2.w * quat1.w - quat2.x * quat1.x - quat2.y * quat1.y - quat2.z * quat1.z};
}

float quaternionDot(Quaternion quat1, Quaternion quat2)
{
  return (quat1.x * quat2.x) + (quat1.y * quat2.y) + (quat1.z * quat2.z) + (quat1.w * quat2.w);
}

float quaternionLength(Quaternion quat)
{
  return sqrtf(quaternionDot(quat, quat));
}

float quaternionAngle(Quaternion quat)
{
  return acosf(quat.w) * 2.0f;
}

Vector quaternionAxis(Quaternion quat)
{
  float const temp1 = 1.0f - quat.w * quat.w;
  if (temp1 < vectorEpsilon)
    return vectorAxisY;

  float const temp2 = 1.0f / sqrtf(temp1);

  return (Vector){.x = quat.x * temp2, .y = quat.y * temp2, .z = quat.z * temp2};
}

Quaternion quaternionNormalize(Quaternion quat)
{
  float const amp = quaternionLength(quat);
  if (!nearlyZeroFloat(amp))
  {
    float const invAmp = 1.0f / amp;
    return (Quaternion){
      .x = quat.x * invAmp, .y = quat.y * invAmp, .z = quat.z * invAmp, .w = quat.w * invAmp};
  }
  else
    return quat;
}

Quaternion quaternionConjugate(Quaternion quat)
{
  return (Quaternion){.x = -quat.x, .y = -quat.y, .z = -quat.z, .w = quat.w};
}

Quaternion quaternionExponentiate(Quaternion quat, float exponent)
{
  if (nearlyEqualFloats(quat.w, 1.0f))
    return quat;

  float const alpha = acosf(quat.w);
  float const newAlpha = alpha * exponent;
  float const compMult = sinf(newAlpha) / sinf(alpha);

  return (Quaternion){
    .x = quat.x * compMult, .y = quat.y * compMult, .z = quat.z * compMult, .w = cosf(newAlpha)};
}

Quaternion quaternionSlerp(Quaternion quat1, Quaternion quat2, float theta)
{
  if (theta <= 0.0f)
    return quat1;
  else if (theta >= 1.0f)
    return quat2;

  float cosOmega = quaternionDot(quat1, quat2);
  Quaternion tempQuat = quat2;

  if (cosOmega < 0.0f)
  {
    tempQuat.x = -tempQuat.x;
    tempQuat.y = -tempQuat.y;
    tempQuat.z = -tempQuat.z;
    tempQuat.w = -tempQuat.w;
    cosOmega = -cosOmega;
  }

  float kappa1, kappa2;

  if (cosOmega > 1.0f - vectorEpsilon)
  {
    kappa1 = 1.0f - theta;
    kappa2 = theta;
  }
  else
  {
    float const sinOmega = sqrtf(1.0f - cosOmega * cosOmega);
    float const omega = atan2f(sinOmega, cosOmega);
    float const oneOverSinOmega = 1.0f / sinOmega;

    kappa1 = sinf((1.0f - theta) * omega) * oneOverSinOmega;
    kappa2 = sinf(theta * omega) * oneOverSinOmega;
  }

  return (Quaternion){
    .x = kappa1 * quat1.x + kappa2 * tempQuat.x,
    .y = kappa1 * quat1.y + kappa2 * tempQuat.y,
    .z = kappa1 * quat1.z + kappa2 * tempQuat.z,
    .w = kappa1 * quat1.w + kappa2 * tempQuat.w};
}

Matrix quaternionMatrix(Quaternion quat)
{
  Matrix temp = matrixIdentity;

  temp.data[0][0] = 1.0f - 2.0f * quat.y * quat.y - 2.0f * quat.z * quat.z;
  temp.data[0][1] = 2.0f * quat.x * quat.y + 2.0f * quat.w * quat.z;
  temp.data[0][2] = 2.0f * quat.x * quat.z - 2.0f * quat.w * quat.y;
  temp.data[1][0] = 2.0f * quat.x * quat.y - 2.0f * quat.w * quat.z;
  temp.data[1][1] = 1.0f - 2.0f * quat.x * quat.x - 2.0f * quat.z * quat.z;
  temp.data[1][2] = 2.0f * quat.y * quat.z + 2.0f * quat.w * quat.x;
  temp.data[2][0] = 2.0f * quat.x * quat.z + 2.0f * quat.w * quat.y;
  temp.data[2][1] = 2.0f * quat.y * quat.z - 2.0f * quat.w * quat.x;
  temp.data[2][2] = 1.0f - 2.0f * quat.x * quat.x - 2.0f * quat.y * quat.y;

  return temp;
}

Quaternion matrixQuaternion(Matrix matrix)
{
  Quaternion const temp = {
    matrix.data[0][0] - matrix.data[1][1] - matrix.data[2][2],
    matrix.data[1][1] - matrix.data[0][0] - matrix.data[2][2],
    matrix.data[2][2] - matrix.data[0][0] - matrix.data[1][1],
    matrix.data[0][0] + matrix.data[1][1] + matrix.data[2][2]};

  int_fast8_t index = 0;
  float maxValue = temp.w;

  if (temp.x > maxValue)
  {
    maxValue = temp.x;
    index = 1;
  }
  else if (temp.y > maxValue)
  {
    maxValue = temp.y;
    index = 2;
  }
  else if (temp.z > maxValue)
  {
    maxValue = temp.z;
    index = 3;
  }

  float const highValue = sqrtf(maxValue + 1.0f) * 0.5f;
  float const multValue = 0.25f / highValue;

  Quaternion res;

  switch (index)
  {
  case 1:
    res.x = highValue;
    res.y = (matrix.data[0][1] + matrix.data[1][0]) * multValue;
    res.z = (matrix.data[2][0] + matrix.data[0][2]) * multValue;
    res.w = (matrix.data[1][2] - matrix.data[2][1]) * multValue;
    break;
  case 2:
    res.x = (matrix.data[0][1] + matrix.data[1][0]) * multValue;
    res.y = highValue;
    res.z = (matrix.data[1][2] + matrix.data[2][1]) * multValue;
    res.w = (matrix.data[2][0] - matrix.data[0][2]) * multValue;
    break;
  case 3:
    res.x = (matrix.data[2][0] + matrix.data[0][2]) * multValue;
    res.y = (matrix.data[1][2] + matrix.data[2][1]) * multValue;
    res.z = highValue;
    res.w = (matrix.data[0][1] - matrix.data[1][0]) * multValue;
    break;
  default:
    res.x = (matrix.data[1][2] - matrix.data[2][1]) * multValue;
    res.y = (matrix.data[2][0] - matrix.data[0][2]) * multValue;
    res.z = (matrix.data[0][1] - matrix.data[1][0]) * multValue;
    res.w = highValue;
    break;
  }
  return res;
}

Quaternion quaternionRotateX(float angle)
{
  float const newAngle = angle * 0.5f;
  return (Quaternion){.x = sinf(newAngle), .y = 0.0f, .z = 0.0f, .w = cosf(newAngle)};
}

Quaternion quaternionRotateY(float angle)
{
  float const newAngle = angle * 0.5f;
  return (Quaternion){.x = 0.0f, .y = sinf(newAngle), .z = 0.0f, .w = cosf(newAngle)};
}

Quaternion quaternionRotateZ(float angle)
{
  float const newAngle = angle * 0.5f;
  return (Quaternion){.x = 0.0f, .y = 0.0f, .z = sinf(newAngle), .w = cosf(newAngle)};
}

Quaternion quaternionRotate(Vector axis, float angle)
{
  float const amp = vectorLength(axis);
  if (!nearlyZeroFloat(amp))
  {
    float const newAngle = angle * 0.5f;
    float const sinDivAmp = sinf(newAngle) / amp;

    return (Quaternion){
      .x = axis.x * sinDivAmp, .y = axis.y * sinDivAmp, .z = axis.z * sinDivAmp, .w = cosf(newAngle)};
  }
  else
    return quaternionIdentity;
}

Quaternion quaternionRotateObjectToInertial(Vector angles)
{
  float const sinHeading = sinf(angles.y);
  float const cosHeading = cosf(angles.y);
  float const sinPitch = sinf(angles.x);
  float const cosPitch = cosf(angles.x);
  float const sinBank = sinf(angles.z);
  float const cosBank = cosf(angles.z);

  return (Quaternion){
    .x = cosHeading * sinPitch * cosBank + sinHeading * cosPitch * sinBank,
    .y = -cosHeading * sinPitch * sinBank + sinHeading * cosPitch * cosBank,
    .z = -sinHeading * sinPitch * cosBank + cosHeading * cosPitch * sinBank,
    .w = cosHeading * cosPitch * cosBank + sinHeading * sinPitch * sinBank};
}

Quaternion quaternionRotateInertialToObject(Vector angles)
{
  float const sinHeading = sinf(angles.y);
  float const cosHeading = cosf(angles.y);
  float const sinPitch = sinf(angles.x);
  float const cosPitch = cosf(angles.x);
  float const sinBank = sinf(angles.z);
  float const cosBank = cosf(angles.z);

  return (Quaternion){
    .x = -cosHeading * sinPitch * cosBank - sinHeading * cosPitch * sinBank,
    .y = cosHeading * sinPitch * sinBank - sinHeading * cosPitch * cosBank,
    .z = sinHeading * sinPitch * cosBank - cosHeading * cosPitch * sinBank,
    .w = cosHeading * cosPitch * cosBank + sinHeading * sinPitch * sinBank};
}

// Rect constants and functions.

Rect const rectZero = {.left = 0, .top = 0, .width = 0, .height = 0};

Rect makeRect(int32_t left, int32_t top, int32_t width, int32_t height)
{
  return (Rect){.left = left, .top = top, .width = width, .height = height};
}

Rect makeBounds(int32_t left, int32_t top, int32_t right, int32_t bottom)
{
  return (Rect){.left = left, .top = top, .width = right - left, .height = bottom - top};
}

Point rectGetSize(Rect rect)
{
  return (Point){.x = rect.width, .y = rect.height};
}

Rect rectSetSize(Rect rect, Point size)
{
  rect.width = size.x;
  rect.height = size.y;
  return rect;
}

int32_t rectGetRight(Rect rect)
{
  return rect.left + rect.width;
}

Rect rectSetRight(Rect rect, int32_t right)
{
  rect.width = right - rect.left;
  return rect;
}

int32_t rectGetBottom(Rect rect)
{
  return rect.top + rect.height;
}

Rect rectSetBottom(Rect rect, int32_t bottom)
{
  rect.height = bottom - rect.top;
  return rect;
}

Point rectGetTopLeft(Rect rect)
{
  return (Point){.x = rect.left, .y = rect.top};
}

Rect rectSetTopLeft(Rect rect, Point topLeft)
{
  rect.left = topLeft.x;
  rect.top = topLeft.y;
  return rect;
}

Point rectGetBottomRight(Rect rect)
{
  return (Point){.x = rect.left + rect.width, .y = rect.top + rect.height};
}

Rect rectSetBottomRight(Rect rect, Point bottomRight)
{
  rect.width = bottomRight.x - rect.left;
  rect.height = bottomRight.y - rect.top;
  return rect;
}

Point rectGetTopRight(Rect rect)
{
  return (Point){.x = rectGetRight(rect), .y = rect.top};
}

Point rectGetBottomLeft(Rect rect)
{
  return (Point){rect.left, rectGetBottom(rect)};
}

bool rectEquals(Rect rect1, Rect rect2)
{
  return rect1.left == rect2.left && rect1.top == rect2.top && rect1.width == rect2.width &&
    rect1.height == rect2.height;
}

bool rectEmpty(Rect rect)
{
  return rect.width <= 0 || rect.height <= 0;
}

bool pointInRect(Point point, Rect rect)
{
  return point.x >= rect.left && point.x < rectGetRight(rect) && point.y >= rect.top &&
    point.y < rectGetBottom(rect);
}

bool rectInRect(Rect rect, Rect largerRect)
{
  return largerRect.left >= rect.left && rectGetRight(largerRect) <= rectGetRight(rect) &&
    largerRect.top >= rect.top && rectGetBottom(largerRect) <= rectGetBottom(rect);
}

bool rectOverlap(Rect rect1, Rect rect2)
{
  return rect2.left < rectGetRight(rect1) && rectGetRight(rect2) > rect1.left &&
    rect2.top < rectGetBottom(rect1) && rectGetBottom(rect2) > rect1.top;
}

Rect rectIntersect(Rect rect1, Rect rect2)
{
  if (rectOverlap(rect1, rect2))
  {
    return makeBounds(
      MAX(rect1.left, rect2.left),
      MAX(rect1.top, rect2.top),
      MIN(rectGetRight(rect1), rectGetRight(rect2)),
      MIN(rectGetBottom(rect1), rectGetBottom(rect2)));
  }
  else
    return rectZero;
}

Rect rectJoin(Rect rect1, Rect rect2)
{
  return makeBounds(
    MIN(rect1.left, rect2.left),
    MIN(rect1.top, rect2.top),
    MAX(rectGetRight(rect1), rectGetRight(rect2)),
    MAX(rectGetBottom(rect1), rectGetBottom(rect2)));
}

Rect rectOffset(Rect rect, Point delta)
{
  return (Rect){
    .left = rect.left + delta.x, .top = rect.top + delta.y, .width = rect.width, .height = rect.height};
}

Rect rectInflate(Rect rect, Point delta)
{
  return (Rect){.left = rect.left - delta.x, .top = rect.top - delta.y,
    .width = rect.width + delta.x * 2, .height = rect.height + delta.y * 2};
}

// RectF constants and functions.

RectF const rectZeroF = {.left = 0.0f, .top = 0.0f, .width = 0.0f, .height = 0.0f};

RectF makeRectF(float left, float top, float width, float height)
{
  return (RectF){.left = left, .top = top, .width = width, .height = height};
}

RectF makeBoundsF(float left, float top, float right, float bottom)
{
  return (RectF){.left = left, .top = top, .width = right - left, .height = bottom - top};
}

RectF rectF(Rect rect)
{
  return (RectF){.left = (float)rect.left, .top = (float)rect.top,
    .width = (float)rect.width, .height = (float)rect.height};
}

Rect rectToIntRectF(RectF rect)
{
  return (Rect){.left = (int32_t)roundf(rect.left), .top = (int32_t)roundf(rect.top),
    .width = (int32_t)roundf(rect.width), .height = (int32_t)roundf(rect.height)};
}

PointF rectGetSizeF(RectF rect)
{
  return (PointF){.x = rect.width, .y = rect.height};
}

RectF rectSetSizeF(RectF rect, PointF size)
{
  rect.width = size.x;
  rect.height = size.y;
  return rect;
}

float rectGetRightF(RectF rect)
{
  return rect.left + rect.width;
}

RectF rectSetRightF(RectF rect, float right)
{
  rect.width = right - rect.left;
  return rect;
}

float rectGetBottomF(RectF rect)
{
  return rect.top + rect.height;
}

RectF rectSetBottomF(RectF rect, float bottom)
{
  rect.height = bottom - rect.top;
  return rect;
}

PointF rectGetTopLeftF(RectF rect)
{
  return (PointF){.x = rect.left, .y = rect.top};
}

RectF rectSetTopLeftF(RectF rect, PointF topLeft)
{
  rect.left = topLeft.x;
  rect.top = topLeft.y;
  return rect;
}

PointF rectGetBottomRightF(RectF rect)
{
  return (PointF){.x = rect.left + rect.width, .y = rect.top + rect.height};
}

RectF rectSetBottomRightF(RectF rect, PointF bottomRight)
{
  rect.width = bottomRight.x - rect.left;
  rect.height = bottomRight.y - rect.top;
  return rect;
}

PointF rectGetTopRightF(RectF rect)
{
  return (PointF){.x = rectGetRightF(rect), .y = rect.top};
}

PointF rectGetBottomLeftF(RectF rect)
{
  return (PointF){rect.left, rectGetBottomF(rect)};
}

bool rectEqualsF(RectF rect1, RectF rect2)
{
  return nearlyEqualFloats(rect1.left, rect2.left) && nearlyEqualFloats(rect1.top, rect2.top) &&
    nearlyEqualFloats(rect1.width, rect2.width) && nearlyEqualFloats(rect1.height, rect2.height);
}

bool rectEmptyF(RectF rect)
{
  return rect.width < vectorEpsilon || rect.height < vectorEpsilon;
}

bool pointInRectF(PointF point, RectF rect)
{
  return point.x >= rect.left && point.x < rectGetRightF(rect) && point.y >= rect.top &&
    point.y < rectGetBottomF(rect);
}

bool rectInRectF(RectF rect, RectF largerRect)
{
  return largerRect.left >= rect.left && rectGetRightF(largerRect) <= rectGetRightF(rect) &&
    largerRect.top >= rect.top && rectGetBottomF(largerRect) <= rectGetBottomF(rect);
}

bool rectOverlapF(RectF rect1, RectF rect2)
{
  return rect2.left < rectGetRightF(rect1) && rectGetRightF(rect2) > rect1.left &&
    rect2.top < rectGetBottomF(rect1) && rectGetBottomF(rect2) > rect1.top;
}

RectF rectIntersectF(RectF rect1, RectF rect2)
{
  if (rectOverlapF(rect1, rect2))
  {
    return makeBoundsF(
      MAX(rect1.left, rect2.left),
      MAX(rect1.top, rect2.top),
      MIN(rectGetRightF(rect1), rectGetRightF(rect2)),
      MIN(rectGetBottomF(rect1), rectGetBottomF(rect2)));
  }
  else
    return rectZeroF;
}

RectF rectJoinF(RectF rect1, RectF rect2)
{
  return makeBoundsF(
    MIN(rect1.left, rect2.left),
    MIN(rect1.top, rect2.top),
    MAX(rectGetRightF(rect1), rectGetRightF(rect2)),
    MAX(rectGetBottomF(rect1), rectGetBottomF(rect2)));
}

RectF rectOffsetF(RectF rect, PointF delta)
{
  return (RectF){
    .left = rect.left + delta.x, .top = rect.top + delta.y, .width = rect.width, .height = rect.height};
}

RectF rectInflateF(RectF rect, PointF delta)
{
  return (RectF){.left = rect.left - delta.x, .top = rect.top - delta.y,
    .width = rect.width + delta.x * 2.0f, .height = rect.height + delta.y * 2.0f};
}

// Quad constants and functions.

Quad const quadZero = {{0.0f, 0.0f}, {0.0f, 0.0f}, {0.0f, 0.0f}, {0.0f, 0.0f}};
Quad const quadUnity = {{0.0f, 0.0f}, {1.0f, 0.0f}, {1.0f, 1.0f}, {0.0f, 1.0f}};

Quad makeQuadWith(PointF topLeft, PointF topRight, PointF bottomRight, PointF bottomLeft)
{
  return (Quad){
    .topLeft = topLeft, .topRight = topRight, .bottomRight = bottomRight, .bottomLeft = bottomLeft};
}

Quad makeQuad(float left, float top, float width, float height)
{
  return (Quad){
    .topLeft = (PointF){.x = left, .y = top},
    .topRight = (PointF){.x = left + width, .y = top},
    .bottomRight = (PointF){.x = left + width, .y = top + height},
    .bottomLeft = (PointF){.x = left, .y = top + height}};
}

Quad makeQuadScaled(float left, float top, float width, float height, float scale, bool centered)
{
  if (nearlyEqualFloats(scale, 1.0f))
    return makeQuad(left, top, width, height);

  if (centered)
  {
    float const newWidth = width * scale;
    float const newHeight = height * scale;

    return makeQuad(
      left + (width - newWidth) * 0.5f, top + (height - newHeight) * 0.5f, newWidth, newHeight);
  }
  else
    return makeQuad(left, top, width * scale, height * scale);
}

Quad makeQuadRotatedAt(PointF origin, PointF size, PointF center, float angle, float scale)
{
  float const sinAngle = sinf(angle);
  float const cosAngle = cosf(angle);

  bool const scaled = !nearlyEqualFloats(scale, 1.0f);

  Quad rect = makeQuad(-center.x, -center.y, size.x, size.y);
  PointF* destPos = &rect.topLeft;

  for (int_fast8_t i = 0; i < 4; ++i)
  {
    if (scaled)
      *destPos = pointRescaleF(*destPos, scale);

    PointF const newPos = makePointF(
      destPos->x * cosAngle - destPos->y * sinAngle,
      destPos->y * cosAngle + destPos->x * sinAngle);

    *destPos = pointAddF(newPos, origin);
    destPos++;
  }
  return rect;
}

Quad makeQuadRotated(PointF origin, PointF size, float angle, float scale)
{
  return makeQuadRotatedAt(origin, size, pointRescaleF(size, 0.5f), angle, scale);
}

Quad makeQuadRotatedTL(PointF topLeft, PointF size, PointF center, float angle, float scale)
{
  return quadOffset(makeQuadRotatedAt(topLeft, size, center, angle, scale), center);
}

Quad makeQuadSkewedHoriz(PointF origin, PointF size, float angle)
{
  float const offset = cosf(angle) * size.y * 0.5f;

  Quad quad = makeQuad(origin.x - size.x * 0.5f, origin.y - size.y * 0.5f, size.x, size.y);
  quad.topLeft.x -= offset;
  quad.topRight.x -= offset;
  quad.bottomRight.x += offset;
  quad.bottomLeft.x += offset;

  return quad;
}

Quad makeQuadSkewedVert(PointF origin, PointF size, float angle)
{
  float const offset = cosf(angle) * size.x * 0.5f;

  Quad quad = makeQuad(origin.x - size.x * 0.5f, origin.y - size.y * 0.5f, size.x, size.y);
  quad.topLeft.y -= offset;
  quad.bottomLeft.y -= offset;
  quad.topRight.y += offset;
  quad.bottomRight.y += offset;

  return quad;
}

Quad quadFromRect(Rect rect)
{
  return (Quad){
    .topLeft = pointF(rectGetTopLeft(rect)),
    .topRight = pointF(rectGetTopRight(rect)),
    .bottomRight = pointF(rectGetBottomRight(rect)),
    .bottomLeft = pointF(rectGetBottomLeft(rect))};
}

Quad quadFromRectF(RectF rect)
{
  return (Quad){
    .topLeft = rectGetTopLeftF(rect),
    .topRight = rectGetTopRightF(rect),
    .bottomRight = rectGetBottomRightF(rect),
    .bottomLeft = rectGetBottomLeftF(rect)};
}

Quad quadScale(Quad quad, float scale, bool centered)
{
  if (nearlyEqualFloats(scale, 1.0f))
    return quad;

  if (centered)
  {
    PointF const center = pointRescaleF(pointAddF(pointAddF(quad.topLeft, quad.topRight),
      pointAddF(quad.bottomRight, quad.bottomLeft)), 0.25f);

    return (Quad){
      .topLeft = pointLerpF(center, quad.topLeft, scale),
      .topRight = pointLerpF(center, quad.topRight, scale),
      .bottomRight = pointLerpF(center, quad.bottomRight, scale),
      .bottomLeft = pointLerpF(center, quad.bottomLeft, scale)};
  }
  else
    return (Quad){
      .topLeft = pointRescaleF(quad.topLeft, scale),
      .topRight = pointRescaleF(quad.topRight, scale),
      .bottomRight = pointRescaleF(quad.bottomRight, scale),
      .bottomLeft = pointRescaleF(quad.bottomLeft, scale)};
}

Quad quadMirror(Quad quad)
{
  return (Quad){.topLeft = quad.topRight, .topRight = quad.topLeft,
    .bottomRight = quad.bottomLeft, .bottomLeft = quad.bottomRight};
}

Quad quadFlip(Quad quad)
{
  return (Quad){.topLeft = quad.bottomLeft, .topRight = quad.bottomRight,
    .bottomRight = quad.topRight, .bottomLeft = quad.topLeft};
}

Quad quadTransform(Quad quad, Matrix3x2 matrix)
{
  return (Quad){
    .topLeft = pointTransformF(quad.topLeft, matrix),
    .topRight = pointTransformF(quad.topRight, matrix),
    .bottomRight = pointTransformF(quad.bottomRight, matrix),
    .bottomLeft = pointTransformF(quad.bottomLeft, matrix)};
}

Quad quadOffset(Quad quad, PointF delta)
{
  return (Quad){
    .topLeft = pointAddF(quad.topLeft, delta),
    .topRight = pointAddF(quad.topRight, delta),
    .bottomRight = pointAddF(quad.bottomRight, delta),
    .bottomLeft = pointAddF(quad.bottomLeft, delta)};
}

// Service functions.

bool nearlyEqualFloats(float left, float right)
{
  return
    left == right ||
    fabsf(left - right) <= fmaxf(fabsf(left), fabsf(right)) * vectorEpsilon;
}

bool nearlyZeroFloat(float value)
{
  return fabsf(value) <= vectorEpsilon;
}

float saturateFloat(float value, float min, float max)
{
  float const tempValue = value > max ? max : value;
  return tempValue < min ? min : tempValue;
}

float lerpFloat(float value1, float value2, float theta)
{
  return value1 + (value2 - value1) * theta;
}

#endif