From 8a91e0bb19c013d0b088c90f3d1a09560dffac5c Mon Sep 17 00:00:00 2001 From: Sebastian Pahnke Date: Sat, 28 Dec 2024 08:51:54 +0100 Subject: [PATCH 1/7] Add regression tests reproducing the issue --- tests/issues/run.bat | 1 + tests/issues/run.sh | 1 + tests/issues/test_issue_4584.odin | 56 +++++++++++++++++++++++++++++++ 3 files changed, 58 insertions(+) create mode 100644 tests/issues/test_issue_4584.odin diff --git a/tests/issues/run.bat b/tests/issues/run.bat index dcea3d483..7ed43205d 100644 --- a/tests/issues/run.bat +++ b/tests/issues/run.bat @@ -16,6 +16,7 @@ set COMMON=-define:ODIN_TEST_FANCY=false -file -vet -strict-style ..\..\..\odin test ..\test_issue_2637.odin %COMMON% || exit /b ..\..\..\odin test ..\test_issue_2666.odin %COMMON% || exit /b ..\..\..\odin test ..\test_issue_4210.odin %COMMON% || exit /b +..\..\..\odin test ..\test_issue_4584.odin %COMMON% || exit /b @echo off diff --git a/tests/issues/run.sh b/tests/issues/run.sh index c3bc00e24..54543980e 100755 --- a/tests/issues/run.sh +++ b/tests/issues/run.sh @@ -17,6 +17,7 @@ $ODIN test ../test_issue_2615.odin $COMMON $ODIN test ../test_issue_2637.odin $COMMON $ODIN test ../test_issue_2666.odin $COMMON $ODIN test ../test_issue_4210.odin $COMMON +$ODIN test ../test_issue_4584.odin $COMMON if [[ $($ODIN build ../test_issue_2395.odin $COMMON 2>&1 >/dev/null | grep -c "Error:") -eq 2 ]] ; then echo "SUCCESSFUL 1/1" else diff --git a/tests/issues/test_issue_4584.odin b/tests/issues/test_issue_4584.odin new file mode 100644 index 000000000..b80ecb3d9 --- /dev/null +++ b/tests/issues/test_issue_4584.odin @@ -0,0 +1,56 @@ +// Tests issue #4584 https://github.com/odin-lang/Odin/issues/4584 +package test_issues + +import "core:testing" +import "core:math/linalg" + +@test +test_adjugate_2x2 :: proc(t: ^testing.T) { + m := matrix[2,2]int { + -3, 2, + -1, 0, + } + expected := matrix[2,2]int { + 0, -2, + 1, -3, + } + testing.expect_value(t, linalg.adjugate(m), expected) + testing.expect_value(t, linalg.determinant(m), 2) + testing.expect_value(t, linalg.adjugate(m) * m, 2 * linalg.identity(matrix[2,2]int)) +} + +@test +test_adjugate_3x3 :: proc(t: ^testing.T) { + m := matrix[3,3]int { + -3, 2, -5, + -1, 0, -2, + 3, -4, 1, + } + expected := matrix[3,3]int { + -8, 18, -4, + -5, 12, -1, + 4, -6, 2, + } + testing.expect_value(t, linalg.adjugate(m), expected) + testing.expect_value(t, linalg.determinant(m), -6) + testing.expect_value(t, linalg.adjugate(m) * m, -6 * linalg.identity(matrix[3,3]int)) +} + +@test +test_adjugate_4x4 :: proc(t: ^testing.T) { + m := matrix[4,4]int { + -3, 2, -5, 1, + -1, 0, -2, 2, + 3, -4, 1, 3, + 4, 5, 6, 7, + } + expected := matrix[4,4]int { + -144, 266, -92, -16, + 57, 92, -5, -16, + 105, -142, 55, 2, + 33, -96, 9, -6, + } + testing.expect_value(t, linalg.adjugate(m), expected) + testing.expect_value(t, linalg.determinant(m), -174) + testing.expect_value(t, linalg.adjugate(m) * m, -174 * linalg.identity(matrix[4,4]int)) +} \ No newline at end of file From e8a202f0a27256162f8fe20f4d92cdb63ec1e967 Mon Sep 17 00:00:00 2001 From: Sebastian Pahnke Date: Sat, 28 Dec 2024 08:56:09 +0100 Subject: [PATCH 2/7] Add tests for glsl and hlsl variants --- tests/issues/test_issue_4584.odin | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/tests/issues/test_issue_4584.odin b/tests/issues/test_issue_4584.odin index b80ecb3d9..176daf6af 100644 --- a/tests/issues/test_issue_4584.odin +++ b/tests/issues/test_issue_4584.odin @@ -3,6 +3,8 @@ package test_issues import "core:testing" import "core:math/linalg" +import glm "core:math/linalg/glsl" +import hlm "core:math/linalg/hlsl" @test test_adjugate_2x2 :: proc(t: ^testing.T) { @@ -17,6 +19,12 @@ test_adjugate_2x2 :: proc(t: ^testing.T) { testing.expect_value(t, linalg.adjugate(m), expected) testing.expect_value(t, linalg.determinant(m), 2) testing.expect_value(t, linalg.adjugate(m) * m, 2 * linalg.identity(matrix[2,2]int)) + + testing.expect_value(t, glm.adjugate(m), expected) + testing.expect_value(t, glm.adjugate(m) * m, 2 * linalg.identity(matrix[2,2]int)) + + testing.expect_value(t, hlm.adjugate(m), expected) + testing.expect_value(t, hlm.adjugate(m) * m, 2 * linalg.identity(matrix[2,2]int)) } @test @@ -34,6 +42,12 @@ test_adjugate_3x3 :: proc(t: ^testing.T) { testing.expect_value(t, linalg.adjugate(m), expected) testing.expect_value(t, linalg.determinant(m), -6) testing.expect_value(t, linalg.adjugate(m) * m, -6 * linalg.identity(matrix[3,3]int)) + + testing.expect_value(t, glm.adjugate(m), expected) + testing.expect_value(t, glm.adjugate(m) * m, -6 * linalg.identity(matrix[3,3]int)) + + testing.expect_value(t, hlm.adjugate(m), expected) + testing.expect_value(t, hlm.adjugate(m) * m, -6 * linalg.identity(matrix[3,3]int)) } @test @@ -53,4 +67,10 @@ test_adjugate_4x4 :: proc(t: ^testing.T) { testing.expect_value(t, linalg.adjugate(m), expected) testing.expect_value(t, linalg.determinant(m), -174) testing.expect_value(t, linalg.adjugate(m) * m, -174 * linalg.identity(matrix[4,4]int)) + + testing.expect_value(t, glm.adjugate(m), expected) + testing.expect_value(t, glm.adjugate(m) * m, -174 * linalg.identity(matrix[4,4]int)) + + testing.expect_value(t, hlm.adjugate(m), expected) + testing.expect_value(t, hlm.adjugate(m) * m, -174 * linalg.identity(matrix[4,4]int)) } \ No newline at end of file From 0d955e55dbe9aaefebddde222d8bd62c59236eae Mon Sep 17 00:00:00 2001 From: Sebastian Pahnke Date: Sat, 28 Dec 2024 09:05:26 +0100 Subject: [PATCH 3/7] Add tests for determinants because their calculation depends on the adjugate --- tests/issues/test_issue_4584.odin | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/tests/issues/test_issue_4584.odin b/tests/issues/test_issue_4584.odin index 176daf6af..f14f1f1d3 100644 --- a/tests/issues/test_issue_4584.odin +++ b/tests/issues/test_issue_4584.odin @@ -21,9 +21,11 @@ test_adjugate_2x2 :: proc(t: ^testing.T) { testing.expect_value(t, linalg.adjugate(m) * m, 2 * linalg.identity(matrix[2,2]int)) testing.expect_value(t, glm.adjugate(m), expected) + testing.expect_value(t, glm.determinant(m), 2) testing.expect_value(t, glm.adjugate(m) * m, 2 * linalg.identity(matrix[2,2]int)) testing.expect_value(t, hlm.adjugate(m), expected) + testing.expect_value(t, hlm.determinant(m), 2) testing.expect_value(t, hlm.adjugate(m) * m, 2 * linalg.identity(matrix[2,2]int)) } @@ -44,9 +46,11 @@ test_adjugate_3x3 :: proc(t: ^testing.T) { testing.expect_value(t, linalg.adjugate(m) * m, -6 * linalg.identity(matrix[3,3]int)) testing.expect_value(t, glm.adjugate(m), expected) + testing.expect_value(t, glm.determinant(m), -6) testing.expect_value(t, glm.adjugate(m) * m, -6 * linalg.identity(matrix[3,3]int)) testing.expect_value(t, hlm.adjugate(m), expected) + testing.expect_value(t, hlm.determinant(m), -6) testing.expect_value(t, hlm.adjugate(m) * m, -6 * linalg.identity(matrix[3,3]int)) } @@ -60,7 +64,7 @@ test_adjugate_4x4 :: proc(t: ^testing.T) { } expected := matrix[4,4]int { -144, 266, -92, -16, - 57, 92, -5, -16, + -57, 92, -5, -16, 105, -142, 55, 2, 33, -96, 9, -6, } @@ -69,8 +73,10 @@ test_adjugate_4x4 :: proc(t: ^testing.T) { testing.expect_value(t, linalg.adjugate(m) * m, -174 * linalg.identity(matrix[4,4]int)) testing.expect_value(t, glm.adjugate(m), expected) + testing.expect_value(t, glm.determinant(m), -174) testing.expect_value(t, glm.adjugate(m) * m, -174 * linalg.identity(matrix[4,4]int)) testing.expect_value(t, hlm.adjugate(m), expected) + testing.expect_value(t, hlm.determinant(m), -174) testing.expect_value(t, hlm.adjugate(m) * m, -174 * linalg.identity(matrix[4,4]int)) } \ No newline at end of file From f23e226854c6777a27f13f70759fb022acafd532 Mon Sep 17 00:00:00 2001 From: Sebastian Pahnke Date: Sat, 28 Dec 2024 09:23:43 +0100 Subject: [PATCH 4/7] Rename adjugate to cofactor to keep existing usages for inverse and determinant correct and add new adjugate procedures --- core/math/linalg/general.odin | 84 ++++++++++++++++++++------ core/math/linalg/glsl/linalg_glsl.odin | 84 ++++++++++++++++++++------ core/math/linalg/hlsl/linalg_hlsl.odin | 84 ++++++++++++++++++++------ 3 files changed, 198 insertions(+), 54 deletions(-) diff --git a/core/math/linalg/general.odin b/core/math/linalg/general.odin index 8c4f2954a..4a0150972 100644 --- a/core/math/linalg/general.odin +++ b/core/math/linalg/general.odin @@ -417,6 +417,13 @@ adjugate :: proc{ matrix4x4_adjugate, } +cofactor :: proc{ + matrix1x1_cofactor, + matrix2x2_cofactor, + matrix3x3_cofactor, + matrix4x4_cofactor, +} + inverse_transpose :: proc{ matrix1x1_inverse_transpose, matrix2x2_inverse_transpose, @@ -479,9 +486,9 @@ matrix3x3_determinant :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (det: T) # } @(require_results) matrix4x4_determinant :: proc "contextless" (m: $M/matrix[4, 4]$T) -> (det: T) #no_bounds_check { - a := adjugate(m) + c := cofactor(m) for i in 0..<4 { - det += m[0, i] * a[0, i] + det += m[0, i] * c[0, i] } return } @@ -497,6 +504,47 @@ matrix1x1_adjugate :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) #no_bo @(require_results) matrix2x2_adjugate :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) #no_bounds_check { + y[0, 0] = +x[1, 1] + y[0, 1] = -x[0, 1] + y[1, 0] = -x[1, 0] + y[1, 1] = +x[0, 0] + return +} + +@(require_results) +matrix3x3_adjugate :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { + y[0, 0] = +(m[1, 1] * m[2, 2] - m[2, 1] * m[1, 2]) + y[1, 0] = -(m[1, 0] * m[2, 2] - m[2, 0] * m[1, 2]) + y[2, 0] = +(m[1, 0] * m[2, 1] - m[2, 0] * m[1, 1]) + y[0, 1] = -(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[2, 1] = -(m[0, 0] * m[2, 1] - m[2, 0] * m[0, 1]) + y[0, 2] = +(m[0, 1] * m[1, 2] - m[1, 1] * m[0, 2]) + y[1, 2] = -(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) #no_bounds_check { + 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, j, i) + } + } + return +} + + +@(require_results) +matrix1x1_cofactor :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) #no_bounds_check { + y = x + return +} + +@(require_results) +matrix2x2_cofactor :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) #no_bounds_check { y[0, 0] = +x[1, 1] y[0, 1] = -x[1, 0] y[1, 0] = -x[0, 1] @@ -505,7 +553,7 @@ matrix2x2_adjugate :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) #no_bo } @(require_results) -matrix3x3_adjugate :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { +matrix3x3_cofactor :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { 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]) @@ -520,7 +568,7 @@ matrix3x3_adjugate :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) #no_bo @(require_results) -matrix4x4_adjugate :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { +matrix4x4_cofactor :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { for i in 0..<4 { for j in 0..<4 { sign: T = 1 if (i + j) % 2 == 0 else -1 @@ -556,19 +604,19 @@ matrix2x2_inverse_transpose :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: @(require_results) matrix3x3_inverse_transpose :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) + c := cofactor(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 + y[i, j] = c[i, j] / d } } } else { id := 1/d for i in 0..<3 { for j in 0..<3 { - y[i, j] = a[i, j] * id + y[i, j] = c[i, j] * id } } } @@ -577,22 +625,22 @@ matrix3x3_inverse_transpose :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: @(require_results) matrix4x4_inverse_transpose :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) + c := cofactor(x) d: T for i in 0..<4 { - d += x[0, i] * a[0, i] + d += x[0, i] * c[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 + y[i, j] = c[i, j] / d } } } else { id := 1/d for i in 0..<4 { for j in 0..<4 { - y[i, j] = a[i, j] * id + y[i, j] = c[i, j] * id } } } @@ -625,19 +673,19 @@ matrix2x2_inverse :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) #no_bou @(require_results) matrix3x3_inverse :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) + c := cofactor(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 + y[i, j] = c[j, i] / d } } } else { id := 1/d for i in 0..<3 { for j in 0..<3 { - y[i, j] = a[j, i] * id + y[i, j] = c[j, i] * id } } } @@ -646,22 +694,22 @@ matrix3x3_inverse :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bou @(require_results) matrix4x4_inverse :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) + c := cofactor(x) d: T for i in 0..<4 { - d += x[0, i] * a[0, i] + d += x[0, i] * c[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 + y[i, j] = c[j, i] / d } } } else { id := 1/d for i in 0..<4 { for j in 0..<4 { - y[i, j] = a[j, i] * id + y[i, j] = c[j, i] * id } } } diff --git a/core/math/linalg/glsl/linalg_glsl.odin b/core/math/linalg/glsl/linalg_glsl.odin index ca61891cb..bd2cf416a 100644 --- a/core/math/linalg/glsl/linalg_glsl.odin +++ b/core/math/linalg/glsl/linalg_glsl.odin @@ -1882,6 +1882,13 @@ adjugate :: proc{ adjugate_matrix4x4, } +cofactor :: proc{ + cofactor_matrix1x1, + cofactor_matrix2x2, + cofactor_matrix3x3, + cofactor_matrix4x4, +} + inverse_transpose :: proc{ inverse_transpose_matrix1x1, inverse_transpose_matrix2x2, @@ -1944,9 +1951,9 @@ determinant_matrix3x3 :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (det: T) { } @(require_results) determinant_matrix4x4 :: proc "contextless" (m: $M/matrix[4, 4]$T) -> (det: T) { - a := adjugate(m) + c := cofactor(m) #no_bounds_check for i in 0..<4 { - det += m[0, i] * a[0, i] + det += m[0, i] * c[0, i] } return } @@ -1962,6 +1969,47 @@ adjugate_matrix1x1 :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { @(require_results) adjugate_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { + y[0, 0] = +x[1, 1] + y[0, 1] = -x[0, 1] + y[1, 0] = -x[1, 0] + 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[1, 0] = -(m[1, 0] * m[2, 2] - m[2, 0] * m[1, 2]) + y[2, 0] = +(m[1, 0] * m[2, 1] - m[2, 0] * m[1, 1]) + y[0, 1] = -(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[2, 1] = -(m[0, 0] * m[2, 1] - m[2, 0] * m[0, 1]) + y[0, 2] = +(m[0, 1] * m[1, 2] - m[1, 1] * m[0, 2]) + y[1, 2] = -(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, j, i) + } + } + return +} + + +@(require_results) +cofactor_matrix1x1 :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { + y = x + return +} + +@(require_results) +cofactor_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] @@ -1970,7 +2018,7 @@ adjugate_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { } @(require_results) -adjugate_matrix3x3 :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) { +cofactor_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]) @@ -1985,7 +2033,7 @@ adjugate_matrix3x3 :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) { @(require_results) -adjugate_matrix4x4 :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) { +cofactor_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 @@ -2021,19 +2069,19 @@ inverse_transpose_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: @(require_results) inverse_transpose_matrix3x3 :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) + c := cofactor(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 + y[i, j] = c[i, j] / d } } } else { id := 1/d for i in 0..<3 { for j in 0..<3 { - y[i, j] = a[i, j] * id + y[i, j] = c[i, j] * id } } } @@ -2042,22 +2090,22 @@ inverse_transpose_matrix3x3 :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: @(require_results) inverse_transpose_matrix4x4 :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) + c := cofactor(x) d: T for i in 0..<4 { - d += x[0, i] * a[0, i] + d += x[0, i] * c[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 + y[i, j] = c[i, j] / d } } } else { id := 1/d for i in 0..<4 { for j in 0..<4 { - y[i, j] = a[i, j] * id + y[i, j] = c[i, j] * id } } } @@ -2090,19 +2138,19 @@ inverse_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { @(require_results) inverse_matrix3x3 :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) + c := cofactor(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 + y[i, j] = c[j, i] / d } } } else { id := 1/d for i in 0..<3 { for j in 0..<3 { - y[i, j] = a[j, i] * id + y[i, j] = c[j, i] * id } } } @@ -2111,22 +2159,22 @@ inverse_matrix3x3 :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bou @(require_results) inverse_matrix4x4 :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) + c := cofactor(x) d: T for i in 0..<4 { - d += x[0, i] * a[0, i] + d += x[0, i] * c[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 + y[i, j] = c[j, i] / d } } } else { id := 1/d for i in 0..<4 { for j in 0..<4 { - y[i, j] = a[j, i] * id + y[i, j] = c[j, i] * id } } } diff --git a/core/math/linalg/hlsl/linalg_hlsl.odin b/core/math/linalg/hlsl/linalg_hlsl.odin index a89fdddd3..cca70f9c8 100644 --- a/core/math/linalg/hlsl/linalg_hlsl.odin +++ b/core/math/linalg/hlsl/linalg_hlsl.odin @@ -1514,6 +1514,13 @@ adjugate :: proc{ adjugate_matrix4x4, } +cofactor :: proc{ + cofactor_matrix1x1, + cofactor_matrix2x2, + cofactor_matrix3x3, + cofactor_matrix4x4, +} + inverse_transpose :: proc{ inverse_transpose_matrix1x1, inverse_transpose_matrix2x2, @@ -1568,9 +1575,9 @@ determinant_matrix3x3 :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (det: T) { } @(require_results) determinant_matrix4x4 :: proc "contextless" (m: $M/matrix[4, 4]$T) -> (det: T) { - a := adjugate(m) + c := cofactor(m) #no_bounds_check for i in 0..<4 { - det += m[0, i] * a[0, i] + det += m[0, i] * c[0, i] } return } @@ -1586,6 +1593,47 @@ adjugate_matrix1x1 :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { @(require_results) adjugate_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { + y[0, 0] = +x[1, 1] + y[0, 1] = -x[0, 1] + y[1, 0] = -x[1, 0] + 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[1, 0] = -(m[1, 0] * m[2, 2] - m[2, 0] * m[1, 2]) + y[2, 0] = +(m[1, 0] * m[2, 1] - m[2, 0] * m[1, 1]) + y[0, 1] = -(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[2, 1] = -(m[0, 0] * m[2, 1] - m[2, 0] * m[0, 1]) + y[0, 2] = +(m[0, 1] * m[1, 2] - m[1, 1] * m[0, 2]) + y[1, 2] = -(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, j, i) + } + } + return +} + + +@(require_results) +cofactor_matrix1x1 :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) { + y = x + return +} + +@(require_results) +cofactor_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] @@ -1594,7 +1642,7 @@ adjugate_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { } @(require_results) -adjugate_matrix3x3 :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) { +cofactor_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]) @@ -1609,7 +1657,7 @@ adjugate_matrix3x3 :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) { @(require_results) -adjugate_matrix4x4 :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) { +cofactor_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 @@ -1645,19 +1693,19 @@ inverse_transpose_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: @(require_results) inverse_transpose_matrix3x3 :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) + c := cofactor(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 + y[i, j] = c[i, j] / d } } } else { id := 1/d for i in 0..<3 { for j in 0..<3 { - y[i, j] = a[i, j] * id + y[i, j] = c[i, j] * id } } } @@ -1666,22 +1714,22 @@ inverse_transpose_matrix3x3 :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: @(require_results) inverse_transpose_matrix4x4 :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) + c := cofactor(x) d: T for i in 0..<4 { - d += x[0, i] * a[0, i] + d += x[0, i] * c[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 + y[i, j] = c[i, j] / d } } } else { id := 1/d for i in 0..<4 { for j in 0..<4 { - y[i, j] = a[i, j] * id + y[i, j] = c[i, j] * id } } } @@ -1714,19 +1762,19 @@ inverse_matrix2x2 :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) { @(require_results) inverse_matrix3x3 :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) + c := cofactor(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 + y[i, j] = c[j, i] / d } } } else { id := 1/d for i in 0..<3 { for j in 0..<3 { - y[i, j] = a[j, i] * id + y[i, j] = c[j, i] * id } } } @@ -1735,22 +1783,22 @@ inverse_matrix3x3 :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bou @(require_results) inverse_matrix4x4 :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check { - a := adjugate(x) + c := cofactor(x) d: T for i in 0..<4 { - d += x[0, i] * a[0, i] + d += x[0, i] * c[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 + y[i, j] = c[j, i] / d } } } else { id := 1/d for i in 0..<4 { for j in 0..<4 { - y[i, j] = a[j, i] * id + y[i, j] = c[j, i] * id } } } From 02a9d8560fb82471ddc89154dee1e9e30afb88ab Mon Sep 17 00:00:00 2001 From: Sebastian Pahnke Date: Sat, 28 Dec 2024 09:33:58 +0100 Subject: [PATCH 5/7] Test symmetry --- tests/issues/test_issue_4584.odin | 30 +++++++++++++++++++++--------- 1 file changed, 21 insertions(+), 9 deletions(-) diff --git a/tests/issues/test_issue_4584.odin b/tests/issues/test_issue_4584.odin index f14f1f1d3..27e5cd789 100644 --- a/tests/issues/test_issue_4584.odin +++ b/tests/issues/test_issue_4584.odin @@ -8,6 +8,7 @@ import hlm "core:math/linalg/hlsl" @test test_adjugate_2x2 :: proc(t: ^testing.T) { + I := linalg.identity(matrix[2,2]int) m := matrix[2,2]int { -3, 2, -1, 0, @@ -18,19 +19,23 @@ test_adjugate_2x2 :: proc(t: ^testing.T) { } testing.expect_value(t, linalg.adjugate(m), expected) testing.expect_value(t, linalg.determinant(m), 2) - testing.expect_value(t, linalg.adjugate(m) * m, 2 * linalg.identity(matrix[2,2]int)) + testing.expect_value(t, linalg.adjugate(m) * m, 2 * I) + testing.expect_value(t, m * linalg.adjugate(m), 2 * I) testing.expect_value(t, glm.adjugate(m), expected) testing.expect_value(t, glm.determinant(m), 2) - testing.expect_value(t, glm.adjugate(m) * m, 2 * linalg.identity(matrix[2,2]int)) + testing.expect_value(t, glm.adjugate(m) * m, 2 * I) + testing.expect_value(t, m * glm.adjugate(m), 2 * I) testing.expect_value(t, hlm.adjugate(m), expected) testing.expect_value(t, hlm.determinant(m), 2) - testing.expect_value(t, hlm.adjugate(m) * m, 2 * linalg.identity(matrix[2,2]int)) + testing.expect_value(t, hlm.adjugate(m) * m, 2 * I) + testing.expect_value(t, m * hlm.adjugate(m), 2 * I) } @test test_adjugate_3x3 :: proc(t: ^testing.T) { + I := linalg.identity(matrix[3,3]int) m := matrix[3,3]int { -3, 2, -5, -1, 0, -2, @@ -43,19 +48,23 @@ test_adjugate_3x3 :: proc(t: ^testing.T) { } testing.expect_value(t, linalg.adjugate(m), expected) testing.expect_value(t, linalg.determinant(m), -6) - testing.expect_value(t, linalg.adjugate(m) * m, -6 * linalg.identity(matrix[3,3]int)) + testing.expect_value(t, linalg.adjugate(m) * m, -6 * I) + testing.expect_value(t, m * linalg.adjugate(m), -6 * I) testing.expect_value(t, glm.adjugate(m), expected) testing.expect_value(t, glm.determinant(m), -6) - testing.expect_value(t, glm.adjugate(m) * m, -6 * linalg.identity(matrix[3,3]int)) + testing.expect_value(t, glm.adjugate(m) * m, -6 * I) + testing.expect_value(t, m * glm.adjugate(m), -6 * I) testing.expect_value(t, hlm.adjugate(m), expected) testing.expect_value(t, hlm.determinant(m), -6) - testing.expect_value(t, hlm.adjugate(m) * m, -6 * linalg.identity(matrix[3,3]int)) + testing.expect_value(t, hlm.adjugate(m) * m, -6 * I) + testing.expect_value(t, m * hlm.adjugate(m), -6 * I) } @test test_adjugate_4x4 :: proc(t: ^testing.T) { + I := linalg.identity(matrix[4,4]int) m := matrix[4,4]int { -3, 2, -5, 1, -1, 0, -2, 2, @@ -70,13 +79,16 @@ test_adjugate_4x4 :: proc(t: ^testing.T) { } testing.expect_value(t, linalg.adjugate(m), expected) testing.expect_value(t, linalg.determinant(m), -174) - testing.expect_value(t, linalg.adjugate(m) * m, -174 * linalg.identity(matrix[4,4]int)) + testing.expect_value(t, linalg.adjugate(m) * m, -174 * I) + testing.expect_value(t, m * linalg.adjugate(m), -174 * I) testing.expect_value(t, glm.adjugate(m), expected) testing.expect_value(t, glm.determinant(m), -174) - testing.expect_value(t, glm.adjugate(m) * m, -174 * linalg.identity(matrix[4,4]int)) + testing.expect_value(t, glm.adjugate(m) * m, -174 * I) + testing.expect_value(t, m * glm.adjugate(m), -174 * I) testing.expect_value(t, hlm.adjugate(m), expected) testing.expect_value(t, hlm.determinant(m), -174) - testing.expect_value(t, hlm.adjugate(m) * m, -174 * linalg.identity(matrix[4,4]int)) + testing.expect_value(t, hlm.adjugate(m) * m, -174 * I) + testing.expect_value(t, m * hlm.adjugate(m), -174 * I) } \ No newline at end of file From ec5ee19c01c30e8fd01c9f04a2249c4d3d18b50d Mon Sep 17 00:00:00 2001 From: Sebastian Pahnke Date: Sat, 28 Dec 2024 10:24:37 +0100 Subject: [PATCH 6/7] Add regression tests for matrix inverse --- tests/issues/test_issue_4584.odin | 95 +++++++++++++++++++++++++++++++ 1 file changed, 95 insertions(+) diff --git a/tests/issues/test_issue_4584.odin b/tests/issues/test_issue_4584.odin index 27e5cd789..a709926d2 100644 --- a/tests/issues/test_issue_4584.odin +++ b/tests/issues/test_issue_4584.odin @@ -2,6 +2,7 @@ package test_issues import "core:testing" +import "core:log" import "core:math/linalg" import glm "core:math/linalg/glsl" import hlm "core:math/linalg/hlsl" @@ -91,4 +92,98 @@ test_adjugate_4x4 :: proc(t: ^testing.T) { testing.expect_value(t, hlm.determinant(m), -174) testing.expect_value(t, hlm.adjugate(m) * m, -174 * I) testing.expect_value(t, m * hlm.adjugate(m), -174 * I) +} + +@test +test_inverse_regression_2x2 :: proc(t: ^testing.T) { + I := linalg.identity(matrix[2,2]f32) + m := matrix[2,2]f32 { + -3, 2, + -1, 0, + } + expected := matrix[2,2]f32 { + 0.0, -1.0, + 1.0/2.0, -3.0/2.0, + } + expect_float_matrix_value(t, linalg.inverse(m), expected) + expect_float_matrix_value(t, linalg.inverse(m) * m, I) + expect_float_matrix_value(t, m * linalg.inverse(m), I) + + expect_float_matrix_value(t, glm.inverse(m), expected) + expect_float_matrix_value(t, glm.inverse(m) * m, I) + expect_float_matrix_value(t, m * glm.inverse(m), I) + + expect_float_matrix_value(t, hlm.inverse(m), expected) + expect_float_matrix_value(t, hlm.inverse(m) * m, I) + expect_float_matrix_value(t, m * hlm.inverse(m), I) +} + +@test +test_inverse_regression_3x3 :: proc(t: ^testing.T) { + I := linalg.identity(matrix[3,3]f32) + m := matrix[3,3]f32 { + -3, 2, -5, + -1, 0, -2, + 3, -4, 1, + } + expected := matrix[3,3]f32 { + 4.0/3.0, -3.0, 2.0/3.0, + 5.0/6.0, -2.0, 1.0/6.0, + -2.0/3.0, 1.0, -1.0/3.0, + } + expect_float_matrix_value(t, linalg.inverse(m), expected) + expect_float_matrix_value(t, linalg.inverse(m) * m, I) + expect_float_matrix_value(t, m * linalg.inverse(m), I) + + expect_float_matrix_value(t, glm.inverse(m), expected) + expect_float_matrix_value(t, glm.inverse(m) * m, I) + expect_float_matrix_value(t, m * glm.inverse(m), I) + + expect_float_matrix_value(t, hlm.inverse(m), expected) + expect_float_matrix_value(t, hlm.inverse(m) * m, I) + expect_float_matrix_value(t, m * hlm.inverse(m), I) +} + +@test +test_inverse_regression_4x4 :: proc(t: ^testing.T) { + I := linalg.identity(matrix[4,4]f32) + m := matrix[4,4]f32 { + -3, 2, -5, 1, + -1, 0, -2, 2, + 3, -4, 1, 3, + 4, 5, 6, 7, + } + expected := matrix[4,4]f32 { + 24.0/29.0, -133.0/87.0, 46.0/87.0, 8.0/87.0, + 19.0/58.0, -46.0/87.0, 5.0/174.0, 8.0/87.0, + -35.0/58.0, 71.0/87.0, -55.0/174.0, -1.0/87.0, + -11.0/58.0, 16.0/29.0, -3.0/58.0, 1.0/29.0, + } + expect_float_matrix_value(t, linalg.inverse(m), expected) + expect_float_matrix_value(t, linalg.inverse(m) * m, I) + expect_float_matrix_value(t, m * linalg.inverse(m), I) + + expect_float_matrix_value(t, glm.inverse(m), expected) + expect_float_matrix_value(t, glm.inverse(m) * m, I) + expect_float_matrix_value(t, m * glm.inverse(m), I) + + expect_float_matrix_value(t, hlm.inverse(m), expected) + expect_float_matrix_value(t, hlm.inverse(m) * m, I) + expect_float_matrix_value(t, m * hlm.inverse(m), I) +} + +@(private="file") +expect_float_matrix_value :: proc(t: ^testing.T, value, expected: $M/matrix[$N, N]f32, loc := #caller_location, value_expr := #caller_expression(value)) -> bool { + ok := true + outer: for i in 0.. 1e-6 { + ok = false + break outer + } + } + } + if !ok do log.errorf("expected %v to be %v, got %v", value_expr, expected, value, location=loc) + return ok } \ No newline at end of file From b21fc1923307602500c696e418efa951a5658ed8 Mon Sep 17 00:00:00 2001 From: Sebastian Pahnke Date: Sat, 28 Dec 2024 10:29:05 +0100 Subject: [PATCH 7/7] Add regression tests for inverse_transpose --- tests/issues/test_issue_4584.odin | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/tests/issues/test_issue_4584.odin b/tests/issues/test_issue_4584.odin index a709926d2..9eea23420 100644 --- a/tests/issues/test_issue_4584.odin +++ b/tests/issues/test_issue_4584.odin @@ -106,14 +106,17 @@ test_inverse_regression_2x2 :: proc(t: ^testing.T) { 1.0/2.0, -3.0/2.0, } expect_float_matrix_value(t, linalg.inverse(m), expected) + expect_float_matrix_value(t, linalg.inverse_transpose(m), linalg.transpose(expected)) expect_float_matrix_value(t, linalg.inverse(m) * m, I) expect_float_matrix_value(t, m * linalg.inverse(m), I) expect_float_matrix_value(t, glm.inverse(m), expected) + expect_float_matrix_value(t, glm.inverse_transpose(m), glm.transpose(expected)) expect_float_matrix_value(t, glm.inverse(m) * m, I) expect_float_matrix_value(t, m * glm.inverse(m), I) expect_float_matrix_value(t, hlm.inverse(m), expected) + expect_float_matrix_value(t, hlm.inverse_transpose(m), hlm.transpose(expected)) expect_float_matrix_value(t, hlm.inverse(m) * m, I) expect_float_matrix_value(t, m * hlm.inverse(m), I) } @@ -132,14 +135,17 @@ test_inverse_regression_3x3 :: proc(t: ^testing.T) { -2.0/3.0, 1.0, -1.0/3.0, } expect_float_matrix_value(t, linalg.inverse(m), expected) + expect_float_matrix_value(t, linalg.inverse_transpose(m), linalg.transpose(expected)) expect_float_matrix_value(t, linalg.inverse(m) * m, I) expect_float_matrix_value(t, m * linalg.inverse(m), I) expect_float_matrix_value(t, glm.inverse(m), expected) + expect_float_matrix_value(t, glm.inverse_transpose(m), glm.transpose(expected)) expect_float_matrix_value(t, glm.inverse(m) * m, I) expect_float_matrix_value(t, m * glm.inverse(m), I) expect_float_matrix_value(t, hlm.inverse(m), expected) + expect_float_matrix_value(t, hlm.inverse_transpose(m), hlm.transpose(expected)) expect_float_matrix_value(t, hlm.inverse(m) * m, I) expect_float_matrix_value(t, m * hlm.inverse(m), I) } @@ -160,14 +166,17 @@ test_inverse_regression_4x4 :: proc(t: ^testing.T) { -11.0/58.0, 16.0/29.0, -3.0/58.0, 1.0/29.0, } expect_float_matrix_value(t, linalg.inverse(m), expected) + expect_float_matrix_value(t, linalg.inverse_transpose(m), linalg.transpose(expected)) expect_float_matrix_value(t, linalg.inverse(m) * m, I) expect_float_matrix_value(t, m * linalg.inverse(m), I) expect_float_matrix_value(t, glm.inverse(m), expected) + expect_float_matrix_value(t, glm.inverse_transpose(m), glm.transpose(expected)) expect_float_matrix_value(t, glm.inverse(m) * m, I) expect_float_matrix_value(t, m * glm.inverse(m), I) expect_float_matrix_value(t, hlm.inverse(m), expected) + expect_float_matrix_value(t, hlm.inverse_transpose(m), hlm.transpose(expected)) expect_float_matrix_value(t, hlm.inverse(m) * m, I) expect_float_matrix_value(t, m * hlm.inverse(m), I) }