/**
 * @file kmath.h
 * @author Travis Vroman (travis@kohiengine.com)
 * @brief This file contains definitions for various important constant values
 * as well as functions for many common math types. Note that this math library
 * is all written to be right-handed (-z forward, +y up) and in column-major format.
 * @version 2.0
 * @date 2025-01-26
 *
 * @copyright Kohi Game Engine is Copyright (c) Travis Vroman 2021-2025
 *
 */

#pragma once

#include "defines.h"
#include "math_types.h"
#include "memory/kmemory.h"

#define KMATH_USE_CUSTOM_RAND 1

/** @brief An approximate representation of PI. */
#define K_PI 3.14159265358979323846f

/** @brief An approximate representation of PI multiplied by 2. */
#define K_2PI (2.0f * K_PI)

/** @brief An approximate representation of PI multiplied by 4. */
#define K_4PI (4.0f * K_PI)

/** @brief An approximate representation of PI divided by 2. */
#define K_HALF_PI (0.5f * K_PI)

/** @brief An approximate representation of PI divided by 4. */
#define K_QUARTER_PI (0.25f * K_PI)

/** @brief One divided by an approximate representation of PI. */
#define K_ONE_OVER_PI (1.0f / K_PI)

/** @brief One divided by half of an approximate representation of PI. */
#define K_ONE_OVER_TWO_PI (1.0f / K_2PI)

/** @brief An approximation of the square root of 2. */
#define K_SQRT_TWO 1.41421356237309504880f

/** @brief An approximation of the square root of 3. */
#define K_SQRT_THREE 1.73205080756887729352f

/** @brief One divided by an approximation of the square root of 2. */
#define K_SQRT_ONE_OVER_TWO 0.70710678118654752440f

/** @brief One divided by an approximation of the square root of 3. */
#define K_SQRT_ONE_OVER_THREE 0.57735026918962576450f

/** @brief A multiplier used to convert degrees to radians. */
#define K_DEG2RAD_MULTIPLIER (K_PI / 180.0f)

/** @brief A multiplier used to convert radians to degrees. */
#define K_RAD2DEG_MULTIPLIER (180.0f / K_PI)

/** @brief The multiplier to convert seconds to microseconds. */
#define K_SEC_TO_US_MULTIPLIER (1000.0f * 1000.0f)

/** @brief The multiplier to convert seconds to milliseconds. */
#define K_SEC_TO_MS_MULTIPLIER 1000.0f

/** @brief The multiplier to convert milliseconds to seconds. */
#define K_MS_TO_SEC_MULTIPLIER 0.001f

/** @brief A huge number that should be larger than any valid number used. */
#define K_INFINITY (1e30f * 1e30f)

/** @brief Smallest positive number where 1.0 + FLOAT_EPSILON != 0 */
#define K_FLOAT_EPSILON 1.192092896e-07f

#define K_FLOAT_MAX 3.40282e+38F
#define K_FLOAT_MIN -K_FLOAT_MAX

// Used to test if a float is close to zero, specifically for sweep tests.
#define SWEEP_EPSILON 1e-5f

#if KMATH_USE_CUSTOM_RAND
#	define KRAND_MAX 0x7FFFFFFF
#else
#	define KRAND_MAX RAND_MAX
#endif

// ------------------------------------------
// General math functions
// ------------------------------------------

/**
 * Swaps the values in the given float pointers.
 * @param a A pointer to the first float.
 * @param b A pointer to the second float.
 */
KINLINE void kswapf (f32 *a, f32 *b) {
	f32 temp = *a;
	*a = *b;
	*b = temp;
}

/** @brief Returns 0.0f if x == 0.0f, -1.0f if negative, otherwise 1.0f. */
KINLINE f32 ksign (f32 x) {
	if (x >= 0.0f) {
		return 1.0f;
	}
	return -1.0f;
}

/** @brief Compares x to edge, returning 0 if x < edge; otherwise 1.0f; */
KINLINE f32 kstep (f32 edge, f32 x) {
	return x < edge ? 0.0f : 1.0f;
}

/**
 * @brief Calculates the sine of x.
 *
 * @param x The number to calculate the sine of.
 * @return The sine of x.
 */
KAPI f32 ksin (f32 x);

/**
 * @brief Calculates the cosine of x.
 *
 * @param x The number to calculate the cosine of.
 * @return The cosine of x.
 */
KAPI f32 kcos (f32 x);

/**
 * @brief Calculates the tangent of x.
 *
 * @param x The number to calculate the tangent of.
 * @return The tangent of x.
 */
KAPI f32 ktan (f32 x);

/**
 * @brief Calculates the arctangent of x.
 *
 * @param x The number to calculate the arctangent of.
 * @return The arctangent of x.
 */
KAPI f32 katan (f32 x);

KAPI f32 katan2 (f32 x, f32 y);

KAPI f32 kasin (f32 x);

/**
 * @brief Calculates the arc cosine of x.
 *
 * @param x The number to calculate the arc cosine of.
 * @return The arc cosine of x.
 */
KAPI f32 kacos (f32 x);

/**
 * @brief Calculates the square root of x.
 *

	f32 half_width = size.x * 0.5f;
	f32 half_height = size.y * 0.5f;
	f32 half_depth = size.z * 0.5f;

	kgeometry out_geometry = {0};
	out_geometry.name = name;
	out_geometry.type = KGEOMETRY_TYPE_3D_STATIC_COLOUR_ONLY;
	out_geometry.generation = INVALID_ID_U16; // NOTE: generation is 0 because this is technically the first "update"
	out_geometry.extents.min = (vec3){-half_width, -half_height, -half_depth};
	out_geometry.extents.max = (vec3){half_width, half_height, half_depth};
	// Always 0 since min/max of each axis are -/+ half of the size.
	out_geometry.center = vec3_zero();
	out_geometry.vertex_element_size = sizeof(colour_vertex_3d);
	out_geometry.vertex_count = 2 * 12; // 12 lines to make a cube.
	out_geometry.vertices = KALLOC_TYPE_CARRAY(colour_vertex_3d, out_geometry.vertex_count);
	out_geometry.vertex_buffer_offset = INVALID_ID_U64;
	// NOTE: line-based boxes do not have/need indices.
	out_geometry.index_count = 0;
	out_geometry.index_element_size = sizeof(u32);
	out_geometry.indices = 0;
	out_geometry.index_buffer_offset = INVALID_ID_U64;

	extents_3d extents = {0};
	extents.min.x = -half_width;
	extents.min.y = -half_height;
	extents.min.z = -half_depth;
	extents.max.x = half_width;
	extents.max.y = half_height;
	extents.max.z = half_depth;

	geometry_recalculate_line_box3d_by_extents(&out_geometry, extents);

	// Set the default colour.
	colour_vertex_3d* verts = (colour_vertex_3d*)out_geometry.vertices;
	for (u32 i = 0; i < out_geometry.vertex_count; ++i) {
		verts[i].colour = vec4_one();
	}

	return out_geometry;
 * @param x The number to calculate the square root of.
 * @return The square root of x.
 */
KAPI f32 ksqrt (f32 x);

/**
 * @brief Calculates the absolute value of x.
 *
 * @param x The number to get the absolute value of.
 * @return The absolute value of x.
 */
KINLINE f32 kabs (f32 x) {
	return x < 0.0f ? -x : x;
}

/**
 * @brief Returns the largest integer value less than or equal to x.
 *
 * @param x The value to be examined.
 * @return the largest integer value less than or equal to x.
 */
KINLINE f32 kfloor (f32 x) {
	i32 i = (i32)x;
	if ((f32)i > x) {
		--i;
	}
	return (f32)i;
}

/**
 * @brief Returns the smallest integer value greater than or equal to x.
 *
 * @param x The value to be examined.
 * @return the smallest integer value greater than or equal to x.
 */
KINLINE f32 kceil (f32 x) {
	i32 i = (i32)x;
	if ((f32)i < x) {
		++i;
	}
	return (f32)i;
}

/**
 * @brief Computes the logarithm of x.
 *
 * @param x The value to be examined.
 * @return The logarithm of x.
 */
KAPI f32 klog (f32 x);

/**
 * @brief Computes the base-2 logarithm of x (i.e. how many times x can be divided by 2).
 *
 * @param x The value to be examined.
 * @return The base-2 logarithm of x.
 */
KAPI f32 klog2 (f32 x);

KAPI f32 kpow (f32 x, f32 y);

KINLINE f32 kmod (f32 x, f32 y) {
	if (y == 0.0f) {
		return 0.0f;
	}
	return x - (f32)((i32)(x / y)) * y;
}

KINLINE f32 kwrap_pi (f32 x) {
	while (x > K_PI) {
		x -= K_2PI;
	}
	while (x < -K_PI) {
		x += K_2PI;
	}
	return x;
}

KINLINE f32 klerp (f32 from, f32 to, f32 t) {
	return from + t * (to - from);
}

KINLINE f32 klerp_inv (f32 from, f32 to, f32 t) {
	return (t - from) / (to - from);
}

KINLINE f32 kfalloff_smooth (f32 t) {
	if (t >= 1.0f) {
		return 0.0f;
	}
	if (t <= 0.0f) {
		return 1.0f;
	}

	return 1.0f - (t * t * (3.0f - 2.0f * t));
}

/**
 * @brief Indicates if the value is a power of 2. 0 is considered _not_ a power
 * of 2.
 * @param value The value to be interpreted.
 * @returns True if a power of 2, otherwise false.
 */
KINLINE b8 is_power_of_2 (u64 value) {
	return (value != 0) && ((value & (value - 1)) == 0);
}

KINLINE b8 kis_finite (f32 x) {
	union {
		f32 f;
		u32 i;
	} u;

	u.f = x;

	/* Exponent bits are bits 23-30.
	   All ones means infinity or NaN. */
	return ((u.i & 0x7F800000u) != 0x7F800000u);
}

/**
 * @brief Decode a hexadecimal digit to it's decimal value (i.e. 0-15).
 *
 * @param c The hex digit to be decoded.
 * @return The decoded hex digit value if valid; otherwise -1.
 */
i32 hex_value (char c);

/**
 * @brief Returns a random integer.
 *
 * @return A random integer.
 */
KAPI i32 krandom (void);

/**
 * @brief Returns a random integer that is within the given range (inclusive).
 *
 * @param min The minimum of the range.
 * @param max The maximum of the range.
 * @return A random integer.
 */
KAPI i32 krandom_in_range (i32 min, i32 max);

/**
 * @brief Returns a random unsigned 64-bit integer.
 * @return A random unsigned 64-bit integer.
 */
KAPI u64 krandom_u64 (void);

/**
 * @brief Returns a random floating-point number.
 *
 * @return A random floating-point number.
 */
KAPI f32 kfrandom (void);

/**
 * @brief Returns a random floating-point number that is within the given range
 * (inclusive).
 *
 * @param min The minimum of the range.
 * @param max The maximum of the range.
 * @return A random floating-point number.
 */
KAPI f32 kfrandom_in_range (f32 min, f32 max);

/**
 * @brief Perform Hermite interpolation between two values.
 *
 * @param edge_0 The lower edge of the Hermite function.
 * @param edge_1 The upper edge of the Hermite function.
 * @param x The value to interpolate.
 * @return The interpolated value.
 */
KINLINE f32 ksmoothstep (f32 edge_0, f32 edge_1, f32 x) {
	f32 t = KCLAMP((x - edge_0) / (edge_1 - edge_0), 0.0f, 1.0f);
	return t * t * (3.0 - 2.0 * t);
}

/**
 * @brief Returns the attenuation of x based off distance from the midpoint of min and max.
 *
 * @param min The minimum value.
 * @param max The maximum value.
 * @param x The value to attenuate.
 * @return The attenuation of x based on distance of the midpoint of min and max.
 */
KAPI f32 kattenuation_min_max (f32 min, f32 max, f32 x);

/**
 * @brief Compares the two floats and returns true if both are less
 * than K_FLOAT_EPSILON apart; otherwise false.
 */
KINLINE b8 kfloat_compare (f32 f_0, f32 f_1) {
	return kabs(f_0 - f_1) < K_FLOAT_EPSILON;
}

// ------------------------------------------
// Vector 2
// ------------------------------------------

/**
 * @brief Creates and returns a new 2-element vector using the supplied values.
 *
 * @param x The x value.
 * @param y The y value.
 * @return A new 2-element vector.
 */
KINLINE vec2 vec2_create (f32 x, f32 y) {
	return (vec2){x, y};
}

/**
 * @brief Creates and returns a new 2-element vector using the supplied scalar for all components.
 *
 * @param scalar The scalar value.
 * @return A new 2-element vector.
 */
KINLINE vec2 vec2_from_scalar (f32 scalar) { return (vec2){scalar, scalar}; }

/**
 * @brief Creates and returns a 2-component vector with all components set to
 * 0.0f.
 */
KINLINE vec2 vec2_zero (void) { return (vec2){0.0f, 0.0f}; }

/**
 * @brief Creates and returns a 2-component vector with all components set
 * to 1.0f.
 */
