Update math and math/linalg; add "pure_none" calling convention

core/math/linalg/extended.odin

@@ -0,0 +1,553 @@
+package linalg
+
+import "builtin"
+import "core:math"
+
+radians :: proc(degrees: $T) -> (out: T) where IS_NUMERIC(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = degrees * RAD_PER_DEG;
+		}
+	} else {
+		out = degrees * RAD_PER_DEG;
+	}
+	return;
+}
+degrees :: proc(radians: $T) -> (out: T) where IS_NUMERIC(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = radians * DEG_PER_RAD;
+		}
+	} else {
+		out = radians * DEG_PER_RAD;
+	}
+	return;
+}
+
+min_double :: proc(a, b: $T) -> (out: T) where IS_NUMERIC(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = builtin.min(a[i], b[i]);
+		}
+	} else {
+		out = builtin.min(a, b);
+	}
+	return;
+}
+
+min_single :: proc(a: $T) -> (out: ELEM_TYPE(T)) where IS_NUMERIC(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		N :: len(T);
+
+		when N == 1 {
+			out = a[0];
+		} else when N == 2 {
+			out = builtin.min(a[0], a[1]);
+		} else {
+			out = builtin.min(a[0], a[1]);
+			for i in 2..<N {
+				out = builtin.min(out, a[i]);
+			}
+		}
+	} else {
+		out = a;
+	}
+	return;
+}
+
+min_triple :: proc(a, b, c: $T) -> T where IS_NUMERIC(ELEM_TYPE(T)) {
+	return min_double(a, min_double(b, c));
+}
+
+min :: proc{min_single, min_double, min_triple};
+
+max_double :: proc(a, b: $T) -> (out: T) where IS_NUMERIC(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = builtin.max(a[i], b[i]);
+		}
+	} else {
+		out = builtin.max(a, b);
+	}
+	return;
+}
+
+max_single :: proc(a: $T) -> (out: ELEM_TYPE(T)) where IS_NUMERIC(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		N :: len(T);
+
+		when N == 1 {
+			out = a[0];
+		} else when N == 2 {
+			out = builtin.max(a[0], a[1]);
+		} else when N == 3 {
+			out = builtin.max(a[0], a[1], a[3]);
+		}else {
+			out = builtin.max(a[0], a[1]);
+			for i in 2..<N {
+				out = builtin.max(out, a[i]);
+			}
+		}
+	} else {
+		out = a;
+	}
+	return;
+}
+
+max_triple :: proc(a, b, c: $T) -> T where IS_NUMERIC(ELEM_TYPE(T)) {
+	return max_double(a, max_double(b, c));
+}
+
+max :: proc{max_single, max_double, max_triple};
+
+abs :: proc(a: $T) -> (out: T) where IS_NUMERIC(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = builtin.abs(a[i]);
+		}
+	} else {
+		out = builtin.abs(a);
+	}
+	return;
+}
+
+sign :: proc(a: $T) -> (out: T) where IS_NUMERIC(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = inline math.sign(a[i]);
+		}
+	} else {
+		out = inline math.sign(a);
+	}
+	return;
+}
+
+clamp :: proc(x, a, b: $T) -> (out: T) where IS_NUMERIC(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = builtin.clamp(x[i], a[i], b[i]);
+		}
+	} else {
+		out = builtin.clamp(x, a, b);
+	}
+	return;
+}
+
+
+saturate :: proc(x: $T) -> T where IS_FLOAT(ELEM_TYPE(T)) {
+	return clamp(x, 0.0, 1.0);
+}
+
+lerp :: proc(a, b, t: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = a[i]*(1-t[i]) + b[i]*t[i];
+		}
+	} else {
+		out = a * (1.0 - t) + b * t;
+	}
+	return;
+}
+mix :: proc(a, b, t: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = a[i]*(1-t[i]) + b[i]*t[i];
+		}
+	} else {
+		out = a * (1.0 - t) + b * t;
+	}
+	return;
+}
+
+unlerp :: proc(a, b, x: $T) -> T where IS_FLOAT(ELEM_TYPE(T)) {
+	return (x - a) / (b - a);
+}
+
+step :: proc(e, x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = x[i] < e[i] ? 0.0 : 1.0;
+		}
+	} else {
+		out = x < e ? 0.0 : 1.0;
+	}
+	return;
+}
+
+smoothstep :: proc(e0, e1, x: $T) -> T where IS_FLOAT(ELEM_TYPE(T)) {
+	t := saturate(unlerp(e0, e1, x));
+	return t * t * (3.0 - 2.0 * t);
+}
+
+smootherstep :: proc(e0, e1, x: $T) -> T where IS_FLOAT(ELEM_TYPE(T)) {
+	t := saturate(unlerp(e0, e1, x));
+	return t * t * t * (t * (6*t - 15) + 10);
+}
+
+
+sqrt :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.sqrt(x[i]);
+		}
+	} else {
+		out = math.sqrt(x);
+	}
+	return;
+}
+
+inverse_sqrt :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = 1.0/math.sqrt(x[i]);
+		}
+	} else {
+		out = 1.0/math.sqrt(x);
+	}
+	return;
+}
+
+cos :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.cos(x[i]);
+		}
+	} else {
+		out = math.cos(x);
+	}
+	return;
+}
+
+sin :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.sin(x[i]);
+		}
+	} else {
+		out = math.sin(x);
+	}
+	return;
+}
+
+tan :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.tan(x[i]);
+		}
+	} else {
+		out = math.tan(x);
+	}
+	return;
+}
+
+acos :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.acos(x[i]);
+		}
+	} else {
+		out = math.acos(x);
+	}
+	return;
+}
+
+asin :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.asin(x[i]);
+		}
+	} else {
+		out = math.asin(x);
+	}
+	return;
+}
+
+atan :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.atan(x[i]);
+		}
+	} else {
+		out = math.atan(x);
+	}
+	return;
+}
+atan2 :: proc(y, x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.atan2(y[i], x[i]);
+		}
+	} else {
+		out = math.atan2(y, x);
+	}
+	return;
+}
+
+
+ln :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.ln(x[i]);
+		}
+	} else {
+		out = math.ln(x);
+	}
+	return;
+}
+
+log2 :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = INVLN2 * math.ln(x[i]);
+		}
+	} else {
+		out = INVLN2 * math.ln(x);
+	}
+	return;
+}
+
+log10 :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = INVLN10 * math.ln(x[i]);
+		}
+	} else {
+		out = INVLN10 * math.ln(x);
+	}
+	return;
+}
+
+log :: proc(x, b: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.ln(x[i]) / math.ln(cast(ELEM_TYPE(T))b[i]);
+		}
+	} else {
+		out = INVLN10 * math.ln(x) / math.ln(cast(ELEM_TYPE(T))b);
+	}
+	return;
+}
+
+exp :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.exp(x[i]);
+		}
+	} else {
+		out = math.exp(x);
+	}
+	return;
+}
+
+exp2 :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.exp(LN2 * x[i]);
+		}
+	} else {
+		out = math.exp(LN2 * x);
+	}
+	return;
+}
+
+exp10 :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.exp(LN10 * x[i]);
+		}
+	} else {
+		out = math.exp(LN10 * x);
+	}
+	return;
+}
+
+pow :: proc(x, e: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = math.pow(x[i], e[i]);
+		}
+	} else {
+		out = math.pow(x, e);
+	}
+	return;
+}
+
+
+ceil :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = inline math.ceil(x[i]);
+		}
+	} else {
+		out = inline math.ceil(x);
+	}
+	return;
+}
+
+floor :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = inline math.floor(x[i]);
+		}
+	} else {
+		out = inline math.floor(x);
+	}
+	return;
+}
+
+round :: proc(x: $T) -> (out: T) where IS_FLOAT(ELEM_TYPE(T)) {
+	when IS_ARRAY(T) {
+		for i in 0..<len(T) {
+			out[i] = inline math.round(x[i]);
+		}
+	} else {
+		out = inline math.round(x);
+	}
+	return;
+}
+
+fract :: proc(x: $T) -> T where IS_FLOAT(ELEM_TYPE(T)) {
+	f := inline floor(x);
+	return x - f;
+}
+
+mod :: proc(x, m: $T) -> T where IS_FLOAT(ELEM_TYPE(T)) {
+	f := inline floor(x / m);
+	return x - f * m;
+}
+
+
+face_forward :: proc(N, I, N_ref: $T) -> (out: T) where IS_ARRAY(T), IS_FLOAT(ELEM_TYPE(T)) {
+	return dot(N_ref, I) < 0 ? N : -N;
+}
+
+distance :: proc(p0, p1: $V/[$N]$E) -> V where IS_NUMERIC(E) {
+	return length(p1 - p0);
+}
+
+reflect :: proc(I, N: $T) -> (out: T) where IS_ARRAY(T), IS_FLOAT(ELEM_TYPE(T)) {
+	b := n * (2 * dot(n, i));
+	return i - b;
+}
+refract :: proc(I, N: $T) -> (out: T) where IS_ARRAY(T), IS_FLOAT(ELEM_TYPE(T)) {
+	dv := dot(n, i);
+	k := 1 - eta*eta - (1 - dv*dv);
+	a := i * eta;
+	b := n * eta*dv*math.sqrt(k);
+	return (a - b) * E(int(k >= 0));
+}
+
+
+
+
+is_nan_single :: proc(x: $T) -> bool where IS_FLOAT(T) {
+	return inline math.is_nan(x);
+}
+
+is_nan_array :: proc(x: $A/[$N]$T) -> (out: [N]bool) where IS_FLOAT(T) {
+	for i in 0..<N {
+		out[i] = inline is_nan(x[i]);
+	}
+	return;
+}
+
+is_inf_single :: proc(x: $T) -> bool where IS_FLOAT(T) {
+	return inline math.is_inf(x);
+}
+
+is_inf_array :: proc(x: $A/[$N]$T) -> (out: [N]bool) where IS_FLOAT(T) {
+	for i in 0..<N {
+		out[i] = inline is_inf(x[i]);
+	}
+	return;
+}
+
+classify_single :: proc(x: $T) -> math.Float_Class where IS_FLOAT(T) {
+	return inline math.classify(x);
+}
+
+classify_array :: proc(x: $A/[$N]$T) -> (out: [N]math.Float_Class) where IS_FLOAT(T) {
+	for i in 0..<N {
+		out[i] = inline classify_single(x[i]);
+	}
+	return;
+}
+
+is_nan :: proc{is_nan_single, is_nan_array};
+is_inf :: proc{is_inf_single, is_inf_array};
+classify :: proc{classify_single, classify_array};
+
+
+less_than_single          :: proc(x, y: $T) -> (out: bool) where IS_FLOAT(T) { return x < y; }
+less_than_equal_single    :: proc(x, y: $T) -> (out: bool) where IS_FLOAT(T) { return x <= y; }
+greater_than_single       :: proc(x, y: $T) -> (out: bool) where IS_FLOAT(T) { return x > y; }
+greater_than_equal_single :: proc(x, y: $T) -> (out: bool) where IS_FLOAT(T) { return x >= y; }
+equal_single              :: proc(x, y: $T) -> (out: bool) where IS_FLOAT(T) { return x == y; }
+not_equal_single          :: proc(x, y: $T) -> (out: bool) where IS_FLOAT(T) { return x != y; }
+
+less_than_array :: proc(x, y: $A/[$N]$T) -> (out: [N]bool) where IS_ARRAY(T), IS_FLOAT(ELEM_TYPE(T)) {
+	for i in 0..<N {
+		out[i] = x[i] < y[i];
+	}
+	return;
+}
+less_than_equal_array :: proc(x, y: $A/[$N]$T) -> (out: [N]bool) where IS_ARRAY(T), IS_FLOAT(ELEM_TYPE(T)) {
+	for i in 0..<N {
+		out[i] = x[i] <= y[i];
+	}
+	return;
+}
+greater_than_array :: proc(x, y: $A/[$N]$T) -> (out: [N]bool) where IS_ARRAY(T), IS_FLOAT(ELEM_TYPE(T)) {
+	for i in 0..<N {
+		out[i] = x[i] > y[i];
+	}
+	return;
+}
+greater_than_equal_array :: proc(x, y: $A/[$N]$T) -> (out: [N]bool) where IS_ARRAY(T), IS_FLOAT(ELEM_TYPE(T)) {
+	for i in 0..<N {
+		out[i] = x[i] >= y[i];
+	}
+	return;
+}
+equal_array :: proc(x, y: $A/[$N]$T) -> (out: [N]bool) where IS_ARRAY(T), IS_FLOAT(ELEM_TYPE(T)) {
+	for i in 0..<N {
+		out[i] = x[i] == y[i];
+	}
+	return;
+}
+not_equal_array :: proc(x, y: $A/[$N]$T) -> (out: [N]bool) where IS_ARRAY(T), IS_FLOAT(ELEM_TYPE(T)) {
+	for i in 0..<N {
+		out[i] = x[i] != y[i];
+	}
+	return;
+}
+
+less_than          :: proc{less_than_single, less_than_array};
+less_than_equal    :: proc{less_than_equal_single, less_than_equal_array};
+greater_than       :: proc{greater_than_single, greater_than_array};
+greater_than_equal :: proc{greater_than_equal_single, greater_than_equal_array};
+equal              :: proc{equal_single, equal_array};
+not_equal          :: proc{not_equal_single, not_equal_array};
+
+any :: proc(x: $A/[$N]bool) -> (out: bool) {
+	for e in x {
+		if x {
+			return true;
+		}
+	}
+	return false;
+}
+all :: proc(x: $A/[$N]bool) -> (out: bool) {
+	for e in x {
+		if !x {
+			return false;
+		}
+	}
+	return true;
+}
+not :: proc(x: $A/[$N]bool) -> (out: A) {
+	for e, i in x {
+		out[i] = !e;
+	}
+	return;
+}

