Skip to content

Commit

Permalink
Merge pull request #4633 from spahnke/fix-matrix-adjugate
Browse files Browse the repository at this point in the history
Fix matrix adjugate
  • Loading branch information
gingerBill authored Jan 6, 2025
2 parents c7739de + b21fc19 commit a3b4280
Show file tree
Hide file tree
Showing 6 changed files with 398 additions and 54 deletions.
84 changes: 66 additions & 18 deletions core/math/linalg/general.odin
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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
}
Expand All @@ -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]
Expand All @@ -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])
Expand All @@ -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
Expand Down Expand Up @@ -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
}
}
}
Expand All @@ -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
}
}
}
Expand Down Expand Up @@ -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
}
}
}
Expand All @@ -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
}
}
}
Expand Down
84 changes: 66 additions & 18 deletions core/math/linalg/glsl/linalg_glsl.odin
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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
}
Expand All @@ -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]
Expand All @@ -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])
Expand All @@ -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
Expand Down Expand Up @@ -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
}
}
}
Expand All @@ -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
}
}
}
Expand Down Expand Up @@ -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
}
}
}
Expand All @@ -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
}
}
}
Expand Down
Loading

0 comments on commit a3b4280

Please sign in to comment.