KINLINE vec2 vec2_one (void) { return (vec2){1.0f, 1.0f}; }

/**
 * @brief Creates and returns a 2-component vector pointing up (0, 1).
 */
KINLINE vec2 vec2_up (void) { return (vec2){0.0f, 1.0f}; }

/**
 * @brief Creates and returns a 2-component vector pointing down (0, -1).
 */
KINLINE vec2 vec2_down (void) { return (vec2){0.0f, -1.0f}; }

/**
 * @brief Creates and returns a 2-component vector pointing left (-1, 0).
 */
KINLINE vec2 vec2_left (void) { return (vec2){-1.0f, 0.0f}; }

/**
 * @brief Creates and returns a 2-component vector pointing right (1, 0).
 */
KINLINE vec2 vec2_right (void) { return (vec2){1.0f, 0.0f}; }

/**
 * @brief Adds vector_1 to vector_0 and returns a copy of the result.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The resulting vector.
 */
KINLINE vec2 vec2_add (vec2 vector_0, vec2 vector_1) {
	return (vec2){vector_0.x + vector_1.x, vector_0.y + vector_1.y};
}

/**
 * @brief Subtracts vector_1 from vector_0 and returns a copy of the result.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The resulting vector.
 */
KINLINE vec2 vec2_sub (vec2 vector_0, vec2 vector_1) {
	return (vec2){vector_0.x - vector_1.x, vector_0.y - vector_1.y};
}

/**
 * @brief Multiplies vector_0 by vector_1 and returns a copy of the result.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The resulting vector.
 */
KINLINE vec2 vec2_mul (vec2 vector_0, vec2 vector_1) {
	return (vec2){vector_0.x * vector_1.x, vector_0.y * vector_1.y};
}

/**
 * @brief Multiplies all elements of vector_0 by scalar and returns a copy of
 * the result.
 *
 * @param vector_0 The vector to be multiplied.
 * @param scalar The scalar value.
 * @return A copy of the resulting vector.
 */
KINLINE vec2 vec2_mul_scalar (vec2 vector_0, f32 scalar) {
	return (vec2){vector_0.x * scalar, vector_0.y * scalar};
}

/**
 * @brief Multiplies vector_0 by vector_1, then adds the result to vector_2.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @param vector_2 The third vector.
 * @return The resulting vector.
 */
KINLINE vec2 vec2_mul_add (vec2 vector_0, vec2 vector_1, vec2 vector_2) {
	return (vec2){
		vector_0.x * vector_1.x + vector_2.x,
		vector_0.y * vector_1.y + vector_2.y};
}

/**
 * @brief Divides vector_0 by vector_1 and returns a copy of the result.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The resulting vector.
 */
KINLINE vec2 vec2_div (vec2 vector_0, vec2 vector_1) {
	return (vec2){vector_0.x / vector_1.x, vector_0.y / vector_1.y};
}

/**
 * @brief Returns the squared length of the provided vector.
 *
 * @param vector The vector to retrieve the squared length of.
 * @return The squared length.
 */
KINLINE f32 vec2_length_squared (vec2 vector) {
	return vector.x * vector.x + vector.y * vector.y;
}

/**
 * @brief Returns the length of the provided vector.
 *
 * @param vector The vector to retrieve the length of.
 * @return The length.
 */
KINLINE f32 vec2_length (vec2 vector) {
	return ksqrt(vec2_length_squared(vector));
}

/**
 * @brief Returns the square root of the sum of the vector's components.
 * Often used to find the hypotenuse of a right triangle. This function
 * has additional guards in place to safely prevent numeric under/overflow
 * as opposed to just using length().
 *
 * @param v The vector whose components are to be used.
 * @return The square root of the sum of the vector's components.
 */
KAPI f32 vec2_hypot (vec2 v);

/**
 * @brief Normalizes the provided vector in place to a unit vector.
 *
 * @param vector A pointer to the vector to be normalized.
 */
KINLINE void vec2_normalize (vec2 *vector) {
	const f32 length = vec2_length(*vector);
	vector->x /= length;
	vector->y /= length;
}

/**
 * @brief Returns a normalized copy of the supplied vector.
 *
 * @param vector The vector to be normalized.
 * @return A normalized copy of the supplied vector
 */
KINLINE vec2 vec2_normalized (vec2 vector) {
	vec2_normalize(&vector);
	return vector;
}

/**
 * @brief Compares all elements of vector_0 and vector_1 and ensures the
 * difference is less than tolerance.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @param tolerance The difference tolerance. Typically K_FLOAT_EPSILON or
 * similar.
 * @return True if within tolerance; otherwise false.
 */
KINLINE b8 vec2_compare (vec2 vector_0, vec2 vector_1, f32 tolerance) {
	if (kabs(vector_0.x - vector_1.x) > tolerance) {
		return false;
	}

	if (kabs(vector_0.y - vector_1.y) > tolerance) {
		return false;
	}

	return true;
}

/**
 * @brief Clamps the provided vector in-place to the given min/max values.
 *
 * @param vector A pointer to the vector to be clamped.
 * @param min The minimum value.
 * @param max The maximum value.
 */
KINLINE void vec2_clamp (vec2 *vector, vec2 min, vec2 max) {
	if (vector) {
		for (u8 i = 0; i < 2; ++i) {
			vector->elements[i] = KCLAMP(vector->elements[i], min.elements[i], max.elements[i]);
		}
	}
}

/**
 * @brief Returns a clamped copy of the provided vector.
 *
 * @param vector The vector to clamp.
 * @param min The minimum value.
 * @param max The maximum value.
 * @return A clamped copy of the provided vector.
 */
KINLINE vec2 vec2_clamped (vec2 vector, vec2 min, vec2 max) {
	vec2_clamp(&vector, min, max);
	return vector;
}

/**
 * @brief Clamps the provided vector in-place to the given min/max values.
 *
 * @param vector A pointer to the vector to be clamped.
 * @param min The minimum value.
 * @param max The maximum value.
 */
KINLINE void vec2_clamp_scalar (vec2 *vector, f32 min, f32 max) {
	if (vector) {
		for (u8 i = 0; i < 2; ++i) {
			vector->elements[i] = KCLAMP(vector->elements[i], min, max);
		}
	}
}

/**
 * @brief Returns a clamped copy of the provided vector.
 *
 * @param vector The vector to clamp.
 * @param min The minimum value.
 * @param max The maximum value.
 * @return A clamped copy of the provided vector.
 */
KINLINE vec2 vec2_clamped_scalar (vec2 vector, f32 min, f32 max) {
	vec2_clamp_scalar(&vector, min, max);
	return vector;
}

/**
 * @brief Returns the distance between vector_0 and vector_1.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The distance between vector_0 and vector_1.
 */
KINLINE f32 vec2_distance (vec2 vector_0, vec2 vector_1) {
	vec2 d = (vec2){vector_0.x - vector_1.x, vector_0.y - vector_1.y};
	return vec2_length(d);
}

/**
 * @brief Returns the squared distance between vector_0 and vector_1.
 * NOTE: If purely for comparison purposes, prefer this over non-squared version to avoid a sqrt call.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The distance between vector_0 and vector_1.
 */
KINLINE f32 vec2_distance_squared (vec2 vector_0, vec2 vector_1) {
	vec2 d = (vec2){vector_0.x - vector_1.x, vector_0.y - vector_1.y};
	return vec2_length_squared(d);
}

KINLINE vec2 vec2_min (vec2 vector_0, vec2 vector_1) {
	return vec2_create(
		KMIN(vector_0.x, vector_1.x),
		KMIN(vector_0.y, vector_1.y));
}

KINLINE vec2 vec2_max (vec2 vector_0, vec2 vector_1) {
	return vec2_create(
		KMAX(vector_0.x, vector_1.x),
		KMAX(vector_0.y, vector_1.y));
}

// ------------------------------------------
// Vector 3
// ------------------------------------------

/**
 * @brief Creates and returns a new 3-element vector using the supplied values.
 *
 * @param x The x value.
 * @param y The y value.
 * @param z The z value.
 * @return A new 3-element vector.
 */
#define vec3_create(x, y, z) \
	(vec3) { x, y, z }

/**
 * @brief Creates and returns a new 3-element vector using the supplied scalar for all components.
 *
 * @param scalar The scalar value.
 * @return A new 3-element vector.
 */
KINLINE vec3 vec3_from_scalar (f32 scalar) { return (vec3){scalar, scalar, scalar}; }

/*
 * @brief Returns a new vec3 containing the x, y and z components of the
 * supplied vec4, essentially dropping the w component.
 *
 * @param vector The 4-component vector to extract from.
 * @return A new vec3
 */
KINLINE vec3 vec3_from_vec4 (vec4 vector) {
	return (vec3){vector.x, vector.y, vector.z};
}

/*
 * @brief Returns a new vec3 containing the x and y components of the
 * supplied vec2, with a z component specified.
 *
 * @param vector The 2-component vector to extract from.
 * @param z The value to use for the z element.
 * @return A new vec3
 */
KINLINE vec3 vec3_from_vec2 (vec2 vector, f32 z) {
	return (vec3){vector.x, vector.y, z};
}

/**
 * @brief Returns a new vec4 using vector as the x, y and z components and w for
 * w.
 *
 * @param vector The 3-component vector.
 * @param w The w component.
 * @return A new vec4
 */
KINLINE vec4 vec3_to_vec4 (vec3 vector, f32 w) {
	return (vec4){vector.x, vector.y, vector.z, w};
}

/**
 * @brief Creates and returns a 3-component vector with all components set to
 * 0.0f.
 */
#define vec3_zero() \
	(vec3){0.0f, 0.0f, 0.0f}

/**
 * @brief Creates and returns a 3-component vector with all components set
 * to 1.0f.
 */
KINLINE vec3 vec3_one (void) { return (vec3){1.0f, 1.0f, 1.0f}; }

/**
 * @brief Creates and returns a 3-component vector pointing up (0, 1, 0).
 */
KINLINE vec3 vec3_up (void) { return (vec3){0.0f, 1.0f, 0.0f}; }

/**
 * @brief Creates and returns a 3-component vector pointing down (0, -1, 0).
 */
KINLINE vec3 vec3_down (void) { return (vec3){0.0f, -1.0f, 0.0f}; }

/**
 * @brief Creates and returns a 3-component vector pointing left (-1, 0, 0).
 */
KINLINE vec3 vec3_left (void) { return (vec3){-1.0f, 0.0f, 0.0f}; }

/**
 * @brief Creates and returns a 3-component vector pointing right (1, 0, 0).
 */
KINLINE vec3 vec3_right (void) { return (vec3){1.0f, 0.0f, 0.0f}; }

/**
 * @brief Creates and returns a 3-component vector pointing forward (0, 0, -1).
 */
KINLINE vec3 vec3_forward (void) { return (vec3){0.0f, 0.0f, -1.0f}; }

/**
 * @brief Creates and returns a 3-component vector pointing backward (0, 0, 1).
 */
KINLINE vec3 vec3_backward (void) { return (vec3){0.0f, 0.0f, 1.0f}; }

/**
 * @brief Adds vector_1 to vector_0 and returns a copy of the result.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The resulting vector.
 */
KINLINE vec3 vec3_add (vec3 vector_0, vec3 vector_1) {
	return (vec3){vector_0.x + vector_1.x, vector_0.y + vector_1.y,
				  vector_0.z + vector_1.z};
}

/**
 * @brief Subtracts vector_1 from vector_0 and returns a copy of the result.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The resulting vector.
 */
KINLINE vec3 vec3_sub (vec3 vector_0, vec3 vector_1) {
	return (vec3){vector_0.x - vector_1.x, vector_0.y - vector_1.y,
				  vector_0.z - vector_1.z};
}

/**
 * @brief Multiplies vector_0 by vector_1 and returns a copy of the result.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The resulting vector.
 */
KINLINE vec3 vec3_mul (vec3 vector_0, vec3 vector_1) {
	return (vec3){vector_0.x * vector_1.x, vector_0.y * vector_1.y,
				  vector_0.z * vector_1.z};
}

/**
 * @brief Multiplies all elements of vector_0 by scalar and returns a copy of
 * the result.
 *
 * @param vector_0 The vector to be multiplied.
 * @param scalar The scalar value.
 * @return A copy of the resulting vector.
 */
KINLINE vec3 vec3_mul_scalar (vec3 vector_0, f32 scalar) {
	return (vec3){vector_0.x * scalar, vector_0.y * scalar, vector_0.z * scalar};
}

/**
 * @brief Multiplies vector_0 by vector_1, then adds the result to vector_2.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @param vector_2 The third vector.
 * @return The resulting vector.
 */
KINLINE vec3 vec3_mul_add (vec3 vector_0, vec3 vector_1, vec3 vector_2) {
	return (vec3){
		vector_0.x * vector_1.x + vector_2.x,
		vector_0.y * vector_1.y + vector_2.y,
		vector_0.z * vector_1.z + vector_2.z};
}

/**
 * @brief Divides vector_0 by vector_1 and returns a copy of the result.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The resulting vector.
 */
