kmath.h

Go to the documentation of this file.

 
#pragma once
 
#include "defines.h"
#include "math_types.h"
#include "memory/kmemory.h"
#include <float.h>
 
#define K_PI 3.14159265358979323846f
 
#define K_2PI (2.0f * K_PI)
 
#define K_4PI (4.0f * K_PI)
 
#define K_HALF_PI (0.5f * K_PI)
 
#define K_QUARTER_PI (0.25f * K_PI)
 
#define K_ONE_OVER_PI (1.0f / K_PI)
 
#define K_ONE_OVER_TWO_PI (1.0f / K_2PI)
 
#define K_SQRT_TWO 1.41421356237309504880f
 
#define K_SQRT_THREE 1.73205080756887729352f
 
#define K_SQRT_ONE_OVER_TWO 0.70710678118654752440f
 
#define K_SQRT_ONE_OVER_THREE 0.57735026918962576450f
 
#define K_DEG2RAD_MULTIPLIER (K_PI / 180.0f)
 
#define K_RAD2DEG_MULTIPLIER (180.0f / K_PI)
 
#define K_SEC_TO_US_MULTIPLIER (1000.0f * 1000.0f)
 
#define K_SEC_TO_MS_MULTIPLIER 1000.0f
 
#define K_MS_TO_SEC_MULTIPLIER 0.001f
 
#define K_INFINITY (1e30f * 1e30f)
 
#define K_FLOAT_EPSILON 1.192092896e-07f
 
#define K_FLOAT_MAX 3.40282e+38F
#define K_FLOAT_MIN -K_FLOAT_MAX
 
// ------------------------------------------
// General math functions
// ------------------------------------------
 
KINLINE void kswapf(f32* a, f32* b) {
    f32 temp = *a;
    *a = *b;
    *b = temp;
}
 
KINLINE f32 ksign(f32 x) {
    return x == 0.0f ? 0.0f : x < 0.0f ? -1.0f
                                       : 1.0f;
}
 
KINLINE f32 kstep(f32 edge, f32 x) {
    return x < edge ? 0.0f : 1.0f;
}
 
KAPI f32 ksin(f32 x);
 
KAPI f32 kcos(f32 x);
 
KAPI f32 ktan(f32 x);
 
KAPI f32 katan(f32 x);
 
KAPI f32 katan2(f32 x, f32 y);
 
KAPI f32 kasin(f32 x);
 
KAPI f32 kacos(f32 x);
 
KAPI f32 ksqrt(f32 x);
 
KAPI f32 kabs(f32 x);
 
KAPI f32 kfloor(f32 x);
 
KAPI f32 kceil(f32 x);
 
KAPI f32 klog(f32 x);
 
KAPI f32 klog2(f32 x);
 
KAPI f32 kpow(f32 x, f32 y);
 
KAPI f32 kmod(f32 x, f32 y);
 
KINLINE f32 klerp(f32 a, f32 b, f32 t) {
    return a + t * (b - a);
}
 
KINLINE b8 is_power_of_2(u64 value) {
    return (value != 0) && ((value & (value - 1)) == 0);
}
 
KAPI i32 krandom(void);
 
KAPI i32 krandom_in_range(i32 min, i32 max);
 
KAPI u64 krandom_u64(void);
 
KAPI f32 kfrandom(void);
 
KAPI f32 kfrandom_in_range(f32 min, f32 max);
 
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);
}
 
KAPI f32 kattenuation_min_max(f32 min, f32 max, f32 x);
 
KINLINE b8 kfloat_compare(f32 f_0, f32 f_1) {
    return kabs(f_0 - f_1) < K_FLOAT_EPSILON;
}
 
// ------------------------------------------
// Vector 2
// ------------------------------------------
 
KINLINE vec2 vec2_create(f32 x, f32 y) {
    return (vec2){x, y};
}
 
KINLINE vec2 vec2_from_scalar(f32 scalar) { return (vec2){scalar, scalar}; }
 
KINLINE vec2 vec2_zero(void) { return (vec2){0.0f, 0.0f}; }
 
KINLINE vec2 vec2_one(void) { return (vec2){1.0f, 1.0f}; }
 
KINLINE vec2 vec2_up(void) { return (vec2){0.0f, 1.0f}; }
 
KINLINE vec2 vec2_down(void) { return (vec2){0.0f, -1.0f}; }
 
KINLINE vec2 vec2_left(void) { return (vec2){-1.0f, 0.0f}; }
 
