local z={}local U=math.abs local Q=math.floor local M=math.sqrt local v=math.cos local n=math.sin local h={}local e=function(z)local U=setmetatable({},h)rawset(U,"\095\118\097\108\117\101\115",z)return U end h.__index=function(z,U)if type(U)=="\110\117\109\098\101\114"then return z._values[U]elseif U=="\076\101\110\103\116\104"then return#z._values end error("\085\110\114\101\099\111\103\110\105\122\101\100 \107\101\121\058 "..tostring(U))end h.__newindex=function()error("\082\111\119\115 \097\114\101 \105\109\109\117\116\097\098\108\101",2)end h.__tostring=function(z)local U=""for Q=1,z.Length,1 do U=U..z[Q]if Q<z.Length then U=U.." "end end return U end h.__unm=function(z)if z.Length==1 then return-1*z[1]end local U={}for Q=1,z.Length,1 do U[Q]=-1*z[Q]end return e(U)end h.__add=function(z,U)if type(z)=="\110\117\109\098\101\114"or type(U)=="\110\117\109\098\101\114"then local Q=type(z)=="\110\117\109\098\101\114"and U or z if Q.Length==1 then local Q=type(z)=="\110\117\109\098\101\114"and z or z[1]local M=type(U)=="\110\117\109\098\101\114"and U or U[1]return Q+M end error("\065\100\100\105\116\105\111\110 \111\102 \097 \115\099\097\108\097\114 \116\111 \097 \114\111\119 \105\115 \117\110\100\101\102\105\110\101\100",2)end if z.Length~=U.Length then error("\065\100\100\105\116\105\111\110 \111\102 \114\111\119\115 \111\102 \100\105\102\102\101\114\101\110\116 \108\101\110\103\116\104\115 \105\115 \117\110\100\101\102\105\110\101\100",2)end if z.Length==1 then return z[1]+U[1]end local Q={}for M=1,z.Length,1 do Q[M]=z[M]+U[M]end return e(Q)end h.__sub=function(z,U)if type(z)=="\110\117\109\098\101\114"or type(U)=="\110\117\109\098\101\114"then local Q=type(z)=="\110\117\109\098\101\114"and U or z if Q.Length==1 then local Q=type(z)=="\110\117\109\098\101\114"and z or z[1]local M=type(U)=="\110\117\109\098\101\114"and U or U[1]return Q-M end error("\083\117\098\116\114\097\099\116\105\111\110 \111\102 \097 \115\099\097\108\097\114 \116\111 \097 \114\111\119 \105\115 \117\110\100\101\102\105\110\101\100",2)end if z.Length~=U.Length then error("\083\117\098\116\114\097\099\116\105\111\110 \111\102 \114\111\119\115 \111\102 \100\105\102\102\101\114\101\110\116 \108\101\110\103\116\104\115 \105\115 \117\110\100\101\102\105\110\101\100",2)end if z.Length==1 then return z[1]-U[1]end local Q={}for M=1,z.Length,1 do Q[M]=z[M]-U[M]end return e(Q)end h.__mul=function(z,U)if type(z)=="\110\117\109\098\101\114"or type(U)=="\110\117\109\098\101\114"then local Q=type(z)=="\110\117\109\098\101\114"and U or z local M=type(z)=="\110\117\109\098\101\114"and z or U if Q.Length==1 then local Q=type(z)=="\110\117\109\098\101\114"and z or z[1]local M=type(U)=="\110\117\109\098\101\114"and U or U[1]return Q*M else local z={}for U=1,Q.Length,1 do z[U]=M*Q[U]end return e(z)end elseif z.Length==1 or U.Length==1 then if z.Length==1 and U.Length==1 then return z[1]*U[1]elseif z.Length==1 then return z[1]*U else return z*U[1]end end error("\077\117\108\116\105\112\108\105\099\097\116\105\111\110 \111\102 \116\119\111 \114\111\119\115 \105\115 \117\110\100\101\102\105\110\101\100",2)end h.__div=function(z,U)if type(z)=="\110\117\109\098\101\114"then error("\068\105\118\105\115\105\111\110 \111\102 \097 \115\099\097\108\097\114 \098\121 \097 \114\111\119 \105\115 \117\110\100\101\102\105\110\101\100",2)elseif type(U)=="\110\117\109\098\101\114"then if z.Length==1 then return z[1]/U else local Q={}for M=1,z.Length,1 do Q[M]=z[M]/U end return e(Q)end elseif z.Length==1 or U.Length==1 then if z.Length==1 and U.Length==1 then return z[1]/U[1]elseif z.Length==1 then error("\068\105\118\105\115\105\111\110 \111\102 \097 \115\099\097\108\097\114 \098\121 \097 \114\111\119 \105\115 \117\110\100\101\102\105\110\101\100",2)else return z/U[1]end end error("\068\105\118\105\115\105\111\110 \111\102 \116\119\111 \114\111\119\115 \105\115 \117\110\100\101\102\105\110\101\100",2)end h.__mod=function(z,U)if type(z)=="\110\117\109\098\101\114"then error("\077\111\100\117\108\111 \111\102 \097 \115\099\097\108\097\114 \098\121 \097 \114\111\119 \105\115 \117\110\100\101\102\105\110\101\100",2)elseif type(U)=="\110\117\109\098\101\114"then if z.Length==1 then return z[1]%U else local Q={}for M=1,z.Length,1 do Q[M]=z[M]%U end return e(Q)end elseif z.Length==1 or U.Length==1 then if z.Length==1 and U.Length==1 then return z[1]%U[1]elseif z.Length==1 then error("\077\111\100\117\108\111 \111\102 \097 \115\099\097\108\097\114 \098\121 \097 \114\111\119 \105\115 \117\110\100\101\102\105\110\101\100",2)else return z%U[1]end end error("\077\111\100\117\108\111 \111\102 \116\119\111 \114\111\119\115 \105\115 \117\110\100\101\102\105\110\101\100",2)end h.__pow=function(z,U)if type(z)=="\110\117\109\098\101\114"then if U.Length==1 then return z^U[1]else error("\069\120\112\111\110\101\110\116\105\097\116\105\111\110 \111\102 \097 \115\099\097\108\097\114 \098\121 \097 \114\111\119 \105\115 \117\110\100\101\102\105\110\101\100",2)end elseif type(U)=="\110\117\109\098\101\114"then if z.Length==1 then return z[1]^U else error("\069\120\112\111\110\101\110\116\105\097\116\105\111\110 \111\102 \097 \114\111\119 \098\121 \097 \115\099\097\108\097\114 \105\115 \117\110\100\101\102\105\110\101\100",2)end elseif z.Length==1 or U.Length==1 then if z.Length==1 and U.Length==1 then return z[1]^U[1]elseif z.Length==1 then error("\069\120\112\111\110\101\110\116\105\097\116\105\111\110 \111\102 \097 \115\099\097\108\097\114 \098\121 \097 \114\111\119 \105\115 \117\110\100\101\102\105\110\101\100",2)else error("\069\120\112\111\110\101\110\116\105\097\116\105\111\110 \111\102 \097 \114\111\119 \098\121 \097 \115\099\097\108\097\114 \105\115 \117\110\100\101\102\105\110\101\100",2)end end error("\069\120\112\111\110\101\110\116\105\097\116\105\111\110 \111\102 \116\119\111 \114\111\119\115 \105\115 \117\110\100\101\102\105\110\101\100",2)end h.__eq=function(z,U)if z.Length~=U.Length then return false end for Q=1,z.Length,1 do if z[Q]~=U[Q]then return false end end return true end local a={}local b=function(z)local U={}for Q=1,#z,1 do U[Q]=e(z[Q])end local Q=setmetatable({},a)rawset(Q,"\095\114\111\119\115",U)return Q end a.__index=function(U,Q)if type(Q)=="\110\117\109\098\101\114"then return U._rows[Q]elseif Q=="\083\104\097\112\101"then return{#U._rows;U._rows[1].Length}elseif Q=="\084"then return z.matrix.transpose(U)elseif Q=="\085\110\105\116"then if U.Shape[2]==1 then return z.vector.unit(U)elseif U.Shape[1]==U.Shape[2]then return z.matrix.identity(U.Shape[1])else error("\085\110\105\116 \111\102 \097 \110\111\110-\115\113\117\097\114\101 \109\097\116\114\105\120 \105\115 \117\110\100\101\102\105\110\101\100",2)end elseif z.matrix[Q]then if type(z.matrix[Q])=="\102\117\110\099\116\105\111\110"then return function(...)return z.matrix[Q](U,...)end end elseif z.vector[Q]and U.Shape[2]==1 then if type(z.vector[Q])=="\102\117\110\099\116\105\111\110"then return function(...)return z.vector[Q](U,...)end end end error("\085\110\114\101\099\111\103\110\105\122\101\100 \107\101\121\058 "..tostring(Q))end a.__newindex=function()error("\077\097\116\114\105\099\101\115 \097\114\101 \105\109\109\117\116\097\098\108\101",2)end a.__tostring=function(z)local U=""for Q=1,z.Shape[1],1 do U=U..tostring(z[Q])if Q<z.Shape[1]then U=U.."\n"end end return U end a.__unm=function(z)local U={}for Q=1,z.Shape[1],1 do U[Q]={}for M=1,z.Shape[2],1 do U[Q][M]=-1*z[Q][M]end end return b(U)end a.__add=function(z,U)if z.Shape[1]~=U.Shape[1]or z.Shape[2]~=U.Shape[2]then error("\067\097\110\110\111\116 \097\100\100 \109\097\116\114\105\099\101\115 \111\102 \115\105\122\101\115 \040"..(z.Shape[1]..(" \120 "..(z.Shape[2]..("\041 \097\110\100 \040"..(U.Shape[1]..(" \120 "..(U.Shape[2].."\041"))))))),2)end local Q={}for M=1,z.Shape[1],1 do Q[M]={}for v=1,z.Shape[2],1 do Q[M][v]=z[M][v]+U[M][v]end end return b(Q)end a.__sub=function(z,U)if z.Shape[1]~=U.Shape[1]or z.Shape[2]~=U.Shape[2]then error("\067\097\110\110\111\116 \115\117\098\116\114\097\099\116 \109\097\116\114\105\099\101\115 \111\102 \115\105\122\101\115 \040"..(z.Shape[1]..(" \120 "..(z.Shape[2]..("\041 \097\110\100 \040"..(U.Shape[1]..(" \120 "..(U.Shape[2].."\041"))))))),2)end local Q={}for M=1,z.Shape[1],1 do Q[M]={}for v=1,z.Shape[2],1 do Q[M][v]=z[M][v]-U[M][v]end end return b(Q)end a.__mul=function(z,U)local Q={}if type(z)=="\110\117\109\098\101\114"or type(U)=="\110\117\109\098\101\114"then local M=type(z)=="\110\117\109\098\101\114"and U or z local v=type(z)=="\110\117\109\098\101\114"and z or U for z=1,M.Shape[1],1 do Q[z]={}for U=1,M.Shape[2],1 do Q[z][U]=v*M[z][U]end end else if z.Shape[2]~=U.Shape[1]then error("\067\097\110\110\111\116 \109\117\108\116\105\112\108\121 \109\097\116\114\105\099\101\115 \111\102 \115\105\122\101\115 \040"..(z.Shape[1]..(" \120 "..(z.Shape[2]..("\041 \097\110\100 \040"..(U.Shape[1]..(" \120 "..(U.Shape[2].."\041"))))))),2)end for M=1,z.Shape[1],1 do Q[M]={}for v=1,U.Shape[2],1 do Q[M][v]=0 for n=1,z.Shape[2],1 do Q[M][v]=Q[M][v]+z[M][n]*U[n][v]end end end end return b(Q)end a.__div=function(z,U)if type(z)=="\110\117\109\098\101\114"then error("\067\097\110\110\111\116 \100\105\118\105\100\101 \097 \115\099\097\108\097\114 \098\121 \097 \109\097\116\114\105\120",2)elseif type(U)~="\110\117\109\098\101\114"then error("\067\097\110\110\111\116 \100\105\118\105\100\101 \097 \109\097\116\114\105\120 \098\121 \097 \109\097\116\114\105\120",2)end local Q={}for M=1,z.Shape[1],1 do Q[M]={}for v=1,z.Shape[2],1 do Q[M][v]=z[M][v]/U end end return b(Q)end a.__mod=function(z,U)if type(z)=="\110\117\109\098\101\114"then error("\067\097\110\110\111\116 \099\111\109\112\117\116\101 \109\111\100\117\108\117\115 \111\102 \097 \115\099\097\108\097\114 \098\121 \097 \109\097\116\114\105\120",2)elseif type(U)~="\110\117\109\098\101\114"then error("\067\097\110\110\111\116 \099\111\109\112\117\116\101 \109\111\100\117\108\117\115 \111\102 \097 \109\097\116\114\105\120 \098\121 \097 \109\097\116\114\105\120",2)end local Q={}for M=1,z.Shape[1],1 do Q[M]={}for v=1,z.Shape[2],1 do Q[M][v]=z[M][v]%U end end return b(Q)end a.__pow=function(U,M)if type(U)~="\110\117\109\098\101\114"and type(M)=="\110\117\109\098\101\114"then if U.Shape[1]~=U.Shape[2]then error("\067\097\110\110\111\116 \101\120\112\111\110\101\110\116\105\097\116\101 \097 \110\111\110-\115\113\117\097\114\101 \109\097\116\114\105\120",2)end if Q(M)~=M or M<1 then error("\067\097\110\110\111\116 \101\120\112\111\110\101\110\116\105\097\116\101 \097 \109\097\116\114\105\120 \098\121 \097 \110\111\110-\112\111\115\105\116\105\118\101 \110\111\110-\105\110\116\101\103\101\114",2)end local v=z.matrix.identity(U.Shape[1])local n=U while M>0 do if M%2==1 then v=v*n end n=n*n M=Q(M/2)end return v elseif type(U)=="\110\117\109\098\101\114"and type(M)~="\110\117\109\098\101\114"then error("\067\097\110\110\111\116 \114\097\105\115\101 \097\110 \097\114\098\105\116\114\097\114\121 \110\117\109\098\101\114 \116\111 \097 \109\097\116\114\105\120 \112\111\119\101\114",2)else error("\067\097\110\110\111\116 \101\120\112\111\110\101\110\116\105\097\116\101 \097 \109\097\116\114\105\120 \098\121 \097 \109\097\116\114\105\120",2)end end a.__eq=function(z,U)if z.Shape[1]~=U.Shape[1]or z.Shape[2]~=U.Shape[2]then return false end for Q=1,z.Shape[1],1 do for M=1,z.Shape[2],1 do if z[Q][M]~=U[Q][M]then return false end end end return true end z.matrix={}z.matrix.new=function(z)return b(z)end z.matrix.fromColumns=function(z)local U={}local Q=z[1].Shape and z[1].Shape[1]or#z[1]for M=1,#z,1 do for Q=1,Q,1 do if not U[Q]then U[Q]={}end U[Q][M]=z[M][Q]*1 end end return b(U)end z.matrix.identity=function(z)local U={}for Q=1,z,1 do U[Q]={}for z=1,z,1 do U[Q][z]=Q==z and 1 or 0 end end return b(U)end z.matrix.zeros=function(z,U)local Q={}for z=1,z,1 do Q[z]={}for U=1,U,1 do Q[z][U]=0 end end return b(Q)end z.matrix.zerosLike=function(U)return z.matrix.zeros(U.Shape[1],U.Shape[2])end z.matrix.diagonal=function(z)local U={}for Q=1,#z,1 do U[Q]={}for M=1,#z,1 do U[Q][M]=Q==M and z[Q]or 0 end end return b(U)end z.matrix.isDiagonal=function(z)for U=1,z.Shape[1],1 do for Q=1,z.Shape[2],1 do if U~=Q and z[U][Q]~=0 then return false end end end return true end z.matrix.isUpperTriangular=function(z)if z.Shape[1]~=z.Shape[2]then return false end for U=1,z.Shape[1],1 do for Q=1,z.Shape[2],1 do if U>Q and z[U][Q]~=0 then return false end end end return true end z.matrix.isLowerTriangular=function(z)if z.Shape[1]~=z.Shape[2]then return false end for U=1,z.Shape[1],1 do for Q=1,z.Shape[2],1 do if U<Q and z[U][Q]~=0 then return false end end end return true end z.matrix.transpose=function(z)local U={}for Q=1,z.Shape[1],1 do for M=1,z.Shape[2],1 do if not U[M]then U[M]={}end U[M][Q]=z[Q][M]end end return b(U)end z.matrix.expm=function(U,Q)if U.Shape[1]~=U.Shape[2]then error("\067\097\110\110\111\116 \102\105\110\100 \116\104\101 \101\120\112\111\110\101\110\116\105\097\108 \111\102 \097 \110\111\110-\115\113\117\097\114\101 \109\097\116\114\105\120",2)end if not Q then Q=128 end local M=z.matrix.identity(U.Shape[1])local v=M local n=1 for z=1,Q,1 do v=v*U n=n*z M=M+(1/n)*v end return M end z.vector={}z.vector.new=function(z)local U={}for Q=1,#z,1 do U[Q]={z[Q]}end return b(U)end z.vector.e=function(z,U)local Q={}for U=1,U,1 do Q[U]=z==U and{1}or{0}end return b(Q)end z.vector.norm={}z.vector.norm.l1=function(z)local Q=0 for M=1,z.Shape[1],1 do Q=Q+U(z[M][1])end return Q end z.vector.norm.l2=function(z)local U=0 for Q=1,z.Shape[1],1 do U=U+z[Q]^2 end return M(U)end z.vector.norm.linf=function(z)local Q=0 for M=1,z.Shape[1],1 do local v=U(z[M][1])if v>Q then Q=v end end return Q end z.vector.ip={}z.vector.ip.dot=function(z,U)if z.Shape[1]~=U.Shape[1]then error("\067\097\110\110\111\116 \099\111\109\112\117\116\101 \116\104\101 \100\111\116 \112\114\111\100\117\099\116 \102\111\114 \118\101\099\116\111\114\115 \111\102 \117\110\101\113\117\097\108 \100\105\109\101\110\115\105\111\110\115",2)end if z.Shape[2]~=1 or U.Shape[2]~=1 then error("\067\097\110\110\111\116 \099\111\109\112\117\116\101 \116\104\101 \100\111\116 \112\114\111\100\117\099\116 \102\111\114 \110\111\110 \099\111\108\117\109\110 \118\101\099\116\111\114\115",2)end local Q=0 for M=1,z.Shape[1],1 do Q=Q+z[M]*U[M]end return Q end z.vector.project=function(U,Q,M,v)if#Q==0 then Q={Q}end if Q[1].Shape[1]~=U.Shape[1]then error("\067\097\110\110\111\116 \112\114\111\106\101\099\116 \118\101\099\116\111\114\115 \111\102 \100\105\102\102\101\114\101\110\116 \100\105\109\101\110\115\105\111\110\115",2)end if Q[1].Shape[2]~=1 or U.Shape[2]~=1 then error("\069\120\112\101\099\116\101\100 \116\119\111 \118\101\099\116\111\114\115\044 \110\111\116 \109\097\116\114\105\099\101\115",2)end if not M then M=z.vector.ip.dot end if not v then v=z.vector.norm.l2 end local n=z.matrix.zerosLike(U)for z=1,#Q,1 do n=n+(M(U,Q[z])/v(Q[z])^2)*Q[z]end return n end z.vector.unit=function(U,Q)if not Q then Q=z.vector.norm.l2 end return U/Q(U)end z.vector.createArbitraryAxisRotationMatrix=function(z,U)if z.Shape[1]~=3 then error("\067\097\110\110\111\116 \099\114\101\097\116\101 \114\111\116\097\116\105\111\110 \109\097\116\114\105\120 \097\098\111\117\116 \097\114\098\105\116\114\097\114\121 \117\110\105\116 \118\101\099\116\111\114 \111\116\104\101\114 \116\104\097\110 \105\110 \082\094\051",2)end local Q=v(U)local M=n(U)local h=1-Q local e=z.Unit return b({{h*e[1]^2+Q,(h*e[1])*e[2]-M*e[3],(h*e[1])*e[3]+M*e[2]};{(h*e[1])*e[2]+M*e[3],h*e[2]^2+Q,(h*e[2])*e[3]-M*e[1]},{(h*e[1])*e[3]-M*e[2],(h*e[2])*e[3]+M*e[1];h*e[3]^2+Q}})end z.gramSchmidt=function(U,Q,M,v,n)if#U==0 then U={U}end if not Q then Q=.01 end if not M then M=U[1].Shape[1]end if not v then v=z.vector.ip.dot end if not n then n=z.vector.norm.l2 end local h={}for z=1,#U,1 do h[z]=U[z].Unit end U=h for z=1,M,1 do local Q={}for U=1,M,1 do Q[U]={z==U and 1 or 0}end table.insert(U,b(Q))end local e={U[1].Unit}for h=1,#U,1 do local a=U[h]-z.vector.project(U[h],e,v,n)if n(a)>=Q then table.insert(e,a.Unit)end if#e==M then break end end return e end return z
local linalg = {}