KINLINE vec3 vec3_div (vec3 vector_0, vec3 vector_1) {
	return (vec3){vector_0.x / vector_1.x, vector_0.y / vector_1.y,
				  vector_0.z / vector_1.z};
}

KINLINE vec3 vec3_div_scalar (vec3 vector_0, f32 scalar) {
	vec3 result;
	for (u64 i = 0; i < 3; ++i) {
		result.elements[i] = vector_0.elements[i] / scalar;
	}
	return result;
}

/**
 * @brief Returns the squared length of the provided vector.
 *
 * @param vector The vector to retrieve the squared length of.
 * @return The squared length.
 */
KINLINE f32 vec3_length_squared (vec3 vector) {
	return vector.x * vector.x + vector.y * vector.y + vector.z * vector.z;
}

/**
 * @brief Returns the length of the provided vector.
 *
 * @param vector The vector to retrieve the length of.
 * @return The length.
 */
KINLINE f32 vec3_length (vec3 vector) {
	return ksqrt(vec3_length_squared(vector));
}

/**
 * @brief Returns the square root of the sum of the vector's components.
 * Often used to find the hypotenuse of a right triangle. This function
 * has additional guards in place to safely prevent numeric under/overflow
 * as opposed to just using length().
 *
 * @param v The vector whose components are to be used.
 * @return The square root of the sum of the vector's components.
 */
KAPI f32 vec3_hypot (vec3 v);

/**
 * @brief Normalizes the provided vector in place to a unit vector.
 *
 * @param vector A pointer to the vector to be normalized.
 */
KINLINE void vec3_normalize (vec3 *vector) {
	const f32 length = vec3_length(*vector);
	vector->x /= length;
	vector->y /= length;
	vector->z /= length;
}

/**
 * @brief Returns a normalized copy of the supplied vector.
 *
 * @param vector The vector to be normalized.
 * @return A normalized copy of the supplied vector
 */
KINLINE vec3 vec3_normalized (vec3 vector) {
	vec3_normalize(&vector);
	return vector;
}

/**
 * @brief Returns the dot product between the provided vectors. Typically used
 * to calculate the difference in direction.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The dot product.
 */
KINLINE f32 vec3_dot (vec3 vector_0, vec3 vector_1) {
	f32 p = 0;
	p += vector_0.x * vector_1.x;
	p += vector_0.y * vector_1.y;
	p += vector_0.z * vector_1.z;
	return p;
}

/**
 * @brief Calculates and returns the cross product of the supplied vectors.
 * The cross product is a new vector which is orthoganal to both provided
 * vectors.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The cross product.
 */
KINLINE vec3 vec3_cross (vec3 vector_0, vec3 vector_1) {
	return (vec3){vector_0.y * vector_1.z - vector_0.z * vector_1.y,
				  vector_0.z * vector_1.x - vector_0.x * vector_1.z,
				  vector_0.x * vector_1.y - vector_0.y * vector_1.x};
}

/**
 * @brief Compares all elements of vector_0 and vector_1 and ensures the
 * difference is less than tolerance.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @param tolerance The difference tolerance. Typically K_FLOAT_EPSILON or
 * similar.
 * @return True if within tolerance; otherwise false.
 */
KINLINE b8 vec3_compare (vec3 vector_0, vec3 vector_1, f32 tolerance) {
	if (kabs(vector_0.x - vector_1.x) > tolerance) {
		return false;
	}

	if (kabs(vector_0.y - vector_1.y) > tolerance) {
		return false;
	}

	if (kabs(vector_0.z - vector_1.z) > tolerance) {
		return false;
	}

	return true;
}

/**
 * @brief Clamps the provided vector in-place to the given min/max values.
 *
 * @param vector A pointer to the vector to be clamped.
 * @param min The minimum value.
 * @param max The maximum value.
 */
KINLINE void vec3_clamp (vec3 *vector, vec3 min, vec3 max) {
	if (vector) {
		for (u8 i = 0; i < 3; ++i) {
			vector->elements[i] = KCLAMP(vector->elements[i], min.elements[i], max.elements[i]);
		}
	}
}

/**
 * @brief Returns a clamped copy of the provided vector.
 *
 * @param vector The vector to clamp.
 * @param min The minimum value.
 * @param max The maximum value.
 * @return A clamped copy of the provided vector.
 */
KINLINE vec3 vec3_clamped (vec3 vector, vec3 min, vec3 max) {
	vec3_clamp(&vector, min, max);
	return vector;
}

/**
 * @brief Clamps the provided vector in-place to the given min/max values.
 *
 * @param vector A pointer to the vector to be clamped.
 * @param min The minimum value.
 * @param max The maximum value.
 */
KINLINE void vec3_clamp_scalar (vec3 *vector, f32 min, f32 max) {
	if (vector) {
		for (u8 i = 0; i < 3; ++i) {
			vector->elements[i] = KCLAMP(vector->elements[i], min, max);
		}
	}
}

/**
 * @brief Returns a clamped copy of the provided vector.
 *
 * @param vector The vector to clamp.
 * @param min The minimum value.
 * @param max The maximum value.
 * @return A clamped copy of the provided vector.
 */
KINLINE vec3 vec3_clamped_scalar (vec3 vector, f32 min, f32 max) {
	vec3_clamp_scalar(&vector, min, max);
	return vector;
}

/**
 * @brief Returns the distance between vector_0 and vector_1.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The distance between vector_0 and vector_1.
 */
KINLINE f32 vec3_distance (vec3 vector_0, vec3 vector_1) {
	vec3 d = (vec3){vector_0.x - vector_1.x, vector_0.y - vector_1.y,
					vector_0.z - vector_1.z};
	return vec3_length(d);
}

/**
 * @brief Returns the squared distance between vector_0 and vector_1.
 * Less intensive than calling the non-squared version due to sqrt.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The squared distance between vector_0 and vector_1.
 */
KINLINE f32 vec3_distance_squared (vec3 vector_0, vec3 vector_1) {
	vec3 d = (vec3){vector_0.x - vector_1.x, vector_0.y - vector_1.y,
					vector_0.z - vector_1.z};
	return vec3_length_squared(d);
}

/**
 * @brief Projects v_0 onto v_1.
 *
 * @param v_0 The first vector.
 * @param v_1 The second vector.
 * @return The projected vector.
 */
KINLINE vec3 vec3_project (vec3 v_0, vec3 v_1) {
	f32 length_sq = vec3_length_squared(v_1);
	if (length_sq == 0.0f) {
		// NOTE: handle divide-by-zero case (i.e. v_1 is a zero vector).
		return vec3_zero();
	}
	f32 scalar = vec3_dot(v_0, v_1) / length_sq;
	return vec3_mul_scalar(v_1, scalar);
}

/**
 * @brief Transform v by m.
 * @param v The vector to transform.
 * @param w Pass 1.0f for a point, or 0.0f for a direction.
 * @param m The matrix to transform by.
 * @return A transformed copy of v.
 */
KINLINE vec3 vec3_transform (vec3 v, f32 w, mat4 m) {
	vec3 out;
	out.x = v.x * m.data[0 + 0] + v.y * m.data[4 + 0] + v.z * m.data[8 + 0] + w * m.data[12 + 0];
	out.y = v.x * m.data[0 + 1] + v.y * m.data[4 + 1] + v.z * m.data[8 + 1] + w * m.data[12 + 1];
	out.z = v.x * m.data[0 + 2] + v.y * m.data[4 + 2] + v.z * m.data[8 + 2] + w * m.data[12 + 2];
	return out;
}

KINLINE vec3 vec3_min (vec3 vector_0, vec3 vector_1) {
	return vec3_create(
		KMIN(vector_0.x, vector_1.x),
		KMIN(vector_0.y, vector_1.y),
		KMIN(vector_0.z, vector_1.z));
}

KINLINE vec3 vec3_max (vec3 vector_0, vec3 vector_1) {
	return vec3_create(
		KMAX(vector_0.x, vector_1.x),
		KMAX(vector_0.y, vector_1.y),
		KMAX(vector_0.z, vector_1.z));
}

KINLINE vec3 vec3_sign (vec3 v) {
	return vec3_create(ksign(v.x), ksign(v.y), ksign(v.z));
}

KINLINE vec3 vec3_rotate (vec3 v, quat q) {
	vec3 u = vec3_create(q.x, q.y, q.z);
	f32 s = q.w;

	return vec3_add(
		vec3_add(
			vec3_mul_scalar(u, 2.0f * vec3_dot(u, v)),
			vec3_mul_scalar(v, (s * s - vec3_dot(u, u)))),
		vec3_mul_scalar(vec3_cross(u, v), 2.0f * s));
}

// ------------------------------------------
// Vector 4
// ------------------------------------------

/**
 * @brief Creates and returns a new 4-element vector using the supplied values.
 *
 * @param x The x value.
 * @param y The y value.
 * @param z The z value.
 * @param w The w value.
 * @return A new 4-element vector.
 */
KINLINE vec4 vec4_create (f32 x, f32 y, f32 z, f32 w) {
	vec4 out_vector;
#if defined(KUSE_SIMD)
	out_vector.data = _mm_setr_ps(x, y, z, w);
#else
	out_vector.x = x;
	out_vector.y = y;
	out_vector.z = z;
	out_vector.w = w;
#endif
	return out_vector;
}

/**
 * @brief Creates and returns a new 4-element vector using the supplied scalar for all components.
 *
 * @param scalar The scalar value.
 * @return A new 4-element vector.
 */
KINLINE vec4 vec4_from_scalar (f32 scalar) { return (vec4){scalar, scalar, scalar, scalar}; }

/**
 * @brief Returns a new vec3 containing the x, y and z components of the
 * supplied vec4, essentially dropping the w component.
 *
 * @param vector The 4-component vector to extract from.
 * @return A new vec3
 */
KINLINE vec3 vec4_to_vec3 (vec4 vector) {
	return (vec3){vector.x, vector.y, vector.z};
}

/**
 * @brief Returns a new vec4 using vector as the x, y and z components and w for
 * w.
 *
 * @param vector The 3-component vector.
 * @param w The w component.
 * @return A new vec4
 */
KINLINE vec4 vec4_from_vec3 (vec3 vector, f32 w) {
#if defined(KUSE_SIMD)
	vec4 out_vector;
	out_vector.data = _mm_setr_ps(x, y, z, w);
	return out_vector;
#else
	return (vec4){vector.x, vector.y, vector.z, w};
#endif
}

/**
 * @brief Creates and returns a 4-component vector with all components set to
 * 0.0f.
 */
KINLINE vec4 vec4_zero (void) { return (vec4){0.0f, 0.0f, 0.0f, 0.0f}; }

/**
 * @brief Creates and returns a 4-component vector with all components set
 * to 1.0f.
 */
KINLINE vec4 vec4_one (void) { return (vec4){1.0f, 1.0f, 1.0f, 1.0f}; }

/**
 * @brief Adds vector_1 to vector_0 and returns a copy of the result.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The resulting vector.
 */
KINLINE vec4 vec4_add (vec4 vector_0, vec4 vector_1) {
	vec4 result;
	for (u64 i = 0; i < 4; ++i) {
		result.elements[i] = vector_0.elements[i] + vector_1.elements[i];
	}
	return result;
}

/**
 * @brief Subtracts vector_1 from vector_0 and returns a copy of the result.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The resulting vector.
 */
KINLINE vec4 vec4_sub (vec4 vector_0, vec4 vector_1) {
	vec4 result;
	for (u64 i = 0; i < 4; ++i) {
		result.elements[i] = vector_0.elements[i] - vector_1.elements[i];
	}
	return result;
}

/**
 * @brief Multiplies vector_0 by vector_1 and returns a copy of the result.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The resulting vector.
 */
KINLINE vec4 vec4_mul (vec4 vector_0, vec4 vector_1) {
	vec4 result;
	for (u64 i = 0; i < 4; ++i) {
		result.elements[i] = vector_0.elements[i] * vector_1.elements[i];
	}
	return result;
}

/**
 * @brief Multiplies all elements of vector_0 by scalar and returns a copy of
 * the result.
 *
 * @param vector_0 The vector to be multiplied.
 * @param scalar The scalar value.
 * @return A copy of the resulting vector.
 */
KINLINE vec4 vec4_mul_scalar (vec4 vector_0, f32 scalar) {
	return (vec4){vector_0.x * scalar, vector_0.y * scalar, vector_0.z * scalar, vector_0.w * scalar};
}

/**
 * @brief Multiplies vector_0 by vector_1, then adds the result to vector_2.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @param vector_2 The third vector.
 * @return The resulting vector.
 */
KINLINE vec4 vec4_mul_add (vec4 vector_0, vec4 vector_1, vec4 vector_2) {
	return (vec4){
		vector_0.x * vector_1.x + vector_2.x,
		vector_0.y * vector_1.y + vector_2.y,
		vector_0.z * vector_1.z + vector_2.z,
		vector_0.w * vector_1.w + vector_2.w,
	};
}

/**
 * @brief Divides vector_0 by vector_1 and returns a copy of the result.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @return The resulting vector.
 */