KINLINE vec2 vec2_right(void) { return (vec2){1.0f, 0.0f}; }
 
KINLINE vec2 vec2_add(vec2 vector_0, vec2 vector_1) {
    return (vec2){vector_0.x + vector_1.x, vector_0.y + vector_1.y};
}
 
KINLINE vec2 vec2_sub(vec2 vector_0, vec2 vector_1) {
    return (vec2){vector_0.x - vector_1.x, vector_0.y - vector_1.y};
}
 
KINLINE vec2 vec2_mul(vec2 vector_0, vec2 vector_1) {
    return (vec2){vector_0.x * vector_1.x, vector_0.y * vector_1.y};
}
 
KINLINE vec2 vec2_mul_scalar(vec2 vector_0, f32 scalar) {
    return (vec2){vector_0.x * scalar, vector_0.y * scalar};
}
 
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};
}
 
KINLINE vec2 vec2_div(vec2 vector_0, vec2 vector_1) {
    return (vec2){vector_0.x / vector_1.x, vector_0.y / vector_1.y};
}
 
KINLINE f32 vec2_length_squared(vec2 vector) {
    return vector.x * vector.x + vector.y * vector.y;
}
 
KINLINE f32 vec2_length(vec2 vector) {
    return ksqrt(vec2_length_squared(vector));
}
 
KINLINE void vec2_normalize(vec2* vector) {
    const f32 length = vec2_length(*vector);
    vector->x /= length;
    vector->y /= length;
}
 
KINLINE vec2 vec2_normalized(vec2 vector) {
    vec2_normalize(&vector);
    return vector;
}
 
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;
}
 
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]);
        }
    }
}
 
KINLINE vec2 vec2_clamped(vec2 vector, vec2 min, vec2 max) {
    vec2_clamp(&vector, min, max);
    return vector;
}
 
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);
        }
    }
}
 
KINLINE vec2 vec2_clamped_scalar(vec2 vector, f32 min, f32 max) {
    vec2_clamp_scalar(&vector, min, max);
    return vector;
}
 
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);
}
 
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
// ------------------------------------------
 
#define vec3_create(x, y, z) \
    (vec3) { x, y, z }
 
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};
}
 
KINLINE vec4 vec3_to_vec4(vec3 vector, f32 w) {
    return (vec4){vector.x, vector.y, vector.z, w};
}
 
KINLINE vec3 vec3_zero(void) { return (vec3){0.0f, 0.0f, 0.0f}; }
 
KINLINE vec3 vec3_one(void) { return (vec3){1.0f, 1.0f, 1.0f}; }
 
KINLINE vec3 vec3_up(void) { return (vec3){0.0f, 1.0f, 0.0f}; }
 
KINLINE vec3 vec3_down(void) { return (vec3){0.0f, -1.0f, 0.0f}; }
 
KINLINE vec3 vec3_left(void) { return (vec3){-1.0f, 0.0f, 0.0f}; }
 
KINLINE vec3 vec3_right(void) { return (vec3){1.0f, 0.0f, 0.0f}; }
 
KINLINE vec3 vec3_forward(void) { return (vec3){0.0f, 0.0f, -1.0f}; }
 
KINLINE vec3 vec3_backward(void) { return (vec3){0.0f, 0.0f, 1.0f}; }
 
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};
}
 
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};
}
 
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};
}
 
KINLINE vec3 vec3_mul_scalar(vec3 vector_0, f32 scalar) {
    return (vec3){vector_0.x * scalar, vector_0.y * scalar, vector_0.z * scalar};
}
 
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};
}
 
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;
}
 
KINLINE f32 vec3_length_squared(vec3 vector) {
    return vector.x * vector.x + vector.y * vector.y + vector.z * vector.z;
}
 
KINLINE f32 vec3_length(vec3 vector) {
    return ksqrt(vec3_length_squared(vector));
}
 
KINLINE void vec3_normalize(vec3* vector) {
    const f32 length = vec3_length(*vector);
    vector->x /= length;
    vector->y /= length;
    vector->z /= length;
}
 
KINLINE vec3 vec3_normalized(vec3 vector) {
    vec3_normalize(&vector);
    return vector;
}
 
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;
}
 
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};
}
 
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;
}
 
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]);
        }
    }
}
 
KINLINE vec3 vec3_clamped(vec3 vector, vec3 min, vec3 max) {
    vec3_clamp(&vector, min, max);
    return vector;
}
 
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);
        }
    }
}
 
