import itertools
import math
import random
from typing import Any,Generic,Iterable,Self,TypeVar, overload

#"pip install numpy"
try:
	import numpy as np
	_MATHUTIL_HAS_NP = True
except ModuleNotFoundError:
	_MATHUTIL_HAS_NP = False



#https://docs.python.org/3/reference/datamodel.html#special-method-names



NaN = float("nan")
Inf = float("inf")



if _MATHUTIL_HAS_NP:
	Tnum = TypeVar( "Tnum", float,int, np.float32,np.int32, np.float64,np.int64 )
else:
	Tnum = TypeVar( "Tnum", float,int )



class MatrixElems( Generic[Tnum] ):
	__slots__ = ( "_shape", "raw" )
	_shape:tuple[int,int]

	if _MATHUTIL_HAS_NP:
		#Note: for performance, 1D vectors are stored as 1D arrays instead of 2D arrays.  As a
		#	consequence, we also have to still store `._shape`.
		raw:np.ndarray

		@property
		def count( self:Self ): return len(self.raw.flat)

	else:
		raw:list[Tnum]

		@property
		def count( self:Self ): return len(self.raw)

	@property
	def nrows( self:Self ): return self._shape[0]
	@property
	def ncols( self:Self ): return self._shape[1]
	@property
	def shape( self:Self ): return self._shape

	@overload
	def __init__( self:Self, other:Self ): ...
	@overload
	def __init__( self:Self, shape:tuple[int,int] ): ...
	@overload
	def __init__( self:Self, shape:tuple[int,int], all_vals:Tnum ): ...
	@overload
	def __init__( self:Self, shape:tuple[int,int], rowmaj_data:Iterable[Tnum] ): ...
	@overload
	def __init__( self:Self, data2D:Iterable[Iterable[Tnum]] ): ...
	if _MATHUTIL_HAS_NP:
		@overload
		def __init__( self:Self, data2D:np.ndarray[Any,Tnum] ): ...
		@overload
		def __init__( self:Self, shape:tuple[int,int], data1D:np.ndarray[Any,Tnum] ): ...
		def __init__( self:Self, *args ):
			if len(args) == 1:
				arg0 = args[0]
				if   type(arg0) == MatrixElems: #other
					self._shape = tuple(arg0._shape)
					self.raw    = np.array(arg0.raw)
				elif type(arg0) == np.ndarray: #2D data
					data2D = arg0
					assert len(data2D.shape)==2, "Pass shape for 1D `np.ndarray`!"
					self._shape = tuple(data2D.shape)
					self.raw    = np.array(data2D)
				else: #data or shape
					if type(arg0)==tuple and len(arg0)==2 and type(arg0[0])==int: #shape
						shape = arg0
						assert type(shape[1]) == int
						self._shape = tuple(shape)
						self.raw = np.empty(shape)
					else: #2D data
						data2D = arg0
						self._shape = ( len(data2D), len(data2D[0]) )
						self.raw = np.array([ val for row in data2D for val in row ]).reshape(self._shape)
			else:
				assert len(args) == 2
				shape = args[0]
				assert len(shape)==2 and type(shape[0])==int and type(shape[1])==int
				self._shape = tuple(shape)

				if type(args[1]) == np.ndarray:
					data1D = args[1]
					assert len(data1D.shape)==1, "Pass shape for 1D `np.ndarray` only!"
					assert shape[0]*shape[1] == data1D.shape[0]
					self.raw    = np.array(data1D)
					if self._shape[0]!=1 and self._shape[1]!=1:
						self.raw = self.raw.reshape(shape)
				else:
					data = args[1]
					if type(data) in ( tuple, list ): #rowmaj_data
						self.raw = np.array(data)
						if self._shape[0]!=1 and self._shape[1]!=1:
							self.raw = self.raw.reshape(shape)
					else: #all_vals
						self.raw = np.full( shape, data )
	else:
		def __init__( self:Self, *args ):
			assert len(args) >= 1
			match args:
				case [ MatrixElems() as other ]:
					self._shape = tuple(other._shape)
					self.raw    = list(other.raw)

				case [ data_or_shape ]:
					if type(data_or_shape)==tuple and len(data_or_shape)==2 and type(data_or_shape[0])==int:
						assert type(data_or_shape[1]) == int
						self._shape = tuple(data_or_shape)
						self.raw    = [ NaN for k in range( self._shape[0] * self._shape[1] ) ]
					else:
						data2D = data_or_shape
						self._shape = ( len(data2D), len(data2D[0]) )
						self.raw    = [ data for k in range( self._shape[0] * self._shape[1] ) ]

				case [ tuple() as shape, data ]:
					assert len(shape)==2 and type(shape[0])==int and type(shape[1])==int
					self._shape = tuple(shape)
					if type(data) in ( tuple, list ):
						self.raw = list(data)
					else:
						self.raw = [ data for k in range( self._shape[0] * self._shape[1] ) ]

				case _:
					raise TypeError("Invalid argument type(s)!")

	@overload
	def __getitem__( self:Self, flat_index:int ) -> Tnum: ...
	@overload
	def __getitem__( self:Self, indices:tuple[int,int] ) -> Tnum: ...
	if _MATHUTIL_HAS_NP:
		@overload
		def __getitem__( self:Self, flat_indices:slice ) -> np.ndarray: ...
		#@overload
		#def __getitem__( self:Self, indices:tuple[slice,slice] ) -> np.ndarray: ...
		def __getitem__( self:Self, *args ):
			match args:
				case [ int()|slice() as flat_ind ]:
					return self.raw.flat[flat_ind]

				case [ (int() as j, int() as i) ]:
					assert j<self._shape[0] and i<self._shape[1]
					return self.raw.flat[ j*self._shape[1] + i ]

				case [ tuple() as indices ]:
					return self.raw[ *indices ]

			raise TypeError("Invalid argument type(s)!")
	else:
		@overload
		def __getitem__( self:Self, flat_indices:slice ) -> list[Tnum]: ...
		#@overload
		#def __getitem__( self:Self, indices:tuple[slice,slice] ) -> list[list[Tnum]]: ...
		def __getitem__( self:Self, *args ):
			match args:
				case [ int()|slice() as flat_ind ]:
					return self.raw[flat_ind]
				case [ (int() as j, int() as i) ]:
					assert j<self._shape[0] and i<self._shape[1]
					return self.raw[ j*self._shape[1] + i ]
				#case [ tuple() as indices ]:
				#	indices = args[0]
				#	if   type(indices[0]) == int:
				#		return self.raw[ indices[0]*self._shape[0] + indices[1] ]
				#	elif type(indices[1]) == slice:
				#		ret:list[list[Tnum]] = []
				#		slice_rows=indices[0]; slice_cols=indices[1]
				#		for     j in range( *slice_rows.indices(self._shape[0]) ):
				#			row:list[Tnum] = []
				#			for i in range( *slice_cols.indices(self._shape[1]) ):
				#				row.append(self.raw[ j*self._shape[1] + i ])
				#			ret.append(row)
				#		return ret
				#	# Fall out
			raise TypeError("Invalid argument type(s)!")
	#These return `None` for performance
	@overload
	def __setitem__( self:Self, flat_index:int, val:Tnum ) -> None: ...
	@overload
	def __setitem__( self:Self, indices:tuple[int,int], val:Tnum ) -> None: ...
	@overload
	def __setitem__( self:Self, flat_indices:slice, vals:Iterable[Tnum] ) -> None: ...
	#@overload
	#def __setitem__( self:Self, indices:tuple[slice,slice], val:Tnum ) -> None: ...
	#TODO: add support for (int,slice) and (slice,int)
	if _MATHUTIL_HAS_NP:
		def __setitem__( self:Self, *args ):
			match args:
				case [ int() as flat_ind, float()|int() as src ] | [ slice() as flat_ind, src ]:
					self.raw.flat[flat_ind] = src

				case [ (int() as j,int() as i), float()|int() as val ]:
					assert j<self._shape[0] and i<self._shape[1]
					self.raw.flat[ j*self._shape[1] + i ] = val

				case _:
					raise TypeError("Invalid argument type(s)!")
	else:
		def __setitem__( self:Self, *args ):
			match args:
				case [ int() as flat_ind, float()|int() as src ] | [ slice() as flat_ind, src ]:
					self.raw[flat_ind] = src

				case [ (int() as j,int() as i), float()|int() as val ]:
					assert j<self._shape[0] and i<self._shape[1]
					self.raw[ j*self._shape[1] + i ] = val

				#case [ tuple() as indices, float()|int() as val ]:
				#	match indices:
				#		case [ int() as j, int() as i ]:
				#			self.raw[ j*self._shape[1] + i ] = val
				#		case [ slice() as slice_rows, slice() as slice_cols ]:
				#			jread = 0
				#			for     jwrite in range( *slice_rows.indices(self._shape[0]) ):
				#				iread = 0
				#				for iwrite in range( *slice_cols.indices(self._shape[1]) ):
				#					self.raw[ jwrite*self._shape[0] + iwrite ] = vals[jread][iread]
				#					iread += 1
				#				assert iread == len(vals[jread])
				#				jread += 1
				#			assert jread == len(vals)
				#		case _:
				#			raise TypeError("Invalid argument type(s)!")

				case _:
					raise TypeError("Invalid argument type(s)!")

	if _MATHUTIL_HAS_NP:
		def __len__        ( self:Self ): return len(self.raw.flat)
		def __length_hint__( self:Self ): return len(self.raw.flat)

		def __iter__( self:Self ):
			#for datum in self.raw.flat: yield datum
			return iter(self.raw.flat)
		#def __reversed__( self:Self ):
		def __contains__( self:Self, item:Tnum ):
			return item in self.raw.flat
	else:
		def __len__        ( self:Self ): return len(self.raw)
		def __length_hint__( self:Self ): return len(self.raw)

		def __iter__( self:Self ):
			#for datum in self.raw: yield datum
			return iter(self.raw)
		#def __reversed__( self:Self ):
		def __contains__( self:Self, item:Tnum ):
			return item in self.raw

	if _MATHUTIL_HAS_NP:
		def __neg__( self:Self ):
			if len(self.raw.shape) == 1: return MatrixElems( self._shape, -self.raw )
			else                       : return MatrixElems(              -self.raw )
		#def __pos__( self:Self ):
		def __abs__( self:Self ):
			if len(self.raw.shape) == 1: return MatrixElems( self._shape, np.absolute(self.raw) )
			else                       : return MatrixElems(              np.absolute(self.raw) )
		def __invert__( self:Self ):
			if len(self.raw.shape) == 1: return MatrixElems( self._shape, np.bitwise_invert(self.raw) )
			else                       : return MatrixElems(              np.bitwise_invert(self.raw) )

		def __sub__( self:Self, other:Self ):
			assert isinstance(other,MatrixElems) and self._shape==other._shape
			ret = MatrixElems(self)
			ret.raw -= other.raw
			return ret
		def __add__( self:Self, other:Self ):
			assert isinstance(other,MatrixElems) and self._shape==other._shape
			ret = MatrixElems(self)
			ret.raw += other.raw
			return ret

		@overload
		def __truediv__( self:Self, divisor:Tnum ) -> Self: ...
		@overload
		def __truediv__( self:Self, divisor:Self ) -> Self: ...
		def __truediv__( self:Self, divisor:Tnum|Self ):
			ret = MatrixElems(self)
			match divisor:
				case float()|int():
					ret.raw /= divisor
				case MatrixElems():
					assert self._shape == divisor._shape
					ret.raw /= divisor.raw
				case _:
					raise TypeError("Invalid argument type!")
			return ret
		@overload
		def __floordiv__( self:Self, divisor:Tnum ) -> Self: ...
		@overload
		def __floordiv__( self:Self, divisor:Self ) -> Self: ...
		def __floordiv__( self:Self, divisor:Self ):
			ret = MatrixElems(self)
			match divisor:
				case float()|int():
					ret.raw //= divisor
				case MatrixElems():
					assert self._shape == divisor._shape
					ret.raw //= divisor.raw
				case _:
					raise TypeError("Invalid argument type!")
			return ret
		@overload
		def __mod__( self:Self, divisor:Tnum ) -> Self: ...
		@overload
		def __mod__( self:Self, divisor:Self ) -> Self: ...
		def __mod__( self:Self, divisor:Self ):
			ret = MatrixElems(self)
			match divisor:
				case float()|int():
					ret.raw %= divisor
				case MatrixElems():
					assert self._shape == divisor._shape
					ret.raw %= divisor.raw
				case _:
					raise TypeError("Invalid argument type!")
			return ret
		@overload
		def __mul__( self:Self, other:Tnum ) -> Self: ...
		@overload
		def __mul__( self:Self, other:Self ) -> Self: ...
		def __mul__( self:Self, other:Tnum|Self ):
			ret = MatrixElems(self)
			match other:
				case float()|int():
					ret.raw *= other
				case MatrixElems():
					assert self._shape == other._shape
					ret.raw *= other.raw
				case _:
					raise TypeError("Invalid argument type!")
			return ret
	else:
		def __neg__( self:Self ):
			return MatrixElems( self._shape, [ -val for val in self.raw ] )
		#def __pos__( self:Self ):
		def __abs__( self:Self ):
			return MatrixElems( self._shape, [ abs(val) for val in self.raw ] )
		def __invert__( self:Self ):
			return MatrixElems( self._shape, [ ~val for val in self.raw ] )

		def __sub__( self:Self, other:Self ):
			assert isinstance(other,MatrixElems) and self._shape==other._shape
			return MatrixElems( self._shape, [ a-b for a,b in zip(self.raw,other.raw) ] )
		def __add__( self:Self, other:Self ):
			assert isinstance(other,MatrixElems) and self._shape==other._shape
			return MatrixElems( self._shape, [ a+b for a,b in zip(self.raw,other.raw) ] )

		@overload
		def __truediv__( self:Self, divisor:Tnum ) -> Self: ...
		@overload
		def __truediv__( self:Self, divisor:Self ) -> Self: ...
		def __truediv__( self:Self, divisor:Tnum|Self ):
			match divisor:
				case float()|int():
					return MatrixElems( self._shape, [ a/divisor for a in self.raw ] )
				case MatrixElems():
					assert self._shape == divisor._shape
					return MatrixElems( self._shape, [ a/b for a,b in zip(self.raw,divisor.raw) ] )
				case _:
					raise TypeError("Invalid argument type!")
		@overload
		def __floordiv__( self:Self, divisor:Tnum ) -> Self: ...
		@overload
		def __floordiv__( self:Self, divisor:Self ) -> Self: ...
		def __floordiv__( self:Self, divisor:Self ):
			match divisor:
				case float()|int():
					return MatrixElems( self._shape, [ a//divisor for a in self.raw ] )
				case MatrixElems():
					assert self._shape == divisor._shape
					return MatrixElems( self._shape, [ a//b for a,b in zip(self.raw,divisor.raw) ] )
				case _:
					raise TypeError("Invalid argument type!")
		@overload
		def __mod__( self:Self, divisor:Tnum ) -> Self: ...
		@overload
		def __mod__( self:Self, divisor:Self ) -> Self: ...
		def __mod__( self:Self, divisor:Self ):
			match divisor:
				case float()|int():
					return MatrixElems( self._shape, [ a%divisor for a in self.raw ] )
				case MatrixElems():
					assert self._shape == divisor._shape
					return MatrixElems( self._shape, [ a%b for a,b in zip(self.raw,divisor.raw) ] )
				case _:
					raise TypeError("Invalid argument type!")
		@overload
		def __mul__( self:Self, other:Tnum ) -> Self: ...
		@overload
		def __mul__( self:Self, other:Self ) -> Self: ...
		def __mul__( self:Self, other:Tnum|Self ):
			match other:
				case float()|int():
					return MatrixElems( self._shape, [ a*other for a in self.raw ] )
				case MatrixElems():
					assert self._shape == other._shape
					return MatrixElems( self._shape, [ a*b for a,b in zip(self.raw,other.raw) ] )
				case _:
					raise TypeError("Invalid argument type!")
	def __rmul__( self:Self, other:Tnum|Self ):
		return self * other

	def __repr__( self:Self ):
		str_elems = ",".join(map( str, self.raw ))
		return f"MatrixElems(({self._shape[0]},{self._shape[1]}),[{str_elems}])"
	def __str__ ( self:Self ):
		contents = ""
		for j in range( self._shape[0] ):
			srow = ", ".join(map( str, self.raw ))
			contents += srow + "\n"
		return f"MatrElems{self._shape[0]}⨯{self._shape[1]}{{\n{contents}}}"
	#def __bytes__( self:Self ):



class Matrix( Generic[Tnum] ):
	__slots__ = ( "elems", )

	elems:MatrixElems[Tnum]

	@property
	def nrows( self:Self ): return self.elems._shape[0]
	@property
	def ncols( self:Self ): return self.elems._shape[1]
	@property
	def shape( self:Self ): return self.elems._shape
	@property
	def rows( self:Self ) -> list[list[Tnum]]: #TODO: list of vectors?
		nrows_, ncols_ = self.elems._shape
		rows:list[list[Tnum]] = []
		for j in range(nrows_):
			row = [ self[j,i] for i in range(ncols_) ]
			rows.append(row)
		return rows

	@overload
	def __init__( self:Self, other_matr:Self ): ...
	@overload
	def __init__( self:Self, other_elems:MatrixElems ): ...
	@overload
	def __init__( self:Self, shape:tuple[int,int] ): ...
	@overload
	def __init__( self:Self, shape:tuple[int,int], all_vals:Tnum ): ...
	@overload
	def __init__( self:Self, shape:tuple[int,int], rowmaj_data:Iterable[Tnum] ): ...
	@overload
	def __init__( self:Self, data:Iterable[Iterable[Tnum]] ): ...
	if _MATHUTIL_HAS_NP:
		@overload
		def __init__( self:Self, data2D:np.ndarray[Any,Tnum] ): ...
		@overload
		def __init__( self:Self, shape:tuple[int,int], data1D:np.ndarray[Any,Tnum] ): ...
	def __init__( self:Self, *args ):
		match args:
			case [ Matrix() as other_matr ]:
				self.elems = MatrixElems(other_matr.elems)
			case [ MatrixElems() as other_elems ]:
				self.elems = MatrixElems(other_elems)
			case _:
				self.elems = MatrixElems( *args )

	@staticmethod
	def augmented( matrA:Self, matrB:Self ) -> Self:
		assert matrA.elems._shape[0]==matrB.elems._shape[0], "Matrices must have the same number of rows!"
		data = []
		for j in range(matrA.elems._shape[0]):
			for i in range(matrA.elems._shape[1]): data.append( matrA[j,i] )
			for i in range(matrB.elems._shape[1]): data.append( matrB[j,i] )
		return Matrix( (matrA.elems._shape[0],matrA.elems._shape[1]+matrB.elems._shape[1]), data )

	#Note static construction methods here, set later because of needing to defined type annotations

	if _MATHUTIL_HAS_NP:
		@staticmethod
		def inverse( matr:Self ) -> Self:
			assert matr.elems._shape[0]==matr.elems._shape[1], "For simple inverse, matrix must be square!"
			return Matrix( np.linalg.inv(matr.elems.raw) )
	else:
		@staticmethod
		def inverse( matr:Self, eps_smallest_pivot:float=1.0e-9 ) -> Self:
			assert matr.elems._shape[0]==matr.elems._shape[1], "For simple inverse, matrix must be square!"

			aug = Matrix.augmented( matr, Matrix.get_identity(matr.elems._shape) )
			aug.to_rref(eps_smallest_pivot)

			ret = Matrix(matr.elems._shape)
			for     j in range(matr.elems._shape[0]):
				for i in range(matr.elems._shape[1]):
					ret[ j, i ] = aug[ j, i+4 ]

			return ret

	@overload
	def __getitem__( self:Self, flat_index:int ) -> Tnum: ...
	@overload
	def __getitem__( self:Self, indices:tuple[int,int] ) -> Tnum: ...
	if _MATHUTIL_HAS_NP:
		@overload
		def __getitem__( self:Self, indices:slice ) -> np.ndarray: ...
		@overload
		def __getitem__( self:Self, indices:tuple[slice,slice] ) -> np.ndarray: ...
	else:
		@overload
		def __getitem__( self:Self, indices:slice ) -> list[Tnum]: ...
		@overload
		def __getitem__( self:Self, indices:tuple[slice,slice] ) -> list[list[Tnum]]: ...
	def __getitem__( self:Self, *args ): return self.elems.__getitem__( *args )

	@overload
	def __setitem__( self:Self, flat_index:int, val:Tnum ) -> None: ...
	@overload
	def __setitem__( self:Self, indices:tuple[int,int], val:Tnum ) -> None: ...
	@overload
	def __setitem__( self:Self, indices:slice, vals:Iterable[Tnum] ) -> None: ...
	def __setitem__( self:Self, *args ): self.elems.__setitem__( *args )

	def __neg__( self:Self ):
		return self.__class__( -self.elems )
	#def __pos__( self:Self ):
	def __abs__( self:Self ):
		return self.__class__( abs(self.elems) )
	#Matrix transpose
	if _MATHUTIL_HAS_NP:
		def __invert__( self:Self ):
			assert self.elems._shape[0]!=1 or self.elems._shape[1]!=1 #should have been picked up by vector
			return self.__class__( self.elems.raw.T )
	else:
		def __invert__( self:Self ):
			nrows_, ncols_ = self.elems._shape
			result = self.__class__(( ncols_, nrows_ )) #TODO: generate data to pass in instead?
			for     j in range(nrows_):
				for i in range(ncols_):
					result[i,j] = self[j,i]
			return result

	def __sub__( self:Self, other:Self ):
		assert isinstance( other, Matrix ) #Shape checked in elems
		return self.__class__( self.elems - other.elems )
	def __add__( self:Self, other:Self ):
		assert isinstance( other, Matrix ) #Shape checked in elems
		return self.__class__( self.elems + other.elems )

	def __truediv__ ( self:Self, divisor:Tnum ):
		assert isinstance( divisor, (float,int) )
		return self.__class__( self.elems /  divisor )
	def __floordiv__( self:Self, divisor:Tnum ):
		assert isinstance( divisor, (float,int) )
		return self.__class__( self.elems // divisor )
	def __mod__     ( self:Self, divisor:Tnum ):
		assert isinstance( divisor, (float,int) )
		return self.__class__( self.elems %  divisor )
	@overload
	def __mul__( self:Self, other:Tnum ) -> Self: ...
	@overload
	def __mul__( self:Self, other:Self ) -> Self: ...
	if _MATHUTIL_HAS_NP:
		def __mul__( self:Self, other:Tnum|Self ) -> Self:
			match other:
				case float()|int():
					return Matrix( self.elems * other )
				case Matrix():
					assert not ( isinstance(self,Vector) or isinstance(other,Vector) )
					return Matrix( np.matmul( self.elems.raw, other.elems.raw ) )
			return NotImplemented
	else:
		def __mul__( self:Self, other:Tnum|Self ) -> Self:
			match other:
				case float()|int() as sc:
					return Matrix( self.elems * sc )
				case Matrix():
					assert not ( isinstance(self,Vector) or isinstance(other,Vector) )
					if self.elems._shape[1] != other.elems._shape[0]:
						raise ValueError("Matrix dimension mismatch!")
					nrows0=self.elems._shape[0]; shared=self.elems._shape[1]; ncols1=other.elems._shape[1]

					ret = Matrix(( nrows0, ncols1 )) #TODO: generate data to pass in instead?
					for     j in range(nrows0):
						for i in range(ncols1):
							ret[j,i] = sum([ self[j,k]*other[k,i] for k in range(shared) ])
					return ret
			return NotImplemented
	def __rmul__( self:Self, other:Tnum ) -> Self:
		assert isinstance( other, (float,int) ) #Should have been picked up by `__mul__` or subclass
		return Matrix( self.elems * other )

	@staticmethod
	def min( matr0:Self, matr1:Self ) -> Self:
		assert matr0.elems._shape==matr1.elems._shape and type(matr0)==type(matr1)
		dat = [ min(a,b) for a,b, in zip(matr0.elems,matr1.elems) ]
		return matr0.__class__( matr0.elems._shape, dat )
	@staticmethod
	def max( matr0:Self, matr1:Self ) -> Self:
		assert matr0.elems._shape==matr1.elems._shape and type(matr0)==type(matr1)
		dat = [ max(a,b) for a,b, in zip(matr0.elems,matr1.elems) ]
		return matr0.__class__( matr0.elems._shape, dat )

	def row_swap( self:Self, j0:int,j1:int ):
		for i in range(self.elems._shape[1]):
			self[j0,i], self[j1,i] = self[j1,i], self[j0,i]
	def row_sc( self:Self, j:int, sc:float ):
		for i in range(self.elems._shape[1]):
			self[j,i] *= sc
	def row_scadd( self:Self, jdst:int, scsrc:float,jsrc:int ):
		for i in range(self.elems._shape[1]):
			self[jdst,i] += scsrc * self[jsrc,i]

	#Compute row echelon form or reduced row echelon form.  See also:
	#	https://en.wikipedia.org/wiki/Gaussian_elimination#Pseudocode
	def to_ref( self:Self, eps_smallest_pivot:float=1.0e-9 ):
		h=0; k=0 #Pivot row and column
		while h<self.elems._shape[0] and k<self.elems._shape[1]:
			#Look for a good pivot as the largest magnitude value in the rest of the column.
			max_j:int
			max_absmag = 0.0
			for j in range( h, self.elems._shape[0], 1 ):
				absmag = abs( self[j,k] )
				assert not math.isnan(absmag), "Invalid matrix element!"
				if absmag > max_absmag:
					max_absmag=absmag; max_j=j

			if max_absmag > eps_smallest_pivot: #Found a good pivot.
				#Swap pivot-containing row up to row `h`.
				self.row_swap( h, max_j )

				#For all the rest of the rows, subtract off a factor `sc` of the pivot row that so
				#	that the leading term is canceled.
				pivot = self[h,k]
				for j in range( h+1, self.elems._shape[0], 1 ):
					sc = self[j,k] / pivot #leading term / pivot

					#self.row_scadd( j, -sc,h ) #would subtract zeros, wasteful
					self[j,k] = 0.0 #Do separately instead of in loop for numerical.
					for i in range( k+1, self.elems._shape[1], 1 ):
						self[j,i] -= sc * self[h,i]

				#Increase both pivot row and column.
				h += 1
			else: #No pivot in rest of column (all smaller in magnitude than `eps_smallest_pivot`).
				#Zero them to ensure it's exact.
				for j in range( h, self.elems._shape[0], 1 ): self[j,k] = 0.0
				#Increase pivot column only.
			k += 1
	def to_rref( self:Self, eps_smallest_pivot:float=1.0e-9 ):
		self.to_ref(eps_smallest_pivot)

		h=0; k=0 #Pivot row and column
		while h < self.elems._shape[0]:
			#Find `pivot`, the first nonzero element in current row `h`.
			pivot:Tnum
			while True:
				pivot = self[h,k]
				if pivot != 0.0: break
				k += 1
				if k >= self.elems._shape[1]: return

			#Subtract it from rows above.  Do this first, for better numerics.
			for j in range(h):
				first_affected_term = self[j,k]
				if first_affected_term == 0.0: continue

				sc = first_affected_term / pivot

				#self.row_scadd( j, -sc,h ) #would subtract zeros, wasteful
				self[j,k] = 0.0; #Do separately instead of in loop for numerical.
				for i in range( k+1, self.elems._shape[1], 1 ):
					self[j,i] -= sc * self[h,i]

			#Scale row so leading element is 1.
			#self.row_sc( h, 1.0/pivot ) #would scale zeros, wasteful (also numerical)
			self[h,k] = 1.0
			for i in range( k+1, self.elems._shape[1], 1 ): self[h,i]/=pivot

			h += 1

	def __repr__( self:Self ):
		str_elems = ",".join(map( str, self.elems.raw ))
		return f"Matrix(({self.elems._shape[0]},{self.elems._shape[1]}),[{str_elems}])"
	def __str__ ( self:Self ):
		contents = ""
		for j in range( self.elems._shape[0] ):
			srow = ", ".join([ "%.5g"%self[j,i] for i in range(self.elems._shape[1]) ])
			contents += "  "+srow+"\n"
		return f"""Matr{self.elems._shape[0]}⨯{self.elems._shape[1]}{{\n{contents}}}"""

	def print( self:Self, fmt:str="" ):
		lines = []
		for row in self.rows:
			line = " "
			for elem in row:
				line += ( fmt % elem )+" "
			lines.append(line)
		for j, line in enumerate(lines):
			if   j == 0:                      line = "⎡"+line+"⎤"
			elif j == self.elems._shape[0]-1: line = "⎣"+line+"⎦"
			else:                             line = "⎢"+line+"⎥"
			print(line)



class Vector( Matrix[Generic[Tnum]] ):
	__slots__ = ()

	@overload
	def __init__( self:Self, other_vec:Self ): ...
	@overload
	def __init__( self:Self, other_elems:MatrixElems ): ...
	@overload
	def __init__( self:Self, count:int ): ...
	@overload
	def __init__( self:Self, count:int, all_vals:Tnum ): ...
	@overload
	def __init__( self:Self,                       data:Iterable[Tnum] ): ...
	@overload
	def __init__( self:Self, shape:tuple[int,int], data:Iterable[Tnum] ): ...
	if _MATHUTIL_HAS_NP:
		@overload
		def __init__( self:Self,                       data:np.ndarray[Any,Tnum] ): ...
		@overload
		def __init__( self:Self, shape:tuple[int,int], data:np.ndarray[Any,Tnum] ): ...
	def __init__( self:Self, *args ):
		match args:
			case [ Vector() as other_vec ]:
				Matrix.__init__( self, other_vec )
			case [ MatrixElems() as other_elems ]:
				Matrix.__init__( self, other_elems )
			case [ int() as count ]:
				Matrix.__init__( self, (count,1) )
			case [ data ]:
				Matrix.__init__( self, (len(data),1), data )
			case [ (nrows,ncols), data ]:
				assert ( nrows==len(data) and ncols==1 ) or ( nrows==1 and ncols==len(data) )
				Matrix.__init__( self, (nrows,ncols), data )

	#Transpose (col-vec <-> row-vec)
	if _MATHUTIL_HAS_NP:
		def __invert__( self:Self ):
			nrows_, ncols_ = self.elems._shape
			return self.__class__( (ncols_,nrows_), self.elems.raw )
	else:
		def __invert__( self:Self ):
			result = self.__class__(self)
			nrows_, ncols_ = result.elems._shape
			result.elems._shape = ( ncols_, nrows_ )
			return result

	if _MATHUTIL_HAS_NP:
		def __mul__( self:Self, other:Tnum|Self|Matrix ) -> float|int|Self|Matrix:
			match other:
				case float()|int() as sc:
					return self.__class__( self.elems * sc )
				case Vector(): #Check first
					if self.elems._shape[0]==1 and other.elems._shape[1]==1:
						return np.vdot( self.elems.raw, other.elems.raw ) #return self | other
					else:
						return Matrix( np.outer( self.elems.raw, other.elems.raw ) ) #outer product
				case Matrix():
					if self.elems._shape[1] != other.elems._shape[0]:
						raise ValueError("Matrix dimension mismatch!")
					#Note `np.matmul(⋯)` works even though 1D
					return Vector._make( (1,self.elems._shape[1]),
						np.matmul( self.elems.raw, other.elems.raw )
					)
			return NotImplemented
	else:
		def __mul__( self:Self, other:Tnum|Self|Matrix ) -> float|int|Self|Matrix:
			match other:
				case float()|int() as sc:
					return self.__class__( self.elems * sc )
				case Vector(): #Check first
					if self.elems._shape[0]==1 and other.elems._shape[1]==1:
						return sum([ a*b for a,b in zip(self.elems.raw,other.elems.raw) ]) #return self | other
					else:
						return super().__mul__( self, other ) #outer product
				case Matrix():
					if self.elems._shape[1] != other.elems._shape[0]:
						raise ValueError("Matrix dimension mismatch!")
					return Vector._make( (1,self.elems._shape[1]), [
						sum([ self.elems.raw[k]*other[k,i] for k in range(other.elems._shape[0]) ])
						for i in range(other.elems._shape[1])
					] )
			return NotImplemented
	def __rmul__( self:Self, other:Tnum|Matrix ) -> Self:
		#Note Matrix * Vector will call Vector.__rmul__ instead of Matrix.__mul__, because Python
		#	thinks it's a feature to delegate to subclass.
		assert isinstance( other, ( float,int, Matrix ) )
		if type(other) == Matrix: #matrix * column vector = column vector
			if self.elems._shape[0] != other.elems._shape[1]:
				raise ValueError("Matrix dimension mismatch!")
			return Vector._make( (other.elems._shape[0],1), [
				sum([ other[j,k]*self.elems.raw[k] for k in range(self.elems._shape[0]) ])
				for j in range(other.elems._shape[0])
			] )
		else:
			return self.__class__( self.elems * other )

	if _MATHUTIL_HAS_NP:
		#Dot product
		def __or__( self:Self, other:Self ) -> float|int:
			assert isinstance(other,Vector) and self.elems._shape==other.elems._shape
			return np.vdot( self.elems.raw, other.elems.raw )

		def get_length_sq( self:Self ) -> float|int:
			#return self | self
			return np.vdot( self.elems.raw, self.elems.raw )
	else:
		#Dot product
		def __or__( self:Self, other:Self ) -> float|int:
			assert isinstance(other,Vector) and self.elems._shape==other.elems._shape
			#return sum( self.elems * other.elems )
			return sum([ a*b for a,b in zip(self.elems.raw,other.elems.raw) ])

		def get_length_sq( self:Self ) -> float|int:
			#return self | self
			return sum([ a*a for a in self.elems.raw ])
	def get_length( self:Self ):
		return math.sqrt( self.get_length_sq() )

	def normalized( self:Self, default_if_zero:Self|None=None ) -> Self:
		l = self.get_length()
		if l == 0.0:
			if default_if_zero == None:
				raise ValueError("Vector has zero length!")
			return default_if_zero
		return self / l

	def __repr__( self:Self ):
		str_elems = ",".join(map( str, self.elems.raw ))
		return f"{self.__class__.__name__}([{str_elems}])"
	def __str__ ( self:Self ):
		str_elems = ",".join(map( str, self.elems.raw ))
		return "⟨"+str_elems+"⟩"
	#def __bytes__( self:Self ):

def _vec_swizzler( Tret, *inds ):
	def getter( self:Vector, Tret:type, inds ):
		return Tret([ self[ind] for ind in inds ])
	def setter( self:Vector,            inds, val:Self ):
		for k in range(len(inds)):
			self[k] = val[ inds[k] ]
	if len(set(inds)) == len(inds):
		return property(
			fget = lambda self,    Tret=Tret,inds=inds: getter(self,Tret,inds    ),
			fset = lambda self,val,          inds=inds: setter(self,     inds,val)
		)
	else:
		return property(
			fget = lambda self,Tret=Tret,inds=inds: getter(self,Tret,inds)
		)

class Vector2( Vector[Generic[Tnum]] ):
	__slots__ = ()

	def _x_get( self:Self ): return self.elems.raw[0]
	def _x_set( self:Self, val:Tnum ): self.elems.raw[0] = val
	x = property( fget=_x_get, fset=_x_set )

	def _y_get( self:Self ): return self.elems.raw[1]
	def _y_set( self:Self, val:Tnum ): self.elems.raw[1] = val
	y = property( fget=_y_get, fset=_y_set )

	@overload
	def __init__( self:Self, other_vec:Self ): ...
	@overload
	def __init__( self:Self, other_elems:MatrixElems ): ...
	@overload
	def __init__( self:Self, all_vals:Tnum ): ...
	@overload
	def __init__( self:Self,                       data1D:Iterable[Tnum] ): ...
	@overload
	def __init__( self:Self, shape:tuple[int,int], data1D:Iterable[Tnum] ): ...
	@overload
	def __init__( self:Self, x:Tnum,y:Tnum ): ...
	if _MATHUTIL_HAS_NP:
		@overload
		def __init__( self:Self, data:np.ndarray[Any,Tnum] ): ...
	def __init__( self:Self, *args ):
		match args:
			case [ Vector() as other_vec ]:
				assert len(other_vec.elems.raw) == 2
				Vector.__init__( self, other_vec )
			case [ MatrixElems() as other_elems ]:
				assert len(other_elems.raw) == 2
				Vector.__init__( self, other_elems )
			case [ float()|int() as all_vals ]:
				Vector.__init__( self, [all_vals,all_vals] )
			case [ data1D ]:
				Vector.__init__( self, data1D )
			case [ tuple() as shape, data1D ]:
				Vector.__init__( self, shape, data1D )
			case [ x, y ]:
				Vector.__init__( self, [x,y] )
			case _:
				raise TypeError("Invalid argument type(s)!")

	#https://en.wikipedia.org/wiki/Rotation_matrix#Basic_3D_rotations
	def rotated( self:Self, angle_rad:float ):
		s=math.sin(angle_rad); c=math.cos(angle_rad)
		return Vector2( self.x*c-self.y*s, self.x*s+self.y*c )

for i0, i1 in itertools.permutations(( 0, 1 )):
	setattr( Vector2, "xy"[i0]+"xy"[i1], _vec_swizzler(Vector2,i0,i1) )

class Vector3( Vector[Generic[Tnum]] ):
	__slots__ = ()

	def _x_get( self:Self ): return self.elems.raw[0]
	def _x_set( self:Self, val:Tnum ): self.elems.raw[0] = val
	x = property( fget=_x_get, fset=_x_set )

	def _y_get( self:Self ): return self.elems.raw[1]
	def _y_set( self:Self, val:Tnum ): self.elems.raw[1] = val
	y = property( fget=_y_get, fset=_y_set )

	def _z_get( self:Self ): return self.elems.raw[2]
	def _z_set( self:Self, val:Tnum ): self.elems.raw[2] = val
	z = property( fget=_z_get, fset=_z_set )

	@overload
	def __init__( self:Self, other_vec:Self ): ...
	@overload
	def __init__( self:Self, other_elems:MatrixElems ): ...
	@overload
	def __init__( self:Self, all_vals:Tnum ): ...
	@overload
	def __init__( self:Self,                       data1D:Iterable[Tnum] ): ...
	@overload
	def __init__( self:Self, shape:tuple[int,int], data1D:Iterable[Tnum] ): ...
	@overload
	def __init__( self:Self, x:Tnum,y:Tnum,z:Tnum ): ...
	@overload
	def __init__( self:Self, xy:Vector2,z:Tnum ): ...
	@overload
	def __init__( self:Self, x:Tnum,yz:Vector2 ): ...
	if _MATHUTIL_HAS_NP:
		@overload
		def __init__( self:Self, data:np.ndarray[Any,Tnum] ): ...
	def __init__( self:Self, *args ):
		match args:
			case [ Vector() as other_vec ]:
				assert len(other_vec.elems.raw) == 3
				Vector.__init__( self, other_vec )
			case [ MatrixElems() as other_elems ]:
				assert len(other_elems.raw) == 3
				Vector.__init__( self, other_elems )
			case [ float()|int() as all_vals ]:
				Vector.__init__( self, [all_vals,all_vals,all_vals] )
			case [ data1D ]:
				Vector.__init__( self, data1D )
			case [ tuple() as shape, data1D ]:
				Vector.__init__( self, shape, data1D )
			case [ Vector2() as xy, z ]:
				Vector.__init__( self, [xy.elems.raw[0],xy.elems.raw[1],z] )
			case [ x, Vector2() as yz ]:
				Vector.__init__( self, [x,yz.elems.raw[0],yz.elems.raw[1]] )
			case [ x, y, z ]:
				Vector.__init__( self, [x,y,z] )
			case _:
				raise TypeError("Invalid argument type(s)!")

	#Cross product
	def __xor__( self:Self, other:Self ) -> Self:
		assert isinstance( other, Vector3 )

		#if _MATHUTIL_HAS_NP:
		#	return Vector3( np.cross( self.elems.raw, other.elems.raw ) )

		#tmp0 = ( self.y, self.z, self.x ) #self.yzx
		#tmp1 = ( tmp0[0]*other.elems.raw[0], tmp0[1]*other.elems.raw[1], tmp0[2]*other.elems.raw[2] ) #tmp0 * other
		#tmp2 = ( tmp0[0]*other.elems.raw[2], tmp0[1]*other.elems.raw[0], tmp0[2]*other.elems.raw[1] ) #tmp0 * other.zxy
		#return Vector3( tmp2[0]-tmp1[1], tmp2[1]-tmp1[2], tmp2[2]-tmp1[0] ) #tmp2 - tmp1.yzx

		u0, u1, u2 = self .elems.raw
		v0, v1, v2 = other.elems.raw
		return Vector3( u1*v2-u2*v1, u2*v0-u0*v2, u0*v1-u1*v0 )

	#https://en.wikipedia.org/wiki/Rotation_matrix#Basic_3D_rotations
	def rotated_z( self:Self, angle_rad:float ):
		s=math.sin(angle_rad); c=math.cos(angle_rad)
		return Vector3( self.x*c-self.y*s, self.x*s+self.y*c, self.z )

for i0, i1 in itertools.combinations( (0,1,2), 2 ):
	setattr( Vector3, "xyz"[i0]+"xyz"[i1], _vec_swizzler(Vector2,i0,i1) )
for i0, i1, i2 in itertools.permutations(( 0, 1, 2 )):
	setattr( Vector3, "xyz"[i0]+"xyz"[i1]+"xyz"[i2], _vec_swizzler(Vector3,i0,i1,i2) )

class Vector4( Vector[Generic[Tnum]] ):
	__slots__ = ()

	def _x_get( self:Self ): return self.elems.raw[0]
	def _x_set( self:Self, val:Tnum ): self.elems.raw[0] = val
	x = property( fget=_x_get, fset=_x_set )

	def _y_get( self:Self ): return self.elems.raw[1]
	def _y_set( self:Self, val:Tnum ): self.elems.raw[1] = val
	y = property( fget=_y_get, fset=_y_set )

	def _z_get( self:Self ): return self.elems.raw[2]
	def _z_set( self:Self, val:Tnum ): self.elems.raw[2] = val
	z = property( fget=_z_get, fset=_z_set )

	def _w_get( self:Self ): return self.elems.raw[3]
	def _w_set( self:Self, val:Tnum ): self.elems.raw[3] = val
	w = property( fget=_w_get, fset=_w_set )

	@overload
	def __init__( self:Self, other_vec:Self ): ...
	@overload
	def __init__( self:Self, other_elems:MatrixElems ): ...
	@overload
	def __init__( self:Self, all_vals:Tnum ): ...
	@overload
	def __init__( self:Self,                       data1D:Iterable[Tnum] ): ...
	@overload
	def __init__( self:Self, shape:tuple[int,int], data1D:Iterable[Tnum] ): ...
	@overload
	def __init__( self:Self, x:Tnum,y:Tnum,z:Tnum,w:Tnum ): ...
	@overload
	def __init__( self:Self, xy:Vector2,z:Tnum,w:Tnum ): ...
	@overload
	def __init__( self:Self, x:Tnum,yz:Vector2,w:Tnum ): ...
	@overload
	def __init__( self:Self, x:Tnum,y:Tnum,zw:Vector2 ): ...
	@overload
	def __init__( self:Self, xy:Vector2,zw:Vector2 ): ...
	@overload
	def __init__( self:Self, x:Tnum,yzw:Vector3 ): ...
	@overload
	def __init__( self:Self, xyz:Vector3,w:Tnum ): ...
	if _MATHUTIL_HAS_NP:
		@overload
		def __init__( self:Self, data:np.ndarray[Any,Tnum] ): ...
	def __init__( self:Self, *args ):
		match args:
			case [ Vector() as other_vec ]:
				assert len(other_vec.elems.raw) == 4
				Vector.__init__( self, other_vec )
			case [ MatrixElems() as other_elems ]:
				assert len(other_elems.raw) == 4
				Vector.__init__( self, other_elems )
			case [ float()|int() as all_vals ]:
				Vector.__init__( self, [all_vals,all_vals,all_vals,all_vals] )
			case [ data1D ]:
				Vector.__init__( self, data1D )
			case [ tuple() as shape, data1D ]:
				Vector.__init__( self, shape, data1D )
			case [ Vector2() as xy, z, w ]:
				Vector.__init__( self, [xy.elems.raw[0],xy.elems.raw[1],z,w] )
			case [ x, Vector2() as yz, w ]:
				Vector.__init__( self, [x,yz.elems.raw[0],yz.elems.raw[1],w] )
			case [ x, y, Vector2() as zw ]:
				Vector.__init__( self, [x,y,zw.elems.raw[0],zw.elems.raw[1]] )
			case [ Vector2() as xy, Vector2() as zw ]:
				Vector.__init__( self, [xy.elems.raw[0],xy.elems.raw[1],zw.elems.raw[0],zw.elems.raw[1]] )
			case [ Vector3() as xyz, w ]:
				Vector.__init__( self, [xyz.elems.raw[0],xyz.elems.raw[1],xyz.elems.raw[2],w] )
			case [ x, Vector3() as yzw ]:
				Vector.__init__( self, [x,yzw.elems.raw[0],yzw.elems.raw[1],yzw.elems.raw[2]] )
			case [ x, y, z, w ]:
				Vector.__init__( self, [x,y,z,w] )
			case _:
				raise TypeError("Invalid argument type(s)!")

for i0, i1 in itertools.combinations( (0,1,2,3), 2 ):
	setattr( Vector4, "xyzw"[i0]+"xyzw"[i1], _vec_swizzler(Vector2,i0,i1) )
for i0, i1, i2 in itertools.combinations( (0,1,2,3), 3 ):
	setattr( Vector4, "xyzw"[i0]+"xyzw"[i1]+"xyzw"[i2], _vec_swizzler(Vector3,i0,i1,i2) )
for i0, i1, i2, i3 in itertools.permutations(( 0, 1, 2, 3 )):
	setattr( Vector4, "xyzw"[i0]+"xyzw"[i1]+"xyzw"[i2]+"xyzw"[i3], _vec_swizzler(Vector4,i0,i1,i2,i3) )

def _vec_make( shape:tuple[int,int], data1D:Iterable[Tnum] ):
	l = len(data1D)
	if   l == 2: vec = Vector2(data1D)
	elif l == 3: vec = Vector3(data1D)
	elif l == 4: vec = Vector4(data1D)
	else       : vec = Vector (data1D)
	return vec
Vector._make = _vec_make



if _MATHUTIL_HAS_NP:
	def _matr_get_empty( shape:tuple[int,int] ):
		return Matrix( np.empty(shape) )
	def _matr_get_zeros( shape:tuple[int,int] ):
		return Matrix( np.zeros(shape) )
	def _matr_get_ones( shape:tuple[int,int] ):
		return Matrix( np.ones(shape) )
	def _matr_get_identity( shape:tuple[int,int] ): #TODO: overload with single res
		return Matrix( np.eye( shape[0], shape[1] ) )
else:
	def _matr_get_empty( shape:tuple[int,int] ):
		return Matrix( shape, NaN )
	def _matr_get_zeros( shape:tuple[int,int] ):
		return Matrix( shape, 0.0 )
	def _matr_get_ones( shape:tuple[int,int] ):
		return Matrix( shape, 0.0 )
	def _matr_get_identity( shape:tuple[int,int] ):
		ret = _matr_get_zeros(shape)
		for k in range(min( shape[1], shape[0] )):
			ret[k,k] = 1.0
		return ret

@overload
def _matr_get_rotation_deg( angle_deg:float, normed_axis:Vector3 ) -> Matrix: ...
@overload
def _matr_get_rotation_deg( angle_deg:float, rx:float,ry:float,rz:float ) -> Matrix: ...
def _matr_get_rotation_deg( angle_deg:float, *args ):
	return _matr_get_rotation_rad( math.radians(angle_deg), *args )
@overload
def _matr_get_rotation_rad( angle_rad:float, normed_axis:Vector3 ) -> Matrix: ...
@overload
def _matr_get_rotation_rad( angle_rad:float, rx:float,ry:float,rz:float ) -> Matrix: ...
def _matr_get_rotation_rad( angle_rad:float, *args ):
	#https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
	#https://registry.khronos.org/OpenGL-Refpages/gl2.1/xhtml/glRotate.xml
	match args:
		case [ x, y, z ]: pass
		case [ Vector3() as normed_axis ]:
			x, y, z = normed_axis
		case _:
			raise TypeError("Invalid argument type(s)!")

	s=math.sin(angle_rad); c=math.cos(angle_rad)
	C = 1 - c

	xs =x* s; ys =y* s; zs =z* s
	xC =x* C; yC =y* C; zC =z* C
	xyC=x*yC; yzC=y*zC; zxC=z*xC

	return Matrix( (4,4), [
		x*xC+c, xyC-zs, zxC+ys, 0,
		xyC+zs, y*yC+c, yzC-xs, 0,
		zxC-ys, yzC+xs, z*zC+c, 0,
		0     , 0     , 0     , 1
	] )
@overload
def _matr_get_scaling( scaling:Vector2 ) -> Matrix: ...
@overload
def _matr_get_scaling( scaling:Vector3 ) -> Matrix: ...
@overload
def _matr_get_scaling( sx:float, sy:float ) -> Matrix: ...
@overload
def _matr_get_scaling( sx:float, sy:float, sz:float ) -> Matrix: ...
def _matr_get_scaling( *args ):
	match args:
		case [ float()|int() as sx, float()|int() as sy, float()|int() as sz ]:
			ret = _matr_get_identity( (4,4) )
			ret[0,0]=sx; ret[1,1]=sy; ret[2,2]=sz
		case [ Vector3() as vec ]: return _matr_get_scaling( vec[0], vec[1], vec[2] )

		case [ float()|int() as sx, float()|int() as sy                      ]:
			ret = _matr_get_identity( (3,3) )
			ret[0,0]=sx; ret[1,1]=sy
		case [ Vector2() as vec ]: return _matr_get_scaling( vec[0], vec[1]         )

		case _: raise TypeError("Invalid argument type(s)!")

	return ret
@overload
def _matr_get_translation( translation:Vector2 ) -> Matrix: ...
@overload
def _matr_get_translation( translation:Vector3 ) -> Matrix: ...
@overload
def _matr_get_translation( tx:float ) -> Matrix: ...
@overload
def _matr_get_translation( tx:float, ty:float ) -> Matrix: ...
@overload
def _matr_get_translation( tx:float, ty:float, tz:float ) -> Matrix: ...
def _matr_get_translation( *args ):
	match args:
		case [ float()|int() as tx, float()|int() as ty, float()|int() as tz ]:
			return Matrix( (4,4), [ 1.0,0.0,0.0,tx, 0.0,1.0,0.0,ty, 0.0,0.0,1.0,tz, 0.0,0.0,0.0,1.0 ] )
		case [ Vector3() as vec ]:
			return _matr_get_translation( vec.elems.raw[0], vec.elems.raw[1], vec.elems.raw[2] )

		case [ float()|int() as tx, float()|int() as ty ]:
			return Matrix( (3,3), [ 1.0,0.0,tx, 0.0,1.0,ty, 0.0,0.0,1.0 ] )
		case [ Vector2() as vec ]:
			return _matr_get_translation( vec.elems.raw[0], vec.elems.raw[1] )

		case [ float()|int() as tx ]:
			return Matrix( (2,2), [ 1.0,tx, 0.0,1.0 ] )

		case _: raise TypeError("Invalid argument type(s)!")
def _matr_get_view_lookat( pos_eye:Vector3, pos_center:Vector3, up_vec:Vector3=Vector3(0,0,1) ):
	#Note: the spec. is probably wrong: the `s` vector needs to be normalized when it is used
	#	inside the matrix, also.  Note normalize checks for zero length, because this can
	#	happen, and it is an error.
	f = ( pos_center - pos_eye ).normalized()
	s = ( f ^ up_vec           ).normalized()
	u = s ^ f

	return Matrix( (4,4), [
		 s[0],  s[1],  s[2], 0,
		 u[0],  u[1],  u[2], 0,
		-f[0], -f[1], -f[2], 0,
		 0   ,  0    , 0   , 1
	] ) * _matr_get_translation( -pos_eye )
def _matr_get_proj_frustum(
	left:float,right:float, bottom:float,top:float, clip_near:float,clip_far:float, parallel:bool
):
	l=left;r=right; b=bottom;t=top; n=clip_near;f=clip_far
	dx=r-l; dy=t-b; dz=f-n
	if parallel:
		#`glOrtho( l,r, b,t, n,f );`
		return Matrix( (4,4), [
			2/dx, 0   ,  0   , -(r+l)/dx,
			0   , 2/dy,  0   , -(t+b)/dy,
			0   , 0   , -2/dz, -(f+n)/dz,
			0   , 0   ,  0   , 1
		] )
	else:
		#`glFrustum( l,r, b,t, n,f );`
		#http://nehe.gamedev.net/article/replacement_for_gluperspective/21002/
		return Matrix( (4,4), [
			2*n/dx, 0       ,    (r+l)/dx,  0       ,
			0     ,   2*n/dy,    (t+b)/dy,  0       ,
			0     ,   0     ,   -(f+n)/dz, -2*f*n/dz,
			0     ,   0     ,   -1       ,  0
		] )
def _matr_get_proj_ortho(
	left:float,right:float, bottom:float,top:float, clip_near:float,clip_far:float
):
	return _matr_get_proj_frustum( left,right, bottom,top, clip_near,clip_far, True )
def _matr_get_proj_persp_rad( vertfov_rad:float, aspect_ratio:float, clip_near:float,clip_far:float ):
	rad_h = math.tan(0.5*vertfov_rad) * clip_near
	rad_w = aspect_ratio * rad_h
	return _matr_get_proj_frustum( -rad_w,rad_w, -rad_h,rad_h, clip_near,clip_far, False )
def _matr_get_proj_persp_deg( vertfov_deg:float, aspect_ratio:float, clip_near:float,clip_far:float ):
	return _matr_get_proj_persp_rad( math.radians(vertfov_deg), aspect_ratio, clip_near,clip_far )

if _MATHUTIL_HAS_NP:
	def _matr_scale1_ip( lhs_matr:Matrix, scale:float ):
		assert lhs_matr.shape == ( 4, 4 )
		for k in range(12):
			if k & 0b11 == 3: continue
			lhs_matr.elems.raw.flat[k] *= scale
	def _matr_translate_ip( lhs_matr:Matrix, translation:Vector3 ):
		assert lhs_matr.shape == ( 4, 4 )
		lhs_matr.elems.raw.flat[ 3] += translation.elems.raw[0]
		lhs_matr.elems.raw.flat[ 7] += translation.elems.raw[1]
		lhs_matr.elems.raw.flat[11] += translation.elems.raw[2]
else:
	def _matr_scale1_ip( lhs_matr:Matrix, scale:float ):
		assert lhs_matr.shape == ( 4, 4 )
		for k in range(12):
			if k & 0b11 == 3: continue
			lhs_matr.elems[k] *= scale
	def _matr_translate_ip( lhs_matr:Matrix, translation:Vector3 ):
		assert lhs_matr.shape == ( 4, 4 )
		lhs_matr.elems[ 3] += translation.elems.raw[0]
		lhs_matr.elems[ 7] += translation.elems.raw[1]
		lhs_matr.elems[11] += translation.elems.raw[2]

Matrix.get_empty    = _matr_get_empty
Matrix.get_zeros    = _matr_get_zeros
Matrix.get_ones     = _matr_get_ones
Matrix.get_identity = _matr_get_identity

Matrix.get_rotation_deg   = _matr_get_rotation_deg
Matrix.get_rotation_rad   = _matr_get_rotation_rad
Matrix.get_scaling        = _matr_get_scaling
Matrix.get_translation    = _matr_get_translation
Matrix.get_view_lookat    = _matr_get_view_lookat
Matrix.get_proj_frustum   = _matr_get_proj_frustum
Matrix.get_proj_ortho     = _matr_get_proj_ortho
Matrix.get_proj_persp_rad = _matr_get_proj_persp_rad
Matrix.get_proj_persp_deg = _matr_get_proj_persp_deg

Matrix.scale1_ip    = _matr_scale1_ip
Matrix.translate_ip = _matr_translate_ip



class Quaternion:
	__slots__ = ( "abcd", )

	def _get_a( self:Self            ): return self.abcd.elems.raw[0]
	def _get_b( self:Self            ): return self.abcd.elems.raw[1]
	def _get_c( self:Self            ): return self.abcd.elems.raw[2]
	def _get_d( self:Self            ): return self.abcd.elems.raw[3]
	def _set_a( self:Self, val:float ): self.abcd.elems.raw[0] = val
	def _set_b( self:Self, val:float ): self.abcd.elems.raw[1] = val
	def _set_c( self:Self, val:float ): self.abcd.elems.raw[2] = val
	def _set_d( self:Self, val:float ): self.abcd.elems.raw[3] = val
	a = property( fget=_get_a, fset=_set_a )
	b = property( fget=_get_b, fset=_set_b )
	c = property( fget=_get_c, fset=_set_c )
	d = property( fget=_get_d, fset=_set_d )

	def _get_angle( self:Self                  ): return self.abcd.elems.raw[0]
	def _set_angle( self:Self, angle_rad:float ): self.abcd.elems.raw[0] = angle_rad
	angle = property( fget=_get_angle, fset=_set_angle )

	def _get_axis( self:Self               ): self.abcd.yzw
	def _set_axis( self:Self, axis:Vector3 ): self.abcd.yzw = axis
	axis = property( fget=_get_axis, fset=_set_axis )

	#https://danceswithcode.net/engineeringnotes/quaternions/quaternions.html
	@overload
	def __init__( self:Self, a:float,b:float,c:float,d:float ): ...
	@overload
	def __init__( self:Self, vec:Vector4 ): ...
	@overload
	def __init__( self:Self, axis:Vector3, angle_rad:float ): ...
	#@overload
	#def __init__( self:Self, other:matr3x3f ): ...
	@overload
	def __init__( self:Self, other:Self ): ...
	def __init__( self:Self, *args ):
		match args:
			case [ float() as a, float() as b, float() as c, float() as d ]:
				self.abcd = Vector4( a, b, c, d )
			case [ Vector3() as axis, float() as angle ]:
				#https://wiki.ogre3d.org/Quaternion+and+Rotation+Primer
				#⟨ w, x, y, z ⟩ = ⟨ cos(a/2), sin(a/2)*nx, sin(a/2)*ny, sin(a/2)*nz ⟩
				theta = 0.5 * angle
				s = math.sin(theta)
				c = math.cos(theta)
				svec = s * axis
				self.abcd = Vector4( c, svec )
			case [ Vector4() as vec ]:
				self.abcd = Vector4(vec)
			case [ Quaternion() as other ]:
				self.abcd = Vector4(other.abcd)
			case _:
				raise TypeError("Invalid argument type(s)!")

	def __neg__( self:Self ):
		return Quaternion( -self.abcd )

	def __sub__( self:Self, other:Self ):
		return Quaternion( self.abcd - other.abcd )
	def __add__( self:Self, other:Self ):
		return Quaternion( self.abcd + other.abcd )

	def __mul__( self:Self, other:float|int|Self ) -> Self:
		match other:
			case float()|int() as sc:
				return Quaternion( self.abcd * sc )

			case Quaternion() as other:
				#(a+bi+cj+dk)*(e+fi+gj+hk) = (ae-bf-cg-dh) + (af+be+ch-dg)i + (ag-bh+ce+df)j + (ah+bg-cf+de)k

				#term0 = self.a * other.abcd
				#term1 = self.b * Vector4( -other.b,  other.a, -other.d,  other.c )
				#term2 = self.c * Vector4( -other.c,  other.d,  other.a, -other.b )
				#term3 = self.d * Vector4( -other.d, -other.c,  other.b,  other.a )
				#return Quaternion( term0 + term1 + term2 + term3 )

				a, b, c, d = self .abcd.elems.raw
				e, f, g, h = other.abcd.elems.raw
				return Quaternion(
					a*e - b*f - c*g - d*h,
					a*f + b*e + c*h - d*g,
					a*g - b*h + c*e + d*f,
					a*h + b*g - c*f + d*e
				)

			case _:
				raise TypeError("Invalid argument type!")
	def __rmul__( self:Self, other:float|int ) -> Self:
		return self * other

	def normalized( self:Self ):
		return Quaternion( self.abcd.normalized() )

	def calc_matr3x3( self:Self ) -> Matrix:
		#https://en.wikipedia.org/wiki/Rotation_matrix#Quaternion

		Nq = self.abcd.get_length_sq()
		s = 2.0/Nq if Nq>0.0 else 0.0

		B =self.b*s; C =self.c*s; D =self.d*s
		aB=self.a*B; aC=self.a*C; aD=self.a*D
		bB=self.b*B; bC=self.b*C; bD=self.b*D
		cC=self.c*C; cD=self.c*D; dD=self.d*D
		return Matrix( (3,3), [
			1.0-cC-dD,     bC-aD,     bD+aC,
			    bC+aD, 1.0-bB-dD,     cD-aB,
			    bD-aC,     cD+aB, 1.0-bB-cC
		] )
	def calc_matr4x4( self:Self ) -> Matrix:
		tmp = self.calc_matr3x3()
		return Matrix( (4,4), [
			tmp[0], tmp[1], tmp[2], 0.0,
			tmp[3], tmp[4], tmp[5], 0.0,
			tmp[6], tmp[7], tmp[8], 0.0,
			0.0   , 0.0   , 0.0   , 1.0
		] )

	def __repr__( self:Self ):
		return f"Quaternion({self.a},{self.b},{self.c},{self.d})"
	def __str__ ( self:Self ):
		return f"{self.a}+{self.b}i+{self.c}j+{self.d}k"



def rndint( val:float ): return int(round(val))

def clamp( val:float|int, lo:float|int,hi:float|int ):
	if val < lo: return lo
	if val > hi: return hi
	return val

def lerp( val0:Any,val1:Any, t:float ):
	return val0*(1.0-t) + val1*t



def parabola_extremum_x  ( coeffs:tuple[float,float,float], eps:float=1e-18 ):
	a, b, c = coeffs
	if   abs(a) > eps:
		return -0.5 * b / a
	elif abs(b) < eps:
		return 0.0
	else:
		return NaN
def parabola_extremum_val( coeffs:tuple[float,float,float], eps:float=1e-18 ):
	a, b, c = coeffs
	if   abs(a) > eps:
		return c - 0.25*b*b/a
	elif abs(b) < eps:
		return c
	else:
		return NaN



def random_sphere():
	x = random.gauss()
	y = random.gauss()
	z = random.gauss()
	return Vector3( x, y, z )



def solve_2x2( A:list[float], b:list[float] ):
	assert len(A)==4 and len(b)==2

	rref = [
		Vector3( A[0],A[1], b[0] ),
		Vector3( A[2],A[3], b[1] ),
	]

	rref[0] /= rref[0].x
	rref[1] -= rref[1].x * rref[0]

	rref[1] /= rref[1].y
	rref[0] -= rref[0].y * rref[1]

	return ( rref[0].z, rref[1].z )
def solve_3x3( A:list[float], b:list[float] ):
	assert len(A)==9 and len(b)==3

	rref = [
		Vector4( A[0],A[1],A[2], b[0] ),
		Vector4( A[3],A[4],A[5], b[1] ),
		Vector4( A[6],A[7],A[8], b[2] )
	]

	rref[0] /= rref[0].x
	rref[1] -= rref[1].x * rref[0]
	rref[2] -= rref[2].x * rref[0]

	rref[1] /= rref[1].y
	rref[0] -= rref[0].y * rref[1]
	rref[2] -= rref[2].y * rref[1]

	rref[2] /= rref[2].z
	rref[0] -= rref[0].z * rref[2]
	rref[1] -= rref[1].z * rref[2]

	return ( rref[0].w, rref[1].w, rref[2].w )



def triangle_normal( v0:Vector3, v1:Vector3, v2:Vector3 ):
	e01 = v1 - v0
	e02 = v2 - v0
	return ( e01 ^ e02 ).normalized()



@overload
def compute_frustrum_bound( apex:Vector3, sph_cen:Vector3,sph_rad:float ) -> tuple[float,float,float]: ...
@overload
def compute_frustrum_bound( apex:Vector3, aabb_lo:Vector3,aabb_hi:Vector3 ) -> tuple[float,float,float]: ...
def compute_frustrum_bound( apex:Vector3, *args ):
	match args:
		case [ Vector3() as sph_cen, float() as sph_rad ]:
			#Computes a frustum from apex that tightly bounds the given sphere (returns apex angle
			#	and near and far distances).
			vec = sph_cen - apex
			d = vec.get_length()
			apex_angle_rad = 2.0*math.asin( sph_rad / d )
			return apex,sph_cen, apex_angle_rad, d-sph_rad,d+sph_rad

		case [ Vector3() as aabb_lo, Vector3() as aabb_hi ]:
			centroid = 0.5*( aabb_lo + aabb_hi )
			rad = ( aabb_hi - centroid ).get_length()
			return compute_frustrum_bound( apex, centroid,rad )

		case _:
			raise TypeError("Invalid argument type(s)!")