KINLINE vec4 vec4_div (vec4 vector_0, vec4 vector_1) {
	vec4 result;
	for (u64 i = 0; i < 4; ++i) {
		result.elements[i] = vector_0.elements[i] / vector_1.elements[i];
	}
	return result;
}

KINLINE vec4 vec4_div_scalar (vec4 vector_0, f32 scalar) {
	vec4 result;
	for (u64 i = 0; i < 4; ++i) {
		result.elements[i] = vector_0.elements[i] / scalar;
	}
	return result;
}

/**
 * @brief Returns the squared length of the provided vector.
 *
 * @param vector The vector to retrieve the squared length of.
 * @return The squared length.
 */
KINLINE f32 vec4_length_squared (vec4 vector) {
	return vector.x * vector.x + vector.y * vector.y + vector.z * vector.z +
		   vector.w * vector.w;
}

/**
 * @brief Returns the length of the provided vector.
 *
 * @param vector The vector to retrieve the length of.
 * @return The length.
 */
KINLINE f32 vec4_length (vec4 vector) {
	return ksqrt(vec4_length_squared(vector));
}

/**
 * @brief Normalizes the provided vector in place to a unit vector.
 *
 * @param vector A pointer to the vector to be normalized.
 */
KINLINE void vec4_normalize (vec4 *vector) {
	const f32 length = vec4_length(*vector);
	vector->x /= length;
	vector->y /= length;
	vector->z /= length;
	vector->w /= length;
}

/**
 * @brief Returns a normalized copy of the supplied vector.
 *
 * @param vector The vector to be normalized.
 * @return A normalized copy of the supplied vector
 */
KINLINE vec4 vec4_normalized (vec4 vector) {
	vec4_normalize(&vector);
	return vector;
}

/**
 * @brief Calculates the dot product using the elements of vec4s provided in
 * split-out format.
 *
 * @param a0 The first element of the a vector.
 * @param a1 The second element of the a vector.
 * @param a2 The third element of the a vector.
 * @param a3 The fourth element of the a vector.
 * @param b0 The first element of the b vector.
 * @param b1 The second element of the b vector.
 * @param b2 The third element of the b vector.
 * @param b3 The fourth element of the b vector.
 * @return The dot product of vectors and b.
 */
KINLINE f32 vec4_dot_f32 (f32 a0, f32 a1, f32 a2, f32 a3, f32 b0, f32 b1, f32 b2,
						  f32 b3) {
	f32 p;
	p = a0 * b0 + a1 * b1 + a2 * b2 + a3 * b3;
	return p;
}

/**
 * @brief Compares all elements of vector_0 and vector_1 and ensures the
 * difference is less than tolerance.
 *
 * @param vector_0 The first vector.
 * @param vector_1 The second vector.
 * @param tolerance The difference tolerance. Typically K_FLOAT_EPSILON or
 * similar.
 * @return True if within tolerance; otherwise false.
 */
KINLINE b8 vec4_compare (vec4 vector_0, vec4 vector_1, f32 tolerance) {
	if (kabs(vector_0.x - vector_1.x) > tolerance) {
		return false;
	}

	if (kabs(vector_0.y - vector_1.y) > tolerance) {
		return false;
	}

	if (kabs(vector_0.z - vector_1.z) > tolerance) {
		return false;
	}

	if (kabs(vector_0.w - vector_1.w) > tolerance) {
		return false;
	}

	return true;
}

/**
 * @brief Clamps the provided vector in-place to the given min/max values.
 *
 * @param vector A pointer to the vector to be clamped.
 * @param min The minimum value.
 * @param max The maximum value.
 */
KINLINE void vec4_clamp (vec4 *vector, vec4 min, vec4 max) {
	if (vector) {
		for (u8 i = 0; i < 4; ++i) {
			vector->elements[i] = KCLAMP(vector->elements[i], min.elements[i], max.elements[i]);
		}
	}
}

/**
 * @brief Returns a clamped copy of the provided vector.
 *
 * @param vector The vector to clamp.
 * @param min The minimum value.
 * @param max The maximum value.
 * @return A clamped copy of the provided vector.
 */
KINLINE vec4 vec4_clamped (vec4 vector, vec4 min, vec4 max) {
	vec4_clamp(&vector, min, max);
	return vector;
}

/**
 * @brief Clamps the provided vector in-place to the given min/max values.
 *
 * @param vector A pointer to the vector to be clamped.
 * @param min The minimum value.
 * @param max The maximum value.
 */
KINLINE void vec4_clamp_scalar (vec4 *vector, f32 min, f32 max) {
	if (vector) {
		for (u8 i = 0; i < 4; ++i) {
			vector->elements[i] = KCLAMP(vector->elements[i], min, max);
		}
	}
}

/**
 * @brief Returns a clamped copy of the provided vector.
 *
 * @param vector The vector to clamp.
 * @param min The minimum value.
 * @param max The maximum value.
 * @return A clamped copy of the provided vector.
 */
KINLINE vec4 vec4_clamped_scalar (vec4 vector, f32 min, f32 max) {
	vec4_clamp_scalar(&vector, min, max);
	return vector;
}

/**
 * @brief Creates and returns an identity matrix:
 *
 * {
 *   {1, 0, 0, 0},
 *   {0, 1, 0, 0},
 *   {0, 0, 1, 0},
 *   {0, 0, 0, 1}
 * }
 *
 * @return A new identity matrix
 */
KINLINE mat4 mat4_identity (void) {
	mat4 out_matrix;
	kzero_memory(out_matrix.data, sizeof(f32) * 16);
	out_matrix.data[0] = 1.0f;
	out_matrix.data[5] = 1.0f;
	out_matrix.data[10] = 1.0f;
	out_matrix.data[15] = 1.0f;
	return out_matrix;
}

/**
 * @brief Returns the result of multiplying matrix_0 and matrix_1.
 *
 * @param matrix_0 The first matrix to be multiplied.
 * @param matrix_1 The second matrix to be multiplied.
 * @return The result of the matrix multiplication.
 */
KINLINE mat4 mat4_mul (mat4 matrix_0, mat4 matrix_1) {
	mat4 out_matrix = mat4_identity();

	const f32 *m1_ptr = matrix_0.data;
	const f32 *m2_ptr = matrix_1.data;
	f32 *dst_ptr = out_matrix.data;

	for (i32 i = 0; i < 4; ++i) {
		for (i32 j = 0; j < 4; ++j) {
			*dst_ptr = m1_ptr[0] * m2_ptr[0 + j] + m1_ptr[1] * m2_ptr[4 + j] +
					   m1_ptr[2] * m2_ptr[8 + j] + m1_ptr[3] * m2_ptr[12 + j];
			dst_ptr++;
		}
		m1_ptr += 4;
	}
	return out_matrix;
}

/**
 * @brief Creates and returns an orthographic projection matrix. Typically used
 * to render flat or 2D scenes.
 *
 * @param left The left side of the view frustum.
 * @param right The right side of the view frustum.
 * @param bottom The bottom side of the view frustum.
 * @param top The top side of the view frustum.
 * @param near_clip The near clipping plane distance.
 * @param far_clip The far clipping plane distance.
 * @return A new orthographic projection matrix.
 */
KINLINE mat4 mat4_orthographic (f32 left, f32 right, f32 bottom, f32 top,
								f32 near_clip, f32 far_clip) {
	mat4 out_matrix = mat4_identity();

	f32 lr = 1.0f / (left - right);
	f32 bt = 1.0f / (bottom - top);
	f32 nf = 1.0f / (near_clip - far_clip);

	out_matrix.data[0] = -2.0f * lr;
	out_matrix.data[5] = -2.0f * bt;
	out_matrix.data[10] = nf;

	out_matrix.data[12] = (left + right) * lr;
	out_matrix.data[13] = (top + bottom) * bt;
	out_matrix.data[14] = -near_clip * nf;

	return out_matrix;
}

/**
 * @brief Creates and returns a perspective matrix. Typically used to render 3d
 * scenes.
 *
 * @param fov_radians The field of view in radians.
 * @param aspect_ratio The aspect ratio.
 * @param near_clip The near clipping plane distance.
 * @param far_clip The far clipping plane distance.
 * @return A new perspective matrix.
 */
KINLINE mat4 mat4_perspective (f32 fov_radians, f32 aspect_ratio, f32 near_clip, f32 far_clip) {
	f32 half_tan_fov = ktan(fov_radians * 0.5f);
	mat4 out_matrix;
	kzero_memory(out_matrix.data, sizeof(f32) * 16);
	out_matrix.data[0] = 1.0f / (aspect_ratio * half_tan_fov);
	out_matrix.data[5] = 1.0f / half_tan_fov;
	out_matrix.data[10] = far_clip / (near_clip - far_clip);
	out_matrix.data[11] = -1.0f;
	out_matrix.data[14] = (far_clip * near_clip) / (near_clip - far_clip);
	return out_matrix;
}

/**
 * @brief Creates and returns a look-at matrix, or a matrix looking
 * at target from the perspective of position.
 *
 * @param position The position of the matrix.
 * @param target The position to "look at".
 * @param up The up vector.
 * @return A matrix looking at target from the perspective of position.
 */
KINLINE mat4 mat4_look_at (vec3 position, vec3 target, vec3 up) {
	mat4 out_matrix;
	vec3 z_axis = vec3_normalized(vec3_sub(target, position));
	vec3 x_axis = vec3_normalized(vec3_cross(z_axis, up));
	vec3 y_axis = vec3_cross(x_axis, z_axis);

	out_matrix.data[0] = x_axis.x;
	out_matrix.data[1] = y_axis.x;
	out_matrix.data[2] = -z_axis.x;
	out_matrix.data[3] = 0;
	out_matrix.data[4] = x_axis.y;
	out_matrix.data[5] = y_axis.y;
	out_matrix.data[6] = -z_axis.y;
	out_matrix.data[7] = 0;
	out_matrix.data[8] = x_axis.z;
	out_matrix.data[9] = y_axis.z;
	out_matrix.data[10] = -z_axis.z;
	out_matrix.data[11] = 0;
	out_matrix.data[12] = -vec3_dot(x_axis, position);
	out_matrix.data[13] = -vec3_dot(y_axis, position);
	out_matrix.data[14] = vec3_dot(z_axis, position);
	out_matrix.data[15] = 1.0f;
	return out_matrix;
}

/**
 * @brief Returns a transposed copy of the provided matrix (rows->colums)
 *
 * @param matrix The matrix to be transposed.
 * @return A transposed copy of of the provided matrix.
 */
KINLINE mat4 mat4_transposed (mat4 matrix) {
	mat4 out_matrix;
	out_matrix.data[0] = matrix.data[0];
	out_matrix.data[1] = matrix.data[4];
	out_matrix.data[2] = matrix.data[8];
	out_matrix.data[3] = matrix.data[12];
	out_matrix.data[4] = matrix.data[1];
	out_matrix.data[5] = matrix.data[5];
	out_matrix.data[6] = matrix.data[9];
	out_matrix.data[7] = matrix.data[13];
	out_matrix.data[8] = matrix.data[2];
	out_matrix.data[9] = matrix.data[6];
	out_matrix.data[10] = matrix.data[10];
	out_matrix.data[11] = matrix.data[14];
	out_matrix.data[12] = matrix.data[3];
	out_matrix.data[13] = matrix.data[7];
	out_matrix.data[14] = matrix.data[11];
	out_matrix.data[15] = matrix.data[15];
	return out_matrix;
}

/**
 * @brief Calculates the determinant of the given matrix.
 *
 * @param matrix The matrix to calculate the determinant of.
 * @return The determinant of the given matrix.
 */
KINLINE f32 mat4_determinant (mat4 matrix) {
	const f32 *m = matrix.data;

	f32 t0 = m[10] * m[15];
	f32 t1 = m[14] * m[11];
	f32 t2 = m[6] * m[15];
	f32 t3 = m[14] * m[7];
	f32 t4 = m[6] * m[11];
	f32 t5 = m[10] * m[7];
	f32 t6 = m[2] * m[15];
	f32 t7 = m[14] * m[3];
	f32 t8 = m[2] * m[11];
	f32 t9 = m[10] * m[3];
	f32 t10 = m[2] * m[7];
	f32 t11 = m[6] * m[3];

	mat3 temp_mat;
	f32 *o = temp_mat.data;

	o[0] = (t0 * m[5] + t3 * m[9] + t4 * m[13]) -
		   (t1 * m[5] + t2 * m[9] + t5 * m[13]);
	o[1] = (t1 * m[1] + t6 * m[9] + t9 * m[13]) -
		   (t0 * m[1] + t7 * m[9] + t8 * m[13]);
	o[2] = (t2 * m[1] + t7 * m[5] + t10 * m[13]) -
		   (t3 * m[1] + t6 * m[5] + t11 * m[13]);
	o[3] = (t5 * m[1] + t8 * m[5] + t11 * m[9]) -
		   (t4 * m[1] + t9 * m[5] + t10 * m[9]);

	f32 determinant = 1.0f / (m[0] * o[0] + m[4] * o[1] + m[8] * o[2] + m[12] * o[3]);
	return determinant;
}