core/math/linalg/general.odin

@@ -5,9 +5,36 @@ import "intrinsics"
 
 // Generic
 
+TAU          :: 6.28318530717958647692528676655900576;
+PI           :: 3.14159265358979323846264338327950288;
+
+E            :: 2.71828182845904523536;
+
+τ :: TAU;
+π :: PI;
+e :: E;
+
+SQRT_TWO     :: 1.41421356237309504880168872420969808;
+SQRT_THREE   :: 1.73205080756887729352744634150587236;
+SQRT_FIVE    :: 2.23606797749978969640917366873127623;
+
+LN2          :: 0.693147180559945309417232121458176568;
+LN10         :: 2.30258509299404568401799145468436421;
+
+MAX_F64_PRECISION :: 16; // Maximum number of meaningful digits after the decimal point for 'f64'
+MAX_F32_PRECISION ::  8; // Maximum number of meaningful digits after the decimal point for 'f32'
+
+RAD_PER_DEG :: TAU/360.0;
+DEG_PER_RAD :: 360.0/TAU;
+
+
+
 @private IS_NUMERIC :: intrinsics.type_is_numeric;
 @private IS_QUATERNION :: intrinsics.type_is_quaternion;
 @private IS_ARRAY :: intrinsics.type_is_array;
+@private IS_FLOAT :: intrinsics.type_is_float;
+@private BASE_TYPE :: intrinsics.type_base_type;
+@private ELEM_TYPE :: intrinsics.type_elem_type;
 
 
 vector_dot :: proc(a, b: $T/[$N]$E) -> (c: E) where IS_NUMERIC(E) {
@@ -41,8 +68,16 @@ vector_cross3 :: proc(a, b: $T/[3]$E) -> (c: T) where IS_NUMERIC(E) {
 	return;
 }
 
+quaternion_cross :: proc(q1, q2: $Q) -> (q3: Q) where IS_QUATERNION(Q) {
+	q3.x = q1.w * q2.x + q1.x * q2.w + q1.y * q2.z - q1.z * q2.y;
+	q3.y = q1.w * q2.y + q1.y * q2.w + q1.z * q2.x - q1.x * q2.z;
+	q3.z = q1.w * q2.z + q1.z * q2.w + q1.x * q2.y - q1.y * q2.x;
+	q3.w = q1.w * q2.w - q1.x * q2.x - q1.y * q2.y - q1.z * q2.z;
+	return;
+}
+
 vector_cross :: proc{vector_cross2, vector_cross3};
-cross :: vector_cross;
+cross :: proc{vector_cross2, vector_cross3, quaternion_cross};
 
 vector_normalize :: proc(v: $T/[$N]$E) -> T where IS_NUMERIC(E) {
 	return v / length(v);
@@ -83,225 +118,6 @@ length :: proc{vector_length, quaternion_length};
 length2 :: proc{vector_length2, quaternion_length2};
 
 
-vector_lerp :: proc(x, y, t: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		ti := t[i];
-		s[i] = x[i]*(1-ti) + y[i]*ti;
-	}
-	return s;
-}
-
-vector_unlerp :: proc(a, b, x: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		ai := a[i];
-		s[i] = (x[i]-ai)/(b[i]-ai);
-	}
-	return s;
-}
-
-vector_sin :: proc(angle: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.sin(angle[i]);
-	}
-	return s;
-}
-
-vector_cos :: proc(angle: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.cos(angle[i]);
-	}
-	return s;
-}
-
-vector_tan :: proc(angle: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.tan(angle[i]);
-	}
-	return s;
-}
-
-
-vector_asin :: proc(x: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.asin(x[i]);
-	}
-	return s;
-}
-
-vector_acos :: proc(x: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.acos(x[i]);
-	}
-	return s;
-}
-
-vector_atan :: proc(x: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.atan(x[i]);
-	}
-	return s;
-}
-
-vector_atan2 :: proc(y, x: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.atan(y[i], x[i]);
-	}
-	return s;
-}
-
-vector_pow :: proc(x, y: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.pow(x[i], y[i]);
-	}
-	return s;
-}
-
-vector_expr :: proc(x: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.expr(x[i]);
-	}
-	return s;
-}
-
-vector_sqrt :: proc(x: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.sqrt(x[i]);
-	}
-	return s;
-}
-
-vector_abs :: proc(x: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = abs(x[i]);
-	}
-	return s;
-}
-
-vector_sign :: proc(v: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.sign(v[i]);
-	}
-	return s;
-}
-
-vector_floor :: proc(v: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.floor(v[i]);
-	}
-	return s;
-}
-
-vector_ceil :: proc(v: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.ceil(v[i]);
-	}
-	return s;
-}
-
-
-vector_mod :: proc(x, y: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = math.mod(x[i], y[i]);
-	}
-	return s;
-}
-
-vector_min :: proc(a, b: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = min(a[i], b[i]);
-	}
-	return s;
-}
-
-vector_max :: proc(a, b: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = max(a[i], b[i]);
-	}
-	return s;
-}
-
-vector_clamp :: proc(x, a, b: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = clamp(x[i], a[i], b[i]);
-	}
-	return s;
-}
-
-vector_mix :: proc(x, y, a: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = x[i]*(1-a[i]) + y[i]*a[i];
-	}
-	return s;
-}
-
-vector_step :: proc(edge, x: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		s[i] = 0 if x[i] < edge[i] else 1;
-	}
-	return s;
-}
-
-vector_smoothstep :: proc(edge0, edge1, x: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		e0, e1 := edge0[i], edge1[i];
-		t := clamp((x[i] - e0) / (e1 - e0), 0, 1);
-		s[i] = t * t * (3 - 2*t);
-	}
-	return s;
-}
-
-vector_smootherstep :: proc(edge0, edge1, x: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	s: V;
-	for i in 0..<N {
-		e0, e1 := edge0[i], edge1[i];
-		t := clamp((x[i] - e0) / (e1 - e0), 0, 1);
-		s[i] = t * t * t * (t * (6*t - 15) + 10);
-	}
-	return s;
-}
-
-vector_distance :: proc(p0, p1: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	return length(p1 - p0);
-}
-
-vector_reflect :: proc(i, n: $V/[$N]$E) -> V where IS_NUMERIC(E) {
-	b := n * (2 * dot(n, i));
-	return i - b;
-}
-
-vector_refract :: proc(i, n: $V/[$N]$E, eta: E) -> V where IS_NUMERIC(E) {
-	dv := dot(n, i);
-	k := 1 - eta*eta - (1 - dv*dv);
-	a := i * eta;
-	b := n * eta*dv*math.sqrt(k);
-	return (a - b) * E(int(k >= 0));
-}
-
-
-
 identity :: proc($T: typeid/[$N][N]$E) -> (m: T) {
 	for i in 0..<N do m[i][i] = E(1);
 	return m;
@@ -337,6 +153,17 @@ matrix_mul :: proc(a, b: $M/[$N][N]$E) -> (c: M)
 	return;
 }
 
