Merge pull request #1245 from odin-lang/new-matrix-type

`matrix` type

core/runtime/core.odin

@@ -162,6 +162,14 @@ Type_Info_Relative_Slice :: struct {
 	slice:        ^Type_Info,
 	base_integer: ^Type_Info,
 }
+Type_Info_Matrix :: struct {
+	elem:         ^Type_Info,
+	elem_size:    int,
+	elem_stride:  int, // elem_stride >= row_count
+	row_count:    int,
+	column_count: int,
+	// Total element count = column_count * elem_stride
+}
 
 Type_Info_Flag :: enum u8 {
 	Comparable     = 0,
@@ -202,6 +210,7 @@ Type_Info :: struct {
 		Type_Info_Simd_Vector,
 		Type_Info_Relative_Pointer,
 		Type_Info_Relative_Slice,
+		Type_Info_Matrix,
 	},
 }
 
@@ -233,6 +242,7 @@ Typeid_Kind :: enum u8 {
 	Simd_Vector,
 	Relative_Pointer,
 	Relative_Slice,
+	Matrix,
 }
 #assert(len(Typeid_Kind) < 32)
 

core/runtime/core_builtin_matrix.odin

@@ -0,0 +1,274 @@
+package runtime
+
+import "core:intrinsics"
+_ :: intrinsics
+
+
+@(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)
+hermitian_adjoint :: proc(m: $M/matrix[$N, N]$T) -> M where intrinsics.type_is_complex(T), N >= 1 {
+	return conj(transpose(m))
+}
+
+@(builtin)
+matrix_trace :: proc(m: $M/matrix[$N, N]$T) -> (trace: T) {
+	for i in 0..<N {
+		trace += m[i, i]
+	}
+	return
+}
+
+@(builtin)
+matrix_minor :: proc(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)
+matrix1x1_determinant :: proc(m: $M/matrix[1, 1]$T) -> (det: T) {
+	return m[0, 0]
+}
+
+@(builtin)
+matrix2x2_determinant :: proc(m: $M/matrix[2, 2]$T) -> (det: T) {
+	return m[0, 0]*m[1, 1] - m[0, 1]*m[1, 0]
+}
+@(builtin)
+matrix3x3_determinant :: proc(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)
+matrix4x4_determinant :: proc(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)
+matrix1x1_adjugate :: proc(x: $M/matrix[1, 1]$T) -> (y: M) {
+	y = x
+	return
+}
+
+@(builtin)
+matrix2x2_adjugate :: proc(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)
+matrix3x3_adjugate :: proc(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)
+matrix4x4_adjugate :: proc(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)
+matrix1x1_inverse_transpose :: proc(x: $M/matrix[1, 1]$T) -> (y: M) {
+	y[0, 0] = 1/x[0, 0]
+	return
+}
+
+@(builtin)
+matrix2x2_inverse_transpose :: proc(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[1, 0] / d
+		y[0, 1] = -x[0, 1] / d
+		y[1, 1] = +x[0, 0] / d
+	} else {
+		id := 1 / d
+		y[0, 0] = +x[1, 1] * id
+		y[1, 0] = -x[1, 0] * id
+		y[0, 1] = -x[0, 1] * id
+		y[1, 1] = +x[0, 0] * id
+	}
+	return
+}
+
+@(builtin)
+matrix3x3_inverse_transpose :: proc(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)
+matrix4x4_inverse_transpose :: proc(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)
+matrix1x1_inverse :: proc(x: $M/matrix[1, 1]$T) -> (y: M) {
+	y[0, 0] = 1/x[0, 0]
+	return
+}
+
+@(builtin)
+matrix2x2_inverse :: proc(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[1, 0] / d
+		y[1, 0] = x[0, 1] / d
+		y[1, 1] = x[0, 0] / d
+	} else {
+		id := 1 / d
+		y[0, 0] = x[1, 1] * id
+		y[0, 1] = x[1, 0] * id
+		y[1, 0] = x[0, 1] * id
+		y[1, 1] = x[0, 0] * id
+	}
+	return
+}
+
+@(builtin)
+matrix3x3_inverse :: proc(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)
+matrix4x4_inverse :: proc(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
+}

core/runtime/error_checks.odin

@@ -96,6 +96,29 @@ dynamic_array_expr_error :: proc "contextless" (file: string, line, column: i32,
 }
 
 