/**
 * @brief Creates and returns an inverse of the provided matrix.
 *
 * @param matrix The matrix to be inverted.
 * @return A inverted copy of the provided matrix.
 */
KINLINE mat4 mat4_inverse (mat4 matrix) {
	const f32 *m = matrix.data;

	f32 t0 = m[10] * m[15];
	f32 t1 = m[14] * m[11];
	f32 t2 = m[6] * m[15];
	f32 t3 = m[14] * m[7];
	f32 t4 = m[6] * m[11];
	f32 t5 = m[10] * m[7];
	f32 t6 = m[2] * m[15];
	f32 t7 = m[14] * m[3];
	f32 t8 = m[2] * m[11];
	f32 t9 = m[10] * m[3];
	f32 t10 = m[2] * m[7];
	f32 t11 = m[6] * m[3];
	f32 t12 = m[8] * m[13];
	f32 t13 = m[12] * m[9];
	f32 t14 = m[4] * m[13];
	f32 t15 = m[12] * m[5];
	f32 t16 = m[4] * m[9];
	f32 t17 = m[8] * m[5];
	f32 t18 = m[0] * m[13];
	f32 t19 = m[12] * m[1];
	f32 t20 = m[0] * m[9];
	f32 t21 = m[8] * m[1];
	f32 t22 = m[0] * m[5];
	f32 t23 = m[4] * m[1];

	mat4 out_matrix;
	f32 *o = out_matrix.data;

	o[0] = (t0 * m[5] + t3 * m[9] + t4 * m[13]) -
		   (t1 * m[5] + t2 * m[9] + t5 * m[13]);
	o[1] = (t1 * m[1] + t6 * m[9] + t9 * m[13]) -
		   (t0 * m[1] + t7 * m[9] + t8 * m[13]);
	o[2] = (t2 * m[1] + t7 * m[5] + t10 * m[13]) -
		   (t3 * m[1] + t6 * m[5] + t11 * m[13]);
	o[3] = (t5 * m[1] + t8 * m[5] + t11 * m[9]) -
		   (t4 * m[1] + t9 * m[5] + t10 * m[9]);

	f32 d = 1.0f / (m[0] * o[0] + m[4] * o[1] + m[8] * o[2] + m[12] * o[3]);

	// Check for singular matrix (determinant near zero)
	if (kabs(d) < 1e-6f) {
		// Return identity matrix if the determinant is close to zero (singular matrix)
		return mat4_identity();
	}

	o[0] = d * o[0];
	o[1] = d * o[1];
	o[2] = d * o[2];
	o[3] = d * o[3];
	o[4] = d * ((t1 * m[4] + t2 * m[8] + t5 * m[12]) -
				(t0 * m[4] + t3 * m[8] + t4 * m[12]));
	o[5] = d * ((t0 * m[0] + t7 * m[8] + t8 * m[12]) -
				(t1 * m[0] + t6 * m[8] + t9 * m[12]));
	o[6] = d * ((t3 * m[0] + t6 * m[4] + t11 * m[12]) -
				(t2 * m[0] + t7 * m[4] + t10 * m[12]));
	o[7] = d * ((t4 * m[0] + t9 * m[4] + t10 * m[8]) -
				(t5 * m[0] + t8 * m[4] + t11 * m[8]));
	o[8] = d * ((t12 * m[7] + t15 * m[11] + t16 * m[15]) -
				(t13 * m[7] + t14 * m[11] + t17 * m[15]));
	o[9] = d * ((t13 * m[3] + t18 * m[11] + t21 * m[15]) -
				(t12 * m[3] + t19 * m[11] + t20 * m[15]));
	o[10] = d * ((t14 * m[3] + t19 * m[7] + t22 * m[15]) -
				 (t15 * m[3] + t18 * m[7] + t23 * m[15]));
	o[11] = d * ((t17 * m[3] + t20 * m[7] + t23 * m[11]) -
				 (t16 * m[3] + t21 * m[7] + t22 * m[11]));
	o[12] = d * ((t14 * m[10] + t17 * m[14] + t13 * m[6]) -
				 (t16 * m[14] + t12 * m[6] + t15 * m[10]));
	o[13] = d * ((t20 * m[14] + t12 * m[2] + t19 * m[10]) -
				 (t18 * m[10] + t21 * m[14] + t13 * m[2]));
	o[14] = d * ((t18 * m[6] + t23 * m[14] + t15 * m[2]) -
				 (t22 * m[14] + t14 * m[2] + t19 * m[6]));
	o[15] = d * ((t22 * m[10] + t16 * m[2] + t21 * m[6]) -
				 (t20 * m[6] + t23 * m[10] + t17 * m[2]));

	return out_matrix;
}

/**
 * @brief Creates and returns a translation matrix from the given position.
 *
 * @param position The position to be used to create the matrix.
 * @return A newly created translation matrix.
 */
KINLINE mat4 mat4_translation (vec3 position) {
	mat4 out_matrix = mat4_identity();
	out_matrix.data[12] = position.x;
	out_matrix.data[13] = position.y;
	out_matrix.data[14] = position.z;
	return out_matrix;
}

/**
 * @brief Returns a scale matrix using the provided scale.
 *
 * @param scale The 3-component scale.
 * @return A scale matrix.
 */
KINLINE mat4 mat4_scale (vec3 scale) {
	mat4 out_matrix = mat4_identity();
	out_matrix.data[0] = scale.x;
	out_matrix.data[5] = scale.y;
	out_matrix.data[10] = scale.z;
	return out_matrix;
}

/**
 * @brief Returns a matrix created from the provided translation, rotation and scale (TRS).
 *
 * @param position The position to be used to create the matrix.
 * @param rotation The quaternion rotation to be used to create the matrix.
 * @param scale The 3-component scale to be used to create the matrix.
 * @return A matrix created in TRS order.
 */
KINLINE mat4 mat4_from_translation_rotation_scale (vec3 t, quat r, vec3 s) {
	mat4 out_matrix;

	out_matrix.data[0] = (1.0f - 2.0f * (r.y * r.y + r.z * r.z)) * s.x;
	out_matrix.data[1] = (r.x * r.y + r.z * r.w) * s.x * 2.0f;
	out_matrix.data[2] = (r.x * r.z - r.y * r.w) * s.x * 2.0f;
	out_matrix.data[3] = 0.0f;
	out_matrix.data[4] = (r.x * r.y - r.z * r.w) * s.y * 2.0f;
	out_matrix.data[5] = (1.0f - 2.0f * (r.x * r.x + r.z * r.z)) * s.y;
	out_matrix.data[6] = (r.y * r.z + r.x * r.w) * s.y * 2.0f;
	out_matrix.data[7] = 0.0f;
	out_matrix.data[8] = (r.x * r.z + r.y * r.w) * s.z * 2.0f;
	out_matrix.data[9] = (r.y * r.z - r.x * r.w) * s.z * 2.0f;
	out_matrix.data[10] = (1.0f - 2.0f * (r.x * r.x + r.y * r.y)) * s.z;
	out_matrix.data[11] = 0.0f;
	out_matrix.data[12] = t.x;
	out_matrix.data[13] = t.y;
	out_matrix.data[14] = t.z;
	out_matrix.data[15] = 1.0f;

	return out_matrix;
}

/**
 * @brief Creates a rotation matrix from the provided x angle.
 *
 * @param angle_radians The x angle in radians.
 * @return A rotation matrix.
 */
KINLINE mat4 mat4_euler_x (f32 angle_radians) {
	mat4 out_matrix = mat4_identity();
	f32 c = kcos(angle_radians);
	f32 s = ksin(angle_radians);

	out_matrix.data[5] = c;
	out_matrix.data[6] = s;
	out_matrix.data[9] = -s;
	out_matrix.data[10] = c;
	return out_matrix;
}

/**
 * @brief Creates a rotation matrix from the provided y angle.
 *
 * @param angle_radians The y angle in radians.
 * @return A rotation matrix.
 */
KINLINE mat4 mat4_euler_y (f32 angle_radians) {
	mat4 out_matrix = mat4_identity();
	f32 c = kcos(angle_radians);
	f32 s = ksin(angle_radians);

	out_matrix.data[0] = c;
	out_matrix.data[2] = -s;
	out_matrix.data[8] = s;
	out_matrix.data[10] = c;
	return out_matrix;
}

/**
 * @brief Creates a rotation matrix from the provided z angle.
 *
 * @param angle_radians The z angle in radians.
 * @return A rotation matrix.
 */
KINLINE mat4 mat4_euler_z (f32 angle_radians) {
	mat4 out_matrix = mat4_identity();

	f32 c = kcos(angle_radians);
	f32 s = ksin(angle_radians);

	out_matrix.data[0] = c;
	out_matrix.data[1] = s;
	out_matrix.data[4] = -s;
	out_matrix.data[5] = c;
	return out_matrix;
}

/**
 * @brief Creates a rotation matrix from the provided x, y and z axis rotations.
 *
 * @param x_radians The x rotation.
 * @param y_radians The y rotation.
 * @param z_radians The z rotation.
 * @return A rotation matrix.
 */
KINLINE mat4 mat4_euler_xyz (f32 x_radians, f32 y_radians, f32 z_radians) {
	mat4 rx = mat4_euler_x(x_radians);
	mat4 ry = mat4_euler_y(y_radians);
	mat4 rz = mat4_euler_z(z_radians);
	mat4 out_matrix = mat4_mul(rx, ry);
	out_matrix = mat4_mul(out_matrix, rz);
	return out_matrix;
}

/**
 * @brief Returns a forward vector relative to the provided matrix.
 *
 * @param matrix The matrix from which to base the vector.
 * @return A 3-component directional vector.
 */
KINLINE vec3 mat4_forward (mat4 matrix) {
	vec3 forward;
	forward.x = -matrix.data[8];
	forward.y = -matrix.data[9];
	forward.z = -matrix.data[10];
	vec3_normalize(&forward);
	return forward;
}

/**
 * @brief Returns a backward vector relative to the provided matrix.
 *
 * @param matrix The matrix from which to base the vector.
 * @return A 3-component directional vector.
 */
KINLINE vec3 mat4_backward (mat4 matrix) {
	vec3 backward;
	backward.x = matrix.data[8];
	backward.y = matrix.data[9];
	backward.z = matrix.data[10];
	vec3_normalize(&backward);
	return backward;
}

/**
 * @brief Returns a upward vector relative to the provided matrix.
 *
 * @param matrix The matrix from which to base the vector.
 * @return A 3-component directional vector.
 */
KINLINE vec3 mat4_up (mat4 matrix) {
	vec3 up;
	up.x = matrix.data[1];
	up.y = matrix.data[5];
	up.z = matrix.data[9];
	vec3_normalize(&up);
	return up;
}

/**
 * @brief Returns a downward vector relative to the provided matrix.
 *
 * @param matrix The matrix from which to base the vector.
 * @return A 3-component directional vector.
 */
KINLINE vec3 mat4_down (mat4 matrix) {
	vec3 down;
	down.x = -matrix.data[1];
	down.y = -matrix.data[5];
	down.z = -matrix.data[9];
	vec3_normalize(&down);
	return down;
}

/**
 * @brief Returns a left vector relative to the provided matrix.
 *
 * @param matrix The matrix from which to base the vector.
 * @return A 3-component directional vector.
 */
KINLINE vec3 mat4_left (mat4 matrix) {
	vec3 left;
	left.x = -matrix.data[0];
	left.y = -matrix.data[1];
	left.z = -matrix.data[2];
	vec3_normalize(&left);
	return left;
}

/**
 * @brief Returns a right vector relative to the provided matrix.
 *
 * @param matrix The matrix from which to base the vector.
 * @return A 3-component directional vector.
 */
KINLINE vec3 mat4_right (mat4 matrix) {
	vec3 right;
	right.x = matrix.data[0];
	right.y = matrix.data[1];
	right.z = matrix.data[2];
	vec3_normalize(&right);
	return right;
}

/**
 * @brief Returns the position relative to the provided matrix.
 *
 * @param matrix The matrix from which to base the vector.
 * @return A 3-component positional vector.
 */
KINLINE vec3 mat4_position (mat4 matrix) {
	vec3 pos;
	pos.x = matrix.data[12];
	pos.y = matrix.data[13];
	pos.z = matrix.data[14];
	return pos;
}

/**
 * @brief Returns the scale relative to the provided matrix.
 *
 * @param matrix The matrix from which to base the vector.
 * @return A 3-component scale vector.
 */
KINLINE vec3 mat4_scale_get (mat4 matrix) {
	return (vec3){matrix.data[0], matrix.data[5], matrix.data[10]};
}

/**
 * @brief Performs m * v
 *
 * @param m The matrix to be multiplied.
 * @param v The vector to multiply by.
 * @return The transformed vector.
 */
KINLINE vec3 mat4_mul_vec3 (mat4 m, vec3 v) {
	return (vec3){
		v.x * m.data[0] + v.y * m.data[4] + v.z * m.data[8] + m.data[12],
		v.x * m.data[1] + v.y * m.data[5] + v.z * m.data[9] + m.data[13],
		v.x * m.data[2] + v.y * m.data[6] + v.z * m.data[10] + m.data[14]};
}

