core_builtin_matrix.odin 5.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274
  1. package runtime
  2. import "core:intrinsics"
  3. _ :: intrinsics
  4. @(builtin)
  5. determinant :: proc{
  6. matrix1x1_determinant,
  7. matrix2x2_determinant,
  8. matrix3x3_determinant,
  9. matrix4x4_determinant,
  10. }
  11. @(builtin)
  12. adjugate :: proc{
  13. matrix1x1_adjugate,
  14. matrix2x2_adjugate,
  15. matrix3x3_adjugate,
  16. matrix4x4_adjugate,
  17. }
  18. @(builtin)
  19. inverse_transpose :: proc{
  20. matrix1x1_inverse_transpose,
  21. matrix2x2_inverse_transpose,
  22. matrix3x3_inverse_transpose,
  23. matrix4x4_inverse_transpose,
  24. }
  25. @(builtin)
  26. inverse :: proc{
  27. matrix1x1_inverse,
  28. matrix2x2_inverse,
  29. matrix3x3_inverse,
  30. matrix4x4_inverse,
  31. }
  32. @(builtin, require_results)
  33. hermitian_adjoint :: proc "contextless" (m: $M/matrix[$N, N]$T) -> M where intrinsics.type_is_complex(T), N >= 1 {
  34. return conj(transpose(m))
  35. }
  36. @(builtin, require_results)
  37. matrix_trace :: proc "contextless" (m: $M/matrix[$N, N]$T) -> (trace: T) {
  38. for i in 0..<N {
  39. trace += m[i, i]
  40. }
  41. return
  42. }
  43. @(builtin, require_results)
  44. matrix_minor :: proc "contextless" (m: $M/matrix[$N, N]$T, row, column: int) -> (minor: T) where N > 1 {
  45. K :: N-1
  46. cut_down: matrix[K, K]T
  47. for col_idx in 0..<K {
  48. j := col_idx + int(col_idx >= column)
  49. for row_idx in 0..<K {
  50. i := row_idx + int(row_idx >= row)
  51. cut_down[row_idx, col_idx] = m[i, j]
  52. }
  53. }
  54. return determinant(cut_down)
  55. }
  56. @(builtin, require_results)
  57. matrix1x1_determinant :: proc "contextless" (m: $M/matrix[1, 1]$T) -> (det: T) {
  58. return m[0, 0]
  59. }
  60. @(builtin, require_results)
  61. matrix2x2_determinant :: proc "contextless" (m: $M/matrix[2, 2]$T) -> (det: T) {
  62. return m[0, 0]*m[1, 1] - m[0, 1]*m[1, 0]
  63. }
  64. @(builtin, require_results)
  65. matrix3x3_determinant :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (det: T) {
  66. a := +m[0, 0] * (m[1, 1] * m[2, 2] - m[1, 2] * m[2, 1])
  67. b := -m[0, 1] * (m[1, 0] * m[2, 2] - m[1, 2] * m[2, 0])
  68. c := +m[0, 2] * (m[1, 0] * m[2, 1] - m[1, 1] * m[2, 0])
  69. return a + b + c
  70. }
  71. @(builtin, require_results)
  72. matrix4x4_determinant :: proc "contextless" (m: $M/matrix[4, 4]$T) -> (det: T) {
  73. a := adjugate(m)
  74. #no_bounds_check for i in 0..<4 {
  75. det += m[0, i] * a[0, i]
  76. }
  77. return
  78. }
  79. @(builtin, require_results)
  80. matrix1x1_adjugate :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) {
  81. y = x
  82. return
  83. }
  84. @(builtin, require_results)
  85. matrix2x2_adjugate :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) {
  86. y[0, 0] = +x[1, 1]
  87. y[0, 1] = -x[1, 0]
  88. y[1, 0] = -x[0, 1]
  89. y[1, 1] = +x[0, 0]
  90. return
  91. }
  92. @(builtin, require_results)
  93. matrix3x3_adjugate :: proc "contextless" (m: $M/matrix[3, 3]$T) -> (y: M) {
  94. y[0, 0] = +(m[1, 1] * m[2, 2] - m[2, 1] * m[1, 2])
  95. y[0, 1] = -(m[1, 0] * m[2, 2] - m[2, 0] * m[1, 2])
  96. y[0, 2] = +(m[1, 0] * m[2, 1] - m[2, 0] * m[1, 1])
  97. y[1, 0] = -(m[0, 1] * m[2, 2] - m[2, 1] * m[0, 2])
  98. y[1, 1] = +(m[0, 0] * m[2, 2] - m[2, 0] * m[0, 2])
  99. y[1, 2] = -(m[0, 0] * m[2, 1] - m[2, 0] * m[0, 1])
  100. y[2, 0] = +(m[0, 1] * m[1, 2] - m[1, 1] * m[0, 2])
  101. y[2, 1] = -(m[0, 0] * m[1, 2] - m[1, 0] * m[0, 2])
  102. y[2, 2] = +(m[0, 0] * m[1, 1] - m[1, 0] * m[0, 1])
  103. return
  104. }
  105. @(builtin, require_results)
  106. matrix4x4_adjugate :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) {
  107. for i in 0..<4 {
  108. for j in 0..<4 {
  109. sign: T = 1 if (i + j) % 2 == 0 else -1
  110. y[i, j] = sign * matrix_minor(x, i, j)
  111. }
  112. }
  113. return
  114. }
  115. @(builtin, require_results)
  116. matrix1x1_inverse_transpose :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) {
  117. y[0, 0] = 1/x[0, 0]
  118. return
  119. }
  120. @(builtin, require_results)
  121. matrix2x2_inverse_transpose :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) {
  122. d := x[0, 0]*x[1, 1] - x[0, 1]*x[1, 0]
  123. when intrinsics.type_is_integer(T) {
  124. y[0, 0] = +x[1, 1] / d
  125. y[1, 0] = -x[0, 1] / d
  126. y[0, 1] = -x[1, 0] / d
  127. y[1, 1] = +x[0, 0] / d
  128. } else {
  129. id := 1 / d
  130. y[0, 0] = +x[1, 1] * id
  131. y[1, 0] = -x[0, 1] * id
  132. y[0, 1] = -x[1, 0] * id
  133. y[1, 1] = +x[0, 0] * id
  134. }
  135. return
  136. }
  137. @(builtin, require_results)
  138. matrix3x3_inverse_transpose :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check {
  139. a := adjugate(x)
  140. d := determinant(x)
  141. when intrinsics.type_is_integer(T) {
  142. for i in 0..<3 {
  143. for j in 0..<3 {
  144. y[i, j] = a[i, j] / d
  145. }
  146. }
  147. } else {
  148. id := 1/d
  149. for i in 0..<3 {
  150. for j in 0..<3 {
  151. y[i, j] = a[i, j] * id
  152. }
  153. }
  154. }
  155. return
  156. }
  157. @(builtin, require_results)
  158. matrix4x4_inverse_transpose :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check {
  159. a := adjugate(x)
  160. d: T
  161. for i in 0..<4 {
  162. d += x[0, i] * a[0, i]
  163. }
  164. when intrinsics.type_is_integer(T) {
  165. for i in 0..<4 {
  166. for j in 0..<4 {
  167. y[i, j] = a[i, j] / d
  168. }
  169. }
  170. } else {
  171. id := 1/d
  172. for i in 0..<4 {
  173. for j in 0..<4 {
  174. y[i, j] = a[i, j] * id
  175. }
  176. }
  177. }
  178. return
  179. }
  180. @(builtin, require_results)
  181. matrix1x1_inverse :: proc "contextless" (x: $M/matrix[1, 1]$T) -> (y: M) {
  182. y[0, 0] = 1/x[0, 0]
  183. return
  184. }
  185. @(builtin, require_results)
  186. matrix2x2_inverse :: proc "contextless" (x: $M/matrix[2, 2]$T) -> (y: M) {
  187. d := x[0, 0]*x[1, 1] - x[0, 1]*x[1, 0]
  188. when intrinsics.type_is_integer(T) {
  189. y[0, 0] = +x[1, 1] / d
  190. y[0, 1] = -x[0, 1] / d
  191. y[1, 0] = -x[1, 0] / d
  192. y[1, 1] = +x[0, 0] / d
  193. } else {
  194. id := 1 / d
  195. y[0, 0] = +x[1, 1] * id
  196. y[0, 1] = -x[0, 1] * id
  197. y[1, 0] = -x[1, 0] * id
  198. y[1, 1] = +x[0, 0] * id
  199. }
  200. return
  201. }
  202. @(builtin, require_results)
  203. matrix3x3_inverse :: proc "contextless" (x: $M/matrix[3, 3]$T) -> (y: M) #no_bounds_check {
  204. a := adjugate(x)
  205. d := determinant(x)
  206. when intrinsics.type_is_integer(T) {
  207. for i in 0..<3 {
  208. for j in 0..<3 {
  209. y[i, j] = a[j, i] / d
  210. }
  211. }
  212. } else {
  213. id := 1/d
  214. for i in 0..<3 {
  215. for j in 0..<3 {
  216. y[i, j] = a[j, i] * id
  217. }
  218. }
  219. }
  220. return
  221. }
  222. @(builtin, require_results)
  223. matrix4x4_inverse :: proc "contextless" (x: $M/matrix[4, 4]$T) -> (y: M) #no_bounds_check {
  224. a := adjugate(x)
  225. d: T
  226. for i in 0..<4 {
  227. d += x[0, i] * a[0, i]
  228. }
  229. when intrinsics.type_is_integer(T) {
  230. for i in 0..<4 {
  231. for j in 0..<4 {
  232. y[i, j] = a[j, i] / d
  233. }
  234. }
  235. } else {
  236. id := 1/d
  237. for i in 0..<4 {
  238. for j in 0..<4 {
  239. y[i, j] = a[j, i] * id
  240. }
  241. }
  242. }
  243. return
  244. }