+matrix_comp_mul :: proc(a, b: $M/[$J][$I]$E) -> (c: M)
+	where !IS_ARRAY(E),
+	     IS_NUMERIC(E) {
+	for j in 0..<J {
+		for i in 0..<I {
+			c[j][i] = a[j][i] * b[j][i];
+		}
+	}
+	return;
+}
+
 matrix_mul_differ :: proc(a: $A/[$J][$I]$E, b: $B/[$K][J]E) -> (c: [K][I]E)
 	where !IS_ARRAY(E),
 		  IS_NUMERIC(E),
@@ -363,6 +190,9 @@ matrix_mul_vector :: proc(a: $A/[$I][$J]$E, b: $B/[I]E) -> (c: B)
 	return;
 }
 
+quaternion_mul_quaternion :: proc(q1, q2: $Q) -> Q where IS_QUATERNION(Q) {
+	return q1 * q2;
+}
 quaternion128_mul_vector3 :: proc(q: $Q/quaternion128, v: $V/[3]$F/f32) -> V {
 	Raw_Quaternion :: struct {xyz: [3]f32, r: f32};
 
@@ -390,6 +220,7 @@ mul :: proc{
 	matrix_mul_vector,
 	quaternion128_mul_vector3,
 	quaternion256_mul_vector3,
+	quaternion_mul_quaternion,
 };
 
 vector_to_ptr :: proc(v: ^$V/[$N]$E) -> ^E where IS_NUMERIC(E), N > 0 #no_bounds_check {
@@ -399,3 +230,7 @@ matrix_to_ptr :: proc(m: ^$A/[$I][$J]$E) -> ^E where IS_NUMERIC(E), I > 0, J > 0
 	return &m[0][0];
 }
 
+to_ptr :: proc{vector_to_ptr, matrix_to_ptr};
+
+
+

core/math/linalg/specific.odin

@@ -5,7 +5,7 @@ import "core:math"
 
 // Specific
 
-Float :: f32;
+Float :: f64 when #config(ODIN_MATH_LINALG_USE_F64, false) else f32;
 
 FLOAT_EPSILON :: 1e-7 when size_of(Float) == 4 else 1e-15;
 
@@ -52,25 +52,16 @@ VECTOR3_Y_AXIS :: Vector3{0, 1, 0};
 VECTOR3_Z_AXIS :: Vector3{0, 0, 1};
 
 