/**
 * @brief Performs v * m
 *
 * @param v The vector to bemultiplied.
 * @param m The matrix to be multiply by.
 * @return The transformed vector.
 */
KINLINE vec3 vec3_mul_mat4 (vec3 v, mat4 m) {
	return (vec3){
		v.x * m.data[0] + v.y * m.data[4] + v.z * m.data[8] + m.data[12],
		v.x * m.data[1] + v.y * m.data[5] + v.z * m.data[9] + m.data[13],
		v.x * m.data[2] + v.y * m.data[6] + v.z * m.data[10] + m.data[14]};
}

/**
 * @brief Performs m * v
 *
 * @param m The matrix to be multiplied.
 * @param v The vector to multiply by.
 * @return The transformed vector.
 */
KINLINE vec4 mat4_mul_vec4 (mat4 m, vec4 v) {
	return (vec4){
		v.x * m.data[0] + v.y * m.data[1] + v.z * m.data[2] + v.w * m.data[3],
		v.x * m.data[4] + v.y * m.data[5] + v.z * m.data[6] + v.w * m.data[7],
		v.x * m.data[8] + v.y * m.data[9] + v.z * m.data[10] + v.w * m.data[11],
		v.x * m.data[12] + v.y * m.data[13] + v.z * m.data[14] + v.w * m.data[15]};
}

/**
 * @brief Performs v * m
 *
 * @param v The vector to bemultiplied.
 * @param m The matrix to be multiply by.
 * @return The transformed vector.
 */
KINLINE vec4 vec4_mul_mat4 (vec4 v, mat4 m) {
	return (vec4){
		v.x * m.data[0] + v.y * m.data[4] + v.z * m.data[8] + v.w * m.data[12],
		v.x * m.data[1] + v.y * m.data[5] + v.z * m.data[9] + v.w * m.data[13],
		v.x * m.data[2] + v.y * m.data[6] + v.z * m.data[10] + v.w * m.data[14],
		v.x * m.data[3] + v.y * m.data[7] + v.z * m.data[11] + v.w * m.data[15]};
}

// ------------------------------------------
// Quaternion
// ------------------------------------------

/**
 * @brief Creates an identity quaternion.
 *
 * @return An identity quaternion.
 */
KINLINE quat quat_identity (void) { return (quat){0, 0, 0, 1.0f}; }

/**
 * @brief Indicates if the provided quaterion is identity. (0, 0, 0, 1)
 */
KAPI b8 quat_is_identity (quat q);

/**
 * @brief Returns the normal of the provided quaternion.
 *
 * @param q The quaternion.
 * @return The normal of the provided quaternion.
 */
KINLINE f32 quat_normal (quat q) {
	return ksqrt(q.x * q.x + q.y * q.y + q.z * q.z + q.w * q.w);
}

/**
 * @brief Returns a normalized copy of the provided quaternion.
 *
 * @param q The quaternion to normalize.
 * @return A normalized copy of the provided quaternion.
 */
KINLINE quat quat_normalize (quat q) {
	f32 normal = quat_normal(q);
	return (quat){q.x / normal, q.y / normal, q.z / normal, q.w / normal};
}

/**
 * @brief Returns the conjugate of the provided quaternion. That is,
 * The x, y and z elements are negated, but the w element is untouched.
 *
 * @param q The quaternion to obtain a conjugate of.
 * @return The conjugate quaternion.
 */
KINLINE quat quat_conjugate (quat q) { return (quat){-q.x, -q.y, -q.z, q.w}; }

/**
 * @brief Returns an inverse copy of the provided quaternion.
 *
 * @param q The quaternion to invert.
 * @return An inverse copy of the provided quaternion.
 */
KINLINE quat quat_inverse (quat q) { return quat_normalize(quat_conjugate(q)); }

/**
 * @brief Multiplies the provided quaternions.
 *
 * @param q_0 The first quaternion.
 * @param q_1 The second quaternion.
 * @return The multiplied quaternion.
 */
KINLINE quat quat_mul (quat q_0, quat q_1) {
	quat out_quaternion;

	out_quaternion.x =
		q_0.x * q_1.w + q_0.y * q_1.z - q_0.z * q_1.y + q_0.w * q_1.x;

	out_quaternion.y =
		-q_0.x * q_1.z + q_0.y * q_1.w + q_0.z * q_1.x + q_0.w * q_1.y;

	out_quaternion.z =
		q_0.x * q_1.y - q_0.y * q_1.x + q_0.z * q_1.w + q_0.w * q_1.z;

	out_quaternion.w =
		-q_0.x * q_1.x - q_0.y * q_1.y - q_0.z * q_1.z + q_0.w * q_1.w;

	return out_quaternion;
}

/**
 * @brief Calculates the dot product of the provided quaternions.
 *
 * @param q_0 The first quaternion.
 * @param q_1 The second quaternion.
 * @return The dot product of the provided quaternions.
 */
KINLINE f32 quat_dot (quat q_0, quat q_1) {
	return q_0.x * q_1.x + q_0.y * q_1.y + q_0.z * q_1.z + q_0.w * q_1.w;
}

/**
 * @brief Creates a rotation matrix from the given quaternion.
 *
 * @param q The quaternion to be used.
 * @return A rotation matrix.
 */
KINLINE mat4 quat_to_mat4 (quat q) {
#if 0
	mat4 out_matrix = mat4_identity();

	// https://stackoverflow.com/questions/1556260/convert-quaternion-rotation-to-rotation-matrix

	quat n = quat_normalize(q);

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

	out_matrix.data[4] = 2.0f * n.x * n.y + 2.0f * n.z * n.w;
	out_matrix.data[5] = 1.0f - 2.0f * n.x * n.x - 2.0f * n.z * n.z;
	out_matrix.data[6] = 2.0f * n.y * n.z - 2.0f * n.x * n.w;

	out_matrix.data[8] = 2.0f * n.x * n.z - 2.0f * n.y * n.w;
	out_matrix.data[9] = 2.0f * n.y * n.z + 2.0f * n.x * n.w;
	out_matrix.data[10] = 1.0f - 2.0f * n.x * n.x - 2.0f * n.y * n.y;

	return out_matrix;
#else
	/*mat4 m;
	float x = q.x, y = q.y, z = q.z, w = q.w;
	float xx = x * x, yy = y * y, zz = z * z;
	float xy = x * y, xz = x * z, yz = y * z;
	float wx = w * x, wy = w * y, wz = w * z;

	// Right (X)
	m.data[0] = 1 - 2 * (yy + zz);
	m.data[1] = 2 * (xy + wz);
	m.data[2] = 2 * (xz - wy);
	m.data[3] = 0;

	// Up (Y)
	m.data[4] = 2 * (xy - wz);
	m.data[5] = 1 - 2 * (xx + zz);
	m.data[6] = 2 * (yz + wx);
	m.data[7] = 0;

	// Forward (-Z)
	m.data[8] = -(2 * (xz + wy));
	m.data[9] = -(2 * (yz - wx));
	m.data[10] = -(1 - 2 * (xx + yy));
	m.data[11] = 0;

	// translation = identity for now
	m.data[12] = 0;
	m.data[13] = 0;
	m.data[14] = 0;
	m.data[15] = 1;

	return m;*/

	/*
	quat n = quat_normalize(q);

	mat4 m = {0};

	float x = n.x, y = n.y, z = n.z, w = n.w;

	// Right (X)
	m.data[0] = 1 - 2 * y * y - 2 * z * z;
	m.data[1] = 2 * x * y + 2 * w * z;
	m.data[2] = 2 * x * z - 2 * w * y;
	m.data[3] = 0;

	// Up (Y)
	m.data[4] = 2 * x * y - 2 * w * z;
	m.data[5] = 1 - 2 * x * x - 2 * z * z;
	m.data[6] = 2 * y * z + 2 * w * x;
	m.data[7] = 0;

	// Forward (-Z)
	m.data[8] = 2 * x * z + 2 * w * y;
	m.data[9] = 2 * y * z - 2 * w * x;
	m.data[10] = -(1 - 2 * x * x - 2 * y * y);
	m.data[11] = 0;

	m.data[12] = 0;
	m.data[13] = 0;
	m.data[14] = 0;
	m.data[15] = 1;

	return m;*/

	quat n = quat_normalize(q);
	float x = n.x, y = n.y, z = n.z, w = n.w;

	mat4 m = {0};

	// Right (X)
	m.data[0] = 1 - 2 * y * y - 2 * z * z;
	m.data[1] = 2 * x * y + 2 * w * z;
	m.data[2] = 2 * x * z - 2 * w * y;
	m.data[3] = 0;

	// Up (Y)
	m.data[4] = 2 * x * y - 2 * w * z;
	m.data[5] = 1 - 2 * x * x - 2 * z * z;
	m.data[6] = 2 * y * z + 2 * w * x;
	m.data[7] = 0;

	// Forward (-Z)
	m.data[8] = -(2 * x * z + 2 * w * y);
	m.data[9] = -(2 * y * z - 2 * w * x);
	m.data[10] = -(1 - 2 * x * x - 2 * y * y);
	m.data[11] = 0;

	// Translation identity
	m.data[12] = 0;
	m.data[13] = 0;
	m.data[14] = 0;
	m.data[15] = 1;

	return m;
#endif
}

/**
 * @brief Calculates a rotation matrix based on the quaternion and the passed in
 * center point.
 *
 * @param q The quaternion.
 * @param center The center point.
 * @return A rotation matrix.
 */
KINLINE mat4 quat_to_rotation_matrix (quat q, vec3 center) {
	mat4 out_matrix;

	f32 *o = out_matrix.data;
	o[0] = (q.x * q.x) - (q.y * q.y) - (q.z * q.z) + (q.w * q.w);
	o[1] = 2.0f * ((q.x * q.y) + (q.z * q.w));
	o[2] = 2.0f * ((q.x * q.z) - (q.y * q.w));
	o[3] = center.x - center.x * o[0] - center.y * o[1] - center.z * o[2];

	o[4] = 2.0f * ((q.x * q.y) - (q.z * q.w));
	o[5] = -(q.x * q.x) + (q.y * q.y) - (q.z * q.z) + (q.w * q.w);
	o[6] = 2.0f * ((q.y * q.z) + (q.x * q.w));
	o[7] = center.y - center.x * o[4] - center.y * o[5] - center.z * o[6];

	o[8] = 2.0f * ((q.x * q.z) + (q.y * q.w));
	o[9] = 2.0f * ((q.y * q.z) - (q.x * q.w));
	o[10] = -(q.x * q.x) - (q.y * q.y) + (q.z * q.z) + (q.w * q.w);
	o[11] = center.z - center.x * o[8] - center.y * o[9] - center.z * o[10];

	o[12] = 0.0f;
	o[13] = 0.0f;
	o[14] = 0.0f;
	o[15] = 1.0f;
	return out_matrix;
}

/**
 * @brief Creates a quaternion from the given axis and angle.
 *
 * @param axis The axis of rotation.
 * @param angle The angle of rotation.
 * @param normalize Indicates if the quaternion should be normalized.
 * @return A new quaternion.
 */
KINLINE quat quat_from_axis_angle (vec3 axis, f32 angle, b8 normalize) {
	const f32 half_angle = 0.5f * angle;
	f32 s = ksin(half_angle);
	f32 c = kcos(half_angle);

	quat q = (quat){s * axis.x, s * axis.y, s * axis.z, c};
	if (normalize) {
		return quat_normalize(q);
	}
	return q;
}

/**
 * @brief Calculates spherical linear interpolation of a given percentage
 * between two quaternions.
 *
 * @param q_0 The first quaternion.
 * @param q_1 The second quaternion.
 * @param percentage The percentage of interpolation, typically a value from
 * 0.0f-1.0f.
 * @return An interpolated quaternion.
 */
KINLINE quat quat_slerp (quat q_0, quat q_1, f32 percentage) {
	quat out_quaternion;
	// Source: https://en.wikipedia.org/wiki/Slerp
	// Only unit quaternions are valid rotations.
	// Normalize to avoid undefined behavior.
	quat v0 = quat_normalize(q_0);
	quat v1 = quat_normalize(q_1);

	// Compute the cosine of the angle between the two vectors.
	f32 dot = quat_dot(v0, v1);

	// If the dot product is negative, slerp won't take
	// the shorter path. Note that v1 and -v1 are equivalent when
	// the negation is applied to all four components. Fix by
	// reversing one quaternion.
	if (dot < 0.0f) {
		v1.x = -v1.x;
		v1.y = -v1.y;
		v1.z = -v1.z;
		v1.w = -v1.w;
		dot = -dot;
	}

	const f32 DOT_THRESHOLD = 0.9995f;
	if (dot > DOT_THRESHOLD) {
		// If the inputs are too close for comfort, linearly interpolate
		// and normalize the result.
		out_quaternion = (quat){v0.x + ((v1.x - v0.x) * percentage),
								v0.y + ((v1.y - v0.y) * percentage),
								v0.z + ((v1.z - v0.z) * percentage),
								v0.w + ((v1.w - v0.w) * percentage)};

		return quat_normalize(out_quaternion);
	}

	// Since dot is in range [0, DOT_THRESHOLD], acos is safe
	f32 theta_0 = kacos(dot);		  // theta_0 = angle between input vectors
	f32 theta = theta_0 * percentage; // theta = angle between v0 and result
	f32 sin_theta = ksin(theta);	  // compute this value only once
	f32 sin_theta_0 = ksin(theta_0);  // compute this value only once

	f32 s0 =
		kcos(theta) -
		dot * sin_theta / sin_theta_0; // == sin(theta_0 - theta) / sin(theta_0)
	f32 s1 = sin_theta / sin_theta_0;

	return (quat){(v0.x * s0) + (v1.x * s1), (v0.y * s0) + (v1.y * s1),
				  (v0.z * s0) + (v1.z * s1), (v0.w * s0) + (v1.w * s1)};
}

