diff options
Diffstat (limited to 'core/math/linalg/general.odin')
| -rw-r--r-- | core/math/linalg/general.odin | 279 |
1 files changed, 274 insertions, 5 deletions
diff --git a/core/math/linalg/general.odin b/core/math/linalg/general.odin index 60185d64d..24bc4c7b3 100644 --- a/core/math/linalg/general.odin +++ b/core/math/linalg/general.odin @@ -1,8 +1,8 @@ package linalg import "core:math" -import "core:builtin" -import "core:intrinsics" +import "base:builtin" +import "base:intrinsics" // Generic @@ -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 +} |