-radians :: proc(degrees: Float) -> Float {
-	return math.TAU * degrees / 360.0;
-}
-
-degrees :: proc(radians: Float) -> Float {
-	return 360.0 * radians / math.TAU;
-}
-
-
-vector2_orthogonal :: proc(v: Vector2) -> Vector2 {
+vector2_orthogonal :: proc(v: $V/[2]$E) -> V where !IS_ARRAY(E), IS_FLOAT(E) {
 	return {-v.y, v.x};
 }
 
-vector3_orthogonal :: proc(v: Vector3) -> Vector3 {
+vector3_orthogonal :: proc(v: $V/[3]$E) -> V where !IS_ARRAY(E), IS_FLOAT(E) {
 	x := abs(v.x);
 	y := abs(v.y);
 	z := abs(v.z);
 
-	other: Vector3;
+	other: V;
 	if x < y {
 		if x < z {
 			other = {1, 0, 0};
@@ -87,6 +78,9 @@ vector3_orthogonal :: proc(v: Vector3) -> Vector3 {
 	return normalize(cross(v, other));
 }
 
+orthogonal :: proc{vector2_orthogonal, vector3_orthogonal};
+
+
 
 vector4_srgb_to_linear :: proc(col: Vector4) -> Vector4 {
 	r := math.pow(col.x, 2.2);
@@ -178,17 +172,41 @@ vector4_rgb_to_hsl :: proc(col: Vector4) -> Vector4 {
 
 
 
-quaternion_angle_axis :: proc(angle_radians: Float, axis: Vector3) -> Quaternion {
+quaternion_angle_axis :: proc(angle_radians: Float, axis: Vector3) -> (q: Quaternion) {
 	t := angle_radians*0.5;
-	w := math.cos(t);
 	v := normalize(axis) * math.sin(t);
-	return quaternion(w, v.x, v.y, v.z);
+	q.x = v.x;
+	q.y = v.y;
+	q.z = v.z;
+	q.w = math.cos(t);
+	return;
 }
 
+angle_from_quaternion :: proc(q: Quaternion) -> Float {
+	if abs(q.w) > math.SQRT_THREE*0.5 {
+		return math.asin(q.x*q.x + q.y*q.y + q.z*q.z) * 2;
+	}
 
-quaternion_from_euler_angles :: proc(roll, pitch, yaw: Float) -> Quaternion {
-	x, y, z := roll, pitch, yaw;
-	a, b, c := x, y, z;
+	return math.cos(q.x) * 2;
+}
+
+axis_from_quaternion :: proc(q: Quaternion) -> Vector3 {
+	t1 := 1 - q.w*q.w;
+	if t1 < 0 {
+		return Vector3{0, 0, 1};
+	}
+	t2 := 1.0 / math.sqrt(t1);
+	return Vector3{q.x*t2, q.y*t2, q.z*t2};
+}
+angle_axis_from_quaternion :: proc(q: Quaternion) -> (angle: Float, axis: Vector3) {
+	angle = angle_from_quaternion(q);
+	axis  = axis_from_quaternion(q);
+	return;
+}
+
+
+quaternion_from_euler_angles :: proc(pitch, yaw, roll: Float) -> Quaternion {
+	a, b, c := pitch, yaw, roll;
 
 	ca, sa := math.cos(a*0.5), math.sin(a*0.5);
 	cb, sb := math.cos(b*0.5), math.sin(b*0.5);
@@ -202,25 +220,30 @@ quaternion_from_euler_angles :: proc(roll, pitch, yaw: Float) -> Quaternion {
 	return q;
 }
 
-euler_angles_from_quaternion :: proc(q: Quaternion) -> (roll, pitch, yaw: Float) {
-	// roll, x-axis rotation
-	sinr_cosp: Float = 2 * (q.w * q.x + q.y * q.z);
-	cosr_cosp: Float = 1 - 2 * (q.x * q.x + q.y * q.y);
-	roll = math.atan2(sinr_cosp, cosr_cosp);
+roll_from_quaternion :: proc(q: Quaternion) -> Float {
+	return math.atan2(2 * q.x*q.y + q.w*q.z, q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z);
+}
 
-	// pitch, y-axis rotation
-	sinp: Float = 2 * (q.w * q.y - q.z * q.x);
-	if abs(sinp) >= 1 {
-		pitch = math.copy_sign(math.TAU * 0.25, sinp);
-	} else {
-		pitch = 2 * math.asin(sinp);
+pitch_from_quaternion :: proc(q: Quaternion) -> Float {
+	y := 2 * (q.y*q.z + q.w*q.w);
+	x := q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z;
+
+	if abs(x) <= FLOAT_EPSILON && abs(y) <= FLOAT_EPSILON {
+		return 2 * math.atan2(q.x, q.w);
 	}
 
-	// yaw, z-axis rotation
-	siny_cosp: Float = 2 * (q.w * q.z + q.x * q.y);
-	cosy_cosp: Float = 1 - 2 * (q.y * q.y + q.z * q.z);
-	yaw = math.atan2(siny_cosp, cosy_cosp);
+	return math.atan2(y, x);
+}
+
+yaw_from_quaternion :: proc(q: Quaternion) -> Float {
+	return math.asin(clamp(-2 * (q.x*q.z - q.w*q.y), -1, 1));
+}
+
 
+euler_angles_from_quaternion :: proc(q: Quaternion) -> (pitch, yaw, roll: Float) {
+	pitch = pitch_from_quaternion(q);
+	yaw = yaw_from_quaternion(q);
+	roll = roll_from_quaternion(q);
 	return;
 }
 
@@ -269,18 +292,20 @@ quaternion_from_forward_and_up :: proc(forward, up: Vector3) -> Quaternion {
 }
 
 quaternion_look_at :: proc(eye, centre: Vector3, up: Vector3) -> Quaternion {
-	return quaternion_from_forward_and_up(centre-eye, up);
+	return quaternion_from_matrix3(matrix3_look_at(eye, centre, up));
 }
 
 
-quaternion_nlerp :: proc(a, b: Quaternion, t: Float) -> Quaternion {
-	c := a + (b-a)*quaternion(t, 0, 0, 0);
+quaternion_nlerp :: proc(a, b: Quaternion, t: Float) -> (c: Quaternion) {
+	c.x = a.x + (b.x-a.x)*t;
+	c.y = a.y + (b.y-a.y)*t;
+	c.z = a.z + (b.z-a.z)*t;
+	c.w = a.w + (b.w-a.w)*t;
 	return normalize(c);
 }
 
 
-quaternion_slerp :: proc(x, y: Quaternion, t: Float) -> Quaternion {
-
+quaternion_slerp :: proc(x, y: Quaternion, t: Float) -> (q: Quaternion) {
 	a, b := x, y;
 	cos_angle := dot(a, b);
 	if cos_angle < 0 {
@@ -288,20 +313,33 @@ quaternion_slerp :: proc(x, y: Quaternion, t: Float) -> Quaternion {
 		cos_angle = -cos_angle;
 	}
 	if cos_angle > 1 - FLOAT_EPSILON {
-		return a + (b-a)*quaternion(t, 0, 0, 0);
+		q.x = a.x + (b.x-a.x)*t;
+		q.y = a.y + (b.y-a.y)*t;
+		q.z = a.z + (b.z-a.z)*t;
+		q.w = a.w + (b.w-a.w)*t;
+		return;
 	}
 
 	angle := math.acos(cos_angle);
 	sin_angle := math.sin(angle);
-	factor_a, factor_b: Quaternion;
-	factor_a = quaternion(math.sin((1-t) * angle) / sin_angle, 0, 0, 0);
-	factor_b = quaternion(math.sin(t * angle)     / sin_angle, 0, 0, 0);
+	factor_a := math.sin((1-t) * angle) / sin_angle;
+	factor_b := math.sin(t * angle)     / sin_angle;
+
+
+	q.x = factor_a * a.x + factor_b * b.x;
+	q.y = factor_a * a.y + factor_b * b.y;
+	q.z = factor_a * a.z + factor_b * b.z;
+	q.w = factor_a * a.w + factor_b * b.w;
+	return;
+}
 
-	return factor_a * a + factor_b * b;
+quaternion_squad :: proc(q1, q2, s1, s2: Quaternion, h: Float) -> Quaternion {
+	slerp :: quaternion_slerp;
+	return slerp(slerp(q1, q2, h), slerp(s1, s2, h), 2 * (1 - h) * h);
 }
 
 
-quaternion_from_matrix4 :: proc(m: Matrix4) -> Quaternion {
+quaternion_from_matrix4 :: proc(m: Matrix4) -> (q: Quaternion) {
 	four_x_squared_minus_1, four_y_squared_minus_1,
 	four_z_squared_minus_1, four_w_squared_minus_1,
 	four_biggest_squared_minus_1: Float;
@@ -336,40 +374,32 @@ quaternion_from_matrix4 :: proc(m: Matrix4) -> Quaternion {
 
 	switch biggest_index {
 	case 0:
-		return quaternion(
-			biggest_value,
-			(m[0][1] + m[1][0]) * mult,
-			(m[2][0] + m[0][2]) * mult,
-			(m[1][2] - m[2][1]) * mult,
-		);
+		q.w = biggest_value;
+		q.x = (m[0][1] + m[1][0]) * mult;
+		q.y = (m[2][0] + m[0][2]) * mult;
+		q.z = (m[1][2] - m[2][1]) * mult;
 	case 1:
-		return quaternion(
-			(m[0][1] + m[1][0]) * mult,
-			biggest_value,
-			(m[1][2] + m[2][1]) * mult,
-			(m[2][0] - m[0][2]) * mult,
-		);
+		q.w = (m[0][1] + m[1][0]) * mult;
+		q.x = biggest_value;
+		q.y = (m[1][2] + m[2][1]) * mult;
+		q.z = (m[2][0] - m[0][2]) * mult;
 	case 2:
-		return quaternion(
-			(m[2][0] + m[0][2]) * mult,
-			(m[1][2] + m[2][1]) * mult,
-			biggest_value,
-			(m[0][1] - m[1][0]) * mult,
-		);
+		q.w = (m[2][0] + m[0][2]) * mult;
+		q.x = (m[1][2] + m[2][1]) * mult;
+		q.y = biggest_value;
+		q.z = (m[0][1] - m[1][0]) * mult;
 	case 3:
-		return quaternion(
-			(m[1][2] - m[2][1]) * mult,
-			(m[2][0] - m[0][2]) * mult,
-			(m[0][1] - m[1][0]) * mult,
-			biggest_value,
-		);
+		q.w = (m[1][2] - m[2][1]) * mult;
+		q.x = (m[2][0] - m[0][2]) * mult;
+		q.y = (m[0][1] - m[1][0]) * mult;
+		q.z = biggest_value;
 	}
 
-	return 0;
+	return;
 }
 
 
-quaternion_from_matrix3 :: proc(m: Matrix3) -> Quaternion {
+quaternion_from_matrix3 :: proc(m: Matrix3) -> (q: Quaternion) {
 	four_x_squared_minus_1, four_y_squared_minus_1,
 	four_z_squared_minus_1, four_w_squared_minus_1,
 	four_biggest_squared_minus_1: Float;
@@ -404,55 +434,54 @@ quaternion_from_matrix3 :: proc(m: Matrix3) -> Quaternion {
 
 	switch biggest_index {
 	case 0:
-		return quaternion(
-			biggest_value,
-			(m[0][1] + m[1][0]) * mult,
-			(m[2][0] + m[0][2]) * mult,
-			(m[1][2] - m[2][1]) * mult,
-		);
+		q.w = biggest_value;
+		q.x = (m[0][1] + m[1][0]) * mult;
+		q.y = (m[2][0] + m[0][2]) * mult;
+		q.z = (m[1][2] - m[2][1]) * mult;
 	case 1:
-		return quaternion(
-			(m[0][1] + m[1][0]) * mult,
-			biggest_value,
-			(m[1][2] + m[2][1]) * mult,
-			(m[2][0] - m[0][2]) * mult,
-		);
+		q.w = (m[0][1] + m[1][0]) * mult;
+		q.x = biggest_value;
+		q.y = (m[1][2] + m[2][1]) * mult;
+		q.z = (m[2][0] - m[0][2]) * mult;
 	case 2:
-		return quaternion(
-			(m[2][0] + m[0][2]) * mult,
-			(m[1][2] + m[2][1]) * mult,
-			biggest_value,
-			(m[0][1] - m[1][0]) * mult,
-		);
+		q.w = (m[2][0] + m[0][2]) * mult;
+		q.x = (m[1][2] + m[2][1]) * mult;
+		q.y = biggest_value;
+		q.z = (m[0][1] - m[1][0]) * mult;
 	case 3:
-		return quaternion(
-			(m[1][2] - m[2][1]) * mult,
-			(m[2][0] - m[0][2]) * mult,
-			(m[0][1] - m[1][0]) * mult,
-			biggest_value,
-		);
+		q.w = (m[1][2] - m[2][1]) * mult;
+		q.x = (m[2][0] - m[0][2]) * mult;
+		q.y = (m[0][1] - m[1][0]) * mult;
+		q.z = biggest_value;
 	}
 
-	return 0;
+	return;
 }
 
-quaternion_between_two_vector3 :: proc(from, to: Vector3) -> Quaternion {
+quaternion_between_two_vector3 :: proc(from, to: Vector3) -> (q: Quaternion) {
 	x := normalize(from);
 	y := normalize(to);
 
 	cos_theta := dot(x, y);
 	if abs(cos_theta + 1) < 2*FLOAT_EPSILON {
 		v := vector3_orthogonal(x);
-		return quaternion(0, v.x, v.y, v.z);
+		q.x = v.x;
+		q.y = v.y;
+		q.z = v.z;
+		q.w = 0;
+		return;
 	}
 	v := cross(x, y);
 	w := cos_theta + 1;
-	return Quaternion(normalize(quaternion(w, v.x, v.y, v.z)));
+	q.w = w;
+	q.x = v.x;
+	q.y = v.y;
+	q.z = v.z;
+	return normalize(q);
 }
 
 
-matrix2_inverse_transpose :: proc(m: Matrix2) -> Matrix2 {
-	c: Matrix2;
+matrix2_inverse_transpose :: proc(m: Matrix2) -> (c: Matrix2) {
 	d := m[0][0]*m[1][1] - m[1][0]*m[0][1];
 	id := 1.0/d;
 	c[0][0] = +m[1][1] * id;
@@ -464,8 +493,7 @@ matrix2_inverse_transpose :: proc(m: Matrix2) -> Matrix2 {
 matrix2_determinant :: proc(m: Matrix2) -> Float {
 	return m[0][0]*m[1][1] - m[1][0]*m[0][1];
 }
-matrix2_inverse :: proc(m: Matrix2) -> Matrix2 {
-	c: Matrix2;
+matrix2_inverse :: proc(m: Matrix2) -> (c: Matrix2) {
 	d := m[0][0]*m[1][1] - m[1][0]*m[0][1];
 	id := 1.0/d;
 	c[0][0] = +m[1][1] * id;
@@ -475,8 +503,7 @@ matrix2_inverse :: proc(m: Matrix2) -> Matrix2 {
 	return c;
 }
 
-matrix2_adjoint :: proc(m: Matrix2) -> Matrix2 {
-	c: Matrix2;
+matrix2_adjoint :: proc(m: Matrix2) -> (c: Matrix2) {
 	c[0][0] = +m[1][1];
 	c[0][1] = -m[1][0];
 	c[1][0] = -m[0][1];
@@ -485,7 +512,7 @@ matrix2_adjoint :: proc(m: Matrix2) -> Matrix2 {
 }
 
 
-matrix3_from_quaternion :: proc(q: Quaternion) -> Matrix3 {
+matrix3_from_quaternion :: proc(q: Quaternion) -> (m: Matrix3) {
 	xx := q.x * q.x;
 	xy := q.x * q.y;
 	xz := q.x * q.z;
@@ -496,7 +523,6 @@ matrix3_from_quaternion :: proc(q: Quaternion) -> Matrix3 {
 	zz := q.z * q.z;
 	zw := q.z * q.w;
 
-	m: Matrix3;
 	m[0][0] = 1 - 2 * (yy + zz);
 	m[1][0] = 2 * (xy - zw);
 	m[2][0] = 2 * (xz + yw);
@@ -524,8 +550,7 @@ matrix3_determinant :: proc(m: Matrix3) -> Float {
 	return a + b + c;
 }
 
-matrix3_adjoint :: proc(m: Matrix3) -> Matrix3 {
-	adjoint: Matrix3;
+matrix3_adjoint :: proc(m: Matrix3) -> (adjoint: Matrix3) {
 	adjoint[0][0] = +(m[1][1] * m[2][2] - m[1][2] * m[2][1]);
 	adjoint[1][0] = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]);
 	adjoint[2][0] = +(m[0][1] * m[1][2] - m[0][2] * m[1][1]);
@@ -553,8 +578,7 @@ matrix3_inverse_transpose :: proc(m: Matrix3) -> Matrix3 {
 }
 
 
-matrix3_scale :: proc(s: Vector3) -> Matrix3 {
-	m: Matrix3;
+matrix3_scale :: proc(s: Vector3) -> (m: Matrix3) {
 	m[0][0] = s[0];
 	m[1][1] = s[1];
 	m[2][2] = s[2];
@@ -597,7 +621,7 @@ matrix3_look_at :: proc(eye, centre, up: Vector3) -> Matrix3 {
 }
 
 matrix4_from_quaternion :: proc(q: Quaternion) -> Matrix4 {
-	m := identity(Matrix4);
+	m := MATRIX4_IDENTITY;
 
 	xx := q.x * q.x;
 	xy := q.x * q.y;
@@ -693,7 +717,7 @@ matrix4_inverse_transpose :: proc(m: Matrix4) -> Matrix4 {
 }
 
 matrix4_translate :: proc(v: Vector3) -> Matrix4 {
-	m := identity(Matrix4);
+	m := MATRIX4_IDENTITY;
 	m[3][0] = v[0];
 	m[3][1] = v[1];
 	m[3][2] = v[2];
@@ -708,7 +732,7 @@ matrix4_rotate :: proc(angle_radians: Float, v: Vector3) -> Matrix4 {
 	a := normalize(v);
 	t := a * (1-c);
 
-	rot := identity(Matrix4);
+	rot := MATRIX4_IDENTITY;
 
 	rot[0][0] = c + t[0]*a[0];
 	rot[0][1] = 0 + t[0]*a[1] + s*a[2];
@@ -737,16 +761,20 @@ matrix4_scale :: proc(v: Vector3) -> Matrix4 {
 	return m;
 }
 
-matrix4_look_at :: proc(eye, centre, up: Vector3) -> Matrix4 {
+matrix4_look_at :: proc(eye, centre, up: Vector3, flip_z_axis := true) -> Matrix4 {
 	f := normalize(centre - eye);
 	s := normalize(cross(f, up));
 	u := cross(s, f);
-	return Matrix4{
+
+	fe := dot(f, eye);
+
+	m := Matrix4{
 		{+s.x, +u.x, -f.x, 0},
 		{+s.y, +u.y, -f.y, 0},
 		{+s.z, +u.z, -f.z, 0},
-		{-dot(s, eye), -dot(u, eye), +dot(f, eye), 1},
+		{-dot(s, eye), -dot(u, eye), +fe if flip_z_axis else -fe, 1},
 	};
+	return m;
 }
 
 
@@ -797,3 +825,5 @@ matrix4_infinite_perspective :: proc(fovy, aspect, near: Float, flip_z_axis := t
 
 	return;
 }
+
+

core/math/linalg/specific_euler_angles.odin

@@ -0,0 +1,624 @@
+package linalg
+
+import "core:math"
+
+
+euler_angle_x :: proc(angle_x: Float) -> (m: Matrix4) {
+	cos_x, sin_x := math.cos(angle_x), math.sin(angle_x);
+	m[0][0] = 1;
+	m[1][1] = +cos_x;
+	m[2][1] = +sin_x;
+	m[1][2] = -sin_x;
+	m[2][2] = +cos_x;
+	m[3][3] = 1;
+	return;
+}
+euler_angle_y :: proc(angle_y: Float) -> (m: Matrix4) {
+	cos_y, sin_y := math.cos(angle_y), math.sin(angle_y);
+	m[0][0] = +cos_y;
+	m[2][0] = -sin_y;
+	m[1][1] = 1;
+	m[0][2] = +sin_y;
+	m[2][2] = +cos_y;
+	m[3][3] = 1;
+	return;
+}
+euler_angle_z :: proc(angle_z: Float) -> (m: Matrix4) {
+	cos_z, sin_z := math.cos(angle_z), math.sin(angle_z);
+	m[0][0] = +cos_z;
+	m[1][0] = +sin_z;
+	m[1][1] = +cos_z;
+	m[0][1] = -sin_z;
+	m[2][2] = 1;
+	m[3][3] = 1;
+	return;
+}
+
+
+derived_euler_angle_x :: proc(angle_x: Float, angular_velocity_x: Float) -> (m: Matrix4) {
+	cos_x := math.cos(angle_x) * angular_velocity_x;
+	sin_x := math.sin(angle_x) * angular_velocity_x;
+	m[0][0] = 1;
+	m[1][1] = +cos_x;
+	m[2][1] = +sin_x;
+	m[1][2] = -sin_x;
+	m[2][2] = +cos_x;
+	m[3][3] = 1;
+	return;
+}
+derived_euler_angle_y :: proc(angle_y: Float, angular_velocity_y: Float) -> (m: Matrix4) {
+	cos_y := math.cos(angle_y) * angular_velocity_y;
+	sin_y := math.sin(angle_y) * angular_velocity_y;
+	m[0][0] = +cos_y;
+	m[2][0] = -sin_y;
+	m[1][1] = 1;
+	m[0][2] = +sin_y;
+	m[2][2] = +cos_y;
+	m[3][3] = 1;
+	return;
+}
+derived_euler_angle_z :: proc(angle_z: Float, angular_velocity_z: Float) -> (m: Matrix4) {
+	cos_z := math.cos(angle_z) * angular_velocity_z;
+	sin_z := math.sin(angle_z) * angular_velocity_z;
+	m[0][0] = +cos_z;
+	m[1][0] = +sin_z;
+	m[1][1] = +cos_z;
+	m[0][1] = -sin_z;
+	m[2][2] = 1;
+	m[3][3] = 1;
+	return;
+}
+
+
+euler_angle_xy :: proc(angle_x, angle_y: Float) -> (m: Matrix4) {
+	cos_x, sin_x := math.cos(angle_x), math.sin(angle_x);
+	cos_y, sin_y := math.cos(angle_y), math.sin(angle_y);
+	m[0][0] = cos_y;
+	m[1][0] = -sin_x * - sin_y;
+	m[2][0] = -cos_x * - sin_y;
+	m[1][1] = cos_x;
+	m[2][1] = sin_x;
+	m[0][2] = sin_y;
+	m[1][2] = -sin_x * cos_y;
+	m[2][2] = cos_x * cos_y;
+	m[3][3] = 1;
+	return;
+}
+
+
+euler_angle_yx :: proc(angle_y, angle_x: Float) -> (m: Matrix4) {
+	cos_x, sin_x := math.cos(angle_x), math.sin(angle_x);
+	cos_y, sin_y := math.cos(angle_y), math.sin(angle_y);
+	m[0][0] = cos_y;
+	m[2][0] = -sin_y;
+	m[0][1] = sin_y*sin_x;
+	m[1][1] = cos_x;
+	m[2][1] = cos_y*sin_x;
+	m[0][2] = sin_y*cos_x;
+	m[1][2] = -sin_x;
+	m[2][2] = cos_y*cos_x;
+	m[3][3] = 1;
+	return;
+}
+
+euler_angle_xz :: proc(angle_x, angle_z: Float) -> (m: Matrix4) {
+	return mul(euler_angle_x(angle_x), euler_angle_z(angle_z));
+}
+euler_angle_zx :: proc(angle_z, angle_x: Float) -> (m: Matrix4) {
+	return mul(euler_angle_z(angle_z), euler_angle_x(angle_x));
+}
+euler_angle_yz :: proc(angle_y, angle_z: Float) -> (m: Matrix4) {
+	return mul(euler_angle_y(angle_y), euler_angle_z(angle_z));
+}
+euler_angle_zy :: proc(angle_z, angle_y: Float) -> (m: Matrix4) {
+	return mul(euler_angle_z(angle_z), euler_angle_y(angle_y));
+}
+
+
+euler_angle_xyz :: proc(t1, t2, t3: Float) -> (m: Matrix4) {
+	c1 := math.cos(-t1);
+	c2 := math.cos(-t2);
+	c3 := math.cos(-t3);
+	s1 := math.sin(-t1);
+	s2 := math.sin(-t2);
+	s3 := math.sin(-t3);
+
+	m[0][0] = c2 * c3;
+	m[0][1] =-c1 * s3 + s1 * s2 * c3;
+	m[0][2] = s1 * s3 + c1 * s2 * c3;
+	m[0][3] = 0;
+	m[1][0] = c2 * s3;
+	m[1][1] = c1 * c3 + s1 * s2 * s3;
+	m[1][2] =-s1 * c3 + c1 * s2 * s3;
+	m[1][3] = 0;
+	m[2][0] =-s2;
+	m[2][1] = s1 * c2;
+	m[2][2] = c1 * c2;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return;
+}
+
+euler_angle_yxz :: proc(yaw, pitch, roll: Float) -> (m: Matrix4) {
+	ch := math.cos(yaw);
+	sh := math.sin(yaw);
+	cp := math.cos(pitch);
+	sp := math.sin(pitch);
+	cb := math.cos(roll);
+	sb := math.sin(roll);
+
+	m[0][0] = ch * cb + sh * sp * sb;
+	m[0][1] = sb * cp;
+	m[0][2] = -sh * cb + ch * sp * sb;
+	m[0][3] = 0;
+	m[1][0] = -ch * sb + sh * sp * cb;
+	m[1][1] = cb * cp;
+	m[1][2] = sb * sh + ch * sp * cb;
+	m[1][3] = 0;
+	m[2][0] = sh * cp;
+	m[2][1] = -sp;
+	m[2][2] = ch * cp;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return;
+}
+
+euler_angle_xzx :: proc(t1, t2, t3: Float) -> (m: Matrix4) {
+	c1 := math.cos(t1);
+	s1 := math.sin(t1);
+	c2 := math.cos(t2);
+	s2 := math.sin(t2);
+	c3 := math.cos(t3);
+	s3 := math.sin(t3);
+
+	m[0][0] = c2;
+	m[0][1] = c1 * s2;
+	m[0][2] = s1 * s2;
+	m[0][3] = 0;
+	m[1][0] =-c3 * s2;
+	m[1][1] = c1 * c2 * c3 - s1 * s3;
+	m[1][2] = c1 * s3 + c2 * c3 * s1;
+	m[1][3] = 0;
+	m[2][0] = s2 * s3;
+	m[2][1] =-c3 * s1 - c1 * c2 * s3;
+	m[2][2] = c1 * c3 - c2 * s1 * s3;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return;
+}
+
+euler_angle_xyx :: proc(t1, t2, t3: Float) -> (m: Matrix4) {
+	c1 := math.cos(t1);
+	s1 := math.sin(t1);
+	c2 := math.cos(t2);
+	s2 := math.sin(t2);
+	c3 := math.cos(t3);
+	s3 := math.sin(t3);
+
+	m[0][0] = c2;
+	m[0][1] = s1 * s2;
+	m[0][2] =-c1 * s2;
+	m[0][3] = 0;
+	m[1][0] = s2 * s3;
+	m[1][1] = c1 * c3 - c2 * s1 * s3;
+	m[1][2] = c3 * s1 + c1 * c2 * s3;
+	m[1][3] = 0;
+	m[2][0] = c3 * s2;
+	m[2][1] =-c1 * s3 - c2 * c3 * s1;
+	m[2][2] = c1 * c2 * c3 - s1 * s3;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return;
+}
+
+euler_angle_yxy :: proc(t1, t2, t3: Float) -> (m: Matrix4) {
+	c1 := math.cos(t1);
+	s1 := math.sin(t1);
+	c2 := math.cos(t2);
+	s2 := math.sin(t2);
+	c3 := math.cos(t3);
+	s3 := math.sin(t3);
+
+	m[0][0] = c1 * c3 - c2 * s1 * s3;
+	m[0][1] = s2* s3;
+	m[0][2] =-c3 * s1 - c1 * c2 * s3;
+	m[0][3] = 0;
+	m[1][0] = s1 * s2;
+	m[1][1] = c2;
+	m[1][2] = c1 * s2;
+	m[1][3] = 0;
+	m[2][0] = c1 * s3 + c2 * c3 * s1;
+	m[2][1] =-c3 * s2;
+	m[2][2] = c1 * c2 * c3 - s1 * s3;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return;
+}
+
+euler_angle_yzy :: proc(t1, t2, t3: Float) -> (m: Matrix4) {
+	c1 := math.cos(t1);
+	s1 := math.sin(t1);
+	c2 := math.cos(t2);
+	s2 := math.sin(t2);
+	c3 := math.cos(t3);
+	s3 := math.sin(t3);
+
+	m[0][0] = c1 * c2 * c3 - s1 * s3;
+	m[0][1] = c3 * s2;
+	m[0][2] =-c1 * s3 - c2 * c3 * s1;
+	m[0][3] = 0;
+	m[1][0] =-c1 * s2;
+	m[1][1] = c2;
+	m[1][2] = s1 * s2;
+	m[1][3] = 0;
+	m[2][0] = c3 * s1 + c1 * c2 * s3;
+	m[2][1] = s2 * s3;
+	m[2][2] = c1 * c3 - c2 * s1 * s3;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return;
+}
+
+euler_angle_zyz :: proc(t1, t2, t3: Float) -> (m: Matrix4) {
+	c1 := math.cos(t1);
+	s1 := math.sin(t1);
+	c2 := math.cos(t2);
+	s2 := math.sin(t2);
+	c3 := math.cos(t3);
+	s3 := math.sin(t3);
+
+	m[0][0] = c1 * c2 * c3 - s1 * s3;
+	m[0][1] = c1 * s3 + c2 * c3 * s1;
+	m[0][2] =-c3 * s2;
+	m[0][3] = 0;
+	m[1][0] =-c3 * s1 - c1 * c2 * s3;
+	m[1][1] = c1 * c3 - c2 * s1 * s3;
+	m[1][2] = s2 * s3;
+	m[1][3] = 0;
+	m[2][0] = c1 * s2;
+	m[2][1] = s1 * s2;
+	m[2][2] = c2;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return;
+}
+
+euler_angle_zxz :: proc(t1, t2, t3: Float) -> (m: Matrix4) {
+	c1 := math.cos(t1);
+	s1 := math.sin(t1);
+	c2 := math.cos(t2);
+	s2 := math.sin(t2);
+	c3 := math.cos(t3);
+	s3 := math.sin(t3);
+
+	m[0][0] = c1 * c3 - c2 * s1 * s3;
+	m[0][1] = c3 * s1 + c1 * c2 * s3;
+	m[0][2] = s2 *s3;
+	m[0][3] = 0;
+	m[1][0] =-c1 * s3 - c2 * c3 * s1;
+	m[1][1] = c1 * c2 * c3 - s1 * s3;
+	m[1][2] = c3 * s2;
+	m[1][3] = 0;
+	m[2][0] = s1 * s2;
+	m[2][1] =-c1 * s2;
+	m[2][2] = c2;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return;
+}
+
+
+euler_angle_xzy :: proc(t1, t2, t3: Float) -> (m: Matrix4) {
+	c1 := math.cos(t1);
+	s1 := math.sin(t1);
+	c2 := math.cos(t2);
+	s2 := math.sin(t2);
+	c3 := math.cos(t3);
+	s3 := math.sin(t3);
+
+	m[0][0] = c2 * c3;
+	m[0][1] = s1 * s3 + c1 * c3 * s2;
+	m[0][2] = c3 * s1 * s2 - c1 * s3;
+	m[0][3] = 0;
+	m[1][0] =-s2;
+	m[1][1] = c1 * c2;
+	m[1][2] = c2 * s1;
+	m[1][3] = 0;
+	m[2][0] = c2 * s3;
+	m[2][1] = c1 * s2 * s3 - c3 * s1;
+	m[2][2] = c1 * c3 + s1 * s2 *s3;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return;
+}
+
+euler_angle_yzx :: proc(t1, t2, t3: Float) -> (m: Matrix4) {
+	c1 := math.cos(t1);
+	s1 := math.sin(t1);
+	c2 := math.cos(t2);
+	s2 := math.sin(t2);
+	c3 := math.cos(t3);
+	s3 := math.sin(t3);
+
+	m[0][0] = c1 * c2;
+	m[0][1] = s2;
+	m[0][2] =-c2 * s1;
+	m[0][3] = 0;
+	m[1][0] = s1 * s3 - c1 * c3 * s2;
+	m[1][1] = c2 * c3;
+	m[1][2] = c1 * s3 + c3 * s1 * s2;
+	m[1][3] = 0;
+	m[2][0] = c3 * s1 + c1 * s2 * s3;
+	m[2][1] =-c2 * s3;
+	m[2][2] = c1 * c3 - s1 * s2 * s3;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return;
+}
+
+euler_angle_zyx :: proc(t1, t2, t3: Float) -> (m: Matrix4) {
+	c1 := math.cos(t1);
+	s1 := math.sin(t1);
+	c2 := math.cos(t2);
+	s2 := math.sin(t2);
+	c3 := math.cos(t3);
+	s3 := math.sin(t3);
+
+	m[0][0] = c1 * c2;
+	m[0][1] = c2 * s1;
+	m[0][2] =-s2;
+	m[0][3] = 0;
+	m[1][0] = c1 * s2 * s3 - c3 * s1;
+	m[1][1] = c1 * c3 + s1 * s2 * s3;
+	m[1][2] = c2 * s3;
+	m[1][3] = 0;
+	m[2][0] = s1 * s3 + c1 * c3 * s2;
+	m[2][1] = c3 * s1 * s2 - c1 * s3;
+	m[2][2] = c2 * c3;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return;
+}
+
+euler_angle_zxy :: proc(t1, t2, t3: Float) -> (m: Matrix4) {
+	c1 := math.cos(t1);
+	s1 := math.sin(t1);
+	c2 := math.cos(t2);
+	s2 := math.sin(t2);
+	c3 := math.cos(t3);
+	s3 := math.sin(t3);
+
+	m[0][0] = c1 * c3 - s1 * s2 * s3;
+	m[0][1] = c3 * s1 + c1 * s2 * s3;
+	m[0][2] =-c2 * s3;
+	m[0][3] = 0;
+	m[1][0] =-c2 * s1;
+	m[1][1] = c1 * c2;
+	m[1][2] = s2;
+	m[1][3] = 0;
+	m[2][0] = c1 * s3 + c3 * s1 * s2;
+	m[2][1] = s1 * s3 - c1 * c3 * s2;
+	m[2][2] = c2 * c3;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return;
+}
+
+
+yaw_pitch_roll :: proc(yaw, pitch, roll: Float) -> (m: Matrix4) {
+	ch := math.cos(yaw);
+	sh := math.sin(yaw);
+	cp := math.cos(pitch);
+	sp := math.sin(pitch);
+	cb := math.cos(roll);
+	sb := math.sin(roll);
+
+	m[0][0] = ch * cb + sh * sp * sb;
+	m[0][1] = sb * cp;
+	m[0][2] = -sh * cb + ch * sp * sb;
+	m[0][3] = 0;
+	m[1][0] = -ch * sb + sh * sp * cb;
+	m[1][1] = cb * cp;
+	m[1][2] = sb * sh + ch * sp * cb;
+	m[1][3] = 0;
+	m[2][0] = sh * cp;
+	m[2][1] = -sp;
+	m[2][2] = ch * cp;
+	m[2][3] = 0;
+	m[3][0] = 0;
+	m[3][1] = 0;
+	m[3][2] = 0;
+	m[3][3] = 1;
+	return m;
+}
+
+extract_euler_angle_xyz :: proc(m: Matrix4) -> (t1, t2, t3: Float) {
+	T1 := math.atan2(m[2][1], m[2][2]);
+	C2 := math.sqrt(m[0][0]*m[0][0] + m[1][0]*m[1][0]);
+	T2 := math.atan2(-m[2][0], C2);
+	S1 := math.sin(T1);
+	C1 := math.cos(T1);
+	T3 := math.atan2(S1*m[0][2] - C1*m[0][1], C1*m[1][1] - S1*m[1][2]);
+	t1 = -T1;
+	t2 = -T2;
+	t3 = -T3;
+	return;
+}
+
+extract_euler_angle_yxz :: proc(m: Matrix4) -> (t1, t2, t3: Float) {
+	T1 := math.atan2(m[2][0], m[2][2]);
+	C2 := math.sqrt(m[0][1]*m[0][1] + m[1][1]*m[1][1]);
+	T2 := math.atan2(-m[2][1], C2);
+	S1 := math.sin(T1);
+	C1 := math.cos(T1);
+	T3 := math.atan2(S1*m[1][2] - C1*m[1][0], C1*m[0][0] - S1*m[0][2]);
+	t1 = T1;
+	t2 = T2;
+	t3 = T3;
+	return;
+}
+
+extract_euler_angle_xzx :: proc(m: Matrix4) -> (t1, t2, t3: Float) {
+	T1 := math.atan2(m[0][2], m[0][1]);
+	S2 := math.sqrt(m[1][0]*m[1][0] + m[2][0]*m[2][0]);
+	T2 := math.atan2(S2, m[0][0]);
+	S1 := math.sin(T1);
+	C1 := math.cos(T1);
+	T3 := math.atan2(C1*m[1][2] - S1*m[1][1], C1*m[2][2] - S1*m[2][1]);
+	t1 = T1;
+	t2 = T2;
+	t3 = T3;
+	return;
+}
+
+extract_euler_angle_xyx :: proc(m: Matrix4) -> (t1, t2, t3: Float) {
+	T1 := math.atan2(m[0][1], -m[0][2]);
+	S2 := math.sqrt(m[1][0]*m[1][0] + m[2][0]*m[2][0]);
+	T2 := math.atan2(S2, m[0][0]);
+	S1 := math.sin(T1);
+	C1 := math.cos(T1);
+	T3 := math.atan2(-C1*m[2][1] - S1*m[2][2], C1*m[1][1] + S1*m[1][2]);
+	t1 = T1;
+	t2 = T2;
+	t3 = T3;
+	return;
+}
+
+extract_euler_angle_yxy :: proc(m: Matrix4) -> (t1, t2, t3: Float) {
+	T1 := math.atan2(m[1][0], m[1][2]);
+	S2 := math.sqrt(m[0][1]*m[0][1] + m[2][1]*m[2][1]);
+	T2 := math.atan2(S2, m[1][1]);
+	S1 := math.sin(T1);
+	C1 := math.cos(T1);
+	T3 := math.atan2(C1*m[2][0] - S1*m[2][2], C1*m[0][0] - S1*m[0][2]);
+	t1 = T1;
+	t2 = T2;
+	t3 = T3;
+	return;
+}
+
+extract_euler_angle_yzy :: proc(m: Matrix4) -> (t1, t2, t3: Float) {
+	T1 := math.atan2(m[1][2], -m[1][0]);
+	S2 := math.sqrt(m[0][1]*m[0][1] + m[2][1]*m[2][1]);
+	T2 := math.atan2(S2, m[1][1]);
+	S1 := math.sin(T1);
+	C1 := math.cos(T1);
+	T3 := math.atan2(-S1*m[0][0] - C1*m[0][2], S1*m[2][0] + C1*m[2][2]);
+	t1 = T1;
+	t2 = T2;
+	t3 = T3;
+	return;
+}
+extract_euler_angle_zyz :: proc(m: Matrix4) -> (t1, t2, t3: Float) {
+	T1 := math.atan2(m[2][1], m[2][0]);
+	S2 := math.sqrt(m[0][2]*m[0][2] + m[1][2]*m[1][2]);
+	T2 := math.atan2(S2, m[2][2]);
+	S1 := math.sin(T1);
+	C1 := math.cos(T1);
+	T3 := math.atan2(C1*m[0][1] - S1*m[0][0], C1*m[1][1] - S1*m[1][0]);
+	t1 = T1;
+	t2 = T2;
+	t3 = T3;
+	return;
+}
+
+extract_euler_angle_zxz :: proc(m: Matrix4) -> (t1, t2, t3: Float) {
+	T1 := math.atan2(m[2][0], -m[2][1]);
+	S2 := math.sqrt(m[0][2]*m[0][2] + m[1][2]*m[1][2]);
+	T2 := math.atan2(S2, m[2][2]);
+	S1 := math.sin(T1);
+	C1 := math.cos(T1);
+	T3 := math.atan2(-C1*m[1][0] - S1*m[1][1], C1*m[0][0] + S1*m[0][1]);
+	t1 = T1;
+	t2 = T2;
+	t3 = T3;
+	return;
+}
+
+extract_euler_angle_xzy :: proc(m: Matrix4) -> (t1, t2, t3: Float) {
+	T1 := math.atan2(m[1][2], m[1][1]);
+	C2 := math.sqrt(m[0][0]*m[0][0] + m[2][0]*m[2][0]);
+	T2 := math.atan2(-m[1][0], C2);
+	S1 := math.sin(T1);
+	C1 := math.cos(T1);
+	T3 := math.atan2(S1*m[0][1] - C1*m[0][2], C1*m[2][2] - S1*m[2][1]);
+	t1 = T1;
+	t2 = T2;
+	t3 = T3;
+	return;
+}
+
+extract_euler_angle_yzx :: proc(m: Matrix4) -> (t1, t2, t3: Float) {
+	T1 := math.atan2(-m[0][2], m[0][0]);
+	C2 := math.sqrt(m[1][1]*m[1][1] + m[2][1]*m[2][1]);
+	T2 := math.atan2(m[0][1], C2);
+	S1 := math.sin(T1);
+	C1 := math.cos(T1);
+	T3 := math.atan2(S1*m[1][0] + C1*m[1][2], S1*m[2][0] + C1*m[2][2]);
+	t1 = T1;
+	t2 = T2;
+	t3 = T3;
+	return;
+}
+
+extract_euler_angle_zyx :: proc(m: Matrix4) -> (t1, t2, t3: Float) {
+	T1 := math.atan2(m[0][1], m[0][0]);
+	C2 := math.sqrt(m[1][2]*m[1][2] + m[2][2]*m[2][2]);
+	T2 := math.atan2(-m[0][2], C2);
+	S1 := math.sin(T1);
+	C1 := math.cos(T1);
+	T3 := math.atan2(S1*m[2][0] - C1*m[2][1], C1*m[1][1] - S1*m[1][0]);
+	t1 = T1;
+	t2 = T2;
+	t3 = T3;
+	return;
+}
+
+extract_euler_angle_zxy :: proc(m: Matrix4) -> (t1, t2, t3: Float) {
+	T1 := math.atan2(-m[1][0], m[1][1]);
+	C2 := math.sqrt(m[0][2]*m[0][2] + m[2][2]*m[2][2]);
+	T2 := math.atan2(m[1][2], C2);
+	S1 := math.sin(T1);
+	C1 := math.cos(T1);
+	T3 := math.atan2(C1*m[2][0] + S1*m[2][1], C1*m[0][0] + S1*m[0][1]);
+	t1 = T1;
+	t2 = T2;
+	t3 = T3;
+	return;
+}

core/math/math.odin

@@ -36,7 +36,7 @@ RAD_PER_DEG :: TAU/360.0;
 DEG_PER_RAD :: 360.0/TAU;
 
 
-@(default_calling_convention="none")
+@(default_calling_convention="pure_none")
 foreign _ {
 	@(link_name="llvm.sqrt.f32")
 	sqrt_f32 :: proc(x: f32) -> f32 ---;
@@ -103,11 +103,12 @@ log10_f64 :: proc(x: f64) -> f64 { return ln(x)/LN10; }
 log10     :: proc{log10_f32, log10_f64};
 
 
-tan_f32 :: proc "c" (θ: f32) -> f32 { return sin(θ)/cos(θ); }
-tan_f64 :: proc "c" (θ: f64) -> f64 { return sin(θ)/cos(θ); }
+tan_f32 :: proc(θ: f32) -> f32 { return sin(θ)/cos(θ); }
+tan_f64 :: proc(θ: f64) -> f64 { return sin(θ)/cos(θ); }
 tan     :: proc{tan_f32, tan_f64};
 
 lerp :: proc(a, b: $T, t: $E) -> (x: T) { return a*(1-t) + b*t; }
+saturate :: proc(a: $T) -> (x: T) { return clamp(a, 0, 1); };
 
 unlerp_f32 :: proc(a, b, x: f32) -> (t: f32) { return (x-a)/(b-a); }
 unlerp_f64 :: proc(a, b, x: f64) -> (t: f64) { return (x-a)/(b-a); }

core/runtime/core.odin

@@ -29,9 +29,13 @@ Calling_Convention :: enum u8 {
 	Invalid     = 0,
 	Odin        = 1,
 	Contextless = 2,
-	C           = 3,
-	Std         = 4,
-	Fast        = 5,
+	Pure        = 3,
+	CDecl       = 4,
+	Std_Call    = 5,
+	Fast_Call   = 6,
+
+	None        = 7,
+	Pure_None   = 8,
 }
 
 Type_Info_Enum_Value :: distinct i64;

core/runtime/internal.odin

@@ -2,31 +2,31 @@ package runtime
 
 import "core:os"
 
-bswap_16 :: proc "none" (x: u16) -> u16 {
+bswap_16 :: proc "pure" (x: u16) -> u16 {
 	return x>>8 | x<<8;
 }
 
-bswap_32 :: proc "none" (x: u32) -> u32 {
+bswap_32 :: proc "pure" (x: u32) -> u32 {
 	return x>>24 | (x>>8)&0xff00 | (x<<8)&0xff0000 | x<<24;
 }
 
-bswap_64 :: proc "none" (x: u64) -> u64 {
+bswap_64 :: proc "pure" (x: u64) -> u64 {
 	return u64(bswap_32(u32(x))) | u64(bswap_32(u32(x>>32)));
 }
 
-bswap_128 :: proc "none" (x: u128) -> u128 {
+bswap_128 :: proc "pure" (x: u128) -> u128 {
 	return u128(bswap_64(u64(x))) | u128(bswap_64(u64(x>>64)));
 }
 
 
-bswap_f32 :: proc "none" (f: f32) -> f32 {
+bswap_f32 :: proc "pure" (f: f32) -> f32 {
 	x := transmute(u32)f;
 	z := x>>24 | (x>>8)&0xff00 | (x<<8)&0xff0000 | x<<24;
 	return transmute(f32)z;
 
 }
 
-bswap_f64 :: proc "none" (f: f64) -> f64 {
+bswap_f64 :: proc "pure" (f: f64) -> f64 {
 	x := transmute(u64)f;
 	z := u64(bswap_32(u32(x))) | u64(bswap_32(u32(x>>32)));
 	return transmute(f64)z;

src/check_decl.cpp

@@ -1197,6 +1197,15 @@ void check_proc_body(CheckerContext *ctx_, Token token, DeclInfo *decl, Type *ty
 	ctx->curr_proc_sig  = type;
 	ctx->curr_proc_calling_convention = type->Proc.calling_convention;
 
+	switch (type->Proc.calling_convention) {
+	case ProcCC_None:
+		error(body, "Procedures with the calling convention \"none\" are not allowed a body");
+		break;
+	case ProcCC_PureNone:
+		error(body, "Procedures with the calling convention \"pure_none\" are not allowed a body");
+		break;
+	}
+
 	ast_node(bs, BlockStmt, body);
 
 	Array<ProcUsingVar> using_entities = {};

src/check_expr.cpp

@@ -7655,7 +7655,7 @@ ExprKind check_call_expr(CheckerContext *c, Operand *operand, Ast *call, Ast *pr
 
 	{
 		if (c->curr_proc_calling_convention == ProcCC_Pure) {
-			if (pt->kind == Type_Proc && pt->Proc.calling_convention != ProcCC_Pure) {
+			if (pt->kind == Type_Proc && pt->Proc.calling_convention != ProcCC_Pure && pt->Proc.calling_convention != ProcCC_PureNone) {
 				error(call, "Only \"pure\" procedure calls are allowed within a \"pure\" procedure");
 			}
 		}

src/check_type.cpp

@@ -2200,7 +2200,7 @@ Type *type_to_abi_compat_param_type(gbAllocator a, Type *original_type, ProcCall
 		return new_type;
 	}
 
-	if (cc == ProcCC_None) {
+	if (cc == ProcCC_None || cc == ProcCC_PureNone) {
 		return new_type;
 	}
 
@@ -2335,7 +2335,7 @@ Type *type_to_abi_compat_result_type(gbAllocator a, Type *original_type, ProcCal
 	if (build_context.ODIN_OS == "windows") {
 		if (build_context.ODIN_ARCH == "amd64") {
 			if (is_type_integer_128bit(single_type)) {
-				if (cc == ProcCC_None) {
+				if (cc == ProcCC_None || cc == ProcCC_PureNone) {
 					return original_type;
 				} else {
 					return alloc_type_simd_vector(2, t_u64);
@@ -2401,7 +2401,7 @@ bool abi_compat_return_by_pointer(gbAllocator a, ProcCallingConvention cc, Type
 	if (abi_return_type == nullptr) {
 		return false;
 	}
-	if (cc == ProcCC_None) {
+	if (cc == ProcCC_None || cc == ProcCC_PureNone) {
 		return false;
 	}
 

src/ir_print.cpp

@@ -1440,6 +1440,7 @@ void ir_print_calling_convention(irFileBuffer *f, irModule *m, ProcCallingConven
 	case ProcCC_StdCall:     ir_write_str_lit(f, "cc 64 "); break;
 	case ProcCC_FastCall:    ir_write_str_lit(f, "cc 65 "); break;
 	case ProcCC_None:        ir_write_str_lit(f, "");       break;
+	case ProcCC_PureNone:    ir_write_str_lit(f, "");       break;
 	default: GB_PANIC("unknown calling convention: %d", cc);
 	}
 }

src/parser.cpp

@@ -2991,6 +2991,7 @@ ProcCallingConvention string_to_calling_convention(String s) {
 	if (s == "odin")        return ProcCC_Odin;
 	if (s == "contextless") return ProcCC_Contextless;
 	if (s == "pure")        return ProcCC_Pure;
+	if (s == "pure_none")   return ProcCC_PureNone;
 	if (s == "cdecl")       return ProcCC_CDecl;
 	if (s == "c")           return ProcCC_CDecl;
 	if (s == "stdcall")     return ProcCC_StdCall;

src/parser.hpp

@@ -178,14 +178,15 @@ enum ProcTag {
 
 enum ProcCallingConvention {
 	ProcCC_Invalid = 0,
-	ProcCC_Odin,
-	ProcCC_Contextless,
-	ProcCC_Pure,
-	ProcCC_CDecl,
-	ProcCC_StdCall,
-	ProcCC_FastCall,
-
-	ProcCC_None,
+	ProcCC_Odin = 1,
+	ProcCC_Contextless = 2,
+	ProcCC_Pure = 3,
+	ProcCC_CDecl = 4,
+	ProcCC_StdCall = 5,
+	ProcCC_FastCall = 6,
+
+	ProcCC_None = 7,
+	ProcCC_PureNone = 8,
 
 	ProcCC_MAX,
 

src/types.cpp

@@ -3556,6 +3556,9 @@ gbString write_type_to_string(gbString str, Type *type) {
 		case ProcCC_FastCall:
 			str = gb_string_appendc(str, " \"fastcall\" ");
 			break;
+		case ProcCC_PureNone:
+			str = gb_string_appendc(str, " \"pure_none\" ");
+			break;
 		case ProcCC_None:
 			str = gb_string_appendc(str, " \"none\" ");
 			break;