/**
 * @brief Converts provided degrees to radians.
 *
 * @param degrees The degrees to be converted.
 * @return The amount in radians.
 */
KINLINE f32 deg_to_rad (f32 degrees) { return degrees * K_DEG2RAD_MULTIPLIER; }

/**
 * @brief Converts provided radians to degrees.
 *
 * @param radians The radians to be converted.
 * @return The amount in degrees.
 */
KINLINE f32 rad_to_deg (f32 radians) { return radians * K_RAD2DEG_MULTIPLIER; }

/**
 * @brief Converts value from the "old" range to the "new" range.
 *
 * @param value The value to be converted.
 * @param from_min The minimum value from the old range.
 * @param from_max The maximum value from the old range.
 * @param to_min The minimum value from the new range.
 * @param to_max The maximum value from the new range.
 * @return The converted value.
 */
KINLINE f32 range_convert_f32 (f32 value, f32 old_min, f32 old_max, f32 new_min,
							   f32 new_max) {
	return (((value - old_min) * (new_max - new_min)) / (old_max - old_min)) +
		   new_min;
}

/**
 * @brief Converts rgb int values [0-255] to a single 32-bit integer.
 *
 * @param r The red value [0-255].
 * @param g The green value [0-255].
 * @param b The blue value [0-255].
 * @param out_u32 A pointer to hold the resulting integer.
 */
KINLINE void rgbu_to_u32 (u32 r, u32 g, u32 b, u32 *out_u32) {
	*out_u32 = (((r & 0x0FF) << 16) | ((g & 0x0FF) << 8) | (b & 0x0FF));
}

/**
 * @brief Converts the given 32-bit integer to rgb values [0-255].
 *
 * @param rgbu The integer holding a rgb value.
 * @param out_r A pointer to hold the red value.
 * @param out_g A pointer to hold the green value.
 * @param out_b A pointer to hold the blue value.
 */
KINLINE void u32_to_rgb (u32 rgbu, u32 *out_r, u32 *out_g, u32 *out_b) {
	*out_r = (rgbu >> 16) & 0x0FF;
	*out_g = (rgbu >> 8) & 0x0FF;
	*out_b = (rgbu) & 0x0FF;
}

/**
 * @brief Converts rgb integer values [0-255] to a vec3 of floating-point values
 * [0.0-1.0]
 *
 * @param r The red value [0-255].
 * @param g The green value [0-255].
 * @param b The blue value [0-255].
 * @param out_v A pointer to hold the vector of floating-point values.
 */
KINLINE void rgb_u32_to_vec3 (u32 r, u32 g, u32 b, vec3 *out_v) {
	out_v->r = r / 255.0f;
	out_v->g = g / 255.0f;
	out_v->b = b / 255.0f;
}

/**
 * @brief Converts a vec3 of rgbvalues [0.0-1.0] to integer rgb values [0-255].
 *
 * @param v The vector of rgb values [0.0-1.0] to be converted.
 * @param out_r A pointer to hold the red value.
 * @param out_g A pointer to hold the green value.
 * @param out_b A pointer to hold the blue value.
 */
KINLINE void vec3_to_rgb_u32 (vec3 v, u32 *out_r, u32 *out_g, u32 *out_b) {
	*out_r = v.r * 255;
	*out_g = v.g * 255;
	*out_b = v.b * 255;
}

KAPI plane_3d plane_3d_create (vec3 p1, vec3 norm);

/**
 * @brief Creates and returns a frustum based on the provided position,
 * direction vectors, aspect, field of view, and near/far clipping planes
 * (typically obtained from a camera). This is typically used for frustum
 * culling.
 *
 * @param position The position to be used.
 * @param target The target vector to be used.
 * @param up The up vector to be used.
 * @param aspect The aspect ratio.
 * @param fov The vertical field of view.
 * @param near The near clipping plane distance.
 * @param far The far clipping plane distance.
 * @return A shiny new frustum.
 */
KAPI kfrustum kfrustum_create (vec3 position, vec3 target, vec3 up, f32 aspect, f32 fov, f32 near, f32 far);

KAPI kfrustum kfrustum_from_view_projection (mat4 view_projection);

/**
 * Calculate the corner points of the provided frustum in world space, using
 * the given projection and view matrices.
 *
 * @param projection_view The combined projection/view matrix from the active camera.
 * @param corners An array of 8 vec4s to hold the caluclated points.
 */
KAPI void kfrustum_corner_points_world_space (mat4 projection_view, vec4 *corners);

/**
 * @brief Obtains the signed distance between the plane p and the provided
 * postion.
 *
 * @param p A constant pointer to a plane.
 * @param position A constant pointer to a position.
 * @return The signed distance from the point to the plane.
 */
KAPI f32 plane_signed_distance (const plane_3d *p, const vec3 *position);

/**
 * @brief Indicates if plane p intersects a sphere constructed via center and
 * radius.
 *
 * @param p A constant pointer to a plane.
 * @param center A constant pointer to a position representing the center of a
 * sphere.
 * @param radius The radius of the sphere.
 * @return True if the sphere intersects the plane; otherwise false.
 */
KAPI b8 plane_intersects_sphere (const plane_3d *p, const vec3 *center,
								 f32 radius);

/**
 * @brief Indicates if the frustum intersects (or contains) a sphere constructed
 * via center and radius.
 *
 * @param f A constant pointer to a frustum.
 * @param center A constant pointer to a position representing the center of a sphere.
 * @param radius The radius of the sphere.
 * @return True if the sphere is intersected by or contained within the frustum f; otherwise false.
 */
KAPI b8 kfrustum_intersects_sphere (const kfrustum *f, const vec3 *center, f32 radius);

/**
 * @brief Indicates if the frustum intersects (or contains) a sphere constructed
 * via center and radius.
 *
 * @param f A constant pointer to a frustum.
 * @param sphere A constant pointer to a sphere.
 * @return True if the sphere is intersected by or contained within the frustum f; otherwise false.
 */
KAPI b8 kfrustum_intersects_ksphere (const kfrustum *f, const ksphere *sphere);

/**
 * @brief Indicates if plane p intersects an axis-aligned bounding box
 * constructed via center and extents.
 *
 * @param p A constant pointer to a plane.
 * @param center A constant pointer to a position representing the center of an
 * axis-aligned bounding box.
 * @param extents The half-extents of an axis-aligned bounding box.
 * @return True if the axis-aligned bounding box intersects the plane; otherwise
 * false.
 */
KAPI b8 plane_intersects_aabb (const plane_3d *p, const vec3 *center, const vec3 *extents);

/**
 * @brief Indicates if frustum f intersects an axis-aligned bounding box
 * constructed via center and extents.
 *
 * @param f A constant pointer to a frustum.
 * @param center A constant pointer to a position representing the center of an
 * axis-aligned bounding box.
 * @param extents The half-extents of an axis-aligned bounding box.
 * @return True if the axis-aligned bounding box is intersected by or contained
 * within the frustum f; otherwise false.
 */
KAPI b8 kfrustum_intersects_aabb (const kfrustum *f, const vec3 *center,
								  const vec3 *extents);

KINLINE b8 rect_2d_contains_point (rect_2d rect, vec2 point) {
	return (point.x >= rect.x && point.x <= rect.x + rect.width) && (point.y >= rect.y && point.y <= rect.y + rect.height);
}

KAPI f32 vec3_distance_to_line (vec3 point, vec3 line_start, vec3 line_direction);

// Angle in degrees between 2 vectors.
KINLINE f32 vec3_angle (vec3 v_0, vec3 v_1);

KINLINE vec3 extents_2d_center (extents_2d extents) {
	return (vec3){
		(extents.min.x + extents.max.x) * 0.5f,
		(extents.min.y + extents.max.y) * 0.5f,
	};
}

KINLINE vec3 extents_2d_half (extents_2d extents) {
	return (vec3){
		kabs(extents.min.x - extents.max.x) * 0.5f,
		kabs(extents.min.y - extents.max.y) * 0.5f,
	};
}

KINLINE extents_3d extents_3d_zero (void) {
	return (extents_3d){
		.min = vec3_zero(),
		.max = vec3_zero()};
}

KINLINE extents_3d extents_3d_one (void) {
	return (extents_3d){
		.min = vec3_from_scalar(-1),
		.max = vec3_one()};
}

KINLINE extents_3d extents_3d_from_scalar (f32 scalar) {
	return (extents_3d){
		.min = vec3_from_scalar(-scalar),
		.max = vec3_from_scalar(scalar)};
}

KINLINE extents_3d extents_3d_from_size (vec3 size) {
	return (extents_3d){
		.min = vec3_mul_scalar(size, -1.0f),
		.max = size};
}

KINLINE vec3 extents_3d_center (extents_3d extents) {
	return (vec3){
		(extents.min.x + extents.max.x) * 0.5f,
		(extents.min.y + extents.max.y) * 0.5f,
		(extents.min.z + extents.max.z) * 0.5f,
	};
}

KINLINE vec3 extents_3d_half (extents_3d extents) {
	return (vec3){
		kabs(extents.min.x - extents.max.x) * 0.5f,
		kabs(extents.min.y - extents.max.y) * 0.5f,
		kabs(extents.min.z - extents.max.z) * 0.5f,
	};
}

KINLINE vec3 size_from_extents_3d (extents_3d extents) {
	vec3 size = vec3_sub(extents.max, extents.min);
	return (vec3){
		kabs(size.x),
		kabs(size.y),
		kabs(size.z)};
}

KINLINE extents_3d extents_combine (extents_3d a, extents_3d b) {
	return (extents_3d){
		.min = vec3_min(a.min, b.min),
		.max = vec3_max(a.max, b.max)};
}

KINLINE b8 extents_3d_is_zero (extents_3d extents) {
	return vec3_compare(size_from_extents_3d(extents), vec3_zero(), K_FLOAT_EPSILON);
}

KINLINE vec2 vec2_mid (vec2 v_0, vec2 v_1) {
	return (vec2){
		(v_0.x - v_1.x) * 0.5f,
		(v_0.y - v_1.y) * 0.5f};
}

KINLINE vec3 vec3_mid (vec3 v_0, vec3 v_1) {
	return (vec3){
		(v_0.x - v_1.x) * 0.5f,
		(v_0.y - v_1.y) * 0.5f,
		(v_0.z - v_1.z) * 0.5f};
}

KINLINE vec3 vec3_lerp (vec3 v_0, vec3 v_1, f32 t) {
	return vec3_add(vec3_mul_scalar(v_0, 1.0f - t), vec3_mul_scalar(v_1, t));
}

KINLINE vec3 triangle_get_normal (const triangle *tri) {
	vec3 edge1 = vec3_sub(tri->verts[1], tri->verts[0]);
	vec3 edge2 = vec3_sub(tri->verts[2], tri->verts[0]);

	vec3 normal = vec3_cross(edge1, edge2);
	return vec3_normalized(normal);
}

KINLINE quat quat_from_surface_normal (vec3 normal, vec3 reference_up) {
	normal = vec3_normalized(normal);
	reference_up = vec3_normalized(reference_up);

	// Compute rotation axis as the cross product
	vec3 axis = vec3_cross(reference_up, normal);
	f32 dot = vec3_dot(reference_up, normal);

	// If dot is near 1, the vectors are already aligned, return identity quaternion
	if (dot > 0.9999f) {
		return quat_identity();
	}

	// If dot is near -1, the vectors are opposite, use an arbitrary perpendicular axis
	if (dot < -0.9999f) {
		axis = vec3_normalized(vec3_cross(reference_up, (vec3){1, 0, 0})); // Try X-axis
		if (vec3_length_squared(axis) < K_FLOAT_EPSILON) {
			axis = vec3_normalized(vec3_cross(reference_up, (vec3){0, 0, 1})); // Try Z-axis
		}
	}

	// Compute the quaternion components
	f32 angle = kacos(dot); // Angle between the vectors
	return quat_from_axis_angle(axis, angle, false);
}

KINLINE aabb aabb_create (vec3 min, vec3 max) {
	return (aabb){
		.min = min,
		.max = max};
}

KINLINE aabb aabb_combine (aabb a_0, aabb a_1) {
	return (aabb){
		.min = vec3_min(a_0.min, a_1.min),
		.max = vec3_max(a_0.max, a_1.max)};
}