-- ARITHMETIC FUNCTIONS
local abs = math.abs
local floor = math.floor
local sqrt = math.sqrt

-- TRIG FUNCTIONS
local cos = math.cos
local sin = math.sin

--[[ ROW CLASS ]]--
local _rowClass = {}
local _newRow = function(values)
	local instance = setmetatable({}, _rowClass)

	rawset(instance, "_values", values)

	return instance
end
_rowClass.__index = function(row, key)
	if type(key) == "number" then
		return row._values[key]
	elseif key == "Length" then
		return #row._values
	end

	error("Unrecognized key: " .. tostring(key))
end
_rowClass.__newindex = function()
	error("Rows are immutable", 2)
end
_rowClass.__tostring = function(row)
	local result = ""

	for i = 1, row.Length do
		result = result .. row[i]
		if i < row.Length then
			result = result .. " "
		end
	end

	return result
end
_rowClass.__unm = function (row)
	if row.Length == 1 then
		return -1 * row[1]
	end

	local resultArray = {}

	for i = 1, row.Length do
		resultArray[i] = -1 * row[i]
	end

	return _newRow(resultArray)
end
_rowClass.__add = function (left, right)
	if type(left) == "number" or type(right) == "number" then
		local row = type(left) == "number" and right or left

		if row.Length == 1 then
			local leftScalar = type(left) == "number" and left or left[1]
			local rightScalar = type(right) == "number" and right or right[1]
			return leftScalar + rightScalar
		end

		error("Addition of a scalar to a row is undefined", 2)
	end

	if left.Length ~= right.Length then
		error("Addition of rows of different lengths is undefined", 2)
	end

	if left.Length == 1 then
		return left[1] + right[1]
	end

	local resultArray = {}

	for i = 1, left.Length do
		resultArray[i] = left[i] + right[i]
	end

	return _newRow(resultArray)