+matrix_bounds_check_error :: proc "contextless" (file: string, line, column: i32, row_index, column_index, row_count, column_count: int) {
+	if 0 <= row_index && row_index < row_count && 
+	   0 <= column_index && column_index < column_count {
+		return
+	}
+	handle_error :: proc "contextless" (file: string, line, column: i32, row_index, column_index, row_count, column_count: int) {
+		print_caller_location(Source_Code_Location{file, line, column, ""})
+		print_string(" Matrix indices [")
+		print_i64(i64(row_index))
+		print_string(", ")
+		print_i64(i64(column_index))
+		print_string(" is out of bounds range [0..<")
+		print_i64(i64(row_count))
+		print_string(", 0..<")
+		print_i64(i64(column_count))
+		print_string("]")
+		print_byte('\n')
+		bounds_trap()
+	}
+	handle_error(file, line, column, row_index, column_index, row_count, column_count)
+}
+
+
 type_assertion_check :: proc "contextless" (ok: bool, file: string, line, column: i32, from, to: typeid) {
 	if ok {
 		return

core/runtime/internal.odin

@@ -2,15 +2,15 @@ package runtime
 
 import "core:intrinsics"
 
-bswap_16 :: proc "none" (x: u16) -> u16 {
+bswap_16 :: proc "contextless" (x: u16) -> u16 {
 	return x>>8 | x<<8
 }
 
-bswap_32 :: proc "none" (x: u32) -> u32 {
+bswap_32 :: proc "contextless" (x: u32) -> u32 {
 	return x>>24 | (x>>8)&0xff00 | (x<<8)&0xff0000 | x<<24
 }
 
-bswap_64 :: proc "none" (x: u64) -> u64 {
+bswap_64 :: proc "contextless" (x: u64) -> u64 {
 	z := x
 	z = (z & 0x00000000ffffffff) << 32 | (z & 0xffffffff00000000) >> 32
 	z = (z & 0x0000ffff0000ffff) << 16 | (z & 0xffff0000ffff0000) >> 16
@@ -18,7 +18,7 @@ bswap_64 :: proc "none" (x: u64) -> u64 {
 	return z
 }
 
-bswap_128 :: proc "none" (x: u128) -> u128 {
+bswap_128 :: proc "contextless" (x: u128) -> u128 {
 	z := transmute([4]u32)x
 	z[0] = bswap_32(z[3])
 	z[1] = bswap_32(z[2])
@@ -27,33 +27,27 @@ bswap_128 :: proc "none" (x: u128) -> u128 {
 	return transmute(u128)z
 }
 
-bswap_f16 :: proc "none" (f: f16) -> f16 {
+bswap_f16 :: proc "contextless" (f: f16) -> f16 {
 	x := transmute(u16)f
 	z := bswap_16(x)
 	return transmute(f16)z
 
 }
 
-bswap_f32 :: proc "none" (f: f32) -> f32 {
+bswap_f32 :: proc "contextless" (f: f32) -> f32 {
 	x := transmute(u32)f
 	z := bswap_32(x)
 	return transmute(f32)z
 
 }
 
-bswap_f64 :: proc "none" (f: f64) -> f64 {
+bswap_f64 :: proc "contextless" (f: f64) -> f64 {
 	x := transmute(u64)f
 	z := bswap_64(x)
 	return transmute(f64)z
 }
 
 
-
-ptr_offset :: #force_inline proc "contextless" (ptr: $P/^$T, n: int) -> P {
-	new := int(uintptr(ptr)) + size_of(T)*n
-	return P(uintptr(new))
-}
-
 is_power_of_two_int :: #force_inline proc(x: int) -> bool {
 	if x <= 0 {
 		return false
@@ -828,12 +822,14 @@ floattidf_unsigned :: proc "c" (a: u128) -> f64 {
 
 @(link_name="__fixunsdfti")
 fixunsdfti :: #force_no_inline proc "c" (a: f64) -> u128 {
+	// TODO(bill): implement `fixunsdfti` correctly
 	x := u64(a)
 	return u128(x)
 }
 
 @(link_name="__fixunsdfdi")
 fixunsdfdi :: #force_no_inline proc "c" (a: f64) -> i128 {
+	// TODO(bill): implement `fixunsdfdi` correctly
 	x := i64(a)
 	return i128(x)
 }

core/runtime/print.odin

@@ -370,5 +370,13 @@ print_type :: proc "contextless" (ti: ^Type_Info) {
 		print_type(info.base_integer)
 		print_string(") ")
 		print_type(info.slice)
+		
+	case Type_Info_Matrix:
+		print_string("matrix[")
+		print_u64(u64(info.row_count))
+		print_string(", ")
+		print_u64(u64(info.column_count))
+		print_string("]")
+		print_type(info.elem)
 	}
 }