tinyfm.h

#pragma once

#include <iostream>
#include <iomanip>
#include <cstddef>
#include <cmath>
#include <cfloat>
#include <algorithm>
#include <vector>

#if !defined(TFM_MATRIX_RTL)
#define TFM_MATRIX_RTL 1
#endif

#if defined(_MSC_VER)
#define TFM_INLINE __forceinline
#elif defined(__GNUC__) || defined(__clang__)
#define TFM_INLINE inline __attribute__((always_inline))
#else
#define TFM_INLINE inline
#endif

#if defined(TFM_FORCE_ALIGNED)
#define TFM_ALIGNAS(x) alignas(x)
#else
#define TFM_ALIGNAS(x)
#endif

#if defined(TFM_FORCE_INTRINSICS)
#if defined(__SSE__) || defined(__AVX__) || defined(__AVX2__)
#define TFM_SSE_AVX_INTRINSICS 1
#include <immintrin.h>
// Loads/stores of tinyfm struct members: aligned instructions when TFM_FORCE_ALIGNED
// guarantees 16-byte alignment, unaligned otherwise. External user pointers always
// use the unaligned intrinsics directly.
#if defined(TFM_FORCE_ALIGNED)
#define TFM_MM_LOAD_PS(p) _mm_load_ps(p)
#define TFM_MM_STORE_PS(p, v) _mm_store_ps(p, v)
#define TFM_MM_LOAD_SI128(p) _mm_load_si128(p)
#define TFM_MM_STORE_SI128(p, v) _mm_store_si128(p, v)
#else
#define TFM_MM_LOAD_PS(p) _mm_loadu_ps(p)
#define TFM_MM_STORE_PS(p, v) _mm_storeu_ps(p, v)
#define TFM_MM_LOAD_SI128(p) _mm_loadu_si128(p)
#define TFM_MM_STORE_SI128(p, v) _mm_storeu_si128(p, v)
#endif
TFM_INLINE __m128 _mm_dot_ps(__m128 a, __m128 b) {
	__m128 m = _mm_mul_ps(a, b);
	__m128 shuf = _mm_movehdup_ps(m);
	__m128 sums = _mm_add_ps(m, shuf);
	shuf = _mm_movehl_ps(shuf, sums);
	sums = _mm_add_ss(sums, shuf);
	return _mm_shuffle_ps(sums, sums, 0x00);
}
TFM_INLINE float _mm_length_f32(__m128 a) {
	__m128 dp = _mm_dot_ps(a, a);
	__m128 len = _mm_sqrt_ps(dp);
	return _mm_cvtss_f32(len);
}
TFM_INLINE __m128 _mm_normalize_ps(__m128 a) {
	const __m128 epsilon = _mm_set1_ps(1e-8f);
	__m128 len_sq = _mm_dot_ps(a, a);
	__m128 mask = _mm_cmpgt_ps(len_sq, epsilon);
	__m128 inv_len = _mm_div_ps(_mm_set1_ps(1.0f), _mm_sqrt_ps(len_sq));
	__m128 normalized = _mm_mul_ps(a, inv_len);
	return _mm_and_ps(normalized, mask);
}
#elif defined(__ARM_NEON) || defined(__ARM_NEON__)
#define TFM_NEON_INTRINSICS 1
#include <arm_neon.h>
#endif
#endif

#pragma push_macro("min")
#pragma push_macro("max")
#undef min
#undef max

namespace tfm {
	constexpr TFM_INLINE float pi() noexcept { return 3.14159265358979323846f; }
	constexpr TFM_INLINE float two_pi()  noexcept { return 6.28318530717958647692f; }
	constexpr TFM_INLINE float half_pi() noexcept { return 1.57079632679489661923f; }

	constexpr TFM_INLINE float radians(float degrees) noexcept { return degrees * 0.017453292519943295f; }

	constexpr TFM_INLINE float lerp(float a, float b, float t) noexcept { return a + (b - a) * t; }

	constexpr TFM_INLINE float mix(float a, float b, float t) noexcept { return a + (b - a) * t; }

	struct vec2i;
	struct vec2f;
	struct vec3f;
	struct vec3i;
	struct vec4i;
	struct vec4f;
	struct quatf;
	struct mat3f;
	struct mat4f;

	struct TFM_ALIGNAS(8) vec2i {
		union {
			struct { int c0, c1; };
			int col[2];
			struct { int x, y; };
			struct { int r, g; };
			struct { int u, v; };
		};

		static TFM_INLINE const vec2i zero() noexcept;
		static TFM_INLINE const vec2i one() noexcept;

		TFM_INLINE vec2i() noexcept = default;
		TFM_INLINE vec2i(int x, int y) noexcept;
		TFM_INLINE vec2i(int xyz) noexcept;
		TFM_INLINE vec2i(const int* ptr) noexcept;

		TFM_INLINE int& operator[](std::size_t i) noexcept;
		TFM_INLINE const int& operator[](std::size_t i) const noexcept;
	};

	TFM_INLINE bool operator==(const vec2i& a, const vec2i& b) noexcept;
	TFM_INLINE bool operator!=(const vec2i& a, const vec2i& b) noexcept;
	TFM_INLINE vec2i abs(const vec2i& v) noexcept;
	TFM_INLINE vec2i max(const vec2i& a, const vec2i& b) noexcept;
	TFM_INLINE vec2i min(const vec2i& a, const vec2i& b) noexcept;
	TFM_INLINE vec2i clamp(const vec2i& v, const vec2i& minv, const vec2i& maxv) noexcept;
	TFM_INLINE int dot(const vec2i& a, const vec2i& b) noexcept;
	TFM_INLINE int manhattan_distance(const vec2i& a, const vec2i& b) noexcept;

	struct TFM_ALIGNAS(8) vec2f {
		union {
			struct { float c0, c1; };
			float col[2];
			struct { float x, y; };
			struct { float r, g; };
			struct { float u, v; };
		};

		static TFM_INLINE const vec2f zero() noexcept;
		static TFM_INLINE const vec2f one() noexcept;
		static TFM_INLINE const vec2f unit_x() noexcept;
		static TFM_INLINE const vec2f unit_y() noexcept;

		TFM_INLINE vec2f() noexcept = default;
		TFM_INLINE vec2f(float x, float y) noexcept;
		TFM_INLINE vec2f(float xy) noexcept;
		TFM_INLINE vec2f(const float* ptr) noexcept;
		TFM_INLINE vec2f(const double* ptr) noexcept;
		TFM_INLINE vec2f(const std::vector<float>& v) noexcept;
		TFM_INLINE vec2f(const std::vector<double>& v) noexcept;
		vec2f(const vec2f& v) noexcept = default;
		vec2f& operator=(const vec2f& v) noexcept = default;
		TFM_INLINE vec2f(const vec3f& v) noexcept;
		TFM_INLINE vec2f(const vec4f& v) noexcept;

		TFM_INLINE float length() const noexcept;
		TFM_INLINE float length_squared() const noexcept;
		TFM_INLINE vec2f normalize() const noexcept;
		TFM_INLINE vec2f round() const noexcept;
		TFM_INLINE vec2f ceil() const noexcept;

		TFM_INLINE float& operator[](std::size_t i) noexcept;
		TFM_INLINE const float& operator[](std::size_t i) const noexcept;
	};