end
_rowClass.__sub = function (left, right)
	if type(left) == "number" or type(right) == "number" then
		local row = type(left) == "number" and right or left

		if row.Length == 1 then
			local leftScalar = type(left) == "number" and left or left[1]
			local rightScalar = type(right) == "number" and right or right[1]
			return leftScalar - rightScalar
		end

		error("Subtraction of a scalar to a row is undefined", 2)
	end

	if left.Length ~= right.Length then
		error("Subtraction of rows of different lengths is undefined", 2)
	end

	if left.Length == 1 then
		return left[1] - right[1]
	end

	local resultArray = {}

	for i = 1, left.Length do
		resultArray[i] = left[i] - right[i]
	end

	return _newRow(resultArray)
end
_rowClass.__mul = function (left, right)
	if type(left) == "number" or type(right) == "number" then
		local row = type(left) == "number" and right or left
		local scalar = type(left) == "number" and left or right

		if row.Length == 1 then
			local leftScalar = type(left) == "number" and left or left[1]
			local rightScalar = type(right) == "number" and right or right[1]
			return leftScalar * rightScalar
		else
			local resultArray = {}

			for i = 1, row.Length do
				resultArray[i] = scalar * row[i]
			end

			return _newRow(resultArray)
		end
	elseif left.Length == 1 or right.Length == 1 then
		if left.Length == 1 and right.Length == 1 then
			return left[1] * right[1]
		elseif left.Length == 1 then
			return left[1] * right
		else
			return left * right[1]
		end
	end

	error("Multiplication of two rows is undefined", 2)