KINLINE vec3 vec3_clamped_scalar(vec3 vector, f32 min, f32 max) {
    vec3_clamp_scalar(&vector, min, max);
    return vector;
}
 
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);
}
 
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);
}
 
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);
}
 
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
// ------------------------------------------
 
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;
}
 
KINLINE vec4 vec4_from_scalar(f32 scalar) { return (vec4){scalar, scalar, scalar, scalar}; }
 
KINLINE vec3 vec4_to_vec3(vec4 vector) {
    return (vec3){vector.x, vector.y, vector.z};
}
 
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
}
 
KINLINE vec4 vec4_zero(void) { return (vec4){0.0f, 0.0f, 0.0f, 0.0f}; }
 
KINLINE vec4 vec4_one(void) { return (vec4){1.0f, 1.0f, 1.0f, 1.0f}; }
 
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;
}
 
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;
}
 
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;
}
 
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};
}
 
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,
    };
}
 
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;
}
 
KINLINE f32 vec4_length_squared(vec4 vector) {
    return vector.x * vector.x + vector.y * vector.y + vector.z * vector.z +
           vector.w * vector.w;
}
 
KINLINE f32 vec4_length(vec4 vector) {
    return ksqrt(vec4_length_squared(vector));
}
 
KINLINE void vec4_normalize(vec4* vector) {
    const f32 length = vec4_length(*vector);
    vector->x /= length;
    vector->y /= length;
    vector->z /= length;
    vector->w /= length;
}
 
KINLINE vec4 vec4_normalized(vec4 vector) {
    vec4_normalize(&vector);
    return vector;
}
 
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;
}
 
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;
}
 
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]);
        }
    }
}
 
KINLINE vec4 vec4_clamped(vec4 vector, vec4 min, vec4 max) {
    vec4_clamp(&vector, min, max);
    return vector;
}
 
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);
        }
    }
}
 
KINLINE vec4 vec4_clamped_scalar(vec4 vector, f32 min, f32 max) {
    vec4_clamp_scalar(&vector, min, max);
    return vector;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
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;
}
 
KINLINE vec3 mat4_scale_get(mat4 matrix) {
    return (vec3){matrix.data[0], matrix.data[5], matrix.data[10]};
}
 
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]};
}
 
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]};
}
 
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]};
}
 
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
// ------------------------------------------
 
KINLINE quat quat_identity(void) { return (quat){0, 0, 0, 1.0f}; }
 
KAPI b8 quat_is_identity(quat q);
 
KINLINE f32 quat_normal(quat q) {
    return ksqrt(q.x * q.x + q.y * q.y + q.z * q.z + q.w * q.w);
}
 
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};
}
 
KINLINE quat quat_conjugate(quat q) { return (quat){-q.x, -q.y, -q.z, q.w}; }
 
KINLINE quat quat_inverse(quat q) { return quat_normalize(quat_conjugate(q)); }
 
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;
}
 
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;
}
 
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
}
 
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;
}
 
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;
}
 
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)};
}
 
KINLINE f32 deg_to_rad(f32 degrees) { return degrees * K_DEG2RAD_MULTIPLIER; }
 
KINLINE f32 rad_to_deg(f32 radians) { return radians * K_RAD2DEG_MULTIPLIER; }
 
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;
}
 
KINLINE void rgbu_to_u32(u32 r, u32 g, u32 b, u32* out_u32) {
    *out_u32 = (((r & 0x0FF) << 16) | ((g & 0x0FF) << 8) | (b & 0x0FF));
}
 
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;
}
 
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;
}
 
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);
 
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);
 
KAPI void kfrustum_corner_points_world_space(mat4 projection_view, vec4* corners);
 
KAPI f32 plane_signed_distance(const plane_3d* p, const vec3* position);
 
KAPI b8 plane_intersects_sphere(const plane_3d* p, const vec3* center,
                                f32 radius);
 
KAPI b8 kfrustum_intersects_sphere(const kfrustum* f, const vec3* center, f32 radius);
 
KAPI b8 kfrustum_intersects_ksphere(const kfrustum* f, const ksphere* sphere);
 
KAPI b8 plane_intersects_aabb(const plane_3d* p, const vec3* center, const vec3* extents);
 
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);
 
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);
}
 
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));
}
 
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));
}
 
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);
}
 
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) {
    return (ray){
        .origin = position,
        .direction = direction};
}
 
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, void* vertices, u32 index_count, 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;
}