KINLINE b8 aabbs_intersect (aabb a_0, aabb a_1) {
	return !(
		a_0.max.x < a_1.min.x || a_0.min.x > a_1.max.x ||
		a_0.max.y < a_1.min.y || a_0.min.y > a_1.max.y ||
		a_0.max.z < a_1.min.z || a_0.min.z > a_1.max.z);
}

KINLINE f32 aabb_surface_area (aabb a) {
	f32 dx = a.max.x - a.min.x;
	f32 dy = a.max.y - a.min.y;
	f32 dz = a.max.z - a.min.z;
	return 2.0f * (dx * dy + dy * dz + dz * dx);
}

/**
 * @brief Expand the AABB by amount on all axes.
 */
KINLINE aabb aabb_expand (aabb a, f32 amount) {
	vec3 d = vec3_from_scalar(amount);
	return aabb_create(vec3_sub(a.min, d), vec3_add(a.max, d));
}

/**
 * @brief Constructs a tight, world-space AABB from an OBB defined by a center point (the
 * translation of mat) and the provided half_extents.
 *
 * @param half_extents Half of the extents of the OBB to be created from.
 * @param mat The world matrix to use when creating the AABB.
 * @returns A newly-created AABB.
 */
KINLINE aabb aabb_from_mat4 (vec3 half_extents, mat4 mat) {
	vec3 center = (vec3){mat.data[12], mat.data[13], mat.data[14]}; // The center
	f32 ax = kabs(mat.data[0]) * half_extents.x + kabs(mat.data[4]) * half_extents.y + kabs(mat.data[8]) * half_extents.z;
	f32 ay = kabs(mat.data[1]) * half_extents.x + kabs(mat.data[5]) * half_extents.y + kabs(mat.data[9]) * half_extents.z;
	f32 az = kabs(mat.data[2]) * half_extents.x + kabs(mat.data[6]) * half_extents.y + kabs(mat.data[10]) * half_extents.z;
	vec3 half = (vec3){ax, ay, az};
	return aabb_create(vec3_sub(center, half), vec3_add(center, half));
}

KINLINE aabb aabb_from_mat4_extents (vec3 local_min, vec3 local_max, mat4 mat) {
	// Local-space center and half extents
	vec3 local_center = vec3_mul_scalar(vec3_add(local_min, local_max), 0.5f);
	vec3 half_extents = vec3_mul_scalar(vec3_sub(local_max, local_min), 0.5f);

	// Transform center into world space
	vec3 center = {
		mat.data[0] * local_center.x + mat.data[4] * local_center.y + mat.data[8] * local_center.z + mat.data[12],
		mat.data[1] * local_center.x + mat.data[5] * local_center.y + mat.data[9] * local_center.z + mat.data[13],
		mat.data[2] * local_center.x + mat.data[6] * local_center.y + mat.data[10] * local_center.z + mat.data[14]};

	// Compute world-space half extents
	f32 ax = kabs(mat.data[0]) * half_extents.x +
			 kabs(mat.data[4]) * half_extents.y +
			 kabs(mat.data[8]) * half_extents.z;

	f32 ay = kabs(mat.data[1]) * half_extents.x +
			 kabs(mat.data[5]) * half_extents.y +
			 kabs(mat.data[9]) * half_extents.z;

	f32 az = kabs(mat.data[2]) * half_extents.x +
			 kabs(mat.data[6]) * half_extents.y +
			 kabs(mat.data[10]) * half_extents.z;

	vec3 half = {ax, ay, az};

	return aabb_create(
		vec3_sub(center, half),
		vec3_add(center, half));
}

/**
 * @brief Indicates if point is inside the provided AABB.
 *
 * @param point The point to test.
 * @param box The box to test against.
 * @returns True if inside; otherwise false.
 */
KINLINE b8 aabb_contains_point (vec3 point, aabb box) {
	return (point.x >= box.min.x && point.x <= box.max.x) &&
		   (point.y >= box.min.y && point.y <= box.max.y) &&
		   (point.z >= box.min.z && point.z <= box.max.z);
}

/**
 * @brief Indicates if AABB b is completely inside AABB a.
 *
 * @param a The outer AABB to test.
 * @param b The inner AABB to test.
 * @returns True if inside; otherwise false.
 */
KINLINE b8 aabb_contains_aabb (aabb a, aabb b) {
	return (
		b.min.x >= a.min.x &&
		b.min.y >= a.min.y &&
		b.min.z >= a.min.z &&
		b.max.x <= a.max.x &&
		b.max.y <= a.max.y &&
		b.max.z <= a.max.z);
}

KINLINE ray ray_create (vec3 position, vec3 direction, f32 max_distance, ray_flags flags) {
	return (ray){
		.origin = position,
		.direction = direction,
		.max_distance = max_distance,
		.flags = flags};
}

KAPI ray ray_transformed (const ray *r, mat4 transform);

KAPI ray ray_from_screen (vec2i screen_pos, rect_2di viewport_rect, vec3 origin, mat4 view, mat4 projection);

KAPI b8 ray_intersects_aabb (aabb box, vec3 origin, vec3 direction, f32 max, f32 *out_min, f32 *out_max);

KAPI b8 raycast_plane_3d (const ray *r, const plane_3d *p, vec3 *out_point, f32 *out_distance);

KAPI b8 raycast_disc_3d (const ray *r, vec3 center, vec3 normal, f32 outer_radius, f32 inner_radius, vec3 *out_point, f32 *out_distance);

KAPI b8 ray_intersects_triangle (const ray *r, const triangle *t);

KAPI b8 ray_pick_triangle (const ray *r, b8 backface_cull, u32 vertex_count, u32 vertex_element_size, const void *vertices, u32 index_count, const u32 *indices, triangle *out_triangle, vec3 *out_hit_pos, vec3 *out_hit_normal);

KAPI b8 ray_intersects_sphere (const ray *r, vec3 center, f32 radius, vec3 *out_point, f32 *out_distance);

typedef struct kintersect_result {
	vec3 normal;
	vec3 closest_point_a;
	vec3 closest_point_b;
	f32 depth;
} kintersect_result;

KINLINE f32 obb_project_extents (const struct obb *o, vec3 axis) {
	return kabs(vec3_dot(o->axis[0], axis)) * o->half_extents.x +
		   kabs(vec3_dot(o->axis[1], axis)) * o->half_extents.y +
		   kabs(vec3_dot(o->axis[2], axis)) * o->half_extents.z;
}

KINLINE vec3 obb_closest_point (const struct obb *o, vec3 p) {
	vec3 d = vec3_sub(p, o->center);
	vec3 out = o->center;

	for (u8 i = 0; i < 3; ++i) {
		f32 dist = vec3_dot(d, o->axis[i]);
		f32 clamped = KCLAMP(dist, -o->half_extents.elements[i], o->half_extents.elements[i]);

		out = vec3_add(out, vec3_mul_scalar(o->axis[i], clamped));
	}

	return out;
}

KINLINE b8 sat_overlap (vec3 axis, vec3 t, const struct obb *a, const struct obb *b, f32 *min_overlap, vec3 *best_axis) {
	f32 dist = kabs(vec3_dot(t, axis));
	f32 ra = obb_project_extents(a, axis);
	f32 rb = obb_project_extents(b, axis);
	f32 overlap = (ra + rb) - dist;

	if (overlap < 0.0f) {
		return false;
	}

	if (overlap < *min_overlap) {
		*min_overlap = overlap;
		*best_axis = axis;
	}

	return true;
}

KINLINE b8 obb_intersects_obb (const struct obb *a, const struct obb *b, kintersect_result *out_result) {
	vec3 t = vec3_sub(b->center, a->center);

	f32 min_overlap = K_FLOAT_MAX;
	vec3 best_axis = vec3_zero();

	for (u8 i = 0; i < 3; ++i) {
		if (!sat_overlap(a->axis[i], t, a, b, &min_overlap, &best_axis)) {
			return false;
		}
	}
	for (u8 i = 0; i < 3; ++i) {
		if (!sat_overlap(b->axis[i], t, a, b, &min_overlap, &best_axis)) {
			return false;
		}
	}

	for (u8 i = 0; i < 3; ++i) {
		for (u8 j = 0; j < 3; ++j) {
			vec3 axis = vec3_cross(a->axis[i], b->axis[j]);

			f32 len = vec3_dot(axis, axis);
			if (len < 1e-8f) {
				continue;
			}

			axis = vec3_mul_scalar(axis, 1.0f / ksqrt(len));

			if (!sat_overlap(axis, t, a, b, &min_overlap, &best_axis)) {
				return false;
			}
		}
	}

	if (out_result) {
		// If there's a collision, get the data.
		vec3 n = best_axis;
		if (vec3_dot(n, t) < 0.0f) {
			n = vec3_mul_scalar(n, -1.0f); // Make sure the normal points from a -> b
		}

		out_result->normal = vec3_normalized(n);
		out_result->depth = min_overlap;

		// Get closts points.
		out_result->closest_point_a = obb_closest_point(a, b->center);
		out_result->closest_point_b = obb_closest_point(b, a->center);
	}

	return true;
}

KINLINE b8 obb_intersects_sphere (const struct obb *o, const struct ksphere *s) {
	vec3 d = vec3_sub(s->position, o->center);

	f32 dist_sq = 0.0f;

	for (u8 i = 0; i < 3; ++i) {
		f32 projected = vec3_dot(d, o->axis[i]);
		f32 clamped = KCLAMP(projected, -o->half_extents.elements[i], o->half_extents.elements[i]);
		f32 diff = projected - clamped;
		dist_sq += (diff * diff);
	}

	return dist_sq <= (s->radius * s->radius);
}

KINLINE b8 sphere_intersects_sphere (const ksphere a, const ksphere b) {
	f32 dist_sq = vec3_distance_squared(a.position, b.position);
	f32 combined_radii_sq = (a.radius * a.radius) + (b.radius * b.radius);

	return dist_sq < combined_radii_sq;
}

KINLINE obb aabb_to_obb (const aabb a, mat4 m) {
	obb out;

	// AABBs by default don't have half-extents and a center. Extract them.
	vec3 half_extents = extents_3d_half(a);
	out.center = extents_3d_center(a);

	out.center = vec3_transform(out.center, 1.0f, m);

	out.axis[0] = (vec3){m.data[0], m.data[4], m.data[8]};
	out.axis[1] = (vec3){m.data[1], m.data[5], m.data[9]};
	out.axis[2] = (vec3){m.data[2], m.data[6], m.data[10]};

	out.half_extents = half_extents;

	return out;
}

KAPI vec3 normal_from_triangle (const triangle *tri);

KAPI vec3 closest_point_on_segment (vec3 a, vec3 b, vec3 pos);

KAPI vec3 vec3_project_on_plane (vec3 pos, vec3 normal);

KAPI vec3 closest_point_on_triangle (vec3 pos, const triangle *t);
// Return whether point P is contained inside 3D region delimited by triangle T0,T1,T2 edges.
KAPI b8 point_inside_triangle (vec3 p, const triangle *tri);
// Return whether point P is contained inside 3D region delimited by parallelogram P0,P1,P2 edges.
KAPI b8 point_inside_parallelogram (vec3 p, vec3 p0, vec3 p1, vec3 p2);
// Return whether point P is contained inside a triangular prism A0,A1,A2-B0,B1,B2.
KAPI b8 point_inside_triangular_prism (vec3 p, vec3 a0, vec3 a1, vec3 a2, vec3 b0, vec3 b1, vec3 b2);
// Sweep sphere C,r with velocity Sv against plane N of triangle T0,T1,T2, ignoring edges.
KAPI b8 sweep_sphere_triangle_plane (sweep_result *sweep, vec3 c, f32 r, vec3 v, vec3 t0, vec3 t1, vec3 t2, vec3 n);
// Sweep sphere C,r with velocity V against plane N of parallelogram P0,P1,P2 ignoring edges.
KAPI b8 sweep_sphere_parallelogram_plane (sweep_result *sweep, vec3 c, f32 r, vec3 v, vec3 p0, vec3 p1, vec3 p2, vec3 n);
// Sweep point P with velocity V against sphere S,r.
KAPI b8 sweep_point_sphere (sweep_result *sweep, vec3 p, vec3 v, vec3 s, f32 r, vec3 fallback_normal);
// Sweep point P with velocity V against cylinder C0,C1,r, ignoring the endcaps.
KAPI b8 sweep_point_uncapped_cylinder (sweep_result *sweep, vec3 p, vec3 v, vec3 c0, vec3 c1, f32 r, vec3 fallback_normal);

// Sweep a capsule C0,C1,Cr with velocity Cv against the triangle T0,T1,T2.
//   c0,c1      capsule line segment endpoints
//   r          capsule radius
//   v          capsule velocity
//   t0,t1,t2   3 triangle vertices
//   returns    whether the capsule and triangle intersect
KAPI b8 sweep_capsule_triangle (sweep_result *s, vec3 c0, vec3 c1, f32 r, vec3 v, vec3 t0, vec3 t1, vec3 t2);