end
_rowClass.__div = function (left, right)
	if type(left) == "number" then
		error("Division of a scalar by a row is undefined", 2)
	elseif type(right) == "number" then
		if left.Length == 1 then
			return left[1] / right
		else
			local resultArray = {}

			for i = 1, left.Length do
				resultArray[i] = left[i] / right
			end

			return _newRow(resultArray)
		end
	elseif left.Length == 1 or right.Length == 1 then
		if left.Length == 1 and right.Length == 1 then
			return left[1] / right[1]
		elseif left.Length == 1 then
			error("Division of a scalar by a row is undefined", 2)
		else
			return left / right[1]
		end
	end

	error("Division of two rows is undefined", 2)
end
_rowClass.__mod = function (left, right)
	if type(left) == "number" then
		error("Modulo of a scalar by a row is undefined", 2)
	elseif type(right) == "number" then
		if left.Length == 1 then
			return left[1] % right
		else
			local resultArray = {}

			for i = 1, left.Length do
				resultArray[i] = left[i] % right
			end

			return _newRow(resultArray)
		end
	elseif left.Length == 1 or right.Length == 1 then
		if left.Length == 1 and right.Length == 1 then
			return left[1] % right[1]
		elseif left.Length == 1 then
			error("Modulo of a scalar by a row is undefined", 2)
		else
			return left % right[1]
		end
	end

	error("Modulo of two rows is undefined", 2)