	TFM_INLINE bool operator==(const vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE bool operator!=(const vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE vec2f& operator+=(vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE vec2f& operator+=(vec2f& a, float b) noexcept;
	TFM_INLINE vec2f& operator-=(vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE vec2f& operator-=(vec2f& a, float b) noexcept;
	TFM_INLINE vec2f& operator*=(vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE vec2f& operator*=(vec2f& a, float b) noexcept;
	TFM_INLINE vec2f& operator/=(vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE vec2f& operator/=(vec2f& a, float b) noexcept;
	TFM_INLINE vec2f operator+(const vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE vec2f operator+(const vec2f& a, float b) noexcept;
	TFM_INLINE vec2f operator+(float a, const vec2f& b) noexcept;
	TFM_INLINE vec2f operator-(const vec2f& a) noexcept;
	TFM_INLINE vec2f operator-(const vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE vec2f operator-(const vec2f& a, float b) noexcept;
	TFM_INLINE vec2f operator-(float a, const vec2f& b) noexcept;
	TFM_INLINE vec2f operator*(const vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE vec2f operator*(const vec2f& a, float b) noexcept;
	TFM_INLINE vec2f operator*(float a, const vec2f& b) noexcept;
	TFM_INLINE vec2f operator/(const vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE vec2f operator/(const vec2f& a, float b) noexcept;
	TFM_INLINE vec2f operator/(float a, const vec2f& b) noexcept;
	TFM_INLINE float dot(const vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE float distance(const vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE float distance_squared(const vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE float cross(const vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE vec2f project(const vec2f& a, const vec2f& onto) noexcept;
	TFM_INLINE vec2f reflect(const vec2f& v, const vec2f& normal) noexcept;
	TFM_INLINE vec2f lerp(const vec2f& a, const vec2f& b, float t) noexcept;
	TFM_INLINE vec2f max(const vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE vec2f min(const vec2f& a, const vec2f& b) noexcept;
	TFM_INLINE vec2f clamp(const vec2f& v, const vec2f& minv, const vec2f& maxv) noexcept;
	TFM_INLINE std::ostream& operator<<(std::ostream& os, const vec2f& vec) noexcept;

	struct TFM_ALIGNAS(16) vec3i {
		union {
			struct { int c0, c1, c2; };
			int col[3];
			struct { int x, y, z; };
			struct { int r, g, b; };
		};

		static TFM_INLINE const vec3i zero() noexcept;
		static TFM_INLINE const vec3i one() noexcept;

		TFM_INLINE vec3i() noexcept = default;
		TFM_INLINE vec3i(int x, int y, int z) noexcept;
		TFM_INLINE vec3i(int x, int y) noexcept;
		TFM_INLINE vec3i(int xyz) noexcept;
		TFM_INLINE vec3i(const int* ptr) noexcept;

		TFM_INLINE int& operator[](std::size_t i) noexcept;
		TFM_INLINE const int& operator[](std::size_t i) const noexcept;
	};

	TFM_INLINE bool operator==(const vec3i& a, const vec3i& b) noexcept;
	TFM_INLINE bool operator!=(const vec3i& a, const vec3i& b) noexcept;
	TFM_INLINE vec3i abs(const vec3i& v) noexcept;
	TFM_INLINE vec3i max(const vec3i& a, const vec3i& b) noexcept;
	TFM_INLINE vec3i min(const vec3i& a, const vec3i& b) noexcept;
	TFM_INLINE vec3i clamp(const vec3i& v, const vec3i& minv, const vec3i& maxv) noexcept;
	TFM_INLINE int dot(const vec3i& a, const vec3i& b) noexcept;
	TFM_INLINE int manhattan_distance(const vec3i& a, const vec3i& b) noexcept;

	struct TFM_ALIGNAS(16) vec4i {
		union {
			struct { int c0, c1, c2, c3; };
			int col[4];
			struct { int x, y, z, w; };
			struct { int r, g, b, a; };
		};

		static TFM_INLINE const vec4i zero() noexcept;
		static TFM_INLINE const vec4i one() noexcept;

		TFM_INLINE vec4i() noexcept = default;
		TFM_INLINE vec4i(int x, int y, int z, int w) noexcept;
		TFM_INLINE vec4i(int x, int y, int z) noexcept;
		TFM_INLINE vec4i(int x, int y) noexcept;
		TFM_INLINE vec4i(int xyzw) noexcept;
		TFM_INLINE vec4i(const int* ptr) noexcept;

		TFM_INLINE int& operator[](std::size_t i) noexcept;
		TFM_INLINE const int& operator[](std::size_t i) const noexcept;
	};

	TFM_INLINE bool operator==(const vec4i& a, const vec4i& b) noexcept;
	TFM_INLINE bool operator!=(const vec4i& a, const vec4i& b) noexcept;
	TFM_INLINE vec4i abs(const vec4i& v) noexcept;
	TFM_INLINE vec4i max(const vec4i& a, const vec4i& b) noexcept;
	TFM_INLINE vec4i min(const vec4i& a, const vec4i& b) noexcept;
	TFM_INLINE vec4i clamp(const vec4i& v, const vec4i& minv, const vec4i& maxv) noexcept;
	TFM_INLINE int dot(const vec4i& a, const vec4i& b) noexcept;
	TFM_INLINE int manhattan_distance(const vec4i& a, const vec4i& b) noexcept;

	struct TFM_ALIGNAS(16) vec3f {
		union {
			struct { float c0, c1, c2; };
			float col[3];
			struct { float x, y, z; };
			struct { float r, g, b; };
		};

		static TFM_INLINE const vec3f zero() noexcept;
		static TFM_INLINE const vec3f one() noexcept;
		static TFM_INLINE const vec3f unit_x() noexcept;
		static TFM_INLINE const vec3f unit_y() noexcept;
		static TFM_INLINE const vec3f unit_z() noexcept;
		static TFM_INLINE const vec3f up() noexcept;
		static TFM_INLINE const vec3f down() noexcept;
		static TFM_INLINE const vec3f right() noexcept;
		static TFM_INLINE const vec3f left() noexcept;
		static TFM_INLINE const vec3f forward() noexcept;
		static TFM_INLINE const vec3f backward() noexcept;

		TFM_INLINE vec3f() noexcept = default;
		TFM_INLINE vec3f(float x, float y, float z) noexcept;
		TFM_INLINE vec3f(float x, float y) noexcept;
		TFM_INLINE vec3f(float xyz) noexcept;
		TFM_INLINE vec3f(const float* ptr) noexcept;
		TFM_INLINE vec3f(const double* ptr) noexcept;
		TFM_INLINE vec3f(const std::vector<float>& v) noexcept;
		TFM_INLINE vec3f(const std::vector<double>& v) noexcept;
		TFM_INLINE vec3f(const vec2f& v) noexcept;
		TFM_INLINE vec3f(const vec2f& v, float z) noexcept;
		vec3f(const vec3f& v) noexcept = default;
		vec3f& operator=(const vec3f& v) noexcept = default;
		TFM_INLINE vec3f(const vec4f& v) noexcept;

		TFM_INLINE float length() const noexcept;
		TFM_INLINE float length_squared() const noexcept;
		TFM_INLINE vec3f normalize() const noexcept;
		TFM_INLINE vec3f round() const noexcept;
		TFM_INLINE vec3f ceil() const noexcept;

		static TFM_INLINE vec3f transform(float x, float y, float z, const mat4f& m) noexcept;
		static TFM_INLINE vec3f transform(const float* ptr, const mat4f& m) noexcept;
		static TFM_INLINE vec3f transform(const double* ptr, const mat4f& m) noexcept;
		static TFM_INLINE vec3f transform(const vec3f& v, const mat4f& m) noexcept;
		static TFM_INLINE vec3f transform_normal(float x, float y, float z, const mat4f& m) noexcept;
		static TFM_INLINE vec3f transform_normal(float* ptr, const mat4f& m) noexcept;
		static TFM_INLINE vec3f transform_normal(double* ptr, const mat4f& m) noexcept;
		static TFM_INLINE vec3f transform_normal(const vec3f& v, const mat4f& m) noexcept;

		TFM_INLINE float& operator[](std::size_t i) noexcept;
		TFM_INLINE const float& operator[](std::size_t i) const noexcept;
	};

	TFM_INLINE bool operator==(const vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE bool operator!=(const vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f& operator+=(vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f& operator+=(vec3f& a, float b) noexcept;
	TFM_INLINE vec3f& operator-=(vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f& operator-=(vec3f& a, float b) noexcept;
	TFM_INLINE vec3f& operator*=(vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f& operator*=(vec3f& a, float b) noexcept;
	TFM_INLINE vec3f& operator/=(vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f& operator/=(vec3f& a, float b) noexcept;
	TFM_INLINE vec3f operator+(const vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f operator+(const vec3f& a, float b) noexcept;
	TFM_INLINE vec3f operator+(float a, const vec3f& b) noexcept;
	TFM_INLINE vec3f operator-(const vec3f& a) noexcept;
	TFM_INLINE vec3f operator-(const vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f operator-(const vec3f& a, float b) noexcept;
	TFM_INLINE vec3f operator-(float a, const vec3f& b) noexcept;
	TFM_INLINE vec3f operator*(const vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f operator*(const vec3f& a, float b) noexcept;
	TFM_INLINE vec3f operator*(float a, const vec3f& b) noexcept;
	TFM_INLINE vec3f operator*(const quatf& q, const vec3f& v) noexcept;
	TFM_INLINE vec3f operator/(const vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f operator/(const vec3f& a, float b) noexcept;
	TFM_INLINE vec3f operator/(float a, const vec3f& b) noexcept;
	TFM_INLINE float dot(const vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE float distance(const vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE float distance_squared(const vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f cross(const vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f rotate(const quatf& q, const vec3f& v) noexcept;
	TFM_INLINE vec3f project(const vec3f& a, const vec3f& onto) noexcept;
	TFM_INLINE vec3f reflect(const vec3f& v, const vec3f& normal) noexcept;
	TFM_INLINE vec3f lerp(const vec3f& a, const vec3f& b, float t) noexcept;
	TFM_INLINE vec3f max(const vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f min(const vec3f& a, const vec3f& b) noexcept;
	TFM_INLINE vec3f clamp(const vec3f& v, const vec3f& minv, const vec3f& maxv) noexcept;
	TFM_INLINE std::ostream& operator<<(std::ostream& os, const vec3f& vec) noexcept;

	struct TFM_ALIGNAS(16) vec4f {
		union {
			struct { float c0, c1, c2, c3; };
			float col[4];
			struct { float x, y, z, w; };
			struct { float r, g, b, a; };
		};

		static TFM_INLINE const vec4f zero() noexcept;
		static TFM_INLINE const vec4f one() noexcept;
		static TFM_INLINE const vec4f unit_x() noexcept;
		static TFM_INLINE const vec4f unit_y() noexcept;
		static TFM_INLINE const vec4f unit_z() noexcept;
		static TFM_INLINE const vec4f unit_w() noexcept;

		TFM_INLINE vec4f() noexcept = default;
		TFM_INLINE vec4f(float x, float y, float z, float w) noexcept;
		TFM_INLINE vec4f(float x, float y) noexcept;
		TFM_INLINE vec4f(float x, float y, float z) noexcept;
		TFM_INLINE vec4f(float xyzw) noexcept;
		TFM_INLINE vec4f(const float* ptr) noexcept;
		TFM_INLINE vec4f(const float* ptr, float w) noexcept;
		TFM_INLINE vec4f(const double* ptr) noexcept;
		TFM_INLINE vec4f(const double* ptr, float w) noexcept;
		TFM_INLINE vec4f(const std::vector<float>& v) noexcept;
		TFM_INLINE vec4f(const std::vector<float>& v, float w) noexcept;
		TFM_INLINE vec4f(const std::vector<double>& v) noexcept;
		TFM_INLINE vec4f(const std::vector<double>& v, float w) noexcept;
		TFM_INLINE vec4f(const vec2f& v) noexcept;
		TFM_INLINE vec4f(const vec2f& v, float z, float w) noexcept;
		TFM_INLINE vec4f(const vec3f& v) noexcept;
		TFM_INLINE vec4f(const vec3f& v, float w) noexcept;
		vec4f(const vec4f& v) noexcept = default;
		vec4f& operator=(const vec4f& v) noexcept = default;

		TFM_INLINE float length() const noexcept;
		TFM_INLINE float length_squared() const noexcept;
		TFM_INLINE vec4f normalize() const noexcept;
		TFM_INLINE vec4f round() const noexcept;
		TFM_INLINE vec4f ceil() const noexcept;

		static TFM_INLINE vec4f transform(const vec4f& v, const mat4f& m) noexcept;
		static TFM_INLINE vec4f transform(float x, float y, float z, float w, const mat4f& m) noexcept;
		static TFM_INLINE vec4f transform(const float* ptr, const mat4f& m) noexcept;
		static TFM_INLINE vec4f transform(const double* ptr, const mat4f& m) noexcept;

		TFM_INLINE float& operator[](std::size_t i) noexcept;
		TFM_INLINE const float& operator[](std::size_t i) const noexcept;
	};

	TFM_INLINE bool operator==(const vec4f& a, const vec4f& b) noexcept;
	TFM_INLINE bool operator!=(const vec4f& a, const vec4f& b) noexcept;
	TFM_INLINE vec4f& operator+=(vec4f& a, vec4f b) noexcept;
	TFM_INLINE vec4f& operator+=(vec4f& a, float b) noexcept;
	TFM_INLINE vec4f& operator-=(vec4f& a, vec4f b) noexcept;
	TFM_INLINE vec4f& operator-=(vec4f& a, float b) noexcept;
	TFM_INLINE vec4f& operator*=(vec4f& a, vec4f b) noexcept;
	TFM_INLINE vec4f& operator*=(vec4f& a, float b) noexcept;
	TFM_INLINE vec4f& operator/=(vec4f& a, vec4f b) noexcept;
	TFM_INLINE vec4f& operator/=(vec4f& a, float b) noexcept;
	TFM_INLINE vec4f operator+(const vec4f& a, const vec4f& b) noexcept;
	TFM_INLINE vec4f operator+(const vec4f& a, float b) noexcept;
	TFM_INLINE vec4f operator+(float a, const vec4f& b) noexcept;
	TFM_INLINE vec4f operator-(const vec4f& a) noexcept;
	TFM_INLINE vec4f operator-(const vec4f& a, const vec4f& b) noexcept;
	TFM_INLINE vec4f operator-(const vec4f& a, float b) noexcept;
	TFM_INLINE vec4f operator-(float a, const vec4f& b) noexcept;
	TFM_INLINE vec4f operator*(const vec4f& a, const vec4f& b) noexcept;
	TFM_INLINE vec4f operator*(const vec4f& a, float b) noexcept;
	TFM_INLINE vec4f operator*(float a, const vec4f& b) noexcept;
	TFM_INLINE vec4f operator/(const vec4f& a, const vec4f& b) noexcept;
	TFM_INLINE vec4f operator/(const vec4f& a, float b) noexcept;
	TFM_INLINE vec4f operator/(float a, const vec4f& b) noexcept;
	TFM_INLINE vec4f operator*(const vec4f& vec, const mat4f& m) noexcept;
	TFM_INLINE vec4f max(const vec4f& a, const vec4f& b) noexcept;
	TFM_INLINE vec4f min(const vec4f& a, const vec4f& b) noexcept;
	TFM_INLINE vec4f clamp(const vec4f& v, const vec4f& minv, const vec4f& maxv) noexcept;
	TFM_INLINE float dot(const vec4f& a, const vec4f& b) noexcept;
	TFM_INLINE float distance(const vec4f& a, const vec4f& b) noexcept;
	TFM_INLINE float distance_squared(const vec4f& a, const vec4f& b) noexcept;
	TFM_INLINE vec4f cross(const vec4f& a, const vec4f& b) noexcept;
	TFM_INLINE vec4f project(const vec4f& a, const vec4f& onto) noexcept;
	TFM_INLINE vec4f reflect(const vec4f& v, const vec4f& normal) noexcept;
	TFM_INLINE vec4f lerp(const vec4f& a, const vec4f& b, float t) noexcept;
	TFM_INLINE vec4f mix(const vec4f& a, const vec4f& b, float t) noexcept;
	TFM_INLINE std::ostream& operator<<(std::ostream& os, const vec4f& vec) noexcept;

	struct TFM_ALIGNAS(16) quatf {
		union {
			struct { float c0, c1, c2, c3; };
			float col[4];
			struct { float x, y, z, w; };
			struct { float r, g, b, a; };
		};

		static TFM_INLINE const quatf identity() noexcept;
		static TFM_INLINE const quatf zero() noexcept;
		static TFM_INLINE const quatf one() noexcept;

		TFM_INLINE quatf() noexcept = default;
		TFM_INLINE quatf(float x, float y, float z, float w) noexcept;
		TFM_INLINE quatf(float xyz) noexcept;
		TFM_INLINE quatf(const float* ptr) noexcept;
		TFM_INLINE quatf(const double* ptr) noexcept;
		TFM_INLINE quatf(const vec2f& v, float z, float w) noexcept;
		TFM_INLINE quatf(const vec3f& v, float w) noexcept;
		TFM_INLINE quatf(const vec4f& v) noexcept;
		quatf(const quatf& v) noexcept = default;
		quatf& operator=(const quatf& v) noexcept = default;

		TFM_INLINE float length() const noexcept;
		TFM_INLINE float length_squared() const noexcept;
		TFM_INLINE quatf normalize() const noexcept;
		TFM_INLINE quatf round() const noexcept;
		TFM_INLINE quatf ceil() const noexcept;
		TFM_INLINE quatf conjugate() const noexcept;
		TFM_INLINE quatf inverse() const noexcept;

		static TFM_INLINE const quatf angle_axis(float angle, float x, float y, float z) noexcept;
		static TFM_INLINE const quatf angle_axis(float angle, const vec3f& v) noexcept;
		static TFM_INLINE const quatf angle_axis(float angle, const vec2f& v, float z) noexcept;
		static TFM_INLINE const quatf angle_axis(float angle, const float* ptr) noexcept;
		static TFM_INLINE const quatf angle_axis(float angle, const float* ptr, float z) noexcept;
		static TFM_INLINE const quatf angle_axis(float angle, const double* ptr) noexcept;
		static TFM_INLINE const quatf angle_axis(float angle, const double* ptr, float z) noexcept;
		static TFM_INLINE const quatf from_euler(float pitch, float yaw, float roll) noexcept;
		static TFM_INLINE const quatf from_euler(const vec3f& angles) noexcept;

		TFM_INLINE float& operator[](std::size_t i) noexcept;
		TFM_INLINE const float& operator[](std::size_t i) const noexcept;
	};

	TFM_INLINE bool operator==(const quatf& a, const quatf& b) noexcept;
	TFM_INLINE bool operator!=(const quatf& a, const quatf& b) noexcept;
	TFM_INLINE quatf& operator+=(quatf& a, quatf b) noexcept;
	TFM_INLINE quatf& operator+=(quatf& a, float b) noexcept;
	TFM_INLINE quatf& operator-=(quatf& a, quatf b) noexcept;
	TFM_INLINE quatf& operator-=(quatf& a, float b) noexcept;
	TFM_INLINE quatf& operator*=(quatf& a, quatf b) noexcept;
	TFM_INLINE quatf& operator*=(quatf& a, float b) noexcept;
	TFM_INLINE quatf& operator/=(quatf& a, quatf b) noexcept;
	TFM_INLINE quatf& operator/=(quatf& a, float b) noexcept;
	TFM_INLINE quatf operator+(const quatf& a, const quatf& b) noexcept;
	TFM_INLINE quatf operator+(const quatf& a, float b) noexcept;
	TFM_INLINE quatf operator+(float a, const quatf& b) noexcept;
	TFM_INLINE quatf operator-(const quatf& a) noexcept;
	TFM_INLINE quatf operator-(const quatf& a, const quatf& b) noexcept;
	TFM_INLINE quatf operator-(const quatf& a, float b) noexcept;
	TFM_INLINE quatf operator-(float a, const quatf& b) noexcept;
	TFM_INLINE quatf operator*(const quatf& a, const quatf& b) noexcept;
	TFM_INLINE quatf operator*(const quatf& a, float b) noexcept;
	TFM_INLINE quatf operator*(float a, const quatf& b) noexcept;
	TFM_INLINE quatf operator/(const quatf& a, const quatf& b) noexcept;
	TFM_INLINE quatf operator/(const quatf& a, float b) noexcept;
	TFM_INLINE quatf operator/(float a, const quatf& b) noexcept;
	TFM_INLINE quatf max(const quatf& a, const quatf& b) noexcept;
	TFM_INLINE quatf min(const quatf& a, const quatf& b) noexcept;
	TFM_INLINE quatf clamp(const quatf& v, const quatf& minv, const quatf& maxv) noexcept;
	TFM_INLINE float dot(const quatf& a, const quatf& b) noexcept;
	TFM_INLINE quatf lerp(const quatf& a, const quatf& b, float t) noexcept;
	TFM_INLINE quatf mix(const quatf& a, const quatf& b, float t) noexcept;
	TFM_INLINE quatf slerp(quatf q1, quatf q2, float t);
	TFM_INLINE std::ostream& operator<<(std::ostream& os, const quatf& vec) noexcept;

	struct TFM_ALIGNAS(16) mat3f {
		union {
			struct { float m00, m01, m02, _pad0, m10, m11, m12, _pad1, m20, m21, m22, _pad2; };
			float m[12];
			float col[3][4];
			vec4f row[3];
		};

		static TFM_INLINE const mat3f zero() noexcept;
		static TFM_INLINE const mat3f identity() noexcept;

		TFM_INLINE mat3f() noexcept = default;
		mat3f(const mat3f& m) noexcept = default;
		mat3f& operator=(const mat3f& m) noexcept = default;
		TFM_INLINE mat3f(
			float m00, float m01, float m02,
			float m10, float m11, float m12,
			float m20, float m21, float m22) noexcept;
		TFM_INLINE mat3f(const float* ptr) noexcept;
		TFM_INLINE mat3f(const double* ptr) noexcept;
		TFM_INLINE mat3f(const vec3f& r0, const vec3f& r1, const vec3f& r2) noexcept;
		TFM_INLINE mat3f(const vec4f& r0, const vec4f& r1, const vec4f& r2) noexcept;
		TFM_INLINE mat3f(const quatf& q) noexcept;
		TFM_INLINE mat3f(const mat4f& m) noexcept;

		static TFM_INLINE mat3f rotation(const quatf& q) noexcept;
		TFM_INLINE vec3f transform(const vec3f& v) const noexcept;
		TFM_INLINE mat3f inverse() const noexcept;
		TFM_INLINE mat3f transpose() const noexcept;
		TFM_INLINE mat3f inverse_transpose() const noexcept;

		TFM_INLINE vec4f& operator[](std::size_t i) noexcept;
		TFM_INLINE const vec4f& operator[](std::size_t i) const noexcept;
	};

	struct TFM_ALIGNAS(16) mat4f {
		union {
			struct { float m00, m01, m02, m03, m10, m11, m12, m13, m20, m21, m22, m23, m30, m31, m32, m33; };
			float m[16];
			float col[4][4];
			vec4f row[4];
		};

		static TFM_INLINE const mat4f zero() noexcept;
		static TFM_INLINE const mat4f identity() noexcept;

		TFM_INLINE mat4f() noexcept = default;
		mat4f(const mat4f& m) noexcept = default;
		mat4f& operator=(const mat4f& m) noexcept = default;
		TFM_INLINE mat4f(
			float m00, float m01, float m02, float m03,
			float m10, float m11, float m12, float m13,
			float m20, float m21, float m22, float m23,
			float m30, float m31, float m32, float m33);
		TFM_INLINE mat4f(const float* ptr) noexcept;
		TFM_INLINE mat4f(const double* ptr) noexcept;
		TFM_INLINE mat4f(const quatf& q) noexcept;
		TFM_INLINE mat4f(const vec4f& r0, const vec4f& r1, const vec4f& r2, const vec4f& r3) noexcept;
		TFM_INLINE mat4f(const mat3f& m) noexcept;

		template <typename T> static TFM_INLINE mat4f transformation(const std::vector<T>& matrix, const std::vector<T>& translation, const std::vector<T>& rotation, const std::vector<T>& scaling);

		static TFM_INLINE mat4f translation(float x, float y, float z) noexcept;
		static TFM_INLINE mat4f translation(const vec3f& v) noexcept;

		TFM_INLINE mat4f translate(float x, float y, float z) const noexcept;
		TFM_INLINE mat4f translate(const vec3f& v) const noexcept;

		static TFM_INLINE mat4f scaling(float x, float y, float z) noexcept;
		static TFM_INLINE mat4f scaling(const vec3f& v) noexcept;
		static TFM_INLINE mat4f scaling(float xyz) noexcept;

		TFM_INLINE mat4f scale(float x, float y, float z) const noexcept;
		TFM_INLINE mat4f scale(const vec3f& v) const noexcept;
		TFM_INLINE mat4f scale(float xyz) const noexcept;

		static TFM_INLINE mat4f rotation(float angle, const vec3f& vec) noexcept;
		static TFM_INLINE mat4f rotation(float angle, float rx, float ry, float rz) noexcept;
		static TFM_INLINE mat4f rotation_x(float angle) noexcept;
		static TFM_INLINE mat4f rotation_y(float angle) noexcept;
		static TFM_INLINE mat4f rotation_z(float angle) noexcept;
		static TFM_INLINE mat4f rotation(const quatf& q) noexcept;

		TFM_INLINE mat4f rotate(float angle, const vec3f& axis) const noexcept;
		TFM_INLINE mat4f rotate(float angle, float rx, float ry, float rz) const noexcept;
		TFM_INLINE mat4f rotate_x(float angle) const noexcept;
		TFM_INLINE mat4f rotate_y(float angle) const noexcept;
		TFM_INLINE mat4f rotate_z(float angle) const noexcept;
		TFM_INLINE mat4f rotate(const quatf& q) const noexcept;

		static TFM_INLINE mat4f perspective_fov_rh_zo(float fovy, float aspect, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f perspective_rh_zo(float left, float right, float bottom, float top, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f perspective_rh_zo(float width, float height, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f perspective_fov_rh_no(float fovy, float aspect, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f perspective_rh_no(float left, float right, float bottom, float top, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f perspective_rh_no(float width, float height, float znear, float zfar) noexcept;

		static TFM_INLINE mat4f perspective_fov_lh_zo(float fovy, float aspect, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f perspective_lh_zo(float left, float right, float bottom, float top, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f perspective_lh_zo(float width, float height, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f perspective_fov_lh_no(float fovy, float aspect, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f perspective_lh_no(float left, float right, float bottom, float top, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f perspective_lh_no(float width, float height, float znear, float zfar) noexcept;

		static TFM_INLINE mat4f frustum_zo_rh(float left, float right, float bottom, float top, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f frustum_zo_lh(float left, float right, float bottom, float top, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f frustum_no_rh(float left, float right, float bottom, float top, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f frustum_no_lh(float left, float right, float bottom, float top, float znear, float zfar) noexcept;

		static TFM_INLINE mat4f ortho_rh_zo(float left, float right, float bottom, float top, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f ortho_rh_zo(float width, float height, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f ortho_rh_no(float left, float right, float bottom, float top, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f ortho_rh_no(float width, float height, float znear, float zfar) noexcept;

		static TFM_INLINE mat4f ortho_lh_zo(float left, float right, float bottom, float top, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f ortho_lh_zo(float width, float height, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f ortho_lh_no(float left, float right, float bottom, float top, float znear, float zfar) noexcept;
		static TFM_INLINE mat4f ortho_lh_no(float width, float height, float znear, float zfar) noexcept;

		static TFM_INLINE mat4f look_to_rh(const vec3f& eye, const vec3f& dir, const vec3f& up) noexcept;
		static TFM_INLINE mat4f look_at_rh(const vec3f& eye, const vec3f& target, const vec3f& up) noexcept;

		static TFM_INLINE mat4f look_to_lh(const vec3f& eye, const vec3f& dir, const vec3f& up) noexcept;
		static TFM_INLINE mat4f look_at_lh(const vec3f& eye, const vec3f& target, const vec3f& up) noexcept;

		TFM_INLINE vec3f translation() const noexcept;
		TFM_INLINE vec3f scaling() const noexcept;
		TFM_INLINE mat3f rotation_matrix() const noexcept;
		TFM_INLINE vec3f transform_point(const vec3f& v) const noexcept;
		TFM_INLINE vec3f transform_vector(const vec3f& v) const noexcept;
		TFM_INLINE mat4f inverse() const noexcept;
		TFM_INLINE mat4f transpose() const noexcept;
		TFM_INLINE mat4f inverse_transpose() const noexcept;

		TFM_INLINE vec4f& operator[](std::size_t i) noexcept;
		TFM_INLINE const vec4f& operator[](std::size_t i) const noexcept;
	};

#if TFM_MATRIX_RTL
	TFM_INLINE mat4f operator*(const mat4f& rhs, const mat4f& lhs) noexcept;
#else
	TFM_INLINE mat4f operator*(const mat4f& lhs, const mat4f& rhs) noexcept;
#endif
	TFM_INLINE mat4f& operator*=(mat4f& a, const mat4f& b) noexcept;
	TFM_INLINE bool operator==(const mat4f& a, const mat4f& b) noexcept;
	TFM_INLINE bool operator!=(const mat4f& a, const mat4f& b) noexcept;
	TFM_INLINE std::ostream& operator<<(std::ostream& os, const mat4f& mat) noexcept;

	/*********/
	/* vec2i */
	/*********/

	TFM_INLINE const vec2i vec2i::zero() noexcept { return vec2i{}; }
	TFM_INLINE const vec2i vec2i::one() noexcept { return vec2i(1); }

	TFM_INLINE vec2i::vec2i(int xyz) noexcept : x(xyz), y(xyz) {}
	TFM_INLINE vec2i::vec2i(int x, int y) noexcept : x(x), y(y) {}
	TFM_INLINE vec2i::vec2i(const int* ptr) noexcept : x(ptr[0]), y(ptr[1]) {}

	TFM_INLINE int& vec2i::operator[](std::size_t i) noexcept {
		return col[i];
	}

	TFM_INLINE const int& vec2i::operator[](std::size_t i) const noexcept {
		return col[i];
	}

	TFM_INLINE bool operator==(const vec2i& a, const vec2i& b) noexcept {
		return a.x == b.x && a.y == b.y;
	}

	TFM_INLINE bool operator!=(const vec2i& a, const vec2i& b) noexcept {
		return !(a == b);
	}

	TFM_INLINE vec2i abs(const vec2i& v) noexcept {
		return vec2i(v.x < 0 ? -v.x : v.x, v.y < 0 ? -v.y : v.y);
	}

	TFM_INLINE vec2i max(const vec2i& a, const vec2i& b) noexcept {
		return vec2i(std::max(a.x, b.x), std::max(a.y, b.y));
	}

	TFM_INLINE vec2i min(const vec2i& a, const vec2i& b) noexcept {
		return vec2i(std::min(a.x, b.x), std::min(a.y, b.y));
	}

	TFM_INLINE vec2i clamp(const vec2i& v, const vec2i& minv, const vec2i& maxv) noexcept {
		return min(max(v, minv), maxv);
	}

	TFM_INLINE int dot(const vec2i& a, const vec2i& b) noexcept {
		return a.x * b.x + a.y * b.y;
	}

	TFM_INLINE int manhattan_distance(const vec2i& a, const vec2i& b) noexcept {
		int dx = a.x - b.x;
		int dy = a.y - b.y;
		return (dx < 0 ? -dx : dx) + (dy < 0 ? -dy : dy);
	}

	/*********/
	/* vec2f */
	/*********/

	TFM_INLINE const vec2f vec2f::zero() noexcept { return { 0.f, 0.f }; }
	TFM_INLINE const vec2f vec2f::one() noexcept { return { 1.f, 1.f }; }
	TFM_INLINE const vec2f vec2f::unit_x() noexcept { return { 1.f, 0.f }; }
	TFM_INLINE const vec2f vec2f::unit_y() noexcept { return { 0.f, 1.f }; }

	TFM_INLINE vec2f::vec2f(float xy) noexcept : x(xy), y(xy) {}
	TFM_INLINE vec2f::vec2f(float x, float y) noexcept : x(x), y(y) {}
	TFM_INLINE vec2f::vec2f(const float* ptr) noexcept : x(ptr[0]), y(ptr[1]) {}
	TFM_INLINE vec2f::vec2f(const double* ptr) noexcept : x(ptr[0]), y(ptr[1]) {}
	TFM_INLINE vec2f::vec2f(const std::vector<float>& v) noexcept : x(v[0]), y(v[1]) {}
	TFM_INLINE vec2f::vec2f(const std::vector<double>& v) noexcept : x(v[0]), y(v[1]) {}
	TFM_INLINE vec2f::vec2f(const vec3f& v) noexcept : x(v.x), y(v.y) {}
	TFM_INLINE vec2f::vec2f(const vec4f& v) noexcept : x(v.x), y(v.y) {}

	TFM_INLINE bool operator==(const vec2f& a, const vec2f& b) noexcept {
		for (int i = 0; i < 2; i++) {
			if (std::fabs(a.col[i] - b.col[i]) > 1e-3f) {
				return false;
			}
		}
		return true;
	}

	TFM_INLINE bool operator!=(const vec2f& a, const vec2f& b) noexcept { return !(a == b); }

	TFM_INLINE vec2f& operator+=(vec2f& a, const vec2f& b) noexcept {
		a.x += b.x;
		a.y += b.y;
		return a;
	}

	TFM_INLINE vec2f& operator+=(vec2f& a, float b) noexcept {
		a.x += b;
		a.y += b;
		return a;
	}

	TFM_INLINE vec2f& operator-=(vec2f& a, const vec2f& b) noexcept {
		a.x -= b.x;
		a.y -= b.y;
		return a;
	}

	TFM_INLINE vec2f& operator-=(vec2f& a, float b) noexcept {
		a.x -= b;
		a.y -= b;
		return a;
	}

	TFM_INLINE vec2f& operator*=(vec2f& a, const vec2f& b) noexcept {
		a.x *= b.x;
		a.y *= b.y;
		return a;
	}

	TFM_INLINE vec2f& operator*=(vec2f& a, float b) noexcept {
		a.x *= b;
		a.y *= b;
		return a;
	}

	TFM_INLINE vec2f& operator/=(vec2f& a, const vec2f& b) noexcept {
		a.x /= b.x;
		a.y /= b.y;
		return a;
	}

	TFM_INLINE vec2f& operator/=(vec2f& a, float b) noexcept {
		a.x /= b;
		a.y /= b;
		return a;
	}

	TFM_INLINE vec2f operator+(const vec2f& a, const vec2f& b) noexcept {
		return vec2f(a.x + b.x, a.y + b.y);
	}

	TFM_INLINE vec2f operator+(const vec2f& a, float b) noexcept {
		return vec2f(a.x + b, a.y + b);
	}

	TFM_INLINE vec2f operator+(float a, const vec2f& b) noexcept {
		return vec2f(a + b.x, a + b.y);
	}

	TFM_INLINE vec2f operator-(const vec2f& a) noexcept {
		return vec2f(-a.x, -a.y);
	}

	TFM_INLINE vec2f operator-(const vec2f& a, const vec2f& b) noexcept {
		return vec2f(a.x - b.x, a.y - b.y);
	}

	TFM_INLINE vec2f operator-(const vec2f& a, float b) noexcept {
		return vec2f(a.x - b, a.y - b);
	}

	TFM_INLINE vec2f operator-(float a, const vec2f& b) noexcept {
		return vec2f(a - b.x, a - b.y);
	}

	TFM_INLINE vec2f operator*(const vec2f& a, const vec2f& b) noexcept {
		return vec2f(a.x * b.x, a.y * b.y);
	}

	TFM_INLINE vec2f operator*(const vec2f& a, float b) noexcept {
		return vec2f(a.x * b, a.y * b);
	}

	TFM_INLINE vec2f operator*(float a, const vec2f& b) noexcept {
		return vec2f(a * b.x, a * b.y);
	}

	TFM_INLINE vec2f operator/(const vec2f& a, const vec2f& b) noexcept {
		return vec2f(a.x / b.x, a.y / b.y);
	}

	TFM_INLINE vec2f operator/(const vec2f& a, float b) noexcept {
		float inv = 1.0f / b;
		return vec2f(a.x * inv, a.y * inv);
	}

	TFM_INLINE vec2f operator/(float a, const vec2f& b) noexcept {
		return vec2f(a / b.x, a / b.y);
	}

	TFM_INLINE float dot(const vec2f& a, const vec2f& b) noexcept {
		return a.x * b.x + a.y * b.y;
	}

	TFM_INLINE float vec2f::length_squared() const noexcept {
		return dot(*this, *this);
	}

	TFM_INLINE float vec2f::length() const noexcept {
		return std::sqrt(length_squared());
	}

	TFM_INLINE float& vec2f::operator[](std::size_t i) noexcept {
		return col[i];
	}

	TFM_INLINE const float& vec2f::operator[](std::size_t i) const noexcept {
		return col[i];
	}

	TFM_INLINE vec2f vec2f::normalize() const noexcept {
		float len = length();
		if (len > 0.0f) {
			float invLen = 1.0f / len;
			return *this * invLen;
		}
		return vec2f::zero();
	}

	TFM_INLINE vec2f vec2f::round() const noexcept {
		return vec2f(std::round(x), std::round(y));
	}

	TFM_INLINE vec2f vec2f::ceil() const noexcept {
		return vec2f(std::ceil(x), std::ceil(y));
	}

	TFM_INLINE float distance(const vec2f& a, const vec2f& b) noexcept {
		return (a - b).length();
	}

	TFM_INLINE float distance_squared(const vec2f& a, const vec2f& b) noexcept {
		return (a - b).length_squared();
	}

	TFM_INLINE float cross(const vec2f& a, const vec2f& b) noexcept {
		return a.x * b.y - a.y * b.x;
	}

	TFM_INLINE vec2f project(const vec2f& a, const vec2f& onto) noexcept {
		float lenSq = onto.length_squared();
		return lenSq > 0.0f ? onto * (dot(a, onto) / lenSq) : vec2f::zero();
	}

	TFM_INLINE vec2f reflect(const vec2f& v, const vec2f& normal) noexcept {
		return v - 2.0f * dot(v, normal) * normal;
	}

	TFM_INLINE vec2f lerp(const vec2f& a, const vec2f& b, float t) noexcept {
		return a + (b - a) * vec2f(t);
	}

	TFM_INLINE vec2f max(const vec2f& a, const vec2f& b) noexcept {
		return vec2f(std::max(a.x, b.x), std::max(a.y, b.y));
	}

	TFM_INLINE vec2f min(const vec2f& a, const vec2f& b) noexcept {
		return vec2f(std::min(a.x, b.x), std::min(a.y, b.y));
	}

	TFM_INLINE vec2f clamp(const vec2f& v, const vec2f& minv, const vec2f& maxv) noexcept {
		return vec2f(std::max(minv.x, std::min(v.x, maxv.x)), std::max(minv.y, std::min(v.y, maxv.y)));
	}

	TFM_INLINE std::ostream& operator<<(std::ostream& os, const vec2f& vec) noexcept {
		for (int j = 0; j < 2; j++) {
			os << std::setw(8) << std::setprecision(4) << vec.col[j] << " ";
		}
		return os;
	}

	/*********/
	/* vec3i */
	/*********/

	TFM_INLINE const vec3i vec3i::zero() noexcept { return vec3i{}; }
	TFM_INLINE const vec3i vec3i::one() noexcept { return vec3i(1); }

	TFM_INLINE vec3i::vec3i(int x, int y, int z) noexcept : x(x), y(y), z(z) {}
	TFM_INLINE vec3i::vec3i(int x, int y) noexcept : x(x), y(y), z(0) {}
	TFM_INLINE vec3i::vec3i(int xyz) noexcept : x(xyz), y(xyz), z(xyz) {}
	TFM_INLINE vec3i::vec3i(const int* ptr) noexcept : x(ptr[0]), y(ptr[1]), z(ptr[2]) {}

	TFM_INLINE int& vec3i::operator[](std::size_t i) noexcept {
		return col[i];
	}

	TFM_INLINE const int& vec3i::operator[](std::size_t i) const noexcept {
		return col[i];
	}

	TFM_INLINE bool operator==(const vec3i& a, const vec3i& b) noexcept {
		return a.x == b.x && a.y == b.y && a.z == b.z;
	}

	TFM_INLINE bool operator!=(const vec3i& a, const vec3i& b) noexcept {
		return !(a == b);
	}

	TFM_INLINE vec3i abs(const vec3i& v) noexcept {
		return vec3i(v.x < 0 ? -v.x : v.x, v.y < 0 ? -v.y : v.y, v.z < 0 ? -v.z : v.z);
	}

	TFM_INLINE vec3i max(const vec3i& a, const vec3i& b) noexcept {
		return vec3i(std::max(a.x, b.x), std::max(a.y, b.y), std::max(a.z, b.z));
	}

	TFM_INLINE vec3i min(const vec3i& a, const vec3i& b) noexcept {
		return vec3i(std::min(a.x, b.x), std::min(a.y, b.y), std::min(a.z, b.z));
	}

	TFM_INLINE vec3i clamp(const vec3i& v, const vec3i& minv, const vec3i& maxv) noexcept {
		return min(max(v, minv), maxv);
	}

	TFM_INLINE int dot(const vec3i& a, const vec3i& b) noexcept {
		return a.x * b.x + a.y * b.y + a.z * b.z;
	}

	TFM_INLINE int manhattan_distance(const vec3i& a, const vec3i& b) noexcept {
		int dx = a.x - b.x;
		int dy = a.y - b.y;
		int dz = a.z - b.z;
		return (dx < 0 ? -dx : dx) + (dy < 0 ? -dy : dy) + (dz < 0 ? -dz : dz);
	}

	/*********/
	/* vec4i */
	/*********/

	TFM_INLINE const vec4i vec4i::zero() noexcept { return vec4i{}; }
	TFM_INLINE const vec4i vec4i::one() noexcept { return vec4i(1); }

	TFM_INLINE vec4i::vec4i(int x, int y, int z, int w) noexcept : x(x), y(y), z(z), w(w) {}
	TFM_INLINE vec4i::vec4i(int x, int y, int z) noexcept : x(x), y(y), z(z), w(0) {}
	TFM_INLINE vec4i::vec4i(int x, int y) noexcept : x(x), y(y), z(0), w(0) {}
	TFM_INLINE vec4i::vec4i(int xyzw) noexcept : x(xyzw), y(xyzw), z(xyzw), w(xyzw) {}
	TFM_INLINE vec4i::vec4i(const int* ptr) noexcept : x(ptr[0]), y(ptr[1]), z(ptr[2]), w(ptr[3]) {}

	TFM_INLINE int& vec4i::operator[](std::size_t i) noexcept {
		return col[i];
	}

	TFM_INLINE const int& vec4i::operator[](std::size_t i) const noexcept {
		return col[i];
	}

	TFM_INLINE bool operator==(const vec4i& a, const vec4i& b) noexcept {
		return a.x == b.x && a.y == b.y && a.z == b.z && a.w == b.w;
	}

	TFM_INLINE bool operator!=(const vec4i& a, const vec4i& b) noexcept {
		return !(a == b);
	}

	TFM_INLINE vec4i abs(const vec4i& v) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		vec4i result;
		TFM_MM_STORE_SI128((__m128i*)result.col, _mm_abs_epi32(TFM_MM_LOAD_SI128((const __m128i*)v.col)));
		return result;
#else
		return vec4i(v.x < 0 ? -v.x : v.x, v.y < 0 ? -v.y : v.y, v.z < 0 ? -v.z : v.z, v.w < 0 ? -v.w : v.w);
#endif
	}

	TFM_INLINE vec4i max(const vec4i& a, const vec4i& b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		vec4i result;
		TFM_MM_STORE_SI128((__m128i*)result.col, _mm_max_epi32(TFM_MM_LOAD_SI128((const __m128i*)a.col), TFM_MM_LOAD_SI128((const __m128i*)b.col)));
		return result;
#else
		return vec4i(std::max(a.x, b.x), std::max(a.y, b.y), std::max(a.z, b.z), std::max(a.w, b.w));
#endif
	}

	TFM_INLINE vec4i min(const vec4i& a, const vec4i& b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		vec4i result;
		TFM_MM_STORE_SI128((__m128i*)result.col, _mm_min_epi32(TFM_MM_LOAD_SI128((const __m128i*)a.col), TFM_MM_LOAD_SI128((const __m128i*)b.col)));
		return result;
#else
		return vec4i(std::min(a.x, b.x), std::min(a.y, b.y), std::min(a.z, b.z), std::min(a.w, b.w));
#endif
	}

	TFM_INLINE vec4i clamp(const vec4i& v, const vec4i& minv, const vec4i& maxv) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		vec4i result;
		const __m128i vv = TFM_MM_LOAD_SI128((const __m128i*)v.col);
		const __m128i vmin = TFM_MM_LOAD_SI128((const __m128i*)minv.col);
		const __m128i vmax = TFM_MM_LOAD_SI128((const __m128i*)maxv.col);
		TFM_MM_STORE_SI128((__m128i*)result.col, _mm_min_epi32(_mm_max_epi32(vv, vmin), vmax));
		return result;
#else
		return min(max(v, minv), maxv);
#endif
	}

	TFM_INLINE int dot(const vec4i& a, const vec4i& b) noexcept {
		return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
	}

	TFM_INLINE int manhattan_distance(const vec4i& a, const vec4i& b) noexcept {
		int dx = a.x - b.x;
		int dy = a.y - b.y;
		int dz = a.z - b.z;
		int dw = a.w - b.w;
		return (dx < 0 ? -dx : dx) + (dy < 0 ? -dy : dy) + (dz < 0 ? -dz : dz) + (dw < 0 ? -dw : dw);
	}

	/********/
	/* vec3f */
	/********/

	TFM_INLINE const vec3f vec3f::zero() noexcept { return { 0.f, 0.f, 0.f }; }
	TFM_INLINE const vec3f vec3f::one() noexcept { return { 1.f, 1.f, 1.f }; }
	TFM_INLINE const vec3f vec3f::unit_x() noexcept { return { 1.f, 0.f, 0.f }; }
	TFM_INLINE const vec3f vec3f::unit_y() noexcept { return { 0.f, 1.f, 0.f }; }
	TFM_INLINE const vec3f vec3f::unit_z() noexcept { return { 0.f, 0.f, 1.f }; }
	TFM_INLINE const vec3f vec3f::up() noexcept { return { 0.f, 1.f, 0.f }; }
	TFM_INLINE const vec3f vec3f::down() noexcept { return { 0.f, -1.f, 0.f }; }
	TFM_INLINE const vec3f vec3f::right() noexcept { return { 1.f, 0.f, 0.f }; }
	TFM_INLINE const vec3f vec3f::left() noexcept { return { -1.f, 0.f, 0.f }; }
	TFM_INLINE const vec3f vec3f::forward() noexcept { return { 0.f, 0.f, -1.f }; }
	TFM_INLINE const vec3f vec3f::backward() noexcept { return { 0.f, 0.f, 1.f }; }

	TFM_INLINE vec3f::vec3f(float xyz) noexcept : x(xyz), y(xyz), z(xyz) {}
	TFM_INLINE vec3f::vec3f(float x, float y, float z) noexcept : x(x), y(y), z(z) {}
	TFM_INLINE vec3f::vec3f(float x, float y) noexcept : x(x), y(y), z(0.0f) {}
	TFM_INLINE vec3f::vec3f(const float* ptr) noexcept : x(ptr[0]), y(ptr[1]), z(ptr[2]) {}
	TFM_INLINE vec3f::vec3f(const double* ptr) noexcept : x(ptr[0]), y(ptr[1]), z(ptr[2]) {}
	TFM_INLINE vec3f::vec3f(const std::vector<float>& v) noexcept : x(v[0]), y(v[1]), z(v[2]) {}
	TFM_INLINE vec3f::vec3f(const std::vector<double>& v) noexcept : x(v[0]), y(v[1]), z(v[2]) {}
	TFM_INLINE vec3f::vec3f(const vec2f& v) noexcept : x(v.x), y(v.y), z(0) {}
	TFM_INLINE vec3f::vec3f(const vec2f& v, float z) noexcept : x(v.x), y(v.y), z(z) {}
	TFM_INLINE vec3f::vec3f(const vec4f& v) noexcept : x(v.x), y(v.y), z(v.z) {}

	TFM_INLINE bool operator==(const vec3f& a, const vec3f& b) noexcept {
		for (int i = 0; i < 3; i++) {
			if (std::fabs(a.col[i] - b.col[i]) > 1e-3f) {
				return false;
			}
		}
		return true;
	}

	TFM_INLINE bool operator!=(const vec3f& a, const vec3f& b) noexcept {
		return !(a == b);
	}

	TFM_INLINE vec3f& operator+=(vec3f& a, const vec3f& b) noexcept {
		a.x += b.x;
		a.y += b.y;
		a.z += b.z;
		return a;
	}

	TFM_INLINE vec3f& operator+=(vec3f& a, float b) noexcept {
		a.x += b;
		a.y += b;
		a.z += b;
		return a;
	}

	TFM_INLINE vec3f& operator-=(vec3f& a, const vec3f& b) noexcept {
		a.x -= b.x;
		a.y -= b.y;
		a.z -= b.z;
		return a;
	}

	TFM_INLINE vec3f& operator-=(vec3f& a, float b) noexcept {
		a.x -= b;
		a.y -= b;
		a.z -= b;
		return a;
	}

	TFM_INLINE vec3f& operator*=(vec3f& a, const vec3f& b) noexcept {
		a.x *= b.x;
		a.y *= b.y;
		a.z *= b.z;
		return a;
	}

	TFM_INLINE vec3f& operator*=(vec3f& a, float b) noexcept {
		a.x *= b;
		a.y *= b;
		a.z *= b;
		return a;
	}

	TFM_INLINE vec3f& operator/=(vec3f& a, const vec3f& b) noexcept {
		a.x /= b.x;
		a.y /= b.y;
		a.z /= b.z;
		return a;
	}

	TFM_INLINE vec3f& operator/=(vec3f& a, float b) noexcept {
		a.x /= b;
		a.y /= b;
		a.z /= b;
		return a;
	}

	TFM_INLINE vec3f operator+(const vec3f& a, const vec3f& b) noexcept {
		return vec3f(a.x + b.x, a.y + b.y, a.z + b.z);
	}

	TFM_INLINE vec3f operator+(const vec3f& a, float b) noexcept {
		return vec3f(a.x + b, a.y + b, a.z + b);
	}

	TFM_INLINE vec3f operator+(float a, const vec3f& b) noexcept {
		return vec3f(a + b.x, a + b.y, a + b.z);
	}

	TFM_INLINE vec3f operator-(const vec3f& a) noexcept {
		return vec3f(-a.x, -a.y, -a.z);
	}

	TFM_INLINE vec3f operator-(const vec3f& a, const vec3f& b) noexcept {
		return vec3f(a.x - b.x, a.y - b.y, a.z - b.z);
	}

	TFM_INLINE vec3f operator-(const vec3f& a, float b) noexcept {
		return vec3f(a.x - b, a.y - b, a.z - b);
	}

	TFM_INLINE vec3f operator-(float a, const vec3f& b) noexcept {
		return vec3f(a - b.x, a - b.y, a - b.z);
	}

	TFM_INLINE vec3f operator*(const vec3f& a, const vec3f& b) noexcept {
		return vec3f(a.x * b.x, a.y * b.y, a.z * b.z);
	}

	TFM_INLINE vec3f operator*(const vec3f& a, float b) noexcept {
		return vec3f(a.x * b, a.y * b, a.z * b);
	}

	TFM_INLINE vec3f operator*(float a, const vec3f& b) noexcept {
		return vec3f(a * b.x, a * b.y, a * b.z);
	}

	TFM_INLINE vec3f operator*(const quatf& q, const vec3f& v) noexcept {
		vec3f qv(q.x, q.y, q.z);

		vec3f t = 2.0f * cross(qv, v);
		return v + q.w * t + cross(qv, t);
	}

	TFM_INLINE vec3f rotate(const quatf& q, const vec3f& v) noexcept {
		return q * v;
	}

	TFM_INLINE vec3f operator/(const vec3f& a, const vec3f& b) noexcept {
		return vec3f(a.x / b.x, a.y / b.y, a.z / b.z);
	}

	TFM_INLINE vec3f operator/(const vec3f& a, float b) noexcept {
		float inv = 1.0f / b;
		return vec3f(a.x * inv, a.y * inv, a.z * inv);
	}

	TFM_INLINE vec3f operator/(float a, const vec3f& b) noexcept {
		return vec3f(a / b.x, a / b.y, a / b.z);
	}

	TFM_INLINE vec3f operator*(const vec3f& v, const tfm::mat4f& m) {
		return vec3f::transform(v, m);
	}

	TFM_INLINE float dot(const vec3f& a, const vec3f& b) noexcept {
		return a.x * b.x + a.y * b.y + a.z * b.z;
	}

	TFM_INLINE float vec3f::length_squared() const noexcept {
		return dot(*this, *this);
	}

	TFM_INLINE float vec3f::length() const noexcept {
		return std::sqrt(length_squared());
	}

	TFM_INLINE float& vec3f::operator[](std::size_t i) noexcept {
		return col[i];
	}

	TFM_INLINE const float& vec3f::operator[](std::size_t i) const noexcept {
		return col[i];
	}

	TFM_INLINE vec3f vec3f::normalize() const noexcept {
		float len = length();
		if (len > 0.0f) {
			float invLen = 1.0f / len;
			return *this * invLen;
		}
		return vec3f::zero();
	}

	TFM_INLINE vec3f vec3f::round() const noexcept {
		return vec3f(std::round(x), std::round(y), std::round(z));
	}

	TFM_INLINE vec3f vec3f::ceil() const noexcept {
		return vec3f(std::ceil(x), std::ceil(y), std::ceil(z));
	}

	TFM_INLINE vec3f vec3f::transform(float x, float y, float z, const mat4f& m) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		__m128 r = _mm_mul_ps(TFM_MM_LOAD_PS(m.row[0].col), _mm_set1_ps(x));
		r = _mm_fmadd_ps(TFM_MM_LOAD_PS(m.row[1].col), _mm_set1_ps(y), r);
		r = _mm_fmadd_ps(TFM_MM_LOAD_PS(m.row[2].col), _mm_set1_ps(z), r);
		r = _mm_add_ps(r, TFM_MM_LOAD_PS(m.row[3].col));
		r = _mm_div_ps(r, _mm_shuffle_ps(r, r, _MM_SHUFFLE(3, 3, 3, 3)));
		return vec3f(
			_mm_cvtss_f32(r),
			_mm_cvtss_f32(_mm_shuffle_ps(r, r, _MM_SHUFFLE(1, 1, 1, 1))),
			_mm_cvtss_f32(_mm_movehl_ps(r, r))
		);
#else
		vec4f r = vec4f::transform(x, y, z, 1, m);

		r = r / r.w;

		return vec3f(r.x, r.y, r.z);
#endif
	}

	TFM_INLINE vec3f vec3f::transform(const float* ptr, const mat4f& m) noexcept {
		return transform(ptr[0], ptr[1], ptr[2], m);
	}

	TFM_INLINE vec3f vec3f::transform(const double* ptr, const mat4f& m) noexcept {
		return transform(ptr[0], ptr[1], ptr[2], m);
	}

	TFM_INLINE vec3f vec3f::transform(const vec3f& v, const mat4f& m) noexcept {
		return transform(v.x, v.y, v.z, m);
	}

	TFM_INLINE vec3f vec3f::transform_normal(float x, float y, float z, const mat4f& m) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		__m128 r = _mm_mul_ps(TFM_MM_LOAD_PS(m.row[0].col), _mm_set1_ps(x));
		r = _mm_fmadd_ps(TFM_MM_LOAD_PS(m.row[1].col), _mm_set1_ps(y), r);
		r = _mm_fmadd_ps(TFM_MM_LOAD_PS(m.row[2].col), _mm_set1_ps(z), r);
		return vec3f(
			_mm_cvtss_f32(r),
			_mm_cvtss_f32(_mm_shuffle_ps(r, r, _MM_SHUFFLE(1, 1, 1, 1))),
			_mm_cvtss_f32(_mm_movehl_ps(r, r))
		);
#else
		vec3f result;
		result.x = x * m.m00 + y * m.m10 + z * m.m20;
		result.y = x * m.m01 + y * m.m11 + z * m.m21;
		result.z = x * m.m02 + y * m.m12 + z * m.m22;
		return result;
#endif
	}

	TFM_INLINE vec3f vec3f::transform_normal(float* ptr, const mat4f& m) noexcept {
		return transform_normal(ptr[0], ptr[1], ptr[2], m);
	}

	TFM_INLINE vec3f vec3f::transform_normal(double* ptr, const mat4f& m) noexcept {
		return transform_normal(ptr[0], ptr[1], ptr[2], m);
	}

	TFM_INLINE vec3f vec3f::transform_normal(const vec3f& vec, const mat4f& m) noexcept {
		return transform_normal(vec.x, vec.y, vec.z, m);
	}

	TFM_INLINE float distance(const vec3f& a, const vec3f& b) noexcept {
		return (a - b).length();
	}

	TFM_INLINE float distance_squared(const vec3f& a, const vec3f& b) noexcept {
		return (a - b).length_squared();
	}

	TFM_INLINE vec3f cross(const vec3f& a, const vec3f& b) noexcept {
		return vec3f(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
	}

	TFM_INLINE vec3f project(const vec3f& a, const vec3f& onto) noexcept {
		float lenSq = onto.length_squared();
		return lenSq > 0.0f ? onto * (dot(a, onto) / lenSq) : vec3f::zero();
	}

	TFM_INLINE vec3f reflect(const vec3f& v, const vec3f& normal) noexcept {
		return v - 2.0f * dot(v, normal) * normal;
	}

	TFM_INLINE vec3f lerp(const vec3f& a, const vec3f& b, float t) noexcept {
		return a + (b - a) * vec3f(t);
	}

	TFM_INLINE vec3f max(const vec3f& a, const vec3f& b) noexcept {
		return vec3f(std::max(a.x, b.x), std::max(a.y, b.y), std::max(a.z, b.z));
	}

	TFM_INLINE vec3f min(const vec3f& a, const vec3f& b) noexcept {
		return vec3f(std::min(a.x, b.x), std::min(a.y, b.y), std::min(a.z, b.z));
	}

	TFM_INLINE vec3f clamp(const vec3f& v, const vec3f& minv, const vec3f& maxv) noexcept {
		return vec3f(std::max(minv.x, std::min(v.x, maxv.x)), std::max(minv.y, std::min(v.y, maxv.y)), std::max(minv.z, std::min(v.z, maxv.z)));
	}

	TFM_INLINE std::ostream& operator<<(std::ostream& os, const vec3f& vec) noexcept {
		for (int j = 0; j < 3; j++) {
			os << std::setw(8) << std::setprecision(4) << vec.col[j] << " ";
		}
		return os;
	}

	/********/
	/* vec4f */
	/********/

	TFM_INLINE const vec4f vec4f::zero() noexcept { return { 0.f, 0.f, 0.f, 0.f }; }
	TFM_INLINE const vec4f vec4f::one() noexcept { return { 1.f, 1.f, 1.f, 1.f }; }
	TFM_INLINE const vec4f vec4f::unit_x() noexcept { return { 1.f, 0.f, 0.f, 0.f }; }
	TFM_INLINE const vec4f vec4f::unit_y() noexcept { return { 0.f, 1.f, 0.f, 0.f }; }
	TFM_INLINE const vec4f vec4f::unit_z() noexcept { return { 0.f, 0.f, 1.f, 0.f }; }
	TFM_INLINE const vec4f vec4f::unit_w() noexcept { return { 0.f, 0.f, 0.f, 1.f }; }

	TFM_INLINE vec4f::vec4f(float xyzw) noexcept : x(xyzw), y(xyzw), z(xyzw), w(xyzw) {}
	TFM_INLINE vec4f::vec4f(float x, float y) noexcept : x(x), y(y), z(0.0f), w(0.0f) {}
	TFM_INLINE vec4f::vec4f(float x, float y, float z) noexcept : x(x), y(y), z(z), w(0.0f) {}
	TFM_INLINE vec4f::vec4f(float x, float y, float z, float w) noexcept : x(x), y(y), z(z), w(w) {}

	TFM_INLINE vec4f::vec4f(const float* ptr) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		TFM_MM_STORE_PS(col, _mm_loadu_ps(ptr));
#elif TFM_NEON_INTRINSICS
		vst1q_f32(col, vld1q_f32(ptr));
#else
		x = ptr[0]; y = ptr[1]; z = ptr[2]; w = ptr[3];
#endif
	}

	TFM_INLINE vec4f::vec4f(const float* ptr, float w) noexcept : x(ptr[0]), y(ptr[1]), z(ptr[2]), w(w) {}

	TFM_INLINE vec4f::vec4f(const double* ptr) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 lo = _mm_cvtpd_ps(_mm_loadu_pd(ptr));
		const __m128 hi = _mm_cvtpd_ps(_mm_loadu_pd(ptr + 2));
		TFM_MM_STORE_PS(col, _mm_movelh_ps(lo, hi));
#else
		x = (float)ptr[0]; y = (float)ptr[1]; z = (float)ptr[2]; w = (float)ptr[3];
#endif
	}

	TFM_INLINE vec4f::vec4f(const double* ptr, float w) noexcept : x(ptr[0]), y(ptr[1]), z(ptr[2]), w(w) {}
	TFM_INLINE vec4f::vec4f(const std::vector<float>& v) noexcept : vec4f(v.data()) {}
	TFM_INLINE vec4f::vec4f(const std::vector<float>& v, float w) noexcept : x(v[0]), y(v[1]), z(v[2]), w(w) {}
	TFM_INLINE vec4f::vec4f(const std::vector<double>& v) noexcept : vec4f(v.data()) {}
	TFM_INLINE vec4f::vec4f(const std::vector<double>& v, float w) noexcept : x(v[0]), y(v[1]), z(v[2]), w(w) {}
	TFM_INLINE vec4f::vec4f(const vec2f& v) noexcept : x(v.x), y(v.y), z(0), w(0) {}
	TFM_INLINE vec4f::vec4f(const vec2f& v, float z, float w) noexcept : x(v.x), y(v.y), z(z), w(w) {}
	TFM_INLINE vec4f::vec4f(const vec3f& v) noexcept : x(v.x), y(v.y), z(v.z), w(0) {}
	TFM_INLINE vec4f::vec4f(const vec3f& v, float w) noexcept : x(v.x), y(v.y), z(v.z), w(w) {}

	TFM_INLINE bool operator==(const vec4f& a, const vec4f& b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		__m128 diff = _mm_sub_ps(va, vb);
		__m128 absd = _mm_andnot_ps(_mm_set1_ps(-0.0f), diff);
		__m128 cmp = _mm_cmple_ps(absd, _mm_set1_ps(1e-3f));
		return _mm_movemask_ps(cmp) == 0xF;
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		float32x4_t diff = vsubq_f32(va, vb);
		float32x4_t absd = vabsq_f32(diff);
		uint32x4_t cmp = vcleq_f32(absd, vdupq_n_f32(1e-3f));

		uint64x2_t tmp = vreinterpretq_u64_u32(cmp);
		return (vgetq_lane_u64(tmp, 0) & vgetq_lane_u64(tmp, 1)) == 0xFFFFFFFFFFFFFFFFull;
#else
		for (int i = 0; i < 4; i++) {
			if (std::fabs(a.col[i] - b.col[i]) > 1e-3f) {
				return false;
			}
		}
		return true;
#endif
	}

	TFM_INLINE bool operator!=(const vec4f& a, const vec4f& b) noexcept {
		return !(a == b);
	}

	TFM_INLINE vec4f& operator+=(vec4f& a, vec4f b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(a.col, _mm_add_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(a.col, vaddq_f32(va, vb));
#else
		a.x += b.x;
		a.y += b.y;
		a.z += b.z;
		a.w += b.w;
#endif
		return a;
	}

	TFM_INLINE vec4f& operator+=(vec4f& a, float b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(a.col, _mm_add_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(a.col, vaddq_f32(va, vdupq_n_f32(b)));
#else
		a.x += b;
		a.y += b;
		a.z += b;
		a.w += b;
#endif
		return a;
	}

	TFM_INLINE vec4f& operator-=(vec4f& a, vec4f b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(a.col, _mm_sub_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(a.col, vsubq_f32(va, vb));
#else
		a.x -= b.x;
		a.y -= b.y;
		a.z -= b.z;
		a.w -= b.w;
#endif
		return a;
	}

	TFM_INLINE vec4f& operator-=(vec4f& a, float b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(a.col, _mm_sub_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(a.col, vsubq_f32(va, vdupq_n_f32(b)));
#else
		a.x -= b;
		a.y -= b;
		a.z -= b;
		a.w -= b;
#endif
		return a;
	}

	TFM_INLINE vec4f& operator*=(vec4f& a, vec4f b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(a.col, _mm_mul_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(a.col, vmulq_f32(va, vb));
#else
		a.x *= b.x;
		a.y *= b.y;
		a.z *= b.z;
		a.w *= b.w;
#endif
		return a;
	}

	TFM_INLINE vec4f& operator*=(vec4f& a, float b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(a.col, _mm_mul_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(a.col, vmulq_f32(va, vdupq_n_f32(b)));
#else
		a.x *= b;
		a.y *= b;
		a.z *= b;
		a.w *= b;
#endif
		return a;
	}

	TFM_INLINE vec4f& operator/=(vec4f& a, vec4f b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(a.col, _mm_div_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		float32x4_t vb = vld1q_f32(b.col);
		float32x4_t recip = vrecpeq_f32(vb);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		vst1q_f32(a.col, vmulq_f32(va, recip));
#else
		a.x /= b.x;
		a.y /= b.y;
		a.z /= b.z;
		a.w /= b.w;
#endif
		return a;
	}

	TFM_INLINE vec4f& operator/=(vec4f& a, float b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(a.col, _mm_div_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		float32x4_t vb = vdupq_n_f32(b);
		float32x4_t recip = vrecpeq_f32(vb);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		vst1q_f32(a.col, vmulq_f32(va, recip));
#else
		a.x /= b;
		a.y /= b;
		a.z /= b;
		a.w /= b;
#endif
		return a;
	}

	TFM_INLINE vec4f operator+(const vec4f& a, const vec4f& b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_add_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vaddq_f32(va, vb));
#else
		result = vec4f(a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w);
#endif
		return result;
	}

	TFM_INLINE vec4f operator+(const vec4f& a, float b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(result.col, _mm_add_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(result.col, vaddq_f32(va, vdupq_n_f32(b)));
#else
		result = vec4f(a.x + b, a.y + b, a.z + b, a.w + b);
#endif
		return result;
	}

	TFM_INLINE vec4f operator+(float a, const vec4f& b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_add_ps(_mm_set1_ps(a), vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vaddq_f32(vdupq_n_f32(a), vb));
#else
		result = vec4f(a + b.x, a + b.y, a + b.z, a + b.w);
#endif
		return result;
	}

	TFM_INLINE vec4f operator-(const vec4f& a) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(result.col, _mm_xor_ps(va, _mm_set1_ps(-0.0f)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(result.col, vnegq_f32(va));
#else
		result = vec4f(-a.x, -a.y, -a.z, -a.w);
#endif
		return result;
	}

	TFM_INLINE vec4f operator-(const vec4f& a, const vec4f& b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_sub_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vsubq_f32(va, vb));
#else
		result = vec4f(a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w);
#endif
		return result;
	}

	TFM_INLINE vec4f operator-(const vec4f& a, float b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(result.col, _mm_sub_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(result.col, vsubq_f32(va, vdupq_n_f32(b)));
#else
		result = vec4f(a.x - b, a.y - b, a.z - b, a.w - b);
#endif
		return result;
	}

	TFM_INLINE vec4f operator-(float a, const vec4f& b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_sub_ps(_mm_set1_ps(a), vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vsubq_f32(vdupq_n_f32(a), vb));
#else
		result = vec4f(a - b.x, a - b.y, a - b.z, a - b.w);
#endif
		return result;
	}

	TFM_INLINE vec4f operator*(const vec4f& a, const vec4f& b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_mul_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vmulq_f32(va, vb));
#else
		result = vec4f(a.x * b.x, a.y * b.y, a.z * b.z, a.w * b.w);
#endif
		return result;
	}

	TFM_INLINE vec4f operator*(const vec4f& a, float b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(result.col, _mm_mul_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(result.col, vmulq_f32(va, vdupq_n_f32(b)));
#else
		result = vec4f(a.x * b, a.y * b, a.z * b, a.w * b);
#endif
		return result;
	}

	TFM_INLINE vec4f operator*(float a, const vec4f& b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_mul_ps(_mm_set1_ps(a), vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vmulq_f32(vdupq_n_f32(a), vb));
#else
		result = vec4f(a * b.x, a * b.y, a * b.z, a * b.w);
#endif
		return result;
	}

	TFM_INLINE vec4f operator/(const vec4f& a, const vec4f& b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_div_ps(va, vb));
#else
		result = vec4f(a.x / b.x, a.y / b.y, a.z / b.z, a.w / b.w);
#endif
		return result;
	}

	TFM_INLINE vec4f operator/(const vec4f& a, float b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(result.col, _mm_div_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		float32x4_t vb = vdupq_n_f32(b);
		float32x4_t recip = vrecpeq_f32(vb);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		vst1q_f32(result.col, vmulq_f32(va, recip));
#else
		result = vec4f(a.x / b, a.y / b, a.z / b, a.w / b);
#endif
		return result;
	}

	TFM_INLINE vec4f operator/(float a, const vec4f& b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_div_ps(_mm_set1_ps(a), vb));
#elif TFM_NEON_INTRINSICS
		float32x4_t va = vdupq_n_f32(a);
		float32x4_t vb = vld1q_f32(b.col);
		float32x4_t recip = vrecpeq_f32(vb);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		vst1q_f32(result.col, vmulq_f32(va, recip));
#else
		result = vec4f(a / b.x, a / b.y, a / b.z, a / b.w);
#endif
		return result;
	}

	TFM_INLINE vec4f operator*(const vec4f& vec, const mat4f& m) noexcept {
		return vec4f::transform(vec, m);
	}

	TFM_INLINE float dot(const vec4f& a, const vec4f& b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		return _mm_cvtss_f32(_mm_dot_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		float32x4_t m = vmulq_f32(va, vb);
		float32x2_t sum = vadd_f32(vget_low_f32(m), vget_high_f32(m));
		sum = vpadd_f32(sum, sum);
		return vget_lane_f32(sum, 0);
#else
		return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
#endif
	}

	TFM_INLINE float vec4f::length_squared() const noexcept {
		return dot(*this, *this);
	}

	TFM_INLINE float vec4f::length() const noexcept {
#if TFM_SSE_AVX_INTRINSICS
		__m128 c = TFM_MM_LOAD_PS(col);
		return _mm_cvtss_f32(_mm_sqrt_ps(_mm_dot_ps(c, c)));
#elif TFM_NEON_INTRINSICS
		float32x4_t c = vld1q_f32(col);
		float32x4_t m = vmulq_f32(c, c);
		float32x2_t sum = vadd_f32(vget_low_f32(m), vget_high_f32(m));
		sum = vpadd_f32(sum, sum);
		float len2 = vget_lane_f32(sum, 0);
		return sqrtf(len2);
#else
		return sqrtf(dot(*this, *this));
#endif
	}

	TFM_INLINE float& vec4f::operator[](std::size_t i) noexcept {
		return col[i];
	}

	TFM_INLINE const float& vec4f::operator[](std::size_t i) const noexcept {
		return col[i];
	}

	TFM_INLINE vec4f vec4f::normalize() const noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 c = TFM_MM_LOAD_PS(col);
		TFM_MM_STORE_PS(result.col, _mm_normalize_ps(c));
#elif TFM_NEON_INTRINSICS
		float32x4_t c = vld1q_f32(col);
		float32x4_t m = vmulq_f32(c, c);
		float32x2_t sum = vadd_f32(vget_low_f32(m), vget_high_f32(m));
		sum = vpadd_f32(sum, sum);
		float32x4_t len2 = vdupq_n_f32(vget_lane_f32(sum, 0));
		float32x4_t invLen = vrsqrteq_f32(len2);
		invLen = vmulq_f32(invLen, vrsqrtsq_f32(vmulq_f32(len2, invLen), invLen));
		invLen = vmulq_f32(invLen, vrsqrtsq_f32(vmulq_f32(len2, invLen), invLen));
		vst1q_f32(result.col, vmulq_f32(c, invLen));
#else
		float len = length();
		result = (len > 0.0f) ? *this / len : vec4f::zero();
#endif
		return result;
	}

	TFM_INLINE vec4f vec4f::round() const noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 c = TFM_MM_LOAD_PS(col);
		TFM_MM_STORE_PS(result.col, _mm_round_ps(c, _MM_FROUND_TO_NEAREST_INT | _MM_FROUND_NO_EXC));
#elif TFM_NEON_INTRINSICS
		const float32x4_t c = vld1q_f32(col);
		vst1q_f32(result.col, vrndnq_f32(c));
#else
		result = vec4f(std::round(x), std::round(y), std::round(z), std::round(w));
#endif
		return result;
	}

	TFM_INLINE vec4f vec4f::ceil() const noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 c = TFM_MM_LOAD_PS(col);
		TFM_MM_STORE_PS(result.col, _mm_ceil_ps(c));
#elif TFM_NEON_INTRINSICS
		const float32x4_t c = vld1q_f32(col);
		vst1q_f32(result.col, vrndpq_f32(c));
#else
		result = vec4f(std::ceil(x), std::ceil(y), std::ceil(z), std::ceil(w));
#endif
		return result;
	}

	TFM_INLINE vec4f vec4f::transform(float x, float y, float z, float w, const mat4f& m) noexcept {
		return transform(vec4f(x, y, z, w), m);
	}

	TFM_INLINE vec4f vec4f::transform(const float* ptr, const mat4f& m) noexcept {
		return transform(vec4f(ptr), m);
	}

	TFM_INLINE vec4f vec4f::transform(const double* ptr, const mat4f& m) noexcept {
		return transform(vec4f(ptr), m);
	}

	TFM_INLINE vec4f vec4f::transform(const vec4f& vec, const mat4f& m) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 r0 = TFM_MM_LOAD_PS(m.row[0].col);
		const __m128 r1 = TFM_MM_LOAD_PS(m.row[1].col);
		const __m128 r2 = TFM_MM_LOAD_PS(m.row[2].col);
		const __m128 r3 = TFM_MM_LOAD_PS(m.row[3].col);
		__m128 r = _mm_mul_ps(r0, _mm_set1_ps(vec.x));
		r = _mm_fmadd_ps(r1, _mm_set1_ps(vec.y), r);
		r = _mm_fmadd_ps(r2, _mm_set1_ps(vec.z), r);
		TFM_MM_STORE_PS(result.col, _mm_fmadd_ps(r3, _mm_set1_ps(vec.w), r));
#elif TFM_NEON_INTRINSICS
		float32x4_t vc = vld1q_f32(vec.col);
		const float32x4_t mr0 = vld1q_f32(m.row[0].col);
		const float32x4_t mr1 = vld1q_f32(m.row[1].col);
		const float32x4_t mr2 = vld1q_f32(m.row[2].col);
		const float32x4_t mr3 = vld1q_f32(m.row[3].col);
		float32x4_t r0 = vmulq_laneq_f32(mr0, vc, 0);
		float32x4_t r1 = vmulq_laneq_f32(mr1, vc, 1);
		float32x4_t r2 = vmulq_laneq_f32(mr2, vc, 2);
		float32x4_t r3 = vmulq_laneq_f32(mr3, vc, 3);
		vst1q_f32(result.col, vaddq_f32(vaddq_f32(r0, r1), vaddq_f32(r2, r3)));
#else
		result.x = vec.x * m.m00 + vec.y * m.m10 + vec.z * m.m20 + vec.w * m.m30;
		result.y = vec.x * m.m01 + vec.y * m.m11 + vec.z * m.m21 + vec.w * m.m31;
		result.z = vec.x * m.m02 + vec.y * m.m12 + vec.z * m.m22 + vec.w * m.m32;
		result.w = vec.x * m.m03 + vec.y * m.m13 + vec.z * m.m23 + vec.w * m.m33;
#endif
		return result;
	}

	TFM_INLINE float distance(const vec4f& a, const vec4f& b) noexcept {
		return (a - b).length();
	}

	TFM_INLINE float distance_squared(const vec4f& a, const vec4f& b) noexcept {
		return (a - b).length_squared();
	}

	TFM_INLINE vec4f cross(const vec4f& a, const vec4f& b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		__m128 va = TFM_MM_LOAD_PS(a.col);
		__m128 vb = TFM_MM_LOAD_PS(b.col);
		__m128 shuf1 = _mm_shuffle_ps(va, va, _MM_SHUFFLE(3, 0, 2, 1));
		__m128 shuf2 = _mm_shuffle_ps(vb, vb, _MM_SHUFFLE(3, 0, 2, 1));
		__m128 c = _mm_sub_ps(_mm_mul_ps(va, shuf2), _mm_mul_ps(shuf1, vb));
		c = _mm_shuffle_ps(c, c, _MM_SHUFFLE(3, 0, 2, 1));

		__m128 mask = _mm_castsi128_ps(_mm_set_epi32(-1, -1, -1, 0));

		TFM_MM_STORE_PS(result.col, _mm_and_ps(c, mask));
#elif TFM_NEON_INTRINSICS
		float32x4_t va = vld1q_f32(a.col);
		float32x4_t vb = vld1q_f32(b.col);
		float32x4_t a_yzx = vextq_f32(va, va, 1);
		float32x4_t b_yzx = vextq_f32(vb, vb, 1);
		float32x4_t c = vsubq_f32(vmulq_f32(va, b_yzx), vmulq_f32(a_yzx, vb));
		c = vextq_f32(c, c, 1);
		c = vsetq_lane_f32(0.0f, c, 3);
		vst1q_f32(result.col, c);
#else
		result = vec4f(
			a.y * b.z - a.z * b.y,
			a.z * b.x - a.x * b.z,
			a.x * b.y - a.y * b.x,
			0.0f
		);
#endif
		return result;
	}

	TFM_INLINE vec4f project(const vec4f& a, const vec4f& onto) noexcept {
		float lenSq = onto.length_squared();
		return lenSq > 0.0f ? onto * (dot(a, onto) / lenSq) : vec4f::zero();
	}

	TFM_INLINE vec4f reflect(const vec4f& v, const vec4f& normal) noexcept {
		return v - 2.0f * dot(v, normal) * normal;
	}

	TFM_INLINE vec4f lerp(const vec4f& a, const vec4f& b, float t) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		__m128 va = TFM_MM_LOAD_PS(a.col);
		__m128 vb = TFM_MM_LOAD_PS(b.col);
		__m128 vt = _mm_set1_ps(t);
		__m128 diff = _mm_sub_ps(vb, va);
		TFM_MM_STORE_PS(result.col, _mm_fmadd_ps(diff, vt, va));
#elif TFM_NEON_INTRINSICS
		float32x4_t va = vld1q_f32(a.col);
		float32x4_t vb = vld1q_f32(b.col);
		float32x4_t diff = vsubq_f32(vb, va);
		vst1q_f32(result.col, vmlaq_n_f32(va, diff, t));
#else
		result = a + (b - a) * vec4f(t);
#endif
		return result;
	}

	TFM_INLINE vec4f mix(const vec4f& a, const vec4f& b, float t) noexcept {
		return lerp(a, b, t);
	}

	TFM_INLINE vec4f max(const vec4f& a, const vec4f& b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_max_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vmaxq_f32(va, vb));
#else
		result = vec4f(std::max(a.x, b.x), std::max(a.y, b.y), std::max(a.z, b.z), std::max(a.w, b.w));
#endif
		return result;
	}

	TFM_INLINE vec4f min(const vec4f& a, const vec4f& b) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_min_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vminq_f32(va, vb));
#else
		result = vec4f(std::min(a.x, b.x), std::min(a.y, b.y), std::min(a.z, b.z), std::min(a.w, b.w));
#endif
		return result;
	}

	TFM_INLINE vec4f clamp(const vec4f& v, const vec4f& minv, const vec4f& maxv) noexcept {
		vec4f result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 vv = TFM_MM_LOAD_PS(v.col);
		const __m128 vmin = TFM_MM_LOAD_PS(minv.col);
		const __m128 vmax = TFM_MM_LOAD_PS(maxv.col);
		TFM_MM_STORE_PS(result.col, _mm_min_ps(_mm_max_ps(vv, vmin), vmax));
#elif TFM_NEON_INTRINSICS
		const float32x4_t vv = vld1q_f32(v.col);
		const float32x4_t vmin = vld1q_f32(minv.col);
		const float32x4_t vmax = vld1q_f32(maxv.col);
		vst1q_f32(result.col, vminq_f32(vmaxq_f32(vv, vmin), vmax));
#else
		result = min(max(v, minv), maxv);
#endif
		return result;
	}

	TFM_INLINE std::ostream& operator<<(std::ostream& os, const vec4f& vec) noexcept {
		for (int j = 0; j < 4; j++) {
			os << std::setw(8) << std::setprecision(4) << vec.col[j] << " ";
		}
		return os;
	}

	/********/
	/* quatf */
	/********/

	TFM_INLINE const quatf quatf::identity() noexcept {
		return quatf(0.0f, 0.0f, 0.0f, 1.0f);
	}

	TFM_INLINE const quatf quatf::zero() noexcept {
		return quatf{};
	}

	TFM_INLINE const quatf quatf::one() noexcept {
		return quatf(1);
	}

	TFM_INLINE quatf::quatf(float xyz) noexcept : x(xyz), y(xyz), z(xyz), w(xyz) {}
	TFM_INLINE quatf::quatf(float x, float y, float z, float w) noexcept : x(x), y(y), z(z), w(w) {}

	TFM_INLINE quatf::quatf(const float* ptr) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		TFM_MM_STORE_PS(col, _mm_loadu_ps(ptr));
#elif TFM_NEON_INTRINSICS
		vst1q_f32(col, vld1q_f32(ptr));
#else
		x = ptr[0]; y = ptr[1]; z = ptr[2]; w = ptr[3];
#endif
	}

	TFM_INLINE quatf::quatf(const double* ptr) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 lo = _mm_cvtpd_ps(_mm_loadu_pd(ptr));
		const __m128 hi = _mm_cvtpd_ps(_mm_loadu_pd(ptr + 2));
		TFM_MM_STORE_PS(col, _mm_movelh_ps(lo, hi));
#else
		x = (float)ptr[0]; y = (float)ptr[1]; z = (float)ptr[2]; w = (float)ptr[3];
#endif
	}
	TFM_INLINE quatf::quatf(const vec2f& v, float z, float w) noexcept : x(v.x), y(v.y), z(z), w(w) {}
	TFM_INLINE quatf::quatf(const vec3f& v, float w) noexcept : x(v.x), y(v.y), z(v.z), w(w) {}
	TFM_INLINE quatf::quatf(const vec4f& v) noexcept : x(v.x), y(v.y), z(v.z), w(v.w) {}

	TFM_INLINE bool operator==(const quatf& a, const quatf& b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		__m128 diff = _mm_sub_ps(va, vb);
		__m128 absd = _mm_andnot_ps(_mm_set1_ps(-0.0f), diff);
		__m128 cmp = _mm_cmple_ps(absd, _mm_set1_ps(1e-3f));
		return _mm_movemask_ps(cmp) == 0xF;
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		float32x4_t diff = vsubq_f32(va, vb);
		float32x4_t absd = vabsq_f32(diff);
		uint32x4_t cmp = vcleq_f32(absd, vdupq_n_f32(1e-3f));

		uint64x2_t tmp = vreinterpretq_u64_u32(cmp);
		return (vgetq_lane_u64(tmp, 0) & vgetq_lane_u64(tmp, 1)) == 0xFFFFFFFFFFFFFFFFull;
#else
		for (int i = 0; i < 4; i++) {
			if (std::fabs(a.col[i] - b.col[i]) > 1e-3f) {
				return false;
			}
		}
		return true;
#endif
	}

	TFM_INLINE bool operator!=(const quatf& a, const quatf& b) noexcept {
		return !(a == b);
	}

	TFM_INLINE quatf& operator+=(quatf& a, quatf b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(a.col, _mm_add_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(a.col, vaddq_f32(va, vb));
#else
		a.x += b.x;
		a.y += b.y;
		a.z += b.z;
		a.w += b.w;
#endif
		return a;
	}

	TFM_INLINE quatf& operator+=(quatf& a, float b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(a.col, _mm_add_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(a.col, vaddq_f32(va, vdupq_n_f32(b)));
#else
		a.x += b;
		a.y += b;
		a.z += b;
		a.w += b;
#endif
		return a;
	}

	TFM_INLINE quatf& operator-=(quatf& a, quatf b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(a.col, _mm_sub_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(a.col, vsubq_f32(va, vb));
#else
		a.x -= b.x;
		a.y -= b.y;
		a.z -= b.z;
		a.w -= b.w;
#endif
		return a;
	}

	TFM_INLINE quatf& operator-=(quatf& a, float b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(a.col, _mm_sub_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(a.col, vsubq_f32(va, vdupq_n_f32(b)));
#else
		a.x -= b;
		a.y -= b;
		a.z -= b;
		a.w -= b;
#endif
		return a;
	}

	TFM_INLINE quatf& operator*=(quatf& a, quatf b) noexcept {
		a = a * b;
		return a;
	}

	TFM_INLINE quatf& operator*=(quatf& a, float b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(a.col, _mm_mul_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(a.col, vmulq_f32(va, vdupq_n_f32(b)));
#else
		a.x *= b;
		a.y *= b;
		a.z *= b;
		a.w *= b;
#endif
		return a;
	}

	TFM_INLINE quatf& operator/=(quatf& a, quatf b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(a.col, _mm_div_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		float32x4_t vb = vld1q_f32(b.col);
		float32x4_t recip = vrecpeq_f32(vb);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		vst1q_f32(a.col, vmulq_f32(va, recip));
#else
		a.x /= b.x;
		a.y /= b.y;
		a.z /= b.z;
		a.w /= b.w;
#endif
		return a;
	}

	TFM_INLINE quatf& operator/=(quatf& a, float b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(a.col, _mm_div_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		float32x4_t vb = vdupq_n_f32(b);
		float32x4_t recip = vrecpeq_f32(vb);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		vst1q_f32(a.col, vmulq_f32(va, recip));
#else
		a.x /= b;
		a.y /= b;
		a.z /= b;
		a.w /= b;
#endif
		return a;
	}

	TFM_INLINE quatf operator+(const quatf& a, const quatf& b) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_add_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vaddq_f32(va, vb));
#else
		result = quatf(a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w);
#endif
		return result;
	}

	TFM_INLINE quatf operator+(const quatf& a, float b) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(result.col, _mm_add_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(result.col, vaddq_f32(va, vdupq_n_f32(b)));
#else
		result = quatf(a.x + b, a.y + b, a.z + b, a.w + b);
#endif
		return result;
	}

	TFM_INLINE quatf operator+(float a, const quatf& b) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_add_ps(_mm_set1_ps(a), vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vaddq_f32(vdupq_n_f32(a), vb));
#else
		result = quatf(a + b.x, a + b.y, a + b.z, a + b.w);
#endif
		return result;
	}

	TFM_INLINE quatf operator-(const quatf& a) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(result.col, _mm_xor_ps(va, _mm_set1_ps(-0.0f)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(result.col, vnegq_f32(va));
#else
		result = quatf(-a.x, -a.y, -a.z, -a.w);
#endif
		return result;
	}

	TFM_INLINE quatf operator-(const quatf& a, const quatf& b) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_sub_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vsubq_f32(va, vb));
#else
		result = quatf(a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w);
#endif
		return result;
	}

	TFM_INLINE quatf operator-(const quatf& a, float b) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(result.col, _mm_sub_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		vst1q_f32(result.col, vsubq_f32(va, vdupq_n_f32(b)));
#else
		result = quatf(a.x - b, a.y - b, a.z - b, a.w - b);
#endif
		return result;
	}

	TFM_INLINE quatf operator-(float a, const quatf& b) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_sub_ps(_mm_set1_ps(a), vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vsubq_f32(vdupq_n_f32(a), vb));
#else
		result = quatf(a - b.x, a - b.y, a - b.z, a - b.w);
#endif
		return result;
	}

	TFM_INLINE quatf operator*(const quatf& a, const quatf& b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		const __m128 ax = _mm_shuffle_ps(va, va, _MM_SHUFFLE(0, 0, 0, 0));
		const __m128 ay = _mm_shuffle_ps(va, va, _MM_SHUFFLE(1, 1, 1, 1));
		const __m128 az = _mm_shuffle_ps(va, va, _MM_SHUFFLE(2, 2, 2, 2));
		const __m128 aw = _mm_shuffle_ps(va, va, _MM_SHUFFLE(3, 3, 3, 3));
		const __m128 b_wzyx = _mm_xor_ps(_mm_shuffle_ps(vb, vb, _MM_SHUFFLE(0, 1, 2, 3)), _mm_setr_ps(0.0f, -0.0f, 0.0f, -0.0f));
		const __m128 b_zwxy = _mm_xor_ps(_mm_shuffle_ps(vb, vb, _MM_SHUFFLE(1, 0, 3, 2)), _mm_setr_ps(0.0f, 0.0f, -0.0f, -0.0f));
		const __m128 b_yxwz = _mm_xor_ps(_mm_shuffle_ps(vb, vb, _MM_SHUFFLE(2, 3, 0, 1)), _mm_setr_ps(-0.0f, 0.0f, 0.0f, -0.0f));
		__m128 r = _mm_mul_ps(aw, vb);
		r = _mm_fmadd_ps(ax, b_wzyx, r);
		r = _mm_fmadd_ps(ay, b_zwxy, r);
		r = _mm_fmadd_ps(az, b_yxwz, r);
		quatf result;
		TFM_MM_STORE_PS(result.col, r);
		return result;
#else
		return quatf(
			a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
			a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
			a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w,
			a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z
		);
#endif
	}

	TFM_INLINE quatf operator*(const quatf& a, float b) noexcept {
		return quatf(a.x * b, a.y * b, a.z * b, a.w * b);
	}

	TFM_INLINE quatf operator*(float a, const quatf& b) noexcept {
		return quatf(a * b.x, a * b.y, a * b.z, a * b.w);
	}

	TFM_INLINE quatf operator/(const quatf& a, const quatf& b) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_div_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		float32x4_t vb = vld1q_f32(b.col);
		float32x4_t recip = vrecpeq_f32(vb);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		vst1q_f32(result.col, vmulq_f32(va, recip));
#else
		result = quatf(a.x / b.x, a.y / b.y, a.z / b.z, a.w / b.w);
#endif
		return result;
	}

	TFM_INLINE quatf operator/(const quatf& a, float b) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		TFM_MM_STORE_PS(result.col, _mm_div_ps(va, _mm_set1_ps(b)));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		float32x4_t vb = vdupq_n_f32(b);
		float32x4_t recip = vrecpeq_f32(vb);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		vst1q_f32(result.col, vmulq_f32(va, recip));
#else
		result = quatf(a.x / b, a.y / b, a.z / b, a.w / b);
#endif
		return result;
	}

	TFM_INLINE quatf operator/(float a, const quatf& b) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_div_ps(_mm_set1_ps(a), vb));
#elif TFM_NEON_INTRINSICS
		float32x4_t va = vdupq_n_f32(a);
		float32x4_t vb = vld1q_f32(b.col);
		float32x4_t recip = vrecpeq_f32(vb);
		recip = vmulq_f32(vrecpsq_f32(vb, recip), recip);
		vst1q_f32(result.col, vmulq_f32(va, recip));
#else
		result = quatf(a / b.x, a / b.y, a / b.z, a / b.w);
#endif
		return result;
	}

	TFM_INLINE quatf max(const quatf& a, const quatf& b) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_max_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vmaxq_f32(va, vb));
#else
		result = quatf(std::max(a.x, b.x), std::max(a.y, b.y), std::max(a.z, b.z), std::max(a.w, b.w));
#endif
		return result;
	}

	TFM_INLINE quatf min(const quatf& a, const quatf& b) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		TFM_MM_STORE_PS(result.col, _mm_min_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		vst1q_f32(result.col, vminq_f32(va, vb));
#else
		result = quatf(std::min(a.x, b.x), std::min(a.y, b.y), std::min(a.z, b.z), std::min(a.w, b.w));
#endif
		return result;
	}

	TFM_INLINE quatf clamp(const quatf& v, const quatf& minv, const quatf& maxv) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 vv = TFM_MM_LOAD_PS(v.col);
		const __m128 vmin = TFM_MM_LOAD_PS(minv.col);
		const __m128 vmax = TFM_MM_LOAD_PS(maxv.col);
		TFM_MM_STORE_PS(result.col, _mm_min_ps(_mm_max_ps(vv, vmin), vmax));
#elif TFM_NEON_INTRINSICS
		const float32x4_t vv = vld1q_f32(v.col);
		const float32x4_t vmin = vld1q_f32(minv.col);
		const float32x4_t vmax = vld1q_f32(maxv.col);
		vst1q_f32(result.col, vminq_f32(vmaxq_f32(vv, vmin), vmax));
#else
		result = min(max(v, minv), maxv);
#endif
		return result;
	}

	TFM_INLINE float dot(const quatf& a, const quatf& b) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 va = TFM_MM_LOAD_PS(a.col);
		const __m128 vb = TFM_MM_LOAD_PS(b.col);
		return _mm_cvtss_f32(_mm_dot_ps(va, vb));
#elif TFM_NEON_INTRINSICS
		const float32x4_t va = vld1q_f32(a.col);
		const float32x4_t vb = vld1q_f32(b.col);
		float32x4_t m = vmulq_f32(va, vb);
		float32x2_t sum = vadd_f32(vget_low_f32(m), vget_high_f32(m));
		sum = vpadd_f32(sum, sum);
		return vget_lane_f32(sum, 0);
#else
		return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
#endif
	}

	TFM_INLINE float quatf::length_squared() const noexcept {
		return dot(*this, *this);
	}

	TFM_INLINE float quatf::length() const noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 c = TFM_MM_LOAD_PS(col);
		return _mm_length_f32(c);
#elif TFM_NEON_INTRINSICS
		float32x4_t c = vld1q_f32(col);
		float32x4_t m = vmulq_f32(c, c);
		float32x2_t sum = vadd_f32(vget_low_f32(m), vget_high_f32(m));
		sum = vpadd_f32(sum, sum);
		float len2 = vget_lane_f32(sum, 0);
		return sqrtf(len2);
#else
		return sqrtf(dot(*this, *this));
#endif
	}

	TFM_INLINE quatf quatf::normalize() const noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 c = TFM_MM_LOAD_PS(col);
		TFM_MM_STORE_PS(result.col, _mm_normalize_ps(c));
#elif TFM_NEON_INTRINSICS
		float32x4_t c = vld1q_f32(col);
		float32x4_t m = vmulq_f32(c, c);
		float32x2_t sum = vadd_f32(vget_low_f32(m), vget_high_f32(m));
		sum = vpadd_f32(sum, sum);
		float32x4_t len2 = vdupq_n_f32(vget_lane_f32(sum, 0));
		float32x4_t invLen = vrsqrteq_f32(len2);
		invLen = vmulq_f32(invLen, vrsqrtsq_f32(vmulq_f32(len2, invLen), invLen));
		invLen = vmulq_f32(invLen, vrsqrtsq_f32(vmulq_f32(len2, invLen), invLen));
		vst1q_f32(result.col, vmulq_f32(c, invLen));
#else
		float len = length();
		result = (len > 0.0f) ? *this / len : quatf::zero();
#endif
		return result;
	}

	TFM_INLINE quatf quatf::round() const noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 c = TFM_MM_LOAD_PS(col);
		TFM_MM_STORE_PS(result.col, _mm_round_ps(c, _MM_FROUND_TO_NEAREST_INT | _MM_FROUND_NO_EXC));
#elif TFM_NEON_INTRINSICS
		const float32x4_t c = vld1q_f32(col);
		vst1q_f32(result.col, vrndnq_f32(c));
#else
		result = quatf(std::round(x), std::round(y), std::round(z), std::round(w));
#endif
		return result;
	}

	TFM_INLINE quatf quatf::ceil() const noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 c = TFM_MM_LOAD_PS(col);
		TFM_MM_STORE_PS(result.col, _mm_ceil_ps(c));
#elif TFM_NEON_INTRINSICS
		const float32x4_t c = vld1q_f32(col);
		vst1q_f32(result.col, vrndpq_f32(c));
#else
		result = quatf(std::ceil(x), std::ceil(y), std::ceil(z), std::ceil(w));
#endif
		return result;
	}

	TFM_INLINE quatf lerp(const quatf& a, const quatf& b, float t) noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		__m128 va = TFM_MM_LOAD_PS(a.col);
		__m128 vb = TFM_MM_LOAD_PS(b.col);
		__m128 ta = _mm_set1_ps(t);
		__m128 diff = _mm_sub_ps(vb, va);
		TFM_MM_STORE_PS(result.col, _mm_fmadd_ps(diff, ta, va));
#else
		result = a + (b - a) * t;
#endif
		return result;
	}

	TFM_INLINE quatf mix(const quatf& a, const quatf& b, float t) noexcept {
		return lerp(a, b, t);
	}

	TFM_INLINE quatf slerp(quatf q1, quatf q2, float t) {
		float dot0 = dot(q1, q2);

		quatf q2m = q2;
		if (dot0 < 0.0f) {
			dot0 = -dot0;
			q2m = -q2m;
		}

		if (dot0 > 1.0f - FLT_EPSILON) return lerp(q1, q2m, t);

		float theta = std::acos(dot0);
		float sin_theta = std::sin(theta);
		float s0 = std::sin((1.0f - t) * theta) / sin_theta;
		float s1 = std::sin(t * theta) / sin_theta;

		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 vq1 = TFM_MM_LOAD_PS(q1.col);
		const __m128 vq2m = TFM_MM_LOAD_PS(q2m.col);
		TFM_MM_STORE_PS(result.col, _mm_fmadd_ps(_mm_set1_ps(s1), vq2m, _mm_mul_ps(_mm_set1_ps(s0), vq1)));
#elif TFM_NEON_INTRINSICS
		float32x4_t vs0 = vdupq_n_f32(s0);
		float32x4_t vs1 = vdupq_n_f32(s1);
		const float32x4_t vq1 = vld1q_f32(q1.col);
		const float32x4_t vq2m = vld1q_f32(q2m.col);
		vst1q_f32(result.col, vfmaq_f32(vmulq_f32(vs0, vq1), vs1, vq2m));
#else
		result.x = s0 * q1.x + s1 * q2m.x;
		result.y = s0 * q1.y + s1 * q2m.y;
		result.z = s0 * q1.z + s1 * q2m.z;
		result.w = s0 * q1.w + s1 * q2m.w;
#endif
		return result;
	}

	TFM_INLINE quatf quatf::conjugate() const noexcept {
		quatf result;
#if TFM_SSE_AVX_INTRINSICS
		const __m128 c = TFM_MM_LOAD_PS(col);
		TFM_MM_STORE_PS(result.col, _mm_xor_ps(c, _mm_setr_ps(-0.0f, -0.0f, -0.0f, 0.0f)));
#elif TFM_NEON_INTRINSICS
		const uint32x4_t signMask = { 0x80000000u, 0x80000000u, 0x80000000u, 0x00000000u };
		const float32x4_t c = vld1q_f32(col);
		const float32x4_t conjugated = vreinterpretq_f32_u32(veorq_u32(vreinterpretq_u32_f32(c), signMask));
		vst1q_f32(result.col, conjugated);
#else
		result = quatf(-x, -y, -z, w);
#endif
		return result;
	}

	TFM_INLINE quatf quatf::inverse() const noexcept {
		float lenSq = length_squared();
		return lenSq > 0.0f ? conjugate() / lenSq : quatf::identity();
	}

	TFM_INLINE float& quatf::operator[](std::size_t i) noexcept {
		return col[i];
	}

	TFM_INLINE const float& quatf::operator[](std::size_t i) const noexcept {
		return col[i];
	}

	TFM_INLINE const quatf quatf::angle_axis(float angle, float x, float y, float z) noexcept {
		float half = angle * 0.5f;
		float s = sinf(half);
		float c = cosf(half);

		return quatf(x * s, y * s, z * s, c);
	}

	TFM_INLINE const quatf quatf::angle_axis(float angle, const vec3f& v) noexcept {
		return angle_axis(angle, v.x, v.y, v.z);
	}

	TFM_INLINE const quatf quatf::angle_axis(float angle, const vec2f& v, float z) noexcept {
		return angle_axis(angle, v.x, v.y, z);
	}

	TFM_INLINE const quatf quatf::angle_axis(float angle, const float* ptr) noexcept {
		return angle_axis(angle, ptr[0], ptr[1], ptr[2]);
	}

	TFM_INLINE const quatf quatf::angle_axis(float angle, const float* ptr, float z) noexcept {
		return angle_axis(angle, ptr[0], ptr[1], z);
	}

	TFM_INLINE const quatf quatf::angle_axis(float angle, const double* ptr) noexcept {
		return angle_axis(angle, (float)ptr[0], (float)ptr[1], (float)ptr[2]);
	}

	TFM_INLINE const quatf quatf::angle_axis(float angle, const double* ptr, float z) noexcept {
		return angle_axis(angle, (float)ptr[0], (float)ptr[1], z);
	}

	TFM_INLINE const quatf quatf::from_euler(float pitch, float yaw, float roll) noexcept {
		// Closed form of angle_axis(pitch, X) * angle_axis(yaw, Y) * angle_axis(roll, Z).
		const float sp = sinf(pitch * 0.5f), cp = cosf(pitch * 0.5f);
		const float sy = sinf(yaw * 0.5f), cy = cosf(yaw * 0.5f);
		const float sr = sinf(roll * 0.5f), cr = cosf(roll * 0.5f);

		return quatf(
			sp * cy * cr + cp * sy * sr,
			cp * sy * cr - sp * cy * sr,
			cp * cy * sr + sp * sy * cr,
			cp * cy * cr - sp * sy * sr
		);
	}

	TFM_INLINE const quatf quatf::from_euler(const vec3f& angles) noexcept {
		return from_euler(angles.x, angles.y, angles.z);
	}

	TFM_INLINE std::ostream& operator<<(std::ostream& os, const quatf& vec) noexcept {
		for (int j = 0; j < 4; j++) {
			os << std::setw(8) << std::setprecision(4) << vec.col[j] << " ";
		}
		return os;
	}

	/*********/
	/* mat3f */
	/*********/

	TFM_INLINE mat3f::mat3f(
		float m00, float m01, float m02,
		float m10, float m11, float m12,
		float m20, float m21, float m22) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		TFM_MM_STORE_PS(row[0].col, _mm_setr_ps(m00, m01, m02, 0.0f));
		TFM_MM_STORE_PS(row[1].col, _mm_setr_ps(m10, m11, m12, 0.0f));
		TFM_MM_STORE_PS(row[2].col, _mm_setr_ps(m20, m21, m22, 0.0f));
#else
		this->m00 = m00; this->m01 = m01; this->m02 = m02; _pad0 = 0;
		this->m10 = m10; this->m11 = m11; this->m12 = m12; _pad1 = 0;
		this->m20 = m20; this->m21 = m21; this->m22 = m22; _pad2 = 0;
#endif
	}

	TFM_INLINE mat3f::mat3f(const float* ptr) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 mask3 = _mm_castsi128_ps(_mm_setr_epi32(-1, -1, -1, 0));
		TFM_MM_STORE_PS(row[0].col, _mm_and_ps(_mm_loadu_ps(ptr), mask3));
		TFM_MM_STORE_PS(row[1].col, _mm_and_ps(_mm_loadu_ps(ptr + 3), mask3));
		TFM_MM_STORE_PS(row[2].col, _mm_setr_ps(ptr[6], ptr[7], ptr[8], 0.0f));
#else
		m00 = ptr[0]; m01 = ptr[1]; m02 = ptr[2]; _pad0 = 0;
		m10 = ptr[3]; m11 = ptr[4]; m12 = ptr[5]; _pad1 = 0;
		m20 = ptr[6]; m21 = ptr[7]; m22 = ptr[8]; _pad2 = 0;
#endif
	}

	TFM_INLINE mat3f::mat3f(const double* ptr) noexcept :
		m00(ptr[0]), m01(ptr[1]), m02(ptr[2]), _pad0(0),
		m10(ptr[3]), m11(ptr[4]), m12(ptr[5]), _pad1(0),
		m20(ptr[6]), m21(ptr[7]), m22(ptr[8]), _pad2(0) {}

	TFM_INLINE mat3f::mat3f(const vec3f& r0, const vec3f& r1, const vec3f& r2) noexcept {
		row[0] = vec4f(r0, 0);
		row[1] = vec4f(r1, 0);
		row[2] = vec4f(r2, 0);

		_pad0 = _pad1 = _pad2 = 0;
	}

	TFM_INLINE mat3f::mat3f(const vec4f& r0, const vec4f& r1, const vec4f& r2) noexcept {
		row[0] = r0;
		row[1] = r1;
		row[2] = r2;
	}

	TFM_INLINE mat3f::mat3f(const quatf& q) noexcept {
		const float xx = q.x * q.x;
		const float yy = q.y * q.y;
		const float zz = q.z * q.z;
		const float xy = q.x * q.y;
		const float xz = q.x * q.z;
		const float yz = q.y * q.z;
		const float wx = q.w * q.x;
		const float wy = q.w * q.y;
		const float wz = q.w * q.z;

		m00 = 1.0f - 2.0f * (yy + zz);
		m01 = 2.0f * (xy + wz);
		m02 = 2.0f * (xz - wy);

		m10 = 2.0f * (xy - wz);
		m11 = 1.0f - 2.0f * (xx + zz);
		m12 = 2.0f * (yz + wx);

		m20 = 2.0f * (xz + wy);
		m21 = 2.0f * (yz - wx);
		m22 = 1.0f - 2.0f * (xx + yy);

		_pad0 = _pad1 = _pad2 = 0;
	}

	TFM_INLINE mat3f::mat3f(const mat4f& m) noexcept {
		row[0] = m.row[0];
		row[1] = m.row[1];
		row[2] = m.row[2];

		_pad0 = _pad1 = _pad2 = 0;
	}

#if TFM_MATRIX_RTL
	TFM_INLINE mat3f operator*(const mat3f& rhs, const mat3f& lhs) noexcept {
#else
	TFM_INLINE mat3f operator*(const mat3f& lhs, const mat3f& rhs) noexcept {
#endif
		mat3f result;
#if TFM_SSE_AVX_INTRINSICS
		for (int i = 0; i < 3; i++) {
			const vec4f& row = rhs.row[i];

			const __m128 lhs0 = TFM_MM_LOAD_PS(lhs.row[0].col);
			const __m128 lhs1 = TFM_MM_LOAD_PS(lhs.row[1].col);
			const __m128 lhs2 = TFM_MM_LOAD_PS(lhs.row[2].col);
			__m128 r = _mm_mul_ps(lhs0, _mm_set1_ps(row.x));
			r = _mm_fmadd_ps(lhs1, _mm_set1_ps(row.y), r);
			TFM_MM_STORE_PS(result.row[i].col, _mm_fmadd_ps(lhs2, _mm_set1_ps(row.z), r));
		}
#elif TFM_NEON_INTRINSICS
		for (int i = 0; i < 3; i++) {
			const vec4f& row = rhs.row[i];

			const float32x4_t lhs0 = vld1q_f32(lhs.row[0].col);
			const float32x4_t lhs1 = vld1q_f32(lhs.row[1].col);
			const float32x4_t lhs2 = vld1q_f32(lhs.row[2].col);
			float32x4_t r = vmulq_f32(lhs0, vdupq_n_f32(row.x));
			r = vfmaq_f32(r, lhs1, vdupq_n_f32(row.y));
			vst1q_f32(result.row[i].col, vfmaq_f32(r, lhs2, vdupq_n_f32(row.z)));
		}
#else
		for (int i = 0; i < 3; i++) {
			const vec4f& bc = rhs.row[i];

			result.row[i] = lhs.row[0] * bc.x + lhs.row[1] * bc.y + lhs.row[2] * bc.z;
		}
#endif
		return result;
	}

	TFM_INLINE bool operator==(const mat3f& a, const mat3f& b) noexcept {
		for (int i = 0; i < 3; i++) {
			if (a.row[i] != b.row[i]) {
				return false;
			}
		}
		return true;
	}

	TFM_INLINE bool operator!=(const mat3f& a, const mat3f& b) noexcept {
		return !(a == b);
	}

	TFM_INLINE std::ostream& operator<<(std::ostream & os, const mat3f& mat) noexcept {
		for (int i = 0; i < 3; i++) {
			os << mat.row[i] << std::endl;
		}
		return os;
	}

	TFM_INLINE const mat3f mat3f::zero() noexcept {
		return mat3f{};
	}

	TFM_INLINE const mat3f mat3f::identity() noexcept {
		return mat3f(
			1, 0, 0,
			0, 1, 0,
			0, 0, 1
		);
	}

	TFM_INLINE mat3f mat3f::rotation(const quatf& q) noexcept {
		return mat3f(q);
	}

	TFM_INLINE vec3f mat3f::transform(const vec3f& v) const noexcept {
#if TFM_SSE_AVX_INTRINSICS
		__m128 r = _mm_mul_ps(TFM_MM_LOAD_PS(row[0].col), _mm_set1_ps(v.x));
		r = _mm_fmadd_ps(TFM_MM_LOAD_PS(row[1].col), _mm_set1_ps(v.y), r);
		r = _mm_fmadd_ps(TFM_MM_LOAD_PS(row[2].col), _mm_set1_ps(v.z), r);
		return vec3f(
			_mm_cvtss_f32(r),
			_mm_cvtss_f32(_mm_shuffle_ps(r, r, _MM_SHUFFLE(1, 1, 1, 1))),
			_mm_cvtss_f32(_mm_movehl_ps(r, r))
		);
#else
		return vec3f(
			v.x * m00 + v.y * m10 + v.z * m20,
			v.x * m01 + v.y * m11 + v.z * m21,
			v.x * m02 + v.y * m12 + v.z * m22
		);
#endif
	}

	TFM_INLINE mat3f mat3f::inverse() const noexcept {
#if TFM_SSE_AVX_INTRINSICS
		// Cofactor rows: c0 = cross(r1, r2), c1 = cross(r2, r0), c2 = cross(r0, r1).
		// inverse = transpose(c0, c1, c2) / dot(r0, c0). Cross w lanes are always 0.
		const __m128 r0 = TFM_MM_LOAD_PS(row[0].col);
		const __m128 r1 = TFM_MM_LOAD_PS(row[1].col);
		const __m128 r2 = TFM_MM_LOAD_PS(row[2].col);
		__m128 c0 = _mm_sub_ps(_mm_mul_ps(r1, _mm_shuffle_ps(r2, r2, _MM_SHUFFLE(3, 0, 2, 1))), _mm_mul_ps(_mm_shuffle_ps(r1, r1, _MM_SHUFFLE(3, 0, 2, 1)), r2));
		c0 = _mm_shuffle_ps(c0, c0, _MM_SHUFFLE(3, 0, 2, 1));
		__m128 c1 = _mm_sub_ps(_mm_mul_ps(r2, _mm_shuffle_ps(r0, r0, _MM_SHUFFLE(3, 0, 2, 1))), _mm_mul_ps(_mm_shuffle_ps(r2, r2, _MM_SHUFFLE(3, 0, 2, 1)), r0));
		c1 = _mm_shuffle_ps(c1, c1, _MM_SHUFFLE(3, 0, 2, 1));
		__m128 c2 = _mm_sub_ps(_mm_mul_ps(r0, _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(3, 0, 2, 1))), _mm_mul_ps(_mm_shuffle_ps(r0, r0, _MM_SHUFFLE(3, 0, 2, 1)), r1));
		c2 = _mm_shuffle_ps(c2, c2, _MM_SHUFFLE(3, 0, 2, 1));

		const float det = _mm_cvtss_f32(_mm_dot_ps(r0, c0));
		if (std::fabs(det) < FLT_EPSILON) return identity();
		const __m128 invDet = _mm_set1_ps(1.0f / det);

		const __m128 zero = _mm_setzero_ps();
		const __m128 t0 = _mm_unpacklo_ps(c0, c1);
		const __m128 t1 = _mm_unpackhi_ps(c0, c1);
		const __m128 t2 = _mm_unpacklo_ps(c2, zero);
		const __m128 t3 = _mm_unpackhi_ps(c2, zero);

		mat3f out;
		TFM_MM_STORE_PS(out.row[0].col, _mm_mul_ps(_mm_movelh_ps(t0, t2), invDet));
		TFM_MM_STORE_PS(out.row[1].col, _mm_mul_ps(_mm_movehl_ps(t2, t0), invDet));
		TFM_MM_STORE_PS(out.row[2].col, _mm_mul_ps(_mm_movelh_ps(t1, t3), invDet));
		return out;
#else
		float det = m00 * (m11 * m22 - m12 * m21) - m01 * (m10 * m22 - m12 * m20) + m02 * (m10 * m21 - m11 * m20);

		if (std::fabs(det) < FLT_EPSILON) return identity();

		float invDet = 1.0f / det;

		mat3f out;
		out.m00 = (m11 * m22 - m12 * m21) * invDet;
		out.m01 = -(m01 * m22 - m02 * m21) * invDet;
		out.m02 = (m01 * m12 - m02 * m11) * invDet;
		out._pad0 = 0;

		out.m10 = -(m10 * m22 - m12 * m20) * invDet;
		out.m11 = (m00 * m22 - m02 * m20) * invDet;
		out.m12 = -(m00 * m12 - m02 * m10) * invDet;
		out._pad1 = 0;

		out.m20 = (m10 * m21 - m11 * m20) * invDet;
		out.m21 = -(m00 * m21 - m01 * m20) * invDet;
		out.m22 = (m00 * m11 - m01 * m10) * invDet;
		out._pad2 = 0;

		return out;
#endif
	}

	TFM_INLINE mat3f mat3f::transpose() const noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 r0 = TFM_MM_LOAD_PS(row[0].col);
		const __m128 r1 = TFM_MM_LOAD_PS(row[1].col);
		const __m128 r2 = TFM_MM_LOAD_PS(row[2].col);
		const __m128 zero = _mm_setzero_ps();
		const __m128 t0 = _mm_unpacklo_ps(r0, r1);
		const __m128 t1 = _mm_unpackhi_ps(r0, r1);
		const __m128 t2 = _mm_unpacklo_ps(r2, zero);
		const __m128 t3 = _mm_unpackhi_ps(r2, zero);
		mat3f out;
		TFM_MM_STORE_PS(out.row[0].col, _mm_movelh_ps(t0, t2));
		TFM_MM_STORE_PS(out.row[1].col, _mm_movehl_ps(t2, t0));
		TFM_MM_STORE_PS(out.row[2].col, _mm_movelh_ps(t1, t3));
		return out;
#else
		return mat3f(
			m00, m10, m20,
			m01, m11, m21,
			m02, m12, m22
		);
#endif
	}

	TFM_INLINE mat3f mat3f::inverse_transpose() const noexcept {
#if TFM_SSE_AVX_INTRINSICS
		// transpose(inverse) = cofactor rows / det directly, no transpose needed.
		const __m128 r0 = TFM_MM_LOAD_PS(row[0].col);
		const __m128 r1 = TFM_MM_LOAD_PS(row[1].col);
		const __m128 r2 = TFM_MM_LOAD_PS(row[2].col);
		__m128 c0 = _mm_sub_ps(_mm_mul_ps(r1, _mm_shuffle_ps(r2, r2, _MM_SHUFFLE(3, 0, 2, 1))), _mm_mul_ps(_mm_shuffle_ps(r1, r1, _MM_SHUFFLE(3, 0, 2, 1)), r2));
		c0 = _mm_shuffle_ps(c0, c0, _MM_SHUFFLE(3, 0, 2, 1));
		__m128 c1 = _mm_sub_ps(_mm_mul_ps(r2, _mm_shuffle_ps(r0, r0, _MM_SHUFFLE(3, 0, 2, 1))), _mm_mul_ps(_mm_shuffle_ps(r2, r2, _MM_SHUFFLE(3, 0, 2, 1)), r0));
		c1 = _mm_shuffle_ps(c1, c1, _MM_SHUFFLE(3, 0, 2, 1));
		__m128 c2 = _mm_sub_ps(_mm_mul_ps(r0, _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(3, 0, 2, 1))), _mm_mul_ps(_mm_shuffle_ps(r0, r0, _MM_SHUFFLE(3, 0, 2, 1)), r1));
		c2 = _mm_shuffle_ps(c2, c2, _MM_SHUFFLE(3, 0, 2, 1));

		const float det = _mm_cvtss_f32(_mm_dot_ps(r0, c0));
		if (std::fabs(det) < FLT_EPSILON) return identity();
		const __m128 invDet = _mm_set1_ps(1.0f / det);

		mat3f out;
		TFM_MM_STORE_PS(out.row[0].col, _mm_mul_ps(c0, invDet));
		TFM_MM_STORE_PS(out.row[1].col, _mm_mul_ps(c1, invDet));
		TFM_MM_STORE_PS(out.row[2].col, _mm_mul_ps(c2, invDet));
		return out;
#else
		return inverse().transpose();
#endif
	}

	TFM_INLINE vec4f& mat3f::operator[](std::size_t i) noexcept {
		return row[i];
	}

	TFM_INLINE const vec4f& mat3f::operator[](std::size_t i) const noexcept {
		return row[i];
	}

	/********/
	/* mat4f */
	/********/

	TFM_INLINE mat4f::mat4f(const float* ptr) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		TFM_MM_STORE_PS(row[0].col, _mm_loadu_ps(ptr));
		TFM_MM_STORE_PS(row[1].col, _mm_loadu_ps(ptr + 4));
		TFM_MM_STORE_PS(row[2].col, _mm_loadu_ps(ptr + 8));
		TFM_MM_STORE_PS(row[3].col, _mm_loadu_ps(ptr + 12));
#elif TFM_NEON_INTRINSICS
		vst1q_f32(row[0].col, vld1q_f32(ptr));
		vst1q_f32(row[1].col, vld1q_f32(ptr + 4));
		vst1q_f32(row[2].col, vld1q_f32(ptr + 8));
		vst1q_f32(row[3].col, vld1q_f32(ptr + 12));
#else
		for (int i = 0; i < 16; i++) m[i] = ptr[i];
#endif
	}

	TFM_INLINE mat4f::mat4f(const double* ptr) noexcept {
#if TFM_SSE_AVX_INTRINSICS
		for (int i = 0; i < 4; i++) {
			const __m128 lo = _mm_cvtpd_ps(_mm_loadu_pd(ptr + 4 * i));
			const __m128 hi = _mm_cvtpd_ps(_mm_loadu_pd(ptr + 4 * i + 2));
			TFM_MM_STORE_PS(row[i].col, _mm_movelh_ps(lo, hi));
		}
#else
		for (int i = 0; i < 16; i++) m[i] = (float)ptr[i];
#endif
	}

	TFM_INLINE mat4f::mat4f(const quatf& q) noexcept {
		const float xx = q.x * q.x;
		const float yy = q.y * q.y;
		const float zz = q.z * q.z;
		const float xy = q.x * q.y;
		const float xz = q.x * q.z;
		const float yz = q.y * q.z;
		const float wx = q.w * q.x;
		const float wy = q.w * q.y;
		const float wz = q.w * q.z;

		m00 = 1.0f - 2.0f * (yy + zz);
		m01 = 2.0f * (xy + wz);
		m02 = 2.0f * (xz - wy);
		m03 = 0.0f;

		m10 = 2.0f * (xy - wz);
		m11 = 1.0f - 2.0f * (xx + zz);
		m12 = 2.0f * (yz + wx);
		m13 = 0.0f;

		m20 = 2.0f * (xz + wy);
		m21 = 2.0f * (yz - wx);
		m22 = 1.0f - 2.0f * (xx + yy);
		m23 = 0.0f;

		m30 = 0.0f;
		m31 = 0.0f;
		m32 = 0.0f;
		m33 = 1.0f;
	}

	TFM_INLINE mat4f::mat4f(const vec4f& r0, const vec4f& r1, const vec4f& r2, const vec4f& r3) noexcept {
		row[0] = r0;
		row[1] = r1;
		row[2] = r2;
		row[3] = r3;
	}

	TFM_INLINE mat4f::mat4f(
		float m00, float m01, float m02, float m03,
		float m10, float m11, float m12, float m13,
		float m20, float m21, float m22, float m23,
		float m30, float m31, float m32, float m33) {
#if TFM_SSE_AVX_INTRINSICS
		TFM_MM_STORE_PS(row[0].col, _mm_setr_ps(m00, m01, m02, m03));
		TFM_MM_STORE_PS(row[1].col, _mm_setr_ps(m10, m11, m12, m13));
		TFM_MM_STORE_PS(row[2].col, _mm_setr_ps(m20, m21, m22, m23));
		TFM_MM_STORE_PS(row[3].col, _mm_setr_ps(m30, m31, m32, m33));
#else
		this->m00 = m00; this->m01 = m01; this->m02 = m02; this->m03 = m03;
		this->m10 = m10; this->m11 = m11; this->m12 = m12; this->m13 = m13;
		this->m20 = m20; this->m21 = m21; this->m22 = m22; this->m23 = m23;
		this->m30 = m30; this->m31 = m31; this->m32 = m32; this->m33 = m33;
#endif
	}

	TFM_INLINE mat4f::mat4f(const mat3f& m) noexcept {
		row[0] = m.row[0];
		row[1] = m.row[1];
		row[2] = m.row[2];
		row[3] = vec4f::zero();
	}

#if TFM_MATRIX_RTL
	TFM_INLINE mat4f operator*(const mat4f& rhs, const mat4f& lhs) noexcept {
#else
	TFM_INLINE mat4f operator*(const mat4f& lhs, const mat4f& rhs) noexcept {
#endif
		mat4f result;
#if TFM_SSE_AVX_INTRINSICS
		for (int i = 0; i < 4; i++) {
			const vec4f& row = rhs.row[i];

			const __m128 lhs0 = TFM_MM_LOAD_PS(lhs.row[0].col);
			const __m128 lhs1 = TFM_MM_LOAD_PS(lhs.row[1].col);
			const __m128 lhs2 = TFM_MM_LOAD_PS(lhs.row[2].col);
			const __m128 lhs3 = TFM_MM_LOAD_PS(lhs.row[3].col);
			__m128 r = _mm_mul_ps(lhs0, _mm_set1_ps(row.x));
			r = _mm_fmadd_ps(lhs1, _mm_set1_ps(row.y), r);
			r = _mm_fmadd_ps(lhs2, _mm_set1_ps(row.z), r);
			TFM_MM_STORE_PS(result.row[i].col, _mm_fmadd_ps(lhs3, _mm_set1_ps(row.w), r));
		}
#elif TFM_NEON_INTRINSICS
		for (int i = 0; i < 4; i++) {
			const vec4f& row = rhs.row[i];

			const float32x4_t lhs0 = vld1q_f32(lhs.row[0].col);
			const float32x4_t lhs1 = vld1q_f32(lhs.row[1].col);
			const float32x4_t lhs2 = vld1q_f32(lhs.row[2].col);
			const float32x4_t lhs3 = vld1q_f32(lhs.row[3].col);
			float32x4_t r = vmulq_f32(lhs0, vdupq_n_f32(row.x));
			r = vfmaq_f32(r, lhs1, vdupq_n_f32(row.y));
			r = vfmaq_f32(r, lhs2, vdupq_n_f32(row.z));
			vst1q_f32(result.row[i].col, vfmaq_f32(r, lhs3, vdupq_n_f32(row.w)));
		}
#else
		for (int i = 0; i < 4; i++) {
			result.row[i] = lhs.row[0] * rhs.row[i].x + lhs.row[1] * rhs.row[i].y + lhs.row[2] * rhs.row[i].z + lhs.row[3] * rhs.row[i].w;
		}
#endif
		return result;
	}

	TFM_INLINE mat4f& operator*=(mat4f & a, const mat4f& b) noexcept {
		return a = a * b;
	}

	TFM_INLINE bool operator==(const mat4f& a, const mat4f& b) noexcept {
		for (int i = 0; i < 4; i++) {
			if (a.row[i] != b.row[i]) {
				return false;
			}
		}
		return true;
	}

	TFM_INLINE bool operator!=(const mat4f& a, const mat4f& b) noexcept {
		return !(a == b);
	}

	TFM_INLINE std::ostream& operator<<(std::ostream & os, const mat4f& mat) noexcept {
		for (int i = 0; i < 4; i++) {
			os << mat.row[i] << std::endl;
		}
		return os;
	}

	TFM_INLINE const mat4f mat4f::zero() noexcept {
		return mat4f{};
	}

	TFM_INLINE const mat4f mat4f::identity() noexcept {
		return mat4f(
			1, 0, 0, 0,
			0, 1, 0, 0,
			0, 0, 1, 0,
			0, 0, 0, 1
		);
	}

	TFM_INLINE vec3f mat4f::translation() const noexcept {
		return vec3f(m30, m31, m32);
	}

	TFM_INLINE vec3f mat4f::scaling() const noexcept {
		return vec3f(
			vec3f(m00, m01, m02).length(),
			vec3f(m10, m11, m12).length(),
			vec3f(m20, m21, m22).length()
		);
	}

	TFM_INLINE mat3f mat4f::rotation_matrix() const noexcept {
		return mat3f(*this);
	}

	TFM_INLINE vec3f mat4f::transform_point(const vec3f& v) const noexcept {
		return vec3f::transform(v, *this);
	}

	TFM_INLINE vec3f mat4f::transform_vector(const vec3f& v) const noexcept {
		return vec3f::transform_normal(v, *this);
	}

	template <typename T> TFM_INLINE mat4f mat4f::transformation(const std::vector<T>& matrix, const std::vector<T>& translation, const std::vector<T>& rotation, const std::vector<T>& scaling) {
		if (!matrix.empty()) return mat4f(matrix.data());

		vec3f t = translation.empty() ? vec3f::zero() : vec3f(translation.data());
		quatf r = rotation.empty() ? quatf(0.0f, 0.0f, 0.0f, 1.0f) : quatf(rotation.data());
		vec3f s = scaling.empty() ? vec3f::one() : vec3f(scaling.data());

		// Direct S * R * T composition: rotation rows scaled per axis, translation in the last row.
		const mat3f rm(r);
		mat4f result;
		result.row[0] = rm.row[0] * s.x;
		result.row[1] = rm.row[1] * s.y;
		result.row[2] = rm.row[2] * s.z;
		result.row[3] = vec4f(t, 1.0f);
		return result;
	}

	TFM_INLINE mat4f mat4f::translation(float x, float y, float z) noexcept {
		return mat4f(
			1.0f, 0.0f, 0.0f, 0.0f,
			0.0f, 1.0f, 0.0f, 0.0f,
			0.0f, 0.0f, 1.0f, 0.0f,
			x, y, z, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::translation(const vec3f& v) noexcept {
		return translation(v.x, v.y, v.z);
	}

	TFM_INLINE mat4f mat4f::translate(float tx, float ty, float tz) const noexcept {
		mat4f result;
		result.row[0] = row[0];
		result.row[1] = row[1];
		result.row[2] = row[2];
		result.row[3] = row[0] * tx + row[1] * ty + row[2] * tz + row[3];
		return result;
	}

	TFM_INLINE mat4f mat4f::translate(const vec3f& v) const noexcept {
		return translate(v.x, v.y, v.z);
	}

	TFM_INLINE mat4f mat4f::scaling(float x, float y, float z) noexcept {
		return mat4f(
			x, 0.0f, 0.0f, 0.0f,
			0.0f, y, 0.0f, 0.0f,
			0.0f, 0.0f, z, 0.0f,
			0.0f, 0.0f, 0.0f, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::scaling(const vec3f& v) noexcept {
		return scaling(v.x, v.y, v.z);
	}

	TFM_INLINE mat4f mat4f::scaling(float xyz) noexcept {
		return scaling(xyz, xyz, xyz);
	}

	TFM_INLINE mat4f mat4f::scale(float sx, float sy, float sz) const noexcept {
		mat4f result;
		result.row[0] = row[0] * sx;
		result.row[1] = row[1] * sy;
		result.row[2] = row[2] * sz;
		result.row[3] = row[3];
		return result;
	}

	TFM_INLINE mat4f mat4f::scale(const vec3f& v) const noexcept {
		return scale(v.x, v.y, v.z);
	}

	TFM_INLINE mat4f mat4f::scale(float xyz) const noexcept {
		return scale(xyz, xyz, xyz);
	}

	TFM_INLINE mat4f mat4f::rotation(float angle, const vec3f& vec) noexcept {
		float c = cosf(angle);
		float s = sinf(angle);

		vec3f axis = vec.normalize();
		vec3f temp = (1.0f - c) * axis;

		mat4f result;
		result.row[0] = vec4f(c + temp.x * axis.x, temp.x * axis.y + s * axis.z, temp.x * axis.z - s * axis.y, 0.0f);
		result.row[1] = vec4f(temp.y * axis.x - s * axis.z, c + temp.y * axis.y, temp.y * axis.z + s * axis.x, 0.0f);
		result.row[2] = vec4f(temp.z * axis.x + s * axis.y, temp.z * axis.y - s * axis.x, c + temp.z * axis.z, 0.0f);
		result.row[3] = vec4f(0.0f, 0.0f, 0.0f, 1.0f);
		return result;
	}

	TFM_INLINE mat4f mat4f::rotation(float angle, float rx, float ry, float rz) noexcept {
		return rotation(angle, vec3f(rx, ry, rz));
	}

	TFM_INLINE mat4f mat4f::rotation_x(float angle) noexcept {
		float c = cosf(angle);
		float s = sinf(angle);

		return mat4f(
			1.0f, 0.0f, 0.0f, 0.0f,
			0.0f, c, s, 0.0f,
			0.0f, -s, c, 0.0f,
			0.0f, 0.0f, 0.0f, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::rotation_y(float angle) noexcept {
		float c = cosf(angle);
		float s = sinf(angle);

		return mat4f(
			c, 0.0f, -s, 0.0f,
			0.0f, 1.0f, 0.0f, 0.0f,
			s, 0.0f, c, 0.0f,
			0.0f, 0.0f, 0.0f, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::rotation_z(float angle) noexcept {
		float c = cosf(angle);
		float s = sinf(angle);

		return mat4f(
			c, s, 0.0f, 0.0f,
			-s, c, 0.0f, 0.0f,
			0.0f, 0.0f, 1.0f, 0.0f,
			0.0f, 0.0f, 0.0f, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::rotation(const quatf& q) noexcept {
		return mat4f(q);
	}

	TFM_INLINE mat4f mat4f::rotate(float angle, const vec3f& vec) const noexcept {
		float c = cosf(angle);
		float s = sinf(angle);

		vec3f axis = vec.normalize();
		vec3f temp = (1.0f - c) * axis;

		float m00 = c + temp.x * axis.x;
		float m01 = temp.x * axis.y + s * axis.z;
		float m02 = temp.x * axis.z - s * axis.y;

		float m10 = temp.y * axis.x - s * axis.z;
		float m11 = c + temp.y * axis.y;
		float m12 = temp.y * axis.z + s * axis.x;
		float m13 = 0.0f;

		float m20 = temp.z * axis.x + s * axis.y;
		float m21 = temp.z * axis.y - s * axis.x;
		float m22 = c + temp.z * axis.z;

		mat4f result;
		result.row[0] = row[0] * m00 + row[1] * m01 + row[2] * m02;
		result.row[1] = row[0] * m10 + row[1] * m11 + row[2] * m12;
		result.row[2] = row[0] * m20 + row[1] * m21 + row[2] * m22;
		result.row[3] = row[3];
		return result;
	}

	TFM_INLINE mat4f mat4f::rotate(float angle, float rx, float ry, float rz) const noexcept {
		return rotate(angle, vec3f(rx, ry, rz));
	}

	TFM_INLINE mat4f mat4f::rotate_x(float angle) const noexcept {
		float c = cosf(angle);
		float s = sinf(angle);

		mat4f result;
		result.row[0] = row[0];
		result.row[1] = row[1] * c + row[2] * s;
		result.row[2] = row[2] * c - row[1] * s;
		result.row[3] = row[3];
		return result;
	}

	TFM_INLINE mat4f mat4f::rotate_y(float angle) const noexcept {
		float c = cosf(angle);
		float s = sinf(angle);

		mat4f result;
		result.row[0] = row[0] * c - row[2] * s;
		result.row[1] = row[1];
		result.row[2] = row[0] * s + row[2] * c;
		result.row[3] = row[3];
		return result;
	}

	TFM_INLINE mat4f mat4f::rotate_z(float angle) const noexcept {
		float c = cosf(angle);
		float s = sinf(angle);

		mat4f result;
		result.row[0] = row[0] * c + row[1] * s;
		result.row[1] = row[1] * c - row[0] * s;
		result.row[2] = row[2];
		result.row[3] = row[3];
		return result;
	}

	TFM_INLINE mat4f mat4f::rotate(const quatf& q) const noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const mat3f r(q);
		const __m128 q0 = TFM_MM_LOAD_PS(r.row[0].col);
		const __m128 q1 = TFM_MM_LOAD_PS(r.row[1].col);
		const __m128 q2 = TFM_MM_LOAD_PS(r.row[2].col);
		mat4f result;
		for (int i = 0; i < 4; i++) {
			const __m128 v = TFM_MM_LOAD_PS(row[i].col);
			__m128 t = _mm_mul_ps(q0, _mm_shuffle_ps(v, v, _MM_SHUFFLE(0, 0, 0, 0)));
			t = _mm_fmadd_ps(q1, _mm_shuffle_ps(v, v, _MM_SHUFFLE(1, 1, 1, 1)), t);
			t = _mm_fmadd_ps(q2, _mm_shuffle_ps(v, v, _MM_SHUFFLE(2, 2, 2, 2)), t);
			t = _mm_blend_ps(t, v, 0x8);
			TFM_MM_STORE_PS(result.row[i].col, t);
		}
		return result;
#else
		return (*this) * mat4f(q);
#endif
	}

	TFM_INLINE mat4f mat4f::perspective_fov_rh_zo(float fovy, float aspect, float znear, float zfar) noexcept {
		const float invTan = 1.0f / tanf(fovy * 0.5f);
		const float invZDiff = 1.0f / (znear - zfar);
		const float zTerm = zfar * invZDiff;
		const float wzTerm = znear * zTerm;

		return mat4f(
			invTan / aspect, 0.0f, 0.0f, 0.0f,
			0.0f, invTan, 0.0f, 0.0f,
			0.0f, 0.0f, zTerm, -1.0f,
			0.0f, 0.0f, wzTerm, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::perspective_rh_zo(float left, float right, float bottom, float top, float znear, float zfar) noexcept {
		float invRL = 1.0f / (right - left);
		float invTB = 1.0f / (top - bottom);
		float invFN = 1.0f / (znear - zfar);

		float A = (right + left) * invRL;
		float B = (top + bottom) * invTB;

		return mat4f(
			2.0f * znear * invRL, 0.0f, 0.0f, 0.0f,
			0.0f, 2.0f * znear * invTB, 0.0f, 0.0f,
			A, B, zfar * invFN, -1.0f,
			0.0f, 0.0f, -znear * invFN, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::perspective_rh_zo(float width, float height, float znear, float zfar) noexcept {
		const float invHalfW = znear / (width * 0.5f);
		const float invHalfH = znear / (height * 0.5f);

		const float invZDiff = 1.0f / (znear - zfar);
		const float zTerm = zfar * invZDiff;
		const float wzTerm = znear * zfar * invZDiff;

		return mat4f(
			invHalfW, 0.0f, 0.0f, 0.0f,
			0.0f, invHalfH, 0.0f, 0.0f,
			0.0f, 0.0f, zTerm, -1.0f,
			0.0f, 0.0f, wzTerm, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::perspective_fov_rh_no(float fovy, float aspect, float znear, float zfar) noexcept {
		const float invTan = 1.0f / tanf(fovy * 0.5f);
		const float invZDiff = 1.0f / (znear - zfar);
		const float zTerm = (zfar + znear) * invZDiff;
		const float wzTerm = 2.0f * znear * zfar * invZDiff;

		return mat4f(
			invTan / aspect, 0.0f, 0.0f, 0.0f,
			0.0f, invTan, 0.0f, 0.0f,
			0.0f, 0.0f, zTerm, -1.0f,
			0.0f, 0.0f, wzTerm, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::perspective_rh_no(float left, float right, float bottom, float top, float znear, float zfar) noexcept {
		float invRL = 1.0f / (right - left);
		float invTB = 1.0f / (top - bottom);
		float invFN = 1.0f / (znear - zfar);

		float A = (right + left) * invRL;
		float B = (top + bottom) * invTB;
		float C = (zfar + znear) * invFN;
		float D = 2.0f * znear * zfar * invFN;

		return mat4f(
			2.0f * znear * invRL, 0.0f, 0.0f, 0.0f,
			0.0f, 2.0f * znear * invTB, 0.0f, 0.0f,
			A, B, C, -1.0f,
			0.0f, 0.0f, D, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::perspective_rh_no(float width, float height, float znear, float zfar) noexcept {
		float invHalfW = znear / (width * 0.5f);
		float invHalfH = znear / (height * 0.5f);

		float invZDiff = 1.0f / (znear - zfar);
		float zTerm = (zfar + znear) * invZDiff;
		float wzTerm = 2.0f * znear * zfar * invZDiff;

		return mat4f(
			invHalfW, 0.0f, 0.0f, 0.0f,
			0.0f, invHalfH, 0.0f, 0.0f,
			0.0f, 0.0f, zTerm, -1.0f,
			0.0f, 0.0f, wzTerm, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::perspective_fov_lh_zo(float fovy, float aspect, float znear, float zfar) noexcept {
		const float invTan = 1.0f / tanf(fovy * 0.5f);
		const float invZDiff = 1.0f / (zfar - znear);
		const float zTerm = zfar * invZDiff;
		const float wzTerm = -znear * zTerm;

		return mat4f(
			invTan / aspect, 0.0f, 0.0f, 0.0f,
			0.0f, invTan, 0.0f, 0.0f,
			0.0f, 0.0f, zTerm, 1,
			0.0f, 0.0f, wzTerm, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::perspective_lh_zo(float left, float right, float bottom, float top, float znear, float zfar) noexcept {
		float invRL = 1.0f / (right - left);
		float invTB = 1.0f / (top - bottom);
		float invFN = 1.0f / (zfar - znear);

		float A = (right + left) * invRL;
		float B = (top + bottom) * invTB;

		return mat4f(
			2.0f * znear * invRL, 0.0f, 0.0f, 0.0f,
			0.0f, 2.0f * znear * invTB, 0.0f, 0.0f,
			A, B, zfar * invFN, 1.0f,
			0.0f, 0.0f, -znear * invFN, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::perspective_lh_zo(float width, float height, float znear, float zfar) noexcept {
		float invHalfW = znear / (width * 0.5f);
		float invHalfH = znear / (height * 0.5f);

		float invZDiff = 1.0f / (zfar - znear);

		float zTerm = zfar * invZDiff;
		float wzTerm = -znear * zfar * invZDiff;

		return mat4f(
			invHalfW, 0.0f, 0.0f, 0.0f,
			0.0f, invHalfH, 0.0f, 0.0f,
			0.0f, 0.0f, zTerm, 1.0f,
			0.0f, 0.0f, wzTerm, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::perspective_fov_lh_no(float fovy, float aspect, float znear, float zfar) noexcept {
		const float invTan = 1.0f / tanf(fovy * 0.5f);
		const float invZDiff = 1.0f / (zfar - znear);
		const float zTerm = (zfar + znear) * invZDiff;
		const float wzTerm = -2.0f * znear * zfar * invZDiff;

		return mat4f(
			invTan / aspect, 0.0f, 0.0f, 0.0f,
			0.0f, invTan, 0.0f, 0.0f,
			0.0f, 0.0f, zTerm, 1.0f,
			0.0f, 0.0f, wzTerm, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::perspective_lh_no(float left, float right, float bottom, float top, float znear, float zfar) noexcept {
		float invRL = 1.0f / (right - left);
		float invTB = 1.0f / (top - bottom);
		float invFN = 1.0f / (zfar - znear);

		float A = (right + left) * invRL;
		float B = (top + bottom) * invTB;
		float C = (zfar + znear) * invFN;
		float D = -2.0f * znear * zfar * invFN;

		return mat4f(
			2.0f * znear * invRL, 0.0f, 0.0f, 0.0f,
			0.0f, 2.0f * znear * invTB, 0.0f, 0.0f,
			A, B, C, 1.0f,
			0.0f, 0.0f, D, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::perspective_lh_no(float width, float height, float znear, float zfar) noexcept {
		float invHalfW = znear / (width * 0.5f);
		float invHalfH = znear / (height * 0.5f);

		float invZDiff = 1.0f / (zfar - znear);
		float zTerm = (zfar + znear) * invZDiff;
		float wzTerm = -2.0f * znear * zfar * invZDiff;

		return mat4f(
			invHalfW, 0.0f, 0.0f, 0.0f,
			0.0f, invHalfH, 0.0f, 0.0f,
			0.0f, 0.0f, zTerm, 1.0f,
			0.0f, 0.0f, wzTerm, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::frustum_zo_rh(float left, float right, float bottom, float top, float znear, float zfar) noexcept {
		float rl = right - left;
		float tb = top - bottom;
		float fn = zfar - znear;

		return mat4f(
			(2.0f * znear) / rl, 0.0f, (right + left) / rl, 0.0f,
			0.0f, (2.0f * znear) / tb, (top + bottom) / tb, 0.0f,
			0.0f, 0.0f, zfar / fn, -(zfar * znear) / fn,
			0.0f, 0.0f, -1.0f, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::frustum_zo_lh(float left, float right, float bottom, float top, float znear, float zfar) noexcept {
		float rl = right - left;
		float tb = top - bottom;
		float fn = zfar - znear;

		return mat4f(
			(2.0f * znear) / rl, 0.0f, (right + left) / rl, 0.0f,
			0.0f, (2.0f * znear) / tb, (top + bottom) / tb, 0.0f,
			0.0f, 0.0f, zfar / fn, -(zfar * znear) / fn,
			0.0f, 0.0f, 1.0f, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::frustum_no_lh(float left, float right, float bottom, float top, float znear, float zfar) noexcept {
		float rl = right - left;
		float tb = top - bottom;
		float fn = zfar - znear;

		return mat4f(
			(2.0f * znear) / rl, 0.0f, (right + left) / rl, 0.0f,
			0.0f, (2.0f * znear) / tb, (top + bottom) / tb, 0.0f,
			0.0f, 0.0f, -(zfar + znear) / fn, -(2.0f * zfar * znear) / fn,
			0.0f, 0.0f, 1.0f, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::frustum_no_rh(float left, float right, float bottom, float top, float znear, float zfar) noexcept {
		float rl = right - left;
		float tb = top - bottom;
		float fn = zfar - znear;

		return mat4f(
			(2.0f * znear) / rl, 0.0f, (right + left) / rl, 0.0f,
			0.0f, (2.0f * znear) / tb, (top + bottom) / tb, 0.0f,
			0.0f, 0.0f, -(zfar + znear) / fn, -(2.0f * zfar * znear) / fn,
			0.0f, 0.0f, -1.0f, 0.0f
		);
	}

	TFM_INLINE mat4f mat4f::ortho_rh_zo(float left, float right, float bottom, float top, float znear, float zfar) noexcept {
		float invRL = 1.0f / (right - left);
		float invTB = 1.0f / (top - bottom);
		float invFN = 1.0f / (zfar - znear);

		return mat4f(
			2.0f * invRL, 0.0f, 0.0f, 0.0f,
			0.0f, 2.0f * invTB, 0.0f, 0.0f,
			0.0f, 0.0f, -invFN, 0.0f,
			-(right + left) * invRL, -(top + bottom) * invTB, -znear * invFN, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::ortho_rh_zo(float width, float height, float znear, float zfar) noexcept {
		float halfW = width * 0.5f;
		float halfH = height * 0.5f;
		float invFN = 1.0f / (zfar - znear);

		return mat4f(
			1.0f / halfW, 0.0f, 0.0f, 0.0f,
			0.0f, 1.0f / halfH, 0.0f, 0.0f,
			0.0f, 0.0f, -invFN, 0.0f,
			0.0f, 0.0f, -znear * invFN, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::ortho_lh_zo(float left, float right, float bottom, float top, float znear, float zfar) noexcept {
		float invRL = 1.0f / (right - left);
		float invTB = 1.0f / (top - bottom);
		float invFN = 1.0f / (zfar - znear);

		return mat4f(
			2.0f * invRL, 0.0f, 0.0f, 0.0f,
			0.0f, 2.0f * invTB, 0.0f, 0.0f,
			0.0f, 0.0f, invFN, 0.0f,
			-(right + left) * invRL, -(top + bottom) * invTB, -znear * invFN, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::ortho_lh_zo(float width, float height, float znear, float zfar) noexcept {
		float halfW = width * 0.5f;
		float halfH = height * 0.5f;
		float invFN = 1.0f / (zfar - znear);

		return mat4f(
			1.0f / halfW, 0.0f, 0.0f, 0.0f,
			0.0f, 1.0f / halfH, 0.0f, 0.0f,
			0.0f, 0.0f, invFN, 0.0f,
			0.0f, 0.0f, -znear * invFN, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::ortho_rh_no(float left, float right, float bottom, float top, float znear, float zfar) noexcept {
		const float invRL = 1.0f / (right - left);
		const float invTB = 1.0f / (top - bottom);
		const float invFN = 1.0f / (zfar - znear);

		return mat4f(
			2.0f * invRL, 0.0f, 0.0f, 0.0f,
			0.0f, 2.0f * invTB, 0.0f, 0.0f,
			0.0f, 0.0f, -2.0f * invFN, 0.0f,
			-(right + left) * invRL, -(top + bottom) * invTB, -(zfar + znear) * invFN, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::ortho_rh_no(float width, float height, float znear, float zfar) noexcept {
		float halfW = width * 0.5f;
		float halfH = height * 0.5f;
		float invFN = 1.0f / (zfar - znear);

		return mat4f(
			1.0f / halfW, 0.0f, 0.0f, 0.0f,
			0.0f, 1.0f / halfH, 0.0f, 0.0f,
			0.0f, 0.0f, -2.0f * invFN, 0.0f,
			0.0f, 0.0f, -(zfar + znear) * invFN, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::ortho_lh_no(float left, float right, float bottom, float top, float znear, float zfar) noexcept {
		const float invRL = 1.0f / (right - left);
		const float invTB = 1.0f / (top - bottom);
		const float invFN = 1.0f / (zfar - znear);

		return mat4f(
			2.0f * invRL, 0.0f, 0.0f, 0.0f,
			0.0f, 2.0f * invTB, 0.0f, 0.0f,
			0.0f, 0.0f, 2.0f * invFN, 0.0f,
			-(right + left) * invRL, -(top + bottom) * invTB, -(zfar + znear) * invFN, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::ortho_lh_no(float width, float height, float znear, float zfar) noexcept {
		float halfW = width * 0.5f;
		float halfH = height * 0.5f;
		float invFN = 1.0f / (zfar - znear);

		return mat4f(
			1.0f / halfW, 0.0f, 0.0f, 0.0f,
			0.0f, 1.0f / halfH, 0.0f, 0.0f,
			0.0f, 0.0f, 2.0f * invFN, 0.0f,
			0.0f, 0.0f, -(zfar + znear) * invFN, 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::look_to_rh(const vec3f& eye, const vec3f& dir, const vec3f& up) noexcept {
		vec3f f = dir.normalize();
		vec3f r = cross(f, up).normalize();
		vec3f u = cross(r, f);

		return mat4f(
			r.x, u.x, -f.x, 0.0f,
			r.y, u.y, -f.y, 0.0f,
			r.z, u.z, -f.z, 0.0f,
			-dot(r, eye), -dot(u, eye), dot(f, eye), 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::look_at_rh(const vec3f& eye, const vec3f& target, const vec3f& up) noexcept {
		return look_to_rh(eye, target - eye, up);
	}

	TFM_INLINE mat4f mat4f::look_to_lh(const vec3f& eye, const vec3f& dir, const vec3f& up) noexcept {
		vec3f f = dir.normalize();
		vec3f r = cross(up, f).normalize();
		vec3f u = cross(f, r);

		return mat4f(
			r.x, u.x, f.x, 0.0f,
			r.y, u.y, f.y, 0.0f,
			r.z, u.z, f.z, 0.0f,
			-dot(r, eye), -dot(u, eye), -dot(f, eye), 1.0f
		);
	}

	TFM_INLINE mat4f mat4f::look_at_lh(const vec3f& eye, const vec3f& target, const vec3f& up) noexcept {
		return look_to_lh(eye, target - eye, up);
	}

	TFM_INLINE mat4f mat4f::inverse() const noexcept {
#if TFM_SSE_AVX_INTRINSICS
		__m128 r0 = TFM_MM_LOAD_PS(row[0].col);
		__m128 r1 = TFM_MM_LOAD_PS(row[1].col);
		__m128 r2 = TFM_MM_LOAD_PS(row[2].col);
		__m128 r3 = TFM_MM_LOAD_PS(row[3].col);

		__m128 r2_2222 = _mm_shuffle_ps(r2, r1, _MM_SHUFFLE(2, 2, 2, 2));
		__m128 r2_1111 = _mm_shuffle_ps(r2, r1, _MM_SHUFFLE(1, 1, 1, 1));
		__m128 r2_3333 = _mm_shuffle_ps(r2, r1, _MM_SHUFFLE(3, 3, 3, 3));
		__m128 r2_0000 = _mm_shuffle_ps(r2, r1, _MM_SHUFFLE(0, 0, 0, 0));
		__m128 r3_3333 = _mm_shuffle_ps(r3, r2, _MM_SHUFFLE(3, 3, 3, 3));
		__m128 r3_2222 = _mm_shuffle_ps(r3, r2, _MM_SHUFFLE(2, 2, 2, 2));
		__m128 r3_1111 = _mm_shuffle_ps(r3, r2, _MM_SHUFFLE(1, 1, 1, 1));
		__m128 r3_0000 = _mm_shuffle_ps(r3, r2, _MM_SHUFFLE(0, 0, 0, 0));

		__m128 shuf_r3_3333_20 = _mm_shuffle_ps(r3_3333, r3_3333, _MM_SHUFFLE(2, 0, 0, 0));
		__m128 shuf_r3_2222_20 = _mm_shuffle_ps(r3_2222, r3_2222, _MM_SHUFFLE(2, 0, 0, 0));
		__m128 shuf_r3_1111_20 = _mm_shuffle_ps(r3_1111, r3_1111, _MM_SHUFFLE(2, 0, 0, 0));
		__m128 shuf_r3_0000_20 = _mm_shuffle_ps(r3_0000, r3_0000, _MM_SHUFFLE(2, 0, 0, 0));

		// Facs
		__m128 fac0 = _mm_sub_ps(_mm_mul_ps(r2_2222, shuf_r3_3333_20), _mm_mul_ps(r2_3333, shuf_r3_2222_20));
		__m128 fac1 = _mm_sub_ps(_mm_mul_ps(r2_1111, shuf_r3_3333_20), _mm_mul_ps(r2_3333, shuf_r3_1111_20));
		__m128 fac2 = _mm_sub_ps(_mm_mul_ps(r2_1111, shuf_r3_2222_20), _mm_mul_ps(r2_2222, shuf_r3_1111_20));
		__m128 fac3 = _mm_sub_ps(_mm_mul_ps(r2_0000, shuf_r3_3333_20), _mm_mul_ps(r2_3333, shuf_r3_0000_20));
		__m128 fac4 = _mm_sub_ps(_mm_mul_ps(r2_0000, shuf_r3_2222_20), _mm_mul_ps(r2_2222, shuf_r3_0000_20));
		__m128 fac5 = _mm_sub_ps(_mm_mul_ps(r2_0000, shuf_r3_1111_20), _mm_mul_ps(r2_1111, shuf_r3_0000_20));

		// Vecs
		__m128 vec0 = _mm_shuffle_ps(_mm_shuffle_ps(r1, r0, _MM_SHUFFLE(0, 0, 0, 0)), _mm_shuffle_ps(r1, r0, _MM_SHUFFLE(0, 0, 0, 0)), _MM_SHUFFLE(2, 2, 2, 0));
		__m128 vec1 = _mm_shuffle_ps(_mm_shuffle_ps(r1, r0, _MM_SHUFFLE(1, 1, 1, 1)), _mm_shuffle_ps(r1, r0, _MM_SHUFFLE(1, 1, 1, 1)), _MM_SHUFFLE(2, 2, 2, 0));
		__m128 vec2f = _mm_shuffle_ps(_mm_shuffle_ps(r1, r0, _MM_SHUFFLE(2, 2, 2, 2)), _mm_shuffle_ps(r1, r0, _MM_SHUFFLE(2, 2, 2, 2)), _MM_SHUFFLE(2, 2, 2, 0));
		__m128 vec3f = _mm_shuffle_ps(_mm_shuffle_ps(r1, r0, _MM_SHUFFLE(3, 3, 3, 3)), _mm_shuffle_ps(r1, r0, _MM_SHUFFLE(3, 3, 3, 3)), _MM_SHUFFLE(2, 2, 2, 0));

		// Signs
		__m128 sign_a = _mm_set_ps(1.0f, -1.0f, 1.0f, -1.0f);
		__m128 sign_b = _mm_set_ps(-1.0f, 1.0f, -1.0f, 1.0f);

		// Inverse columns
		__m128 inv0 = _mm_mul_ps(sign_b, _mm_add_ps(_mm_sub_ps(_mm_mul_ps(vec1, fac0), _mm_mul_ps(vec2f, fac1)), _mm_mul_ps(vec3f, fac2)));
		__m128 inv1 = _mm_mul_ps(sign_a, _mm_add_ps(_mm_sub_ps(_mm_mul_ps(vec0, fac0), _mm_mul_ps(vec2f, fac3)), _mm_mul_ps(vec3f, fac4)));
		__m128 inv2 = _mm_mul_ps(sign_b, _mm_add_ps(_mm_sub_ps(_mm_mul_ps(vec0, fac1), _mm_mul_ps(vec1, fac3)), _mm_mul_ps(vec3f, fac5)));
		__m128 inv3 = _mm_mul_ps(sign_a, _mm_add_ps(_mm_sub_ps(_mm_mul_ps(vec0, fac2), _mm_mul_ps(vec1, fac4)), _mm_mul_ps(vec2f, fac5)));

		// Determinant
		__m128 row0 = _mm_shuffle_ps(inv0, inv1, _MM_SHUFFLE(0, 0, 0, 0));
		__m128 row1 = _mm_shuffle_ps(inv2, inv3, _MM_SHUFFLE(0, 0, 0, 0));
		__m128 row2 = _mm_shuffle_ps(row0, row1, _MM_SHUFFLE(2, 0, 2, 0));
		__m128 det0 = _mm_dot_ps(r0, row2);
		__m128 rcp0 = _mm_div_ps(_mm_set1_ps(1.0f), det0);

		mat4f out;
		TFM_MM_STORE_PS(out.row[0].col, _mm_mul_ps(inv0, rcp0));
		TFM_MM_STORE_PS(out.row[1].col, _mm_mul_ps(inv1, rcp0));
		TFM_MM_STORE_PS(out.row[2].col, _mm_mul_ps(inv2, rcp0));
		TFM_MM_STORE_PS(out.row[3].col, _mm_mul_ps(inv3, rcp0));
		return out;
#elif TFM_NEON_INTRINSICS

		float32x4_t r0 = vld1q_f32(row[0].col);
		float32x4_t r1 = vld1q_f32(row[1].col);
		float32x4_t r2 = vld1q_f32(row[2].col);
		float32x4_t r3 = vld1q_f32(row[3].col);

		// Broadcasts (equivalente aos _mm_shuffle_ps usados)
		float32x4_t r2_2222 = vdupq_laneq_f32(r2, 2);
		float32x4_t r2_1111 = vdupq_laneq_f32(r2, 1);
		float32x4_t r2_3333 = vdupq_laneq_f32(r2, 3);
		float32x4_t r2_0000 = vdupq_laneq_f32(r2, 0);

		float32x4_t r3_3333 = vdupq_laneq_f32(r3, 3);
		float32x4_t r3_2222 = vdupq_laneq_f32(r3, 2);
		float32x4_t r3_1111 = vdupq_laneq_f32(r3, 1);
		float32x4_t r3_0000 = vdupq_laneq_f32(r3, 0);

		// Facs
		float32x4_t fac0 = vsubq_f32(vmulq_f32(r2_2222, r3_3333), vmulq_f32(r2_3333, r3_2222));
		float32x4_t fac1 = vsubq_f32(vmulq_f32(r2_1111, r3_3333), vmulq_f32(r2_3333, r3_1111));
		float32x4_t fac2 = vsubq_f32(vmulq_f32(r2_1111, r3_2222), vmulq_f32(r2_2222, r3_1111));
		float32x4_t fac3 = vsubq_f32(vmulq_f32(r2_0000, r3_3333), vmulq_f32(r2_3333, r3_0000));
		float32x4_t fac4 = vsubq_f32(vmulq_f32(r2_0000, r3_2222), vmulq_f32(r2_2222, r3_0000));
		float32x4_t fac5 = vsubq_f32(vmulq_f32(r2_0000, r3_1111), vmulq_f32(r2_1111, r3_0000));

		float32x4_t vec0 = vdupq_laneq_f32(r1, 0);
		float32x4_t vec1 = vdupq_laneq_f32(r1, 1);
		float32x4_t vec2f = vdupq_laneq_f32(r1, 2);
		float32x4_t vec3f = vdupq_laneq_f32(r1, 3);

		float32x4_t signA = { -1.f,  1.f, -1.f,  1.f };
		float32x4_t signB = { 1.f, -1.f,  1.f, -1.f };

		float32x4_t inv0 = vmulq_f32(signB, vaddq_f32(vsubq_f32(vmulq_f32(vec1, fac0), vmulq_f32(vec2f, fac1)), vmulq_f32(vec3f, fac2)));
		float32x4_t inv1 = vmulq_f32(signA, vaddq_f32(vsubq_f32(vmulq_f32(vec0, fac0), vmulq_f32(vec2f, fac3)), vmulq_f32(vec3f, fac4)));
		float32x4_t inv2 = vmulq_f32(signB, vaddq_f32(vsubq_f32(vmulq_f32(vec0, fac1), vmulq_f32(vec1, fac3)), vmulq_f32(vec3f, fac5)));
		float32x4_t inv3 = vmulq_f32(signA, vaddq_f32(vsubq_f32(vmulq_f32(vec0, fac2), vmulq_f32(vec1, fac4)), vmulq_f32(vec2f, fac5)));

		float32x4_t mul = vmulq_f32(r0, inv0);
		float32x2_t sum2 = vpadd_f32(vget_low_f32(mul), vget_high_f32(mul));
		float32x2_t sum1 = vpadd_f32(sum2, sum2);
		float32x4_t det = vdupq_n_f32(vget_lane_f32(sum1, 0));

		float32x4_t rcp = vdivq_f32(vdupq_n_f32(1.0f), det);

		mat4f out;
		vst1q_f32(out.row[0].col, vmulq_f32(inv0, rcp));
		vst1q_f32(out.row[1].col, vmulq_f32(inv1, rcp));
		vst1q_f32(out.row[2].col, vmulq_f32(inv2, rcp));
		vst1q_f32(out.row[3].col, vmulq_f32(inv3, rcp));
		return out;
#else
		mat4f inv;

		inv.m00 = m11 * m22 * m33 - m11 * m23 * m32 - m21 * m12 * m33 + m21 * m13 * m32 + m31 * m12 * m23 - m31 * m13 * m22;
		inv.m01 = -m01 * m22 * m33 + m01 * m23 * m32 + m21 * m02 * m33 - m21 * m03 * m32 - m31 * m02 * m23 + m31 * m03 * m22;
		inv.m02 = m01 * m12 * m33 - m01 * m13 * m32 - m11 * m02 * m33 + m11 * m03 * m32 + m31 * m02 * m13 - m31 * m03 * m12;
		inv.m03 = -m01 * m12 * m23 + m01 * m13 * m22 + m11 * m02 * m23 - m11 * m03 * m22 - m21 * m02 * m13 + m21 * m03 * m12;

		inv.m10 = -m10 * m22 * m33 + m10 * m23 * m32 + m20 * m12 * m33 - m20 * m13 * m32 - m30 * m12 * m23 + m30 * m13 * m22;
		inv.m11 = m00 * m22 * m33 - m00 * m23 * m32 - m20 * m02 * m33 + m20 * m03 * m32 + m30 * m02 * m23 - m30 * m03 * m22;
		inv.m12 = -m00 * m12 * m33 + m00 * m13 * m32 + m10 * m02 * m33 - m10 * m03 * m32 - m30 * m02 * m13 + m30 * m03 * m12;
		inv.m13 = m00 * m12 * m23 - m00 * m13 * m22 - m10 * m02 * m23 + m10 * m03 * m22 + m20 * m02 * m13 - m20 * m03 * m12;

		inv.m20 = m10 * m21 * m33 - m10 * m23 * m31 - m20 * m11 * m33 + m20 * m13 * m31 + m30 * m11 * m23 - m30 * m13 * m21;
		inv.m21 = -m00 * m21 * m33 + m00 * m23 * m31 + m20 * m01 * m33 - m20 * m03 * m31 - m30 * m01 * m23 + m30 * m03 * m21;
		inv.m22 = m00 * m11 * m33 - m00 * m13 * m31 - m10 * m01 * m33 + m10 * m03 * m31 + m30 * m01 * m13 - m30 * m03 * m11;
		inv.m23 = -m00 * m11 * m23 + m00 * m13 * m21 + m10 * m01 * m23 - m10 * m03 * m21 - m20 * m01 * m13 + m20 * m03 * m11;

		inv.m30 = -m10 * m21 * m32 + m10 * m22 * m31 + m20 * m11 * m32 - m20 * m12 * m31 - m30 * m11 * m22 + m30 * m12 * m21;
		inv.m31 = m00 * m21 * m32 - m00 * m22 * m31 - m20 * m01 * m32 + m20 * m02 * m31 + m30 * m01 * m22 - m30 * m02 * m21;
		inv.m32 = -m00 * m11 * m32 + m00 * m12 * m31 + m10 * m01 * m32 - m10 * m02 * m31 - m30 * m01 * m12 + m30 * m02 * m11;
		inv.m33 = m00 * m11 * m22 - m00 * m12 * m21 - m10 * m01 * m22 + m10 * m02 * m21 + m20 * m01 * m12 - m20 * m02 * m11;

		float det = m00 * inv.m00 + m01 * inv.m10 + m02 * inv.m20 + m03 * inv.m30;

		if (std::fabs(det) < FLT_EPSILON) return identity();

		float invDet = 1.0f / det;

		for (int i = 0; i < 4; i++)
			inv.row[i] *= invDet;

		return inv;
#endif
	}

	TFM_INLINE mat4f mat4f::transpose() const noexcept {
#if TFM_SSE_AVX_INTRINSICS
		const __m128 r0 = TFM_MM_LOAD_PS(row[0].col);
		const __m128 r1 = TFM_MM_LOAD_PS(row[1].col);
		const __m128 r2 = TFM_MM_LOAD_PS(row[2].col);
		const __m128 r3 = TFM_MM_LOAD_PS(row[3].col);
		const __m128 t0 = _mm_unpacklo_ps(r0, r1);
		const __m128 t1 = _mm_unpackhi_ps(r0, r1);
		const __m128 t2 = _mm_unpacklo_ps(r2, r3);
		const __m128 t3 = _mm_unpackhi_ps(r2, r3);
		mat4f out;
		TFM_MM_STORE_PS(out.row[0].col, _mm_movelh_ps(t0, t2));
		TFM_MM_STORE_PS(out.row[1].col, _mm_movehl_ps(t2, t0));
		TFM_MM_STORE_PS(out.row[2].col, _mm_movelh_ps(t1, t3));
		TFM_MM_STORE_PS(out.row[3].col, _mm_movehl_ps(t3, t1));
		return out;
#elif TFM_NEON_INTRINSICS
		const float32x4x2_t p01 = vtrnq_f32(vld1q_f32(row[0].col), vld1q_f32(row[1].col));
		const float32x4x2_t p23 = vtrnq_f32(vld1q_f32(row[2].col), vld1q_f32(row[3].col));
		mat4f out;
		vst1q_f32(out.row[0].col, vcombine_f32(vget_low_f32(p01.val[0]), vget_low_f32(p23.val[0])));
		vst1q_f32(out.row[1].col, vcombine_f32(vget_low_f32(p01.val[1]), vget_low_f32(p23.val[1])));
		vst1q_f32(out.row[2].col, vcombine_f32(vget_high_f32(p01.val[0]), vget_high_f32(p23.val[0])));
		vst1q_f32(out.row[3].col, vcombine_f32(vget_high_f32(p01.val[1]), vget_high_f32(p23.val[1])));
		return out;
#else
		return mat4f(
			m00, m10, m20, m30,
			m01, m11, m21, m31,
			m02, m12, m22, m32,
			m03, m13, m23, m33
		);
#endif
	}

	TFM_INLINE mat4f mat4f::inverse_transpose() const noexcept {
		return inverse().transpose();
	}

	TFM_INLINE vec4f& mat4f::operator[](std::size_t i) noexcept {
		return row[i];
	}

	TFM_INLINE const vec4f& mat4f::operator[](std::size_t i) const noexcept {
		return row[i];
	}
}

#pragma pop_macro("min")
#pragma pop_macro("max")