diff options
| author | gingerBill <bill@gingerbill.org> | 2024-01-28 21:28:54 +0000 |
|---|---|---|
| committer | gingerBill <bill@gingerbill.org> | 2024-01-28 21:28:54 +0000 |
| commit | f67691c457f9a369013aac61bcd99890ce252c51 (patch) | |
| tree | efa3d14e4d0640965e2d75fe332123b870af2540 | |
| parent | d04c82e5471bb291cddbef883a36b1caad3b2b99 (diff) | |
Move matrix related procedures to the different `linalg` packages
| -rw-r--r-- | base/runtime/core_builtin_matrix.odin | 283 | ||||
| -rw-r--r-- | core/math/linalg/general.odin | 275 | ||||
| -rw-r--r-- | core/math/linalg/glsl/linalg_glsl.odin | 298 | ||||
| -rw-r--r-- | core/math/linalg/hlsl/linalg_hlsl.odin | 291 | ||||
| -rw-r--r-- | core/math/linalg/specific.odin | 12 |
5 files changed, 829 insertions, 330 deletions
diff --git a/base/runtime/core_builtin_matrix.odin b/base/runtime/core_builtin_matrix.odin deleted file mode 100644 index ed1b5b1e6..000000000 --- a/base/runtime/core_builtin_matrix.odin +++ /dev/null @@ -1,283 +0,0 @@ -package runtime - -import "core:intrinsics" -_ :: intrinsics - -@(builtin) -transpose :: intrinsics.transpose -@(builtin) -outer_product :: intrinsics.outer_product -@(builtin) -hadamard_product :: intrinsics.hadamard_product -@(builtin) -matrix_flatten :: intrinsics.matrix_flatten - - -@(builtin) -determinant :: proc{ - matrix1x1_determinant, - matrix2x2_determinant, - matrix3x3_determinant, - matrix4x4_determinant, -} - -@(builtin) -adjugate :: proc{ - matrix1x1_adjugate, - matrix2x2_adjugate, - matrix3x3_adjugate, - matrix4x4_adjugate, -} - -@(builtin) -inverse_transpose :: proc{ - matrix1x1_inverse_transpose, - matrix2x2_inverse_transpose, - matrix3x3_inverse_transpose, - matrix4x4_inverse_transpose, -} - - -@(builtin) -inverse :: proc{ - matrix1x1_inverse, - matrix2x2_inverse, - matrix3x3_inverse, - matrix4x4_inverse, -} - -@(builtin, require_results) -hermitian_adjoint :: proc "contextless" (m: $M/matrix[$N, N]$T) -> M where intrinsics.type_is_complex(T), N >= 1 { - return conj(transpose(m)) -} - -@(builtin, require_results) -matrix_trace :: proc "contextless" (m: $M/matrix[$N, N]$T) -> (trace: T) { - for i in 0..<N { - trace += m[i, i] - } - return -} - -@(builtin, require_results) -matrix_minor :: proc "contextless" (m: $M/matrix[$N, N]$T, row, column: int) -> (minor: T) where N > 1 { - K :: N-1 - cut_down: matrix[K, K]T - for col_idx in 0..<K { - j := col_idx + int(col_idx >= column) - for row_idx in 0..<K { - i := row_idx + int(row_idx >= row) - cut_down[row_idx, col_idx] = m[i, j] - } - } - return determinant(cut_down) -} - - - -@(builtin, require_results) -matrix1x1_determinant :: proc "contextless" (m: $M/matrix[1, 1]$T) -> (det: T) { - return m[0, 0] -} - -@(builtin, require_results) -matrix2x2_determinant :: proc "contextless" (m: $M/matrix[2, 2]$T) -> (det: T) { - return m[0, 0]*m[1, 1] - m[0, 1]*m[1, 0] -} -@(builtin, require_results) -matrix3x3_determinant :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (det: T) { - a := +m[0, 0] * (m[1, 1] * m[2, 2] - m[1, 2] * m[2, 1]) - b := -m[0, 1] * (m[1, 0] * m[2, 2] - m[1, 2] * m[2, 0]) - c := +m[0, 2] * (m[1, 0] * m[2, 1] - m[1, 1] * m[2, 0]) - return a + b + c -} -@(builtin, require_results) -matrix4x4_determinant :: proc "contextless" (m: $M/matrix[4, 4]$T) -> (det: T) { - a := adjugate(m) - #no_bounds_check for i in 0..<4 { - det += m[0, i] * a[0, i] - } - return -} - - - - -@(builtin, require_results) -matrix1x1_adjugate :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { - y = x - return -} - -@(builtin, require_results) -matrix2x2_adjugate :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { - y[0, 0] = +x[1, 1] - y[0, 1] = -x[1, 0] - y[1, 0] = -x[0, 1] - y[1, 1] = +x[0, 0] - return -} - -@(builtin, require_results) -matrix3x3_adjugate :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) { - y[0, 0] = +(m[1, 1] * m[2, 2] - m[2, 1] * m[1, 2]) - y[0, 1] = -(m[1, 0] * m[2, 2] - m[2, 0] * m[1, 2]) - y[0, 2] = +(m[1, 0] * m[2, 1] - m[2, 0] * m[1, 1]) - y[1, 0] = -(m[0, 1] * m[2, 2] - m[2, 1] * m[0, 2]) - y[1, 1] = +(m[0, 0] * m[2, 2] - m[2, 0] * m[0, 2]) - y[1, 2] = -(m[0, 0] * m[2, 1] - m[2, 0] * m[0, 1]) - y[2, 0] = +(m[0, 1] * m[1, 2] - m[1, 1] * m[0, 2]) - y[2, 1] = -(m[0, 0] * m[1, 2] - m[1, 0] * m[0, 2]) - y[2, 2] = +(m[0, 0] * m[1, 1] - m[1, 0] * m[0, 1]) - return -} - - -@(builtin, require_results) -matrix4x4_adjugate :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) { - for i in 0..<4 { - for j in 0..<4 { - sign: T = 1 if (i + j) % 2 == 0 else -1 - y[i, j] = sign * matrix_minor(x, i, j) - } - } - return -} - -@(builtin, require_results) -matrix1x1_inverse_transpose :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { - y[0, 0] = 1/x[0, 0] - return -} - -@(builtin, require_results) -matrix2x2_inverse_transpose :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { - d := x[0, 0]*x[1, 1] - x[0, 1]*x[1, 0] - when intrinsics.type_is_integer(T) { - y[0, 0] = +x[1, 1] / d - y[1, 0] = -x[0, 1] / d - y[0, 1] = -x[1, 0] / d - y[1, 1] = +x[0, 0] / d - } else { - id := 1 / d - y[0, 0] = +x[1, 1] * id - y[1, 0] = -x[0, 1] * id - y[0, 1] = -x[1, 0] * id - y[1, 1] = +x[0, 0] * id - } - return -} - -@(builtin, require_results) -matrix3x3_inverse_transpose :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) - d := determinant(x) - when intrinsics.type_is_integer(T) { - for i in 0..<3 { - for j in 0..<3 { - y[i, j] = a[i, j] / d - } - } - } else { - id := 1/d - for i in 0..<3 { - for j in 0..<3 { - y[i, j] = a[i, j] * id - } - } - } - return -} - -@(builtin, require_results) -matrix4x4_inverse_transpose :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) - d: T - for i in 0..<4 { - d += x[0, i] * a[0, i] - } - when intrinsics.type_is_integer(T) { - for i in 0..<4 { - for j in 0..<4 { - y[i, j] = a[i, j] / d - } - } - } else { - id := 1/d - for i in 0..<4 { - for j in 0..<4 { - y[i, j] = a[i, j] * id - } - } - } - return -} - -@(builtin, require_results) -matrix1x1_inverse :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { - y[0, 0] = 1/x[0, 0] - return -} - -@(builtin, require_results) -matrix2x2_inverse :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { - d := x[0, 0]*x[1, 1] - x[0, 1]*x[1, 0] - when intrinsics.type_is_integer(T) { - y[0, 0] = +x[1, 1] / d - y[0, 1] = -x[0, 1] / d - y[1, 0] = -x[1, 0] / d - y[1, 1] = +x[0, 0] / d - } else { - id := 1 / d - y[0, 0] = +x[1, 1] * id - y[0, 1] = -x[0, 1] * id - y[1, 0] = -x[1, 0] * id - y[1, 1] = +x[0, 0] * id - } - return -} - -@(builtin, require_results) -matrix3x3_inverse :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) - d := determinant(x) - when intrinsics.type_is_integer(T) { - for i in 0..<3 { - for j in 0..<3 { - y[i, j] = a[j, i] / d - } - } - } else { - id := 1/d - for i in 0..<3 { - for j in 0..<3 { - y[i, j] = a[j, i] * id - } - } - } - return -} - -@(builtin, require_results) -matrix4x4_inverse :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) - d: T - for i in 0..<4 { - d += x[0, i] * a[0, i] - } - when intrinsics.type_is_integer(T) { - for i in 0..<4 { - for j in 0..<4 { - y[i, j] = a[j, i] / d - } - } - } else { - id := 1/d - for i in 0..<4 { - for j in 0..<4 { - y[i, j] = a[j, i] * id - } - } - } - return -} diff --git a/core/math/linalg/general.odin b/core/math/linalg/general.odin index 60185d64d..8743cb2f0 100644 --- a/core/math/linalg/general.odin +++ b/core/math/linalg/general.odin @@ -66,7 +66,7 @@ quaternion256_dot :: proc "contextless" (a, b: $T/quaternion256) -> (c: f64) { dot :: proc{scalar_dot, vector_dot, quaternion64_dot, quaternion128_dot, quaternion256_dot} inner_product :: dot -outer_product :: builtin.outer_product +outer_product :: intrinsics.outer_product @(require_results) quaternion_inverse :: proc "contextless" (q: $Q) -> Q where IS_QUATERNION(Q) { @@ -179,8 +179,7 @@ identity :: proc "contextless" ($T: typeid/[$N][N]$E) -> (m: T) #no_bounds_check return m } -trace :: builtin.matrix_trace -transpose :: builtin.transpose +transpose :: intrinsics.transpose @(require_results) matrix_mul :: proc "contextless" (a, b: $M/matrix[$N, N]$E) -> (c: M) @@ -355,3 +354,273 @@ matrix_cast :: proc "contextless" (v: $A/matrix[$M, $N]$T, $Elem_Type: typeid) - @(require_results) to_quaternion64 :: #force_inline proc(v: $A/[$N]$T) -> [N]quaternion64 { return array_cast(v, quaternion64) } @(require_results) to_quaternion128 :: #force_inline proc(v: $A/[$N]$T) -> [N]quaternion128 { return array_cast(v, quaternion128) } @(require_results) to_quaternion256 :: #force_inline proc(v: $A/[$N]$T) -> [N]quaternion256 { return array_cast(v, quaternion256) } + + +hadamard_product :: intrinsics.hadamard_product +matrix_flatten :: intrinsics.matrix_flatten + + +determinant :: proc{ + matrix1x1_determinant, + matrix2x2_determinant, + matrix3x3_determinant, + matrix4x4_determinant, +} + +adjugate :: proc{ + matrix1x1_adjugate, + matrix2x2_adjugate, + matrix3x3_adjugate, + matrix4x4_adjugate, +} + +inverse_transpose :: proc{ + matrix1x1_inverse_transpose, + matrix2x2_inverse_transpose, + matrix3x3_inverse_transpose, + matrix4x4_inverse_transpose, +} + + +inverse :: proc{ + matrix1x1_inverse, + matrix2x2_inverse, + matrix3x3_inverse, + matrix4x4_inverse, +} + +@(require_results) +hermitian_adjoint :: proc "contextless" (m: $M/matrix[$N, N]$T) -> M where intrinsics.type_is_complex(T), N >= 1 { + return conj(transpose(m)) +} + +@(require_results) +trace :: proc "contextless" (m: $M/matrix[$N, N]$T) -> (trace: T) { + for i in 0..<N { + trace += m[i, i] + } + return +} + +@(require_results) +matrix_minor :: proc "contextless" (m: $M/matrix[$N, N]$T, #any_int row, column: int) -> (minor: T) where N > 1 { + K :: int(N-1) + cut_down: matrix[K, K]T + for col_idx in 0..<K { + j := col_idx + int(col_idx >= column) + for row_idx in 0..<K { + i := row_idx + int(row_idx >= row) + cut_down[row_idx, col_idx] = m[i, j] + } + } + return determinant(cut_down) +} + + + +@(require_results) +matrix1x1_determinant :: proc "contextless" (m: $M/matrix[1, 1]$T) -> (det: T) { + return m[0, 0] +} + +@(require_results) +matrix2x2_determinant :: proc "contextless" (m: $M/matrix[2, 2]$T) -> (det: T) { + return m[0, 0]*m[1, 1] - m[0, 1]*m[1, 0] +} +@(require_results) +matrix3x3_determinant :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (det: T) { + a := +m[0, 0] * (m[1, 1] * m[2, 2] - m[1, 2] * m[2, 1]) + b := -m[0, 1] * (m[1, 0] * m[2, 2] - m[1, 2] * m[2, 0]) + c := +m[0, 2] * (m[1, 0] * m[2, 1] - m[1, 1] * m[2, 0]) + return a + b + c +} +@(require_results) +matrix4x4_determinant :: proc "contextless" (m: $M/matrix[4, 4]$T) -> (det: T) { + a := adjugate(m) + #no_bounds_check for i in 0..<4 { + det += m[0, i] * a[0, i] + } + return +} + + + + +@(require_results) +matrix1x1_adjugate :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { + y = x + return +} + +@(require_results) +matrix2x2_adjugate :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { + y[0, 0] = +x[1, 1] + y[0, 1] = -x[1, 0] + y[1, 0] = -x[0, 1] + y[1, 1] = +x[0, 0] + return +} + +@(require_results) +matrix3x3_adjugate :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) { + y[0, 0] = +(m[1, 1] * m[2, 2] - m[2, 1] * m[1, 2]) + y[0, 1] = -(m[1, 0] * m[2, 2] - m[2, 0] * m[1, 2]) + y[0, 2] = +(m[1, 0] * m[2, 1] - m[2, 0] * m[1, 1]) + y[1, 0] = -(m[0, 1] * m[2, 2] - m[2, 1] * m[0, 2]) + y[1, 1] = +(m[0, 0] * m[2, 2] - m[2, 0] * m[0, 2]) + y[1, 2] = -(m[0, 0] * m[2, 1] - m[2, 0] * m[0, 1]) + y[2, 0] = +(m[0, 1] * m[1, 2] - m[1, 1] * m[0, 2]) + y[2, 1] = -(m[0, 0] * m[1, 2] - m[1, 0] * m[0, 2]) + y[2, 2] = +(m[0, 0] * m[1, 1] - m[1, 0] * m[0, 1]) + return +} + + +@(require_results) +matrix4x4_adjugate :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) { + for i in 0..<4 { + for j in 0..<4 { + sign: T = 1 if (i + j) % 2 == 0 else -1 + y[i, j] = sign * matrix_minor(x, i, j) + } + } + return +} + +@(require_results) +matrix1x1_inverse_transpose :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { + y[0, 0] = 1/x[0, 0] + return +} + +@(require_results) +matrix2x2_inverse_transpose :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { + d := x[0, 0]*x[1, 1] - x[0, 1]*x[1, 0] + when intrinsics.type_is_integer(T) { + y[0, 0] = +x[1, 1] / d + y[1, 0] = -x[0, 1] / d + y[0, 1] = -x[1, 0] / d + y[1, 1] = +x[0, 0] / d + } else { + id := 1 / d + y[0, 0] = +x[1, 1] * id + y[1, 0] = -x[0, 1] * id + y[0, 1] = -x[1, 0] * id + y[1, 1] = +x[0, 0] * id + } + return +} + +@(require_results) +matrix3x3_inverse_transpose :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { + a := adjugate(x) + d := determinant(x) + when intrinsics.type_is_integer(T) { + for i in 0..<3 { + for j in 0..<3 { + y[i, j] = a[i, j] / d + } + } + } else { + id := 1/d + for i in 0..<3 { + for j in 0..<3 { + y[i, j] = a[i, j] * id + } + } + } + return +} + +@(require_results) +matrix4x4_inverse_transpose :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { + a := adjugate(x) + d: T + for i in 0..<4 { + d += x[0, i] * a[0, i] + } + when intrinsics.type_is_integer(T) { + for i in 0..<4 { + for j in 0..<4 { + y[i, j] = a[i, j] / d + } + } + } else { + id := 1/d + for i in 0..<4 { + for j in 0..<4 { + y[i, j] = a[i, j] * id + } + } + } + return +} + +@(require_results) +matrix1x1_inverse :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { + y[0, 0] = 1/x[0, 0] + return +} + +@(require_results) +matrix2x2_inverse :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { + d := x[0, 0]*x[1, 1] - x[0, 1]*x[1, 0] + when intrinsics.type_is_integer(T) { + y[0, 0] = +x[1, 1] / d + y[0, 1] = -x[0, 1] / d + y[1, 0] = -x[1, 0] / d + y[1, 1] = +x[0, 0] / d + } else { + id := 1 / d + y[0, 0] = +x[1, 1] * id + y[0, 1] = -x[0, 1] * id + y[1, 0] = -x[1, 0] * id + y[1, 1] = +x[0, 0] * id + } + return +} + +@(require_results) +matrix3x3_inverse :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { + a := adjugate(x) + d := determinant(x) + when intrinsics.type_is_integer(T) { + for i in 0..<3 { + for j in 0..<3 { + y[i, j] = a[j, i] / d + } + } + } else { + id := 1/d + for i in 0..<3 { + for j in 0..<3 { + y[i, j] = a[j, i] * id + } + } + } + return +} + +@(require_results) +matrix4x4_inverse :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { + a := adjugate(x) + d: T + for i in 0..<4 { + d += x[0, i] * a[0, i] + } + when intrinsics.type_is_integer(T) { + for i in 0..<4 { + for j in 0..<4 { + y[i, j] = a[j, i] / d + } + } + } else { + id := 1/d + for i in 0..<4 { + for j in 0..<4 { + y[i, j] = a[j, i] * id + } + } + } + return +} diff --git a/core/math/linalg/glsl/linalg_glsl.odin b/core/math/linalg/glsl/linalg_glsl.odin index 0d91ad4a3..68a05be53 100644 --- a/core/math/linalg/glsl/linalg_glsl.odin +++ b/core/math/linalg/glsl/linalg_glsl.odin @@ -2,6 +2,7 @@ package math_linalg_glsl import "core:builtin" +import "core:intrinsics" TAU :: 6.28318530717958647692528676655900576 PI :: 3.14159265358979323846264338327950288 @@ -1838,30 +1839,281 @@ dquatMulDvec3 :: proc "c" (q: dquat, v: dvec3) -> dvec3 { -@(require_results) inverse_mat2 :: proc "c" (m: mat2) -> mat2 { return builtin.inverse(m) } -@(require_results) inverse_mat3 :: proc "c" (m: mat3) -> mat3 { return builtin.inverse(m) } -@(require_results) inverse_mat4 :: proc "c" (m: mat4) -> mat4 { return builtin.inverse(m) } -@(require_results) inverse_dmat2 :: proc "c" (m: dmat2) -> dmat2 { return builtin.inverse(m) } -@(require_results) inverse_dmat3 :: proc "c" (m: dmat3) -> dmat3 { return builtin.inverse(m) } -@(require_results) inverse_dmat4 :: proc "c" (m: dmat4) -> dmat4 { return builtin.inverse(m) } +@(require_results) inverse_mat2 :: proc "c" (m: mat2) -> mat2 { return inverse_matrix2x2(m) } +@(require_results) inverse_mat3 :: proc "c" (m: mat3) -> mat3 { return inverse_matrix3x3(m) } +@(require_results) inverse_mat4 :: proc "c" (m: mat4) -> mat4 { return inverse_matrix4x4(m) } +@(require_results) inverse_dmat2 :: proc "c" (m: dmat2) -> dmat2 { return inverse_matrix2x2(m) } +@(require_results) inverse_dmat3 :: proc "c" (m: dmat3) -> dmat3 { return inverse_matrix3x3(m) } +@(require_results) inverse_dmat4 :: proc "c" (m: dmat4) -> dmat4 { return inverse_matrix4x4(m) } @(require_results) inverse_quat :: proc "c" (q: quat) -> quat { return 1/q } @(require_results) inverse_dquat :: proc "c" (q: dquat) -> dquat { return 1/q } + +transpose :: intrinsics.transpose + + +determinant :: proc{ + determinant_matrix1x1, + determinant_matrix2x2, + determinant_matrix3x3, + determinant_matrix4x4, +} + +adjugate :: proc{ + adjugate_matrix1x1, + adjugate_matrix2x2, + adjugate_matrix3x3, + adjugate_matrix4x4, +} + +inverse_transpose :: proc{ + inverse_transpose_matrix1x1, + inverse_transpose_matrix2x2, + inverse_transpose_matrix3x3, + inverse_transpose_matrix4x4, +} + + inverse :: proc{ - inverse_mat2, - inverse_mat3, - inverse_mat4, - inverse_dmat2, - inverse_dmat3, - inverse_dmat4, - inverse_quat, - inverse_dquat, -} - -transpose :: builtin.transpose -inverse_transpose :: builtin.inverse_transpose -adjugate :: builtin.adjugate -hermitian_adjoint :: builtin.hermitian_adjoint -minor :: builtin.matrix_minor -determinant :: builtin.determinant -trace :: builtin.matrix_trace
\ No newline at end of file + inverse_matrix1x1, + inverse_matrix2x2, + inverse_matrix3x3, + inverse_matrix4x4, +} + +@(require_results) +hermitian_adjoint :: proc "contextless" (m: $M/matrix[$N, N]$T) -> M where intrinsics.type_is_complex(T), N >= 1 { + return conj(transpose(m)) +} + +@(require_results) +trace :: proc "contextless" (m: $M/matrix[$N, N]$T) -> (trace: T) { + for i in 0..<N { + trace += m[i, i] + } + return +} + +@(require_results) +matrix_minor :: proc "contextless" (m: $M/matrix[$N, N]$T, #any_int row, column: int) -> (minor: T) where N > 1 { + K :: int(N-1) + cut_down: matrix[K, K]T + for col_idx in 0..<K { + j := col_idx + int(col_idx >= column) + for row_idx in 0..<K { + i := row_idx + int(row_idx >= row) + cut_down[row_idx, col_idx] = m[i, j] + } + } + return determinant(cut_down) +} + + + +@(require_results) +determinant_matrix1x1 :: proc "contextless" (m: $M/matrix[1, 1]$T) -> (det: T) { + return m[0, 0] +} + +@(require_results) +determinant_matrix2x2 :: proc "contextless" (m: $M/matrix[2, 2]$T) -> (det: T) { + return m[0, 0]*m[1, 1] - m[0, 1]*m[1, 0] +} +@(require_results) +determinant_matrix3x3 :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (det: T) { + a := +m[0, 0] * (m[1, 1] * m[2, 2] - m[1, 2] * m[2, 1]) + b := -m[0, 1] * (m[1, 0] * m[2, 2] - m[1, 2] * m[2, 0]) + c := +m[0, 2] * (m[1, 0] * m[2, 1] - m[1, 1] * m[2, 0]) + return a + b + c +} +@(require_results) +determinant_matrix4x4 :: proc "contextless" (m: $M/matrix[4, 4]$T) -> (det: T) { + a := adjugate(m) + #no_bounds_check for i in 0..<4 { + det += m[0, i] * a[0, i] + } + return +} + + + + +@(require_results) +adjugate_matrix1x1 :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { + y = x + return +} + +@(require_results) +adjugate_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { + y[0, 0] = +x[1, 1] + y[0, 1] = -x[1, 0] + y[1, 0] = -x[0, 1] + y[1, 1] = +x[0, 0] + return +} + +@(require_results) +adjugate_matrix3x3 :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) { + y[0, 0] = +(m[1, 1] * m[2, 2] - m[2, 1] * m[1, 2]) + y[0, 1] = -(m[1, 0] * m[2, 2] - m[2, 0] * m[1, 2]) + y[0, 2] = +(m[1, 0] * m[2, 1] - m[2, 0] * m[1, 1]) + y[1, 0] = -(m[0, 1] * m[2, 2] - m[2, 1] * m[0, 2]) + y[1, 1] = +(m[0, 0] * m[2, 2] - m[2, 0] * m[0, 2]) + y[1, 2] = -(m[0, 0] * m[2, 1] - m[2, 0] * m[0, 1]) + y[2, 0] = +(m[0, 1] * m[1, 2] - m[1, 1] * m[0, 2]) + y[2, 1] = -(m[0, 0] * m[1, 2] - m[1, 0] * m[0, 2]) + y[2, 2] = +(m[0, 0] * m[1, 1] - m[1, 0] * m[0, 1]) + return +} + + +@(require_results) +adjugate_matrix4x4 :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) { + for i in 0..<4 { + for j in 0..<4 { + sign: T = 1 if (i + j) % 2 == 0 else -1 + y[i, j] = sign * matrix_minor(x, i, j) + } + } + return +} + +@(require_results) +inverse_transpose_matrix1x1 :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { + y[0, 0] = 1/x[0, 0] + return +} + +@(require_results) +inverse_transpose_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { + d := x[0, 0]*x[1, 1] - x[0, 1]*x[1, 0] + when intrinsics.type_is_integer(T) { + y[0, 0] = +x[1, 1] / d + y[1, 0] = -x[0, 1] / d + y[0, 1] = -x[1, 0] / d + y[1, 1] = +x[0, 0] / d + } else { + id := 1 / d + y[0, 0] = +x[1, 1] * id + y[1, 0] = -x[0, 1] * id + y[0, 1] = -x[1, 0] * id + y[1, 1] = +x[0, 0] * id + } + return +} + +@(require_results) +inverse_transpose_matrix3x3 :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { + a := adjugate(x) + d := determinant(x) + when intrinsics.type_is_integer(T) { + for i in 0..<3 { + for j in 0..<3 { + y[i, j] = a[i, j] / d + } + } + } else { + id := 1/d + for i in 0..<3 { + for j in 0..<3 { + y[i, j] = a[i, j] * id + } + } + } + return +} + +@(require_results) +inverse_transpose_matrix4x4 :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { + a := adjugate(x) + d: T + for i in 0..<4 { + d += x[0, i] * a[0, i] + } + when intrinsics.type_is_integer(T) { + for i in 0..<4 { + for j in 0..<4 { + y[i, j] = a[i, j] / d + } + } + } else { + id := 1/d + for i in 0..<4 { + for j in 0..<4 { + y[i, j] = a[i, j] * id + } + } + } + return +} + +@(require_results) +inverse_matrix1x1 :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { + y[0, 0] = 1/x[0, 0] + return +} + +@(require_results) +inverse_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { + d := x[0, 0]*x[1, 1] - x[0, 1]*x[1, 0] + when intrinsics.type_is_integer(T) { + y[0, 0] = +x[1, 1] / d + y[0, 1] = -x[0, 1] / d + y[1, 0] = -x[1, 0] / d + y[1, 1] = +x[0, 0] / d + } else { + id := 1 / d + y[0, 0] = +x[1, 1] * id + y[0, 1] = -x[0, 1] * id + y[1, 0] = -x[1, 0] * id + y[1, 1] = +x[0, 0] * id + } + return +} + +@(require_results) +inverse_matrix3x3 :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { + a := adjugate(x) + d := determinant(x) + when intrinsics.type_is_integer(T) { + for i in 0..<3 { + for j in 0..<3 { + y[i, j] = a[j, i] / d + } + } + } else { + id := 1/d + for i in 0..<3 { + for j in 0..<3 { + y[i, j] = a[j, i] * id + } + } + } + return +} + +@(require_results) +inverse_matrix4x4 :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { + a := adjugate(x) + d: T + for i in 0..<4 { + d += x[0, i] * a[0, i] + } + when intrinsics.type_is_integer(T) { + for i in 0..<4 { + for j in 0..<4 { + y[i, j] = a[j, i] / d + } + } + } else { + id := 1/d + for i in 0..<4 { + for j in 0..<4 { + y[i, j] = a[j, i] * id + } + } + } + return +} + diff --git a/core/math/linalg/hlsl/linalg_hlsl.odin b/core/math/linalg/hlsl/linalg_hlsl.odin index 351aa7ea3..66278345e 100644 --- a/core/math/linalg/hlsl/linalg_hlsl.odin +++ b/core/math/linalg/hlsl/linalg_hlsl.odin @@ -2,6 +2,7 @@ package math_linalg_hlsl import "core:builtin" +import "core:intrinsics" TAU :: 6.28318530717958647692528676655900576 PI :: 3.14159265358979323846264338327950288 @@ -1471,14 +1472,14 @@ not :: proc{ -@(require_results) inverse_float1x1 :: proc "c" (m: float1x1) -> float1x1 { return builtin.inverse(m) } -@(require_results) inverse_float2x2 :: proc "c" (m: float2x2) -> float2x2 { return builtin.inverse(m) } -@(require_results) inverse_float3x3 :: proc "c" (m: float3x3) -> float3x3 { return builtin.inverse(m) } -@(require_results) inverse_float4x4 :: proc "c" (m: float4x4) -> float4x4 { return builtin.inverse(m) } -@(require_results) inverse_double1x1 :: proc "c" (m: double1x1) -> double1x1 { return builtin.inverse(m) } -@(require_results) inverse_double2x2 :: proc "c" (m: double2x2) -> double2x2 { return builtin.inverse(m) } -@(require_results) inverse_double3x3 :: proc "c" (m: double3x3) -> double3x3 { return builtin.inverse(m) } -@(require_results) inverse_double4x4 :: proc "c" (m: double4x4) -> double4x4 { return builtin.inverse(m) } +@(require_results) inverse_float1x1 :: proc "c" (m: float1x1) -> float1x1 { return inverse_matrix1x1(m) } +@(require_results) inverse_float2x2 :: proc "c" (m: float2x2) -> float2x2 { return inverse_matrix2x2(m) } +@(require_results) inverse_float3x3 :: proc "c" (m: float3x3) -> float3x3 { return inverse_matrix3x3(m) } +@(require_results) inverse_float4x4 :: proc "c" (m: float4x4) -> float4x4 { return inverse_matrix4x4(m) } +@(require_results) inverse_double1x1 :: proc "c" (m: double1x1) -> double1x1 { return inverse_matrix1x1(m) } +@(require_results) inverse_double2x2 :: proc "c" (m: double2x2) -> double2x2 { return inverse_matrix2x2(m) } +@(require_results) inverse_double3x3 :: proc "c" (m: double3x3) -> double3x3 { return inverse_matrix3x3(m) } +@(require_results) inverse_double4x4 :: proc "c" (m: double4x4) -> double4x4 { return inverse_matrix4x4(m) } inverse :: proc{ inverse_float1x1, @@ -1489,15 +1490,275 @@ inverse :: proc{ inverse_double2x2, inverse_double3x3, inverse_double4x4, + + inverse_matrix1x1, + inverse_matrix2x2, + inverse_matrix3x3, + inverse_matrix4x4, +} + +transpose :: intrinsics.transpose + + +determinant :: proc{ + determinant_matrix1x1, + determinant_matrix2x2, + determinant_matrix3x3, + determinant_matrix4x4, +} + +adjugate :: proc{ + adjugate_matrix1x1, + adjugate_matrix2x2, + adjugate_matrix3x3, + adjugate_matrix4x4, +} + +inverse_transpose :: proc{ + inverse_transpose_matrix1x1, + inverse_transpose_matrix2x2, + inverse_transpose_matrix3x3, + inverse_transpose_matrix4x4, +} + +@(require_results) +hermitian_adjoint :: proc "contextless" (m: $M/matrix[$N, N]$T) -> M where intrinsics.type_is_complex(T), N >= 1 { + return conj(transpose(m)) +} + +@(require_results) +trace :: proc "contextless" (m: $M/matrix[$N, N]$T) -> (trace: T) { + for i in 0..<N { + trace += m[i, i] + } + return +} + +@(require_results) +matrix_minor :: proc "contextless" (m: $M/matrix[$N, N]$T, #any_int row, column: int) -> (minor: T) where N > 1 { + K :: int(N-1) + cut_down: matrix[K, K]T + for col_idx in 0..<K { + j := col_idx + int(col_idx >= column) + for row_idx in 0..<K { + i := row_idx + int(row_idx >= row) + cut_down[row_idx, col_idx] = m[i, j] + } + } + return determinant(cut_down) +} + + + +@(require_results) +determinant_matrix1x1 :: proc "contextless" (m: $M/matrix[1, 1]$T) -> (det: T) { + return m[0, 0] +} + +@(require_results) +determinant_matrix2x2 :: proc "contextless" (m: $M/matrix[2, 2]$T) -> (det: T) { + return m[0, 0]*m[1, 1] - m[0, 1]*m[1, 0] +} +@(require_results) +determinant_matrix3x3 :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (det: T) { + a := +m[0, 0] * (m[1, 1] * m[2, 2] - m[1, 2] * m[2, 1]) + b := -m[0, 1] * (m[1, 0] * m[2, 2] - m[1, 2] * m[2, 0]) + c := +m[0, 2] * (m[1, 0] * m[2, 1] - m[1, 1] * m[2, 0]) + return a + b + c +} +@(require_results) +determinant_matrix4x4 :: proc "contextless" (m: $M/matrix[4, 4]$T) -> (det: T) { + a := adjugate(m) + #no_bounds_check for i in 0..<4 { + det += m[0, i] * a[0, i] + } + return +} + + + + +@(require_results) +adjugate_matrix1x1 :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { + y = x + return +} + +@(require_results) +adjugate_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { + y[0, 0] = +x[1, 1] + y[0, 1] = -x[1, 0] + y[1, 0] = -x[0, 1] + y[1, 1] = +x[0, 0] + return } -transpose :: builtin.transpose -inverse_transpose :: builtin.inverse_transpose -adjugate :: builtin.adjugate -hermitian_adjoint :: builtin.hermitian_adjoint -minor :: builtin.matrix_minor -determinant :: builtin.determinant -trace :: builtin.matrix_trace +@(require_results) +adjugate_matrix3x3 :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) { + y[0, 0] = +(m[1, 1] * m[2, 2] - m[2, 1] * m[1, 2]) + y[0, 1] = -(m[1, 0] * m[2, 2] - m[2, 0] * m[1, 2]) + y[0, 2] = +(m[1, 0] * m[2, 1] - m[2, 0] * m[1, 1]) + y[1, 0] = -(m[0, 1] * m[2, 2] - m[2, 1] * m[0, 2]) + y[1, 1] = +(m[0, 0] * m[2, 2] - m[2, 0] * m[0, 2]) + y[1, 2] = -(m[0, 0] * m[2, 1] - m[2, 0] * m[0, 1]) + y[2, 0] = +(m[0, 1] * m[1, 2] - m[1, 1] * m[0, 2]) + y[2, 1] = -(m[0, 0] * m[1, 2] - m[1, 0] * m[0, 2]) + y[2, 2] = +(m[0, 0] * m[1, 1] - m[1, 0] * m[0, 1]) + return +} + + +@(require_results) +adjugate_matrix4x4 :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) { + for i in 0..<4 { + for j in 0..<4 { + sign: T = 1 if (i + j) % 2 == 0 else -1 + y[i, j] = sign * matrix_minor(x, i, j) + } + } + return +} + +@(require_results) +inverse_transpose_matrix1x1 :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { + y[0, 0] = 1/x[0, 0] + return +} + +@(require_results) +inverse_transpose_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { + d := x[0, 0]*x[1, 1] - x[0, 1]*x[1, 0] + when intrinsics.type_is_integer(T) { + y[0, 0] = +x[1, 1] / d + y[1, 0] = -x[0, 1] / d + y[0, 1] = -x[1, 0] / d + y[1, 1] = +x[0, 0] / d + } else { + id := 1 / d + y[0, 0] = +x[1, 1] * id + y[1, 0] = -x[0, 1] * id + y[0, 1] = -x[1, 0] * id + y[1, 1] = +x[0, 0] * id + } + return +} + +@(require_results) +inverse_transpose_matrix3x3 :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { + a := adjugate(x) + d := determinant(x) + when intrinsics.type_is_integer(T) { + for i in 0..<3 { + for j in 0..<3 { + y[i, j] = a[i, j] / d + } + } + } else { + id := 1/d + for i in 0..<3 { + for j in 0..<3 { + y[i, j] = a[i, j] * id + } + } + } + return +} + +@(require_results) +inverse_transpose_matrix4x4 :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { + a := adjugate(x) + d: T + for i in 0..<4 { + d += x[0, i] * a[0, i] + } + when intrinsics.type_is_integer(T) { + for i in 0..<4 { + for j in 0..<4 { + y[i, j] = a[i, j] / d + } + } + } else { + id := 1/d + for i in 0..<4 { + for j in 0..<4 { + y[i, j] = a[i, j] * id + } + } + } + return +} + +@(require_results) +inverse_matrix1x1 :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { + y[0, 0] = 1/x[0, 0] + return +} + +@(require_results) +inverse_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { + d := x[0, 0]*x[1, 1] - x[0, 1]*x[1, 0] + when intrinsics.type_is_integer(T) { + y[0, 0] = +x[1, 1] / d + y[0, 1] = -x[0, 1] / d + y[1, 0] = -x[1, 0] / d + y[1, 1] = +x[0, 0] / d + } else { + id := 1 / d + y[0, 0] = +x[1, 1] * id + y[0, 1] = -x[0, 1] * id + y[1, 0] = -x[1, 0] * id + y[1, 1] = +x[0, 0] * id + } + return +} + +@(require_results) +inverse_matrix3x3 :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { + a := adjugate(x) + d := determinant(x) + when intrinsics.type_is_integer(T) { + for i in 0..<3 { + for j in 0..<3 { + y[i, j] = a[j, i] / d + } + } + } else { + id := 1/d + for i in 0..<3 { + for j in 0..<3 { + y[i, j] = a[j, i] * id + } + } + } + return +} + +@(require_results) +inverse_matrix4x4 :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { + a := adjugate(x) + d: T + for i in 0..<4 { + d += x[0, i] * a[0, i] + } + when intrinsics.type_is_integer(T) { + for i in 0..<4 { + for j in 0..<4 { + y[i, j] = a[j, i] / d + } + } + } else { + id := 1/d + for i in 0..<4 { + for j in 0..<4 { + y[i, j] = a[j, i] * id + } + } + } + return +} + + + asfloat :: proc{ asfloat_float, diff --git a/core/math/linalg/specific.odin b/core/math/linalg/specific.odin index 1f96eb178..656def1cc 100644 --- a/core/math/linalg/specific.odin +++ b/core/math/linalg/specific.odin @@ -1447,16 +1447,16 @@ matrix3_adjoint :: proc{ @(require_results) -matrix3_inverse_transpose_f16 :: proc "contextless" (m: Matrix3f16) -> (inverse_transpose: Matrix3f16) { - return builtin.inverse_transpose(m) +matrix3_inverse_transpose_f16 :: proc "contextless" (m: Matrix3f16) -> (p: Matrix3f16) { + return inverse_transpose(m) } @(require_results) -matrix3_inverse_transpose_f32 :: proc "contextless" (m: Matrix3f32) -> (inverse_transpose: Matrix3f32) { - return builtin.inverse_transpose(m) +matrix3_inverse_transpose_f32 :: proc "contextless" (m: Matrix3f32) -> (p: Matrix3f32) { + return inverse_transpose(m) } @(require_results) -matrix3_inverse_transpose_f64 :: proc "contextless" (m: Matrix3f64) -> (inverse_transpose: Matrix3f64) { - return builtin.inverse_transpose(m) +matrix3_inverse_transpose_f64 :: proc "contextless" (m: Matrix3f64) -> (p: Matrix3f64) { + return inverse_transpose(m) } matrix3_inverse_transpose :: proc{ matrix3_inverse_transpose_f16, |