end
_rowClass.__pow = function (left, right)
	if type(left) == "number" then
		if right.Length == 1 then
			return left ^ right[1]
		else
			error("Exponentiation of a scalar by a row is undefined", 2)
		end
	elseif type(right) == "number" then
		if left.Length == 1 then
			return left[1] ^ right
		else
			error("Exponentiation of a row by a scalar is undefined", 2)
		end
	elseif left.Length == 1 or right.Length == 1 then
		if left.Length == 1 and right.Length == 1 then
			return left[1] ^ right[1]
		elseif left.Length == 1 then
			error("Exponentiation of a scalar by a row is undefined", 2)
		else
			error("Exponentiation of a row by a scalar is undefined", 2)
		end
	end

	error("Exponentiation of two rows is undefined", 2)
end
_rowClass.__eq = function(left, right)
	if left.Length ~= right.Length then
		return false
	end

	for i = 1, left.Length do
		if left[i] ~= right[i] then
			return false
		end
	end

	return true
end

--[[ MATRIX CLASS ]]--
local _matrixClass = {}
local _newMatrix = function(rows)
	local rowInstances = {}
	for i = 1, #rows do
		rowInstances[i] = _newRow(rows[i])
	end

	local instance = setmetatable({}, _matrixClass)
	rawset(instance, "_rows", rowInstances)

	return instance
end
_matrixClass.__index = function(mat, key)
	if type(key) == "number" then
		return mat._rows[key]
	elseif key == "Shape" then
		return { #mat._rows, mat._rows[1].Length }
	elseif key == "T" then
		return linalg.matrix.transpose(mat)
	elseif key == "Unit" then
		if mat.Shape[2] == 1 then
			return linalg.vector.unit(mat)
		elseif mat.Shape[1] == mat.Shape[2] then
			return linalg.matrix.identity(mat.Shape[1])
		else
			error("Unit of a non-square matrix is undefined", 2)
		end
	elseif linalg.matrix[key] then
		if type(linalg.matrix[key]) == "function" then
			return function (...) return linalg.matrix[key](mat, ...) end
		end
	elseif linalg.vector[key] and mat.Shape[2] == 1 then
		if type(linalg.vector[key]) == "function" then
			return function (...) return linalg.vector[key](mat, ...) end
		end
	end

	error("Unrecognized key: " .. tostring(key))
end
_matrixClass.__newindex = function()
	error("Matrices are immutable", 2)
end
_matrixClass.__tostring = function(mat)
	local result = ""

	for i = 1, mat.Shape[1] do
		result = result .. tostring(mat[i])
		if i < mat.Shape[1] then
			result = result .. "\n"
		end
	end

	return result
end
_matrixClass.__unm = function(mat)
	local resultRows = {}

	for i = 1, mat.Shape[1] do
		resultRows[i] = {}
		for j = 1, mat.Shape[2] do
			resultRows[i][j] = -1 * mat[i][j]
		end
	end

	return _newMatrix(resultRows)
end
_matrixClass.__add = function(left, right)
	if left.Shape[1] ~= right.Shape[1] or left.Shape[2] ~= right.Shape[2] then
		error(
			"Cannot add matrices of sizes (" ..
				left.Shape[1] .. " x " .. left.Shape[2] .. ") and (" .. right.Shape[1] .. " x " .. right.Shape[2] .. ")",
			2
		)
	end

	local resultRows = {}

	for i = 1, left.Shape[1] do
		resultRows[i] = {}
		for j = 1, left.Shape[2] do
			resultRows[i][j] = left[i][j] + right[i][j]
		end
	end

	return _newMatrix(resultRows)
end
_matrixClass.__sub = function(left, right)
	if left.Shape[1] ~= right.Shape[1] or left.Shape[2] ~= right.Shape[2] then
		error(
			"Cannot subtract matrices of sizes (" ..
				left.Shape[1] .. " x " .. left.Shape[2] .. ") and (" .. right.Shape[1] .. " x " .. right.Shape[2] .. ")",
			2
		)
	end

	local resultRows = {}

	for i = 1, left.Shape[1] do
		resultRows[i] = {}
		for j = 1, left.Shape[2] do
			resultRows[i][j] = left[i][j] - right[i][j]
		end
	end

	return _newMatrix(resultRows)
end
_matrixClass.__mul = function(left, right)
	local resultRows = {}

	if type(left) == "number" or type(right) == "number" then
		local mat = type(left) == "number" and right or left
		local scalar = type(left) == "number" and left or right

		for i = 1, mat.Shape[1] do
			resultRows[i] = {}
			for j = 1, mat.Shape[2] do
				resultRows[i][j] = scalar * mat[i][j]
			end
		end
	else
		if left.Shape[2] ~= right.Shape[1] then
			error(
				"Cannot multiply matrices of sizes (" ..
					left.Shape[1] .. " x " .. left.Shape[2] .. ") and (" .. right.Shape[1] .. " x " .. right.Shape[2] .. ")",
				2
			)
		end

		for i = 1, left.Shape[1] do
			resultRows[i] = {}
			for j = 1, right.Shape[2] do
				resultRows[i][j] = 0
				for k = 1, left.Shape[2] do
					resultRows[i][j] = resultRows[i][j] + (left[i][k] * right[k][j])
				end
			end
		end
	end

	return _newMatrix(resultRows)
end
_matrixClass.__div = function(mat, scalar)
	if type(mat) == "number" then
		error("Cannot divide a scalar by a matrix", 2)
	elseif type(scalar) ~= "number" then
		error("Cannot divide a matrix by a matrix", 2)
	end

	local resultRows = {}

	for i = 1, mat.Shape[1] do
		resultRows[i] = {}
		for j = 1, mat.Shape[2] do
			resultRows[i][j] = mat[i][j] / scalar
		end
	end

	return _newMatrix(resultRows)
end
_matrixClass.__mod = function(mat, scalar)
	if type(mat) == "number" then
		error("Cannot compute modulus of a scalar by a matrix", 2)
	elseif type(scalar) ~= "number" then
		error("Cannot compute modulus of a matrix by a matrix", 2)
	end

	local resultRows = {}

	for i = 1, mat.Shape[1] do
		resultRows[i] = {}
		for j = 1, mat.Shape[2] do
			resultRows[i][j] = mat[i][j] % scalar
		end
	end

	return _newMatrix(resultRows)
end
_matrixClass.__pow = function(base, power)
	if type(base) ~= "number" and type(power) == "number" then
		if base.Shape[1] ~= base.Shape[2] then
			error("Cannot exponentiate a non-square matrix", 2)
		end
		if floor(power) ~= power or power < 1 then
			error("Cannot exponentiate a matrix by a non-positive non-integer", 2)
		end

		-- Runs in O(n^2 log_2(n)) time
		-- See https://www.hackerearth.com/practice/notes/matrix-exponentiation-1/
		local resultMatrix = linalg.matrix.identity(base.Shape[1])
		local curMat = base

		while power > 0 do
			if power % 2 == 1 then
				resultMatrix = resultMatrix * curMat
			end

			curMat = curMat * curMat
			power = floor(power / 2)
		end

		return resultMatrix
	elseif type(base) == "number" and type(power) ~= "number" then
		error("Cannot raise an arbitrary number to a matrix power", 2)
	else
		error("Cannot exponentiate a matrix by a matrix", 2)
		--TODO (You should be able to do exactly this)
	end
end
_matrixClass.__eq = function(left, right)
	if left.Shape[1] ~= right.Shape[1] or left.Shape[2] ~= right.Shape[2] then
		return false
	end

	for i = 1, left.Shape[1] do
		for j = 1, left.Shape[2] do
			if left[i][j] ~= right[i][j] then
				return false
			end
		end
	end

	return true
end

-- MATRIX FUNCTIONS
linalg.matrix = {}

--[[**
	Creates a new matrix

	@param [t:array<array<number>>] rows A (m x n) array of numbers to fill the matrix with

	@returns [t:(m x n) matrix] The new matrix
**--]]
linalg.matrix.new = function(rows)
	return _newMatrix(rows)
end

linalg.matrix.fromColumns = function(columnVectors)
	local resultRows = {}
	local n = columnVectors[1]["Shape"] and columnVectors[1].Shape[1] or #columnVectors[1]

	for i = 1, #columnVectors do
		for j = 1, n do
			if not resultRows[j] then
				resultRows[j] = {}
			end

			resultRows[j][i] = columnVectors[i][j] * 1 --*1 to ensure scalar functionality
		end
	end

	return _newMatrix(resultRows)
end

--[[**
	Creates an identity matrix of size (n x n)

	@param [t:number] n The size of the matrix

	@returns [t:(n x n) matrix] The identity matrix
**--]]
linalg.matrix.identity = function(n)
	local rows = {}

	for i = 1, n do
		rows[i] = {}
		for j = 1, n do
			rows[i][j] = i == j and 1 or 0
		end
	end

	return _newMatrix(rows)
end

--[[**
	Creates a matrix of all zeros of size (m x n)

	@param [t:number] m
	@param [t:number] n

	@returns [t:(m x n) matrix] The zeros matrix
**--]]
linalg.matrix.zeros = function(m, n)
	local rows = {}

	for i = 1, m do
		rows[i] = {}
		for j = 1, n do
			rows[i][j] = 0
		end
	end

	return _newMatrix(rows)
end

--[[**
	Creates a matrix of all zeros of the same size as the provided matrix

	@param [t:(m x n) matrix] mat The matrix to copy the size of

	@returns [t:(m x n) matrix] The zeros matrix
**--]]
linalg.matrix.zerosLike = function(mat)
	return linalg.matrix.zeros(mat.Shape[1], mat.Shape[2])
end

--[[**
	Creates a diagonal matrix with the values provided as the diagonal entries

	@param [t:array<number>] values An n-length array whose entries will be set as the diagonal entries

	@returns [t:(n x n) matrix] The resulting diagonal matrix
**--]]
linalg.matrix.diagonal = function(values)
	local resultRows = {}

	for i = 1, #values do
		resultRows[i] = {}
		for j = 1, #values do
			resultRows[i][j] = i == j and values[i] or 0
		end
	end

	return _newMatrix(resultRows)
end

--[[**
	Determines whether a matrix is diagonal
	Does not exclusively refer to square matrices
	Refers strictly to whether all non-zero values are on the main diagonal (i.e., a_{ij} = 0 for all i, j where i ~= j)

	@param [t:(m x n) matrix] mat The matrix to check

	@returns [t:boolean] True if the matrix is diagonal, false otherwise
**--]]
linalg.matrix.isDiagonal = function(mat)
	for i = 1, mat.Shape[1] do
		for j = 1, mat.Shape[2] do
			if i ~= j and mat[i][j] ~= 0 then
				return false
			end
		end
	end

	return true
end

--[[**
	Determines whether a matrix is upper triangular
	Note that any non-square matrix will return false

	@param [t:(m x n) matrix] mat The matrix to check

	@returns [t:boolean] True if the matrix is upper triangular, false otherwise
**--]]
linalg.matrix.isUpperTriangular = function(mat)
	if mat.Shape[1] ~= mat.Shape[2] then
		return false
	end

	for i = 1, mat.Shape[1] do
		for j = 1, mat.Shape[2] do
			if i > j and mat[i][j] ~= 0 then
				return false
			end
		end
	end

	return true
end

--[[**
	Determines whether a matrix is lower triangular
	Note that any non-square matrix will return false

	@param [t:(m x n) matrix] mat The matrix to check

	@returns [t:boolean] True if the matrix is lower triangular, false otherwise
**--]]
linalg.matrix.isLowerTriangular = function(mat)
	if mat.Shape[1] ~= mat.Shape[2] then
		return false
	end

	for i = 1, mat.Shape[1] do
		for j = 1, mat.Shape[2] do
			if i < j and mat[i][j] ~= 0 then
				return false
			end
		end
	end

	return true
end

--[[**
	Creates a new matrix that is the transpose of the provided matrix

	@param [t: (m x n) matrix] mat The matrix to create the transpose of

	@returns [t: (n x m) matrix] The transpose of mat
**--]]
linalg.matrix.transpose = function(mat)
	local resultRows = {}

	for i = 1, mat.Shape[1] do
		for j = 1, mat.Shape[2] do
			if not resultRows[j] then
				resultRows[j] = {}
			end

			resultRows[j][i] = mat[i][j]
		end
	end

	return _newMatrix(resultRows)
end

--[[**
	Solves for e^mat
	Defined as: e^A = \sum_{k=0}^{\infinity} \frac{1}{k!} A^k
	Implemented in a naive way to approximate by using iterations
	Runtime of O(n^3)

	@param [t:(n x n) matrix] mat The matrix to use as the exponent
	@param [t:number] numIterations The number of iterations to take the sum of the taylor series to

	@returns The matrix exponential approximation
**--]]
linalg.matrix.expm = function(mat, numIterations)
	if mat.Shape[1] ~= mat.Shape[2] then
		error("Cannot find the exponential of a non-square matrix", 2)
	end

	if not numIterations then
		numIterations = 128
	end

	local resultMatrix = linalg.matrix.identity(mat.Shape[1])
	local curMat = resultMatrix
	local factorial = 1

	for i = 1, numIterations do
		curMat = curMat * mat
		factorial = factorial * i

		resultMatrix = resultMatrix + ((1 / factorial) * curMat)
	end

	return resultMatrix
end

-- VECTOR FUNCTIONS
linalg.vector = {}

--[[**
	Creates a new column vector

	@param [t:array<number>] values The values to have for the column vector

	@returns [t:(n x 1) matrix] The new column vector
**--]]
linalg.vector.new = function(values)
	local rows = {}

	for i = 1, #values do
		rows[i] = {values[i]}
	end

	return _newMatrix(rows)
end

--[[**
	Creates the standard basis vector i for R^n
	That is, creates a vector of length n with all zeros except at index i which will have value 1

	@param [t:number] i The index of e
	@param [t:number] n The dimensionality of the vector

	@returns [t:(n x 1) matrix] The standard basis vector e_i in R^n
**--]]
linalg.vector.e = function(i, n)
	local rows = {}

	for j = 1, n do
		rows[j] = i == j and {1} or {0}
	end

	return _newMatrix(rows)
end

linalg.vector.norm = {}

--[[**
	The L1 norm of a vector
	sum_i{|v_i|}

	@param [t:(n x 1) matrix] v The vector

	@returns [t:number] The resulting value
**--]]
linalg.vector.norm.l1 = function(v)
	local sum = 0

	for i = 1, v.Shape[1] do
		sum = sum + abs(v[i][1])
	end

	return sum
end

--[[**
	The L2 norm of a vector
	sqrt(sum_i{(v_i)^2})

	@param [t:(n x 1) matrix] v The vector

	@returns [t:number] The resulting value
**--]]
linalg.vector.norm.l2 = function(v)
	local sum = 0

	for i = 1, v.Shape[1] do
		sum = sum + v[i]^2
	end

	return sqrt(sum)
end

--[[**
	The L-infinity norm of a vector
	max{v}

	@param [t:(n x 1) matrix] v The vector

	@returns [t:number] The resulting value
**--]]
linalg.vector.norm.linf = function(v)
	local maxComponentValue = 0

	for i = 1, v.Shape[1] do
		local componentValue = abs(v[i][1])
		if componentValue > maxComponentValue then
			maxComponentValue = componentValue
		end
	end

	return maxComponentValue
end

-- Inner product functions
linalg.vector.ip = {}

--[[**
	Computes the standard dot product of two vectors
	Defined as \sum_{i=0}^{n-1} v1[i] * v2[i]

	@param [t:(n x 1) matrix] v1 The first vector
	@param [t:(n x 1) matrix] v2 The second vector

	@returns [t:number] The result
**--]]
linalg.vector.ip.dot = function(v1, v2)
	if v1.Shape[1] ~= v2.Shape[1] then
		error("Cannot compute the dot product for vectors of unequal dimensions", 2)
	end
	if v1.Shape[2] ~= 1 or v2.Shape[2] ~= 1 then
		error("Cannot compute the dot product for non column vectors", 2)
	end

	local result = 0

	for i = 1, v1.Shape[1] do
		result = result + (v1[i] * v2[i])
	end

	return result
end

--[[**
	Projects vector v onto vector space u
	Defined as \sum_{i=0}^{m-1} <v, u[i]>/|u[i]|^2 * u[i]

	@param [t:(n x 1) matrix] v The vector to project onto u
	@param [t:array<(n x 1) matrix>] u The vector space to project v onto (can also be just one vector)
	@param [t:function([(n x 1) matrix], [(n x 1) matrix])?] innerProductFunc The inner product function to use; Defaults to the dot product
	@param [t:function([(n x 1) matrix])?] normFunc The norm function to use; Defaults to the L2 norm

	@returns The vector projection of v onto u
**--]]
linalg.vector.project = function (v, u, innerProductFunc, normFunc)
	if #u == 0 then -- u must be a single vector
		u = { u }
	end

	if u[1].Shape[1] ~= v.Shape[1] then
		error("Cannot project vectors of different dimensions", 2)
	end
	if u[1].Shape[2] ~= 1 or v.Shape[2] ~= 1 then
		error("Expected two vectors, not matrices", 2)
	end

	if not innerProductFunc then
		innerProductFunc = linalg.vector.ip.dot
	end
	if not normFunc then
		normFunc = linalg.vector.norm.l2
	end

	local result = linalg.matrix.zerosLike(v)

	for i = 1, #u do
		result = result + (innerProductFunc(v, u[i]) / normFunc(u[i])^2) * u[i]
	end

	return result
end

--[[**
	Gets the unit vector with the same direction as the provided vectr

	@param [t:(n x 1) matrix] v The vector with the appropriate direction
	@param [t:function?] normFunc The function to use as the norm; Defaults to the L2 norm

	@returns [t:(n x 1) matrix] The unit vector
**--]]
linalg.vector.unit = function(v, normFunc)
	if not normFunc then normFunc = linalg.vector.norm.l2 end
	return v / normFunc(v)
end

--[[**
	Creates a matrix that rotates a vector about an arbitrary vector
	Only works for 3 dimensions

	@param [t:(n x 1) matrix] v The vector to rotate about (should be a unit vector)
	@param [t:number] theta The angle to rotate by (in radians)

	@returns [t:nxn matrix] The resulting linear operator
**--]]
linalg.vector.createArbitraryAxisRotationMatrix = function(v, theta)
	if v.Shape[1] ~= 3 then
		error("Cannot create rotation matrix about arbitrary unit vector other than in R^3", 2)
	end

	local costheta = cos(theta)
	local sintheta = sin(theta)
	local versine = 1 - costheta
	local u = v.Unit

	return _newMatrix({
		{versine * u[1] ^ 2 + costheta, (versine * u[1] * u[2]) - (sintheta * u[3]), (versine * u[1] * u[3]) + (sintheta * u[2])},
		{(versine * u[1] * u[2]) + (sintheta * u[3]), versine * u[2] ^ 2 + costheta, (versine * u[2] * u[3]) - (sintheta * u[1])},
		{(versine * u[1] * u[3]) - (sintheta * u[2]), (versine * u[2] * u[3]) + (sintheta * u[1]), versine * u[3] ^ 2 + costheta}
	})
end

--[[ GENERAL FUNCTIONS ]]--
--[[**
	Creates an orthonormal basis for a dim-dimensional inner product space

	@param [t:array<(n x 1) matrix>] u The list of matrices to add to the basis (will be converted to unit vectors) (can be a single vector instead of an array)
	@param [t:number?] epsilon The minimum norm value for a vector to count to be added to the basis; Defaults to 0.01
	@param [t:function([(n x 1) matrix], [(n x 1) matrix])?] innerProductFunc The inner product function to use; Defaults to the dot product
	@param [t:function([(n x 1) matrix])?] normFunc The norm function to use; Defaults to the L2 norm

	@returns [t:array<(n x 1) matrix>] An orthonormal basis that includes the unit vectors of the original u
**--]]
linalg.gramSchmidt = function (u, epsilon, dim, innerProductFunc, normFunc)
	if #u == 0 then -- u must be just a single vector
		u = { u }
	end
	if not epsilon then
		epsilon = 0.01
	end
	if not dim then
		dim = u[1].Shape[1]
	end
	if not innerProductFunc then
		innerProductFunc = linalg.vector.ip.dot
	end
	if not normFunc then
		normFunc = linalg.vector.norm.l2
	end

	-- Normalize all elements in u
	local newU = {}
	for i = 1, #u do
		newU[i] = u[i].Unit
	end
	u = newU

	for i = 1, dim do
		local eiRows = {}
		for j = 1, dim do
			eiRows[j] = { i == j and 1 or 0 }
		end
		table.insert(u, _newMatrix(eiRows))
	end

	local basisVectors = { u[1].Unit }

	for i = 1, #u do
		local potentialBasisVector = u[i] - linalg.vector.project(u[i], basisVectors, innerProductFunc, normFunc)
		if normFunc(potentialBasisVector) >= epsilon then
			table.insert(basisVectors, potentialBasisVector.Unit)
		end

		if #basisVectors == dim then
			break
		end
	end

	return basisVectors
end

return linalg
