查看“︁Module:Complex Number/Matrix”︁的源代码

local p = {}
local comp_lib = require("Module:Complex_Number")
local sol_lib = require("Module:Complex_Number/Solver")
local cmath = comp_lib.cmath.init()
local toCnumber = cmath.constructor
p.cmath = cmath
p.toCnumber = toCnumber
function p._isNaN(x)
	return (not (x==x)) and  (x~=x)
end

--[[
	p.mmath.matrix(p.mmath.row(1,5),p.mmath.row(7,9))
	p.mmath.matrix(p.mmath.row(1,3),p.mmath.row(2,1))
	p.mmath.matrix(p.mmath.row(1, 2, 3), p.mmath.row(5, 7, 6), p.mmath.row(7, 9, 8))

	p.determinant(p.mmath.matrix(p.mmath.row(1,5),p.mmath.row(7,9)))
	p.determinant(p.mmath.matrix(p.mmath.row(1, 2, 3), p.mmath.row(5, 7, 6), p.mmath.row(7, 9, 8)))
	p.mmath.determinant(p.mmath.matrix(p.mmath.row(1, 2, 3, 0), p.mmath.row(5, 7, 6, 2), p.mmath.row(7, 9, 8, 1), p.mmath.row(2, 1, 3, 9)))
	p.mmath.inverse(p.mmath.matrix(p.mmath.row(1, 2, 3, 0), p.mmath.row(5, 7, 6, 2), p.mmath.row(7, 9, 8, 1), p.mmath.row(2, 1, 3, 9)))
]]
p.mmath={
	rows = function(mat)return #mat end,
	cols = function(mat)
		local max_val = -1
		for i = 1,#mat do if #(mat[i]) > max_val then max_val = #(mat[i]) end end
		return max_val
	end,
	determinant = function(matrix)return p._determinant(matrix, p.mmath.rows(matrix))end,
	adj = function(matrix)return p.mmath.adjoint(matrix)end,
	cof = function(matrix, _c_p, _c_q)return p.mmath.cofactor(matrix, _c_p, _c_q)end,
	det = function(matrix)return p._determinant(matrix, p.mmath.rows(matrix))end,
	mathform = function(mat)
		for i = 1,p.mmath.rows(mat) do mat[i].mathform = true end
		mat.mathform = true
		return mat
	end,
	abs = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.abs(toCnumber(matrix)) end
		local sum = 0
		local rows = p.mmath.rows(matrix)
		local cols = p.mmath.cols(matrix)
		for i = 1,rows do 
			local tmp_row = matrix[i]or{}
			for j = 1,cols do 
				local item = tmp_row[j]or 0
				sum = sum + item * item
			end
		end
		return cmath.sqrt(sum)
	end,

	floor = function(matrix)return p._applyMathFunc(cmath.floor, matrix)end,
	ceil = function(matrix)return p._applyMathFunc(cmath.ceil, matrix)end,
	round = function(op1, op2, op3)return p._applyMathFunc(function(matrix)return cmath.round(matrix, op2, op3)end, op1)end,
	div = function(op1, op2)
		if p.mmath.isMatrix(op1) then
			return op1 / op2
		else
			local scalar = toCnumber(op1)
			local diver = p.mmath.inverse(op2)
			return diver * scalar
		end
	end,

	re = function(matrix)return p._applyMathFunc(cmath.re, matrix)end,
	im = function(matrix)return p._applyMathFunc(cmath.im, matrix)end,
	conjugate = function(matrix)return p.mmath.transpose(p._applyMathFunc(cmath.conjugate, matrix))end,
	inverse = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.inverse(toCnumber(matrix)) end
		local N = p.mmath.rows(matrix)
		local M = p.mmath.cols(matrix)
		local dim = math.max(N, M)
		
		if dim <= 2 then
			if N~=M then error("in inverse() :Matrix Error: only square matrix can inverse") end
			if dim==0 then return p.mmath.matrix(p.mmath.row(1))
			elseif dim==1 then 
				local det = matrix[1][1]or 0
			    if toCnumber(det) == cmath.zero or det == 0 then
			    	error("in inverse() :Matrix Error: Singular matrix, can't find its inverse")
			    end
				return p.mmath.matrix(p.mmath.row(cmath.inverse(det)))
			elseif dim==2 then
				local a, b, c, d = matrix[1][1], matrix[1][2], matrix[2][1], matrix[2][2]
				local det = a*d - b*c
			    if toCnumber(det) == cmath.zero or det == 0 then
			    	error("in inverse() :Matrix Error: Singular matrix, can't find its inverse")
			    end
			    local idet = cmath.inverse(det)
				return p.mmath.matrix(
					p.mmath.row(d*idet, -b*idet),
					p.mmath.row(-c*idet, a*idet))
			end
		else
			local result = {}
		    -- Find determinant of matrix[][]
		    local det = p._determinant(matrix, dim);
		    if toCnumber(det) == cmath.zero or det == 0 then
		    	error("in inverse() :Matrix Error: Singular matrix, can't find its inverse")
		    end
		 
		    -- Find adjoint
		    local adj = p.mmath.adjoint(matrix)
		 
		    -- Find Inverse using formula "inverse(matrix) = adj(matrix)/det(matrix)"
			for i = 1,dim do
				local res_row = {}
				for j = 1,dim do
					res_row[#res_row + 1] = ((adj[i]or{})[j]or 0) / det
				end
				result[#result + 1] = p.mmath.row(unpack(res_row))
			end
			result = p.mmath.matrix(unpack(result)) 
		    return result
		end
	end,
	trunc = function(op1, op2)return p._applyMathFunc(function(matrix)return cmath.trunc(matrix, op2)end, op1)end,
	digits = function(matrix)return p._applyMathFunc(cmath.digits, matrix)end,
	sqrt = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.sqrt(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.sqrt, matrix)
	end,
	sin = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.sin(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.sin, matrix)
	end,
	cos = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.cos(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.cos, matrix)
	end,
	tan = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.tan(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.tan, matrix)
	end,
	cot = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.cot(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.cot, matrix)
	end,

	asin = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.asin(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.asin, matrix)
	end,
	acos = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.acos(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.acos, matrix)
	end,
	atan = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.atan(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.atan, matrix)
	end,
	acot = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.acot(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.acot, matrix)
	end,
	
	sinh = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.sinh(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.sinh, matrix)
	end,
	cosh = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.cosh(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.cosh, matrix)
	end,
	tanh = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.tanh(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.tanh, matrix)
	end,
	coth = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.coth(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.coth, matrix)
	end,
	
	asinh = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.asinh(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.asinh, matrix)
	end,
	acosh = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.acosh(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.acosh, matrix)
	end,
	atanh = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.atanh(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.atanh, matrix)
	end,
	acoth = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.acoth(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.acoth, matrix)
	end,
	
	dot = function(op1, op2)
		local op1_is_mat = p.mmath.isMatrix(op1)
		local op2_is_mat = p.mmath.isMatrix(op2)
		local dot_scalar = false
		local dot_matrix = op1
		if (not op1_is_mat) and (not op2_is_mat) then
			return cmath.dot(toCnumber(op1), toCnumber(op2))
		elseif not (op1_is_mat and op2_is_mat) then
			if op1_is_mat then
				dot_matrix = op1
				dot_scalar = op2
			else
				dot_matrix = op2
				dot_scalar = op1
			end
		end
		local sum = 0
		local rows = p.mmath.rows(dot_matrix)
		local cols = p.mmath.cols(dot_matrix)
		for i = 1,rows do 
			local tmp_row = dot_matrix[i]or{}
			for j = 1,cols do 
				local item = tmp_row[j]or 0
				local dot_num = dot_scalar
				if dot_num == false then dot_num = (op2[i]or{})[j]or 0 end
				sum = sum + item * dot_num
			end
		end
		return sum
	end,
	sgn = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.sgn(toCnumber(matrix)) end
		return matrix / p.mmath.abs(matrix)
	end,
	arg = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.arg(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.arg, matrix)
	end,
	cis = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.cis(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.cis, matrix)
	end,
	exp = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.exp(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.exp, matrix)
	end,
	elog = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.elog(toCnumber(matrix)) end
		return p._applyMatrixFunction(cmath.elog, matrix)
	end,
	log = function(matrix, basez)
		if basez ~= nil then
			if not p.mmath.isMatrix(matrix) and not p.mmath.isMatrix(basez) then return cmath.log(toCnumber(matrix), toCnumber(basez)) end
			if not p.mmath.isMatrix(matrix) and p.mmath.isMatrix(basez) then
				return p._applyMatrixFunction(function(base_z)return cmath.log(toCnumber(matrix), toCnumber(base_z)) end, basez)
			elseif p.mmath.isMatrix(matrix) and not p.mmath.isMatrix(basez) then
				return p._applyMatrixFunction(function(mat)return cmath.log(toCnumber(mat), toCnumber(basez)) end, matrix)
			end
			return p.mmath.elog(basez) * p.mmath.inverse(p.mmath.elog(matrix))
		else
			if not p.mmath.isMatrix(matrix) then return cmath.log(toCnumber(matrix)) end
			return p._applyMatrixFunction(cmath.log, matrix)
		end
	end,
	pow = function(op1, op2)
		if not p.mmath.isMatrix(op1) and not p.mmath.isMatrix(op2) then return cmath.pow(toCnumber(op1), toCnumber(op2)) end
		if not p.mmath.isMatrix(op1) and p.mmath.isMatrix(op2) then
			return p._applyMatrixFunction(function(op_2)return cmath.pow(toCnumber(op1), toCnumber(op_2)) end, op2)
		elseif p.mmath.isMatrix(op1) and not p.mmath.isMatrix(op2) then
			local dim = p.mmath.rows(op1)
			if dim <= 3 then
				return p._applyMatrixFunction(function(op_1)return cmath.pow(toCnumber(op_1), toCnumber(op2)) end, op1)
			else
				local times = tonumber(tostring(op2))or-1
				if times == 0 then return p.mmath.identity(dim)
				elseif times == 1 then return p.mmath.clone(op1)
				elseif times > 1 then
					local result_mat = p.mmath.identity(dim)
					for i=1,times do
						result_mat = result_mat * op1
					end
					return result_mat
				end
			end
		end
		local a, z = op1, op2 --a ^ z
		local log_a = p._applyMatrixFunction(cmath.elog, a)
		return p._applyMatrixFunction(cmath.exp, log_a * z)
	end,
	
	random = function(op1, op2)
		if not p.mmath.isMatrix(op1) and not p.mmath.isMatrix(op2) then return cmath.random(op1, op2) end
		local mat1, mat2 = op1, op2
		if not p.mmath.isMatrix(op1) or not p.mmath.isMatrix(op2) then
			if p.mmath.isMatrix(op1) then
				mat1, mat2 = op1, op2
			else
				mat1, mat2 = op2, op1
			end
		end
		local res_mat = {}
		local rows = p.mmath.rows(mat1)
		local cols = p.mmath.cols(mat1)
		for i = 1,rows do
			local tmp_row = mat1[i]or{}
			local res_row = {}
			for j = 1,cols do
				local mat2_cell = p.mmath.isMatrix(mat2) and (mat2[i][j] or 0) or mat2
				local cell_str = tostring(mat2_cell):lower()
				if cell_str=="nan" or cell_str=="-nan" then 
					res_row[#res_row + 1] = cmath.random()
				else
					res_row[#res_row + 1] = cmath.random(tmp_row[j]or 0, mat2_cell)
				end
			end
			res_mat[#res_mat + 1] = p.mmath.row(unpack(res_row))
		end
		return p.mmath.matrix(unpack(res_mat)) 
	end,
	MatrixFunction = function(func, matrix)
		if not p.mmath.isMatrix(matrix) then return func(toCnumber(matrix)) end
		return p._applyMatrixFunction(func, matrix)
	end,
	--sqrt,sin,cos...etc need to solve eigenvalue and eigenvector. dimension more then 3x3 not support here.
	GaussElimination = function(mat)
		if not p.mmath.isMatrix(mat) then return mat end
		return p._applyGaussElimination(mat) 
	end,
	rank = function(mat)
		if not p.mmath.isMatrix(mat) then return (mat == 0 or mat == cmath.zero) and 0 or 1 end
		local elimination = p._applyGaussElimination(mat)
		local erows = p.mmath.rows(elimination)
		local ecols = p.mmath.cols(elimination)
		local rank_val = 0
		for i = 1,erows do 
			for j = 1,ecols do 
				if elimination[i][j] ~= cmath.zero then
					rank_val = rank_val + 1
					break
				end
			end
		end
		return rank_val
	end,
	transpose = function(mat)
		if not p.mmath.isMatrix(mat) then return mat end
		local res_mat = {}
		for i = 1,p.mmath.cols(mat) do
			local res_row = {}
			for j = 1,p.mmath.rows(mat) do
				res_row[#res_row + 1] = (mat[j]or 0)[i]or 0
			end
			res_mat[#res_mat + 1] = p.mmath.row(unpack(res_row))
			res_mat[#res_mat].mathform = mat.mathform
		end
		local result_mat = p.mmath.matrix(unpack(res_mat)) 
		result_mat.mathform = mat.mathform
		return result_mat
	end,
	clone = function(mat)
		if not p.mmath.isMatrix(mat) then return mat end
		local res_mat = {}
		for i = 1,p.mmath.rows(mat) do
			local res_row = {}
			for j = 1,p.mmath.cols(mat) do
				res_row[#res_row + 1] = (mat[i]or 0)[j]or 0
			end
			res_mat[#res_mat + 1] = p.mmath.row(unpack(res_row))
			res_mat[#res_mat].mathform = mat.mathform
		end
		local result_mat = p.mmath.matrix(unpack(res_mat)) 
		result_mat.mathform = mat.mathform
		return result_mat
	end,
	identity = function(num)
		if p.mmath.isMatrix(num) then return error("in identity() :Argument Exception: dimension can not be matrix") end
		local size = tonumber(tostring(num)) or 1
		local res_mat = {}
		for i = 1,size do
			local res_row = {}
			for j = 1,size do
				res_row[#res_row + 1] = (i==j)and 1 or 0
			end
			res_mat[#res_mat + 1] = p.mmath.row(unpack(res_row))
		end
		return p.mmath.matrix(unpack(res_mat)) 
	end,
	diag = function(...)
		local diag_items = {...}
		local vec_checker = diag_items[1]
		if type(vec_checker) == type({})then
			if #vec_checker == 1 then
				return p._diagonal(unpack(vec_checker[1]))
			elseif #vec_checker > 1 then
				local res_row = {}
				for i = 1,#vec_checker do
					res_row[#res_row + 1] = (type(vec_checker[i])==type({}))and vec_checker[i][1]or vec_checker[i]
				end
				return p._diagonal(unpack(res_row))
			end
		end
		return p._diagonal(...)
	end,
	adjoint = function(matrix)
		if not p.mmath.isMatrix(matrix) then return cmath.one end
		local M = p.mmath.rows(matrix)
		local N = p.mmath.cols(matrix)
		local result = {}
		if M <= 1 and N <= 1 then
		    --result[0][0] = 1;
		    return p.mmath.matrix(p.mmath.row(math.max(M,N)))
		end
		
		for i = 1,M do
			local res_row = {}
			for j = 1,N do
				res_row[#res_row + 1] = 0
			end
			result[#result + 1] = p.mmath.row(unpack(res_row))
		end
		result = p.mmath.matrix(unpack(result)) 
		
		-- temp is used to store cofactors of matrix[][]
		local temp = {}
		local sign = 1
	 
		for i=1,M do
			for j=1,N do
				-- Get cofactor of matrix[i][j]
				temp = p.mmath.cofactor(matrix, i, j);
				
				-- sign of result[j][i] positive if sum of row
				-- and column indexes is even.
				sign = ((i+j)%2==0)and 1 or -1
				
				-- Interchanging rows and columns to get the
				-- transpose of the cofactor matrix
				result[j][i] = (sign)*(p._determinant(temp, M-1));
			end
		end
		return result
	end,
	cofactor = function(matrix, _c_p, _c_q)
		local _p, _q = tonumber(tostring(_c_p)), tonumber(tostring(_c_q))
		if _p==nil or _q==nil then
			return p.mmath.transpose(p.mmath.adjoint(matrix))
		end
		local result = {}
		local m = p.mmath.rows(matrix)
		local n = p.mmath.cols(matrix)
	    local i, j = 1, 1;
		for i = 2,m do
			local res_row = {}
			for j = 2,n do
				res_row[#res_row + 1] = 0
			end
			result[#result + 1] = p.mmath.row(unpack(res_row))
		end
		result = p.mmath.matrix(unpack(result)) 
	    -- Looping for each element of the matrix
		for row = 1,m do
			for col = 1,n do
				--  Copying into temporary matrix only those element
				--  which are not in given row and column
				if row ~= _p and col ~= _q then
					result[i][j] = (matrix[row]or{})[col]or 0
					j = j + 1
					-- Row is filled, so increase row index and
					-- reset col index
					if j == m then
						j = 1;
						i = i + 1
					end
				end
			end
		end
		return result
	end,
	eigenvalue=function(matrix, it)
		local i = tonumber(tostring(it))
		local eigenvalues =  p._eigenvalues(matrix)
		if it~=nil then
			if eigenvalues[i]~=nil then
				return eigenvalues[i]
			end
		end
		return eigenvalues
	end,
	eigenvector=function(matrix, it)
		local i = tonumber(tostring(it))
		local eigenvector =  p._eigenvectors(matrix)
		if it~=nil then
			if eigenvector[i]~=nil then
				return p.mmath.vector(unpack(eigenvector[i]))
			end
		end
		return eigenvector
	end,
	matrixMeta = {
		__add = function (op1, op2)
			local result_obj = {}
			local all_count = math.max(p.mmath.rows(op1), p.mmath.rows(op2))
			for i = 1,all_count do result_obj[i] = (op1[i]or p.mmath.row(0))+(op2[i]or p.mmath.row(0))end

			local result_mat = p.mmath.matrix(unpack(result_obj)) 
			result_mat.mathform = op1.mathform or op2.mathform
			return result_mat
		end,
		__sub = function (op1, op2)
			local result_obj = {}
			local all_count = math.max(p.mmath.rows(op1), p.mmath.rows(op2))
			for i = 1,all_count do result_obj[i] = (op1[i]or p.mmath.row(0))-(op2[i]or p.mmath.row(0))end
			local result_mat = p.mmath.matrix(unpack(result_obj)) 
			result_mat.mathform = op1.mathform or op2.mathform
			return result_mat
		end,
		__mul = function (op1, op2)
			local scalar_flag = false
			local can_mul = true
			if type(op2)==type({})then
				if type(op2.objType)==type("string") then
					if op2.objType~='matrix' then can_mul = false end
				else
					scalar_flag = true
				end
			elseif type(op2)==type(tostring)then
				can_mul = false
			else
				scalar_flag = true
			end
			if not can_mul then error( "in mul() :Argument Exception: Can not multiply" ) end
			local res_mat = {}
			if scalar_flag then
				local rows1 = p.mmath.rows(op1)
				local cols1 = p.mmath.cols(op1)
				for i = 1,rows1 do
					local ri = op1[i]or{}
					local res_row = {}
					for j = 1,cols1 do
						res_row[#res_row + 1] = op2 * (ri[j] or 0)
					end
					res_mat[#res_mat + 1] = p.mmath.row(unpack(res_row))
					res_mat[#res_mat].mathform = op1.mathform or op2.mathform
				end
			else
				local _m = p.mmath.rows(op1)
				local _n = p.mmath.cols(op1)
				local _p = p.mmath.cols(op2)
				if _n ~= p.mmath.rows(op2) then error( "in mul(): Matrix Error: dimension mismatch" ) end

				for _i = 1,_m do 
					local res_row = {}
					for _j = 1,_p do
						local sum = 0
						for _r = 1,_n do 
							sum = sum + ((op1[_i]or{})[_r]or 0) * ((op2[_r]or{})[_j]or 0)
						end
						res_row[#res_row + 1] = sum
					end
					res_mat[#res_mat + 1] = p.mmath.row(unpack(res_row))
					res_mat[#res_mat].mathform = op1.mathform or op2.mathform
				end
			end
			local result_mat = p.mmath.matrix(unpack(res_mat)) 
			result_mat.mathform = op1.mathform or op2.mathform
			return result_mat
		end,
		__div = function (op1, op2)
			local scalar_flag = false
			local can_mul = true
			if type(op2)==type({})then
				if type(op2.objType)==type("string") then
					if op2.objType~='matrix' then can_mul = false end
				else
					scalar_flag = true
				end
			elseif type(op2)==type(tostring)then
				can_mul = false
			else
				scalar_flag = true
			end
			if not can_mul then error( "in div() :Argument Exception: Can not divide" ) end
			local res_mat = {}
			if scalar_flag then
				local op2_inv = cmath.inverse(op2)
				local rows1 = p.mmath.rows(op1)
				local cols1 = p.mmath.cols(op1)
				for i = 1,rows1 do
					local ri = op1[i]or{}
					local res_row = {}
					for j = 1,cols1 do
						res_row[#res_row + 1] = op2_inv * (ri[j] or 0)
					end
					res_mat[#res_mat + 1] = p.mmath.row(unpack(res_row))
					res_mat[#res_mat].mathform = op1.mathform or op2.mathform
				end
				local result_mat = p.mmath.matrix(unpack(res_mat)) 
				result_mat.mathform = op1.mathform or op2.mathform
				return result_mat
			else
				local op2_inv = p.mmath.inverse(op2)
				return op1 * op2_inv
			end
		end,
		__mod = function (op1, op2)
			local scalar_flag = false
			local can_mul = true
			if type(op2)==type({})then
				if type(op2.objType)==type("string") then
					if op2.objType~='matrix' then can_mul = false end
				else
					scalar_flag = true
				end
			elseif type(op2)==type(tostring)then
				can_mul = false
			else
				scalar_flag = true
			end
			if not can_mul then error( "in mod() :Argument Exception: Can not divide" ) end
			local res_mat = {}
			if scalar_flag then
				local diver = toCnumber(op2)
				return p._applyMathFunc(function(_A) return _A % diver end, op1)
			else
				return op1 - p.mmath.floor(op1 / op2) * op2
			end
		end,
		__tostring = function (this)
			local ret = ''
			for i = 1,#this do 
				ret = ret..(ret==''and''or(this.mathform and' \\\\\n'or','))..tostring(this[i])
			end
			if this.mathform then return'\\begin{bmatrix}'..ret..'\\end{bmatrix}'
			else return'{'..ret..'}'end
		end,
		__unm = function (this)return p._applyMathFunc(function(_A) return -_A end, this)end,
		__eq = function (op1, op2)
			local result_flag = true
			local all_count = math.max(p.mmath.rows(op1), p.mmath.rows(op2))
			for i = 1,all_count do result_flag = result_flag and ((op1[i]or 0)==(op2[i]or 0))end
			return result_flag
		end,
	},
	rowMeta = {
		__add = function (op1, op2)
			local result_obj = {}
			local all_count = math.max(#op1, #op2)
			for i = 1,all_count do result_obj[i] = (op1[i]or 0)+(op2[i]or 0)end
			return p.mmath.row(unpack(result_obj)) 
		end,
		__sub = function (op1, op2)
			local result_obj = {}
			local all_count = math.max(#op1, #op2)
			for i = 1,all_count do result_obj[i] = (op1[i]or 0)-(op2[i]or 0)end
			return p.mmath.row(unpack(result_obj)) 
		end,
		__mul = function (op1, op2)
			local result_obj = {}
			local all_count = math.max(#op1, #op2)
			for i = 1,all_count do result_obj[i] = (op1[i]or 0)*(op2[i]or 0)end
			return p.mmath.row(unpack(result_obj)) 
		end,
		__div = function (op1, op2)
			local result_obj = {}
			local all_count = math.max(#op1, #op2)
			for i = 1,all_count do result_obj[i] = (op1[i]or 0)/(op2[i]or 0)end
			return p.mmath.row(unpack(result_obj)) 
		end,
		__mod = function (op1, op2)
			local result_obj = {}
			local all_count = math.max(#op1, #op2)
			for i = 1,all_count do result_obj[i] = (op1[i]or 0)%(op2[i]or 0)end
			return p.mmath.row(unpack(result_obj)) 
		end,
		__unm = function (this)
			local result_obj = {}
			for i = 1,#this do result_obj[i] = -this[i]end
			return p.mmath.row(unpack(result_obj)) 
		end,
		__eq = function (op1, op2)
			local flag = true
			local all_count = math.max(#op1, #op2)
			for i = 1,all_count do flag = flag and ((op1[i]or 0)==(op2[i]or 0))end
			return flag
		end,
		__tostring = function (this)
			local ret = ''
			for i = 1,#this do 
				local to_text = tostring(this[i])
				if this.mathform then to_text='{'..to_text..'}'end
				ret = ret..(ret==''and''or(this.mathform and' & 'or','))..to_text
			end
			if this.mathform then return ret
			else return'{'..ret..'}'end
		end,
	},
	row = function(...)
		local row_obj = {...}
		for i=1,#row_obj do row_obj[i] = toCnumber(row_obj[i]) end
		setmetatable(row_obj,p.mmath.rowMeta)
		row_obj.objType = 'row'
		row_obj.numberType = "matrix"
		row_obj.mathform = false
		return row_obj
	end,
	vector = function(...)
		local vec_items = {...}
		local size = #vec_items
		local res_mat = {}
		for i = 1,size do res_mat[#res_mat + 1] = p.mmath.row(vec_items[i])end
		return p.mmath.matrix(unpack(res_mat))
	end,
	rowvector = function(...)return p.mmath.matrix(p.mmath.row(...))end,
	matrix = function(...)
		local col_obj = {...}
		setmetatable(col_obj, p.mmath.matrixMeta)
		col_obj.objType = 'matrix'
		col_obj.numberType = "matrix"
		col_obj.mathform = false
		return col_obj
	end,
	isMatrix = function(input_data)
		if type(input_data) == type({}) then
			if type(input_data.objType)==type('objType') then
				return input_data
			end
		end
		return false
	end,
	toMatrix = function(input_data)
		local matrix_checker = tostring(input_data)
		local matrix_flag, pre_matrix = xpcall(function()
			if mw.ustring.find(matrix_checker, "%{")then
				if mw.ustring.find(matrix_checker, "%\\begin") then
					local text = ''..matrix_checker
					text = mw.ustring.gsub(text, "[%s\n]+", " "):gsub("\\%w+%{%w+%}",""):gsub("[%{%}]+","")
					local rowdatas = mw.text.split(text, "\\\\")
					local load_matrix = {}
					for i=1,#rowdatas do
						rowdatas[i] = mw.text.split(rowdatas[i] , "&")
						for j=1,#(rowdatas[i]) do
							rowdatas[i][j] = toCnumber(mw.text.trim(rowdatas[i][j]))or cmath.zero
						end
						load_matrix[#load_matrix + 1] = p.mmath.row(unpack(rowdatas[i]))
						load_matrix[#load_matrix].mathform = true
					end
					load_matrix = p.mmath.matrix(unpack(load_matrix)) 
					load_matrix.mathform = true
					return load_matrix
				else
					local text = ''..matrix_checker
					text = "},"..mw.ustring.gsub(text,"%{%s*%{","{"):gsub("%}%s*%}","}")..",{"

					local rowdatas_raw = mw.text.split(text, "%}%s*,%s*{")
					local rowdatas = {}
					for _,sub_text in ipairs(rowdatas_raw) do
						local it_text = mw.text.trim(sub_text)
						if it_text ~= "" then rowdatas[#rowdatas+1]=it_text end
					end
					local load_matrix = {}
					for i=1,#rowdatas do
						rowdatas[i] = mw.text.split(rowdatas[i] , ",")
						for j=1,#(rowdatas[i]) do
							rowdatas[i][j] = toCnumber(mw.text.trim(rowdatas[i][j]))or cmath.zero
						end
						load_matrix[#load_matrix + 1] = p.mmath.row(unpack(rowdatas[i]))
					end
					load_matrix = p.mmath.matrix(unpack(load_matrix))
					return load_matrix
				end
				return nil
			end
		end, function()end)
		if matrix_flag then
			if (pre_matrix or{}).objType == 'matrix' then
				return pre_matrix
			end
		end
		if type(input_data) == type({}) then
			if type(input_data.objType)==type('objType') then
				return input_data
			end
		end
		return toCnumber(input_data)
	end,
	init = function()
		p.mmath.e = cmath.e
		p.mmath.pi = cmath.pi
		p.mmath["π"] = cmath["π"] 
		p.mmath["°"] = cmath["°"]
		p.mmath.nan = cmath.nan
		p.mmath.i = cmath.i
		
		p.mmath.numberType = comp_lib._numberType
		p.mmath.constructor = p.mmath.toMatrix
		return p.mmath
	end
}

-- ================ Eigens ================
--[[
p.mmath.eigenvalue(p.mmath.toMatrix("{{1, 2, 4, 3}, {5, 7, 6, 8}, {10, 9, 12, 11}, {13, 16, 14, 15}}"))
p.mmath.eigenvector(p.mmath.toMatrix("{{1, 2, 4, 3}, {5, 7, 6, 8}, {10, 9, 12, 11}, {13, 16, 14, 15}}"))
{35.4422, 2.72101, -2.05098, -1.11227}
{{0.199913, 0.471587, 0.73557, 1.}, {-0.665669, 1.20938, -1.64109, 1.}, 
{-0.79522, -0.522759, 0.117929, 1.}, {0.230581, -0.767239, -0.488143, 1.}}
]]
function p._eigenvectors(matrix)
	local N = p.mmath.rows(matrix)
	local M = p.mmath.cols(matrix)
	if N~=M then error("only square matrix has eigenvectors") end
	if N==1 then
		return p.mmath.matrix(p.mmath.row(1))
	elseif N==2 then
		local a, b, c, d = matrix[1][1], matrix[1][2], matrix[2][1], matrix[2][2]
		if (toCnumber(b) == cmath.zero or b==0) and (toCnumber(c) == cmath.zero or c==0) then 
			return p.mmath.matrix(p.mmath.row(1,0),p.mmath.row(0,1)) 
		end
		return p.mmath.matrix(p.mmath.row(-(-a+d+cmath.sqrt(a*a + 4*b*c-2*a*d+d*d))/(2*c),1),
							p.mmath.row(-(-a+d-cmath.sqrt(a*a + 4*b*c-2*a*d+d*d))/(2*c),1))
	elseif N==3 then
		local a = matrix
		if	(toCnumber(a[1][2]) == cmath.zero or a[1][2]==0) and (toCnumber(a[1][3]) == cmath.zero or a[1][3]==0) and
			(toCnumber(a[2][1]) == cmath.zero or a[2][1]==0) and (toCnumber(a[2][3]) == cmath.zero or a[2][3]==0) and
			(toCnumber(a[3][1]) == cmath.zero or a[3][1]==0) and (toCnumber(a[3][2]) == cmath.zero or a[3][2]==0) then 
				return p.mmath.matrix(p.mmath.row(1,0,0),p.mmath.row(0,1,0),p.mmath.row(0,0,1)) 
		end
		local cbRoot = sol_lib._cubicRoot(1, -a[1][1] - a[2][2] - a[3][3], 
						-a[1][2]*a[2][1] + a[1][1]*a[2][2] - a[1][3]*a[3][1] - a[2][3]*a[3][2] + a[1][1]*a[3][3] + a[2][2]*a[3][3], 
						a[1][3]*a[2][2]*a[3][1] - a[1][2]*a[2][3]*a[3][1] - a[1][3]*a[2][1]*a[3][2] + a[1][1]*a[2][3]*a[3][2] + a[1][2]*a[2][1]*a[3][3] - a[1][1]*a[2][2]*a[3][3])
		return	p.mmath.matrix(p.mmath.row(-(1/a[3][1])*(a[3][3] - cbRoot[1]) + (a[3][2]*(-a[2][3]*a[3][1] + a[2][1]*a[3][3] - a[2][1]*cbRoot[1]))/(a[3][1]*(-a[2][2]*a[3][1] + a[2][1]*a[3][2] + a[3][1]*cbRoot[1])), -((-a[2][3]*a[3][1] + a[2][1]*a[3][3] - a[2][1]*cbRoot[1])/(-a[2][2]*a[3][1] + a[2][1]*a[3][2] + a[3][1]*cbRoot[1])), 1), 
				p.mmath.row(-(1/a[3][1])*(a[3][3] - cbRoot[2]) + (a[3][2]*(-a[2][3]*a[3][1] + a[2][1]*a[3][3] - a[2][1]*cbRoot[2]))/(a[3][1]*(-a[2][2]*a[3][1] + a[2][1]*a[3][2] + a[3][1]*cbRoot[2])), -((-a[2][3]*a[3][1] + a[2][1]*a[3][3] - a[2][1]*cbRoot[2])/(-a[2][2]*a[3][1] + a[2][1]*a[3][2] + a[3][1]*cbRoot[2])), 1), 
				p.mmath.row(-(1/a[3][1])*(a[3][3] - cbRoot[3]) + (a[3][2]*(-a[2][3]*a[3][1] + a[2][1]*a[3][3] - a[2][1]*cbRoot[3]))/(a[3][1]*(-a[2][2]*a[3][1] + a[2][1]*a[3][2] + a[3][1]*cbRoot[3])), -((-a[2][3]*a[3][1] + a[2][1]*a[3][3] - a[2][1]*cbRoot[3])/(-a[2][2]*a[3][1] + a[2][1]*a[3][2] + a[3][1]*cbRoot[3])), 1))
	else
		local check_diag = true
		local should_break = false
		local diag_eigen = {}
		for i = 1,N do
			local res_row = {}
			for j = 1,M do
				if i~=j then
					res_row[#res_row + 1] = 0
					if toCnumber(matrix[i][j]) ~= cmath.zero and matrix[i][j] ~= 0 then
						check_diag = false
						should_break = true
						break
					end
				else res_row[#res_row + 1] = 1
				end
			end
			diag_eigen[#diag_eigen + 1] = p.mmath.row(unpack(res_row))
			if should_break then break end
		end
		if check_diag then 
			diag_eigen = p.mmath.matrix(unpack(diag_eigen))
			return diag_eigen
		end
		
		local result_matrix = {noMatrix=true}
		if N==4 then
			local eigenVec = sol_lib._quarticEigenData(matrix)
			local result_eigen = {}
			for i = 1,N do
				local res_row = {}
				for j = 1,M do
					res_row[#res_row + 1] = eigenVec[i][j]
				end
				result_eigen[#result_eigen + 1] = p.mmath.row(unpack(res_row))
			end
			result_matrix = p.mmath.matrix(unpack(result_eigen))
		else
			error("in eigenvectors(): Not Implemented Exception : dimension = "..N.."x"..M.." not support.")
		end
		if not result_matrix.noMatrix then
			mw.addWarning("WARNING: calculate eigenvector of matrix with dimension = "..N.."x"..M.." will be large numerical deviations.")
			return result_matrix
		end
	end
end
function p._eigenvalues(matrix)
	local N = p.mmath.rows(matrix)
	local M = p.mmath.cols(matrix)
	if N~=M then error("in eigenvalues(): Matrix Error: only square matrix has eigenvectors") end
	if N==1 then
		return p.mmath.matrix(p.mmath.row(matrix[1][1]))
	elseif N==2 then
		local a, b, c, d = matrix[1][1], matrix[1][2], matrix[2][1], matrix[2][2]
		if b==0 and c==0 then return p.mmath.row(a,d) end
		return p.mmath.row((a+d-cmath.sqrt(a*a+4*b*c-2*a*d+d*d))/2,
							(a+d+cmath.sqrt(a*a+4*b*c-2*a*d+d*d))/2)
	elseif N==3 then
		local a = matrix
		if a[1][2]==0 and a[1][3]==0 and a[2][1]==0 and a[2][3]==0 and a[3][1]==0 and a[3][2]==0 then 
			return p.mmath.row(a[1][1],a[2][2],a[3][3])
		end
		local cbRoot = sol_lib._cubicRoot(1, -a[1][1] - a[2][2] - a[3][3], 
						-a[1][2]*a[2][1] + a[1][1]*a[2][2] - a[1][3]*a[3][1] - a[2][3]*a[3][2] + a[1][1]*a[3][3] + a[2][2]*a[3][3], 
						a[1][3]*a[2][2]*a[3][1] - a[1][2]*a[2][3]*a[3][1] - a[1][3]*a[2][1]*a[3][2] + a[1][1]*a[2][3]*a[3][2] + a[1][2]*a[2][1]*a[3][3] - a[1][1]*a[2][2]*a[3][3])
		return p.mmath.row(cbRoot[1],cbRoot[2],cbRoot[3])
	else
		local check_diag = true
		local should_break = false
		local diag_eigen = {}
		for i = 1,N do
			for j = 1,M do
				if i~=j then
					if toCnumber(matrix[i][j]) ~= cmath.zero and matrix[i][j] ~= 0 then
						check_diag = false
						should_break = true
						break
					end
				else diag_eigen[#diag_eigen + 1] = matrix[i][j]
				end
			end
			if should_break then break end
		end
		if check_diag then 
			diag_eigen = p.mmath.row(unpack(diag_eigen))
			return diag_eigen
		end
		if N==4 then
			local eigenRoot = sol_lib._quarticEigenRoot(matrix)
			return p.mmath.row(eigenRoot[1],eigenRoot[2],eigenRoot[3],eigenRoot[4])
		else
			error("in eigenvalues(): Not Implemented Exception : dimension = "..N.."x"..M.." not support.")
		end
	end
end

-- ================ Matrix Functions ================
function p._applyMatrixFunction(func, matrix)
	local EV = p._eigenvalues(matrix)
	local FD = {}
	local size = #EV
	for i = 1,size do
		local res_row = {}
		for j = 1,size do
			if i==j then
				res_row[#res_row + 1] = func(toCnumber(EV[i]))
			else
				res_row[#res_row + 1] = 0
			end
		end
		FD[#FD + 1] = p.mmath.row(unpack(res_row))
	end
	FD = p.mmath.matrix(unpack(FD))
	local P = p.mmath.transpose(p._eigenvectors(matrix))
	local P_inv = p.mmath.inverse(P)
	return P * FD * P_inv
end

function p._applyMathFunc(func, matrix)
	if not p.mmath.isMatrix(matrix) then return func(toCnumber(matrix)) end
	local res_mat = {}
	local rows = p.mmath.rows(matrix)
	local cols = p.mmath.cols(matrix)
	for i = 1,rows do
		local tmp_row = matrix[i]or{}
		local res_row = {}
		for j = 1,cols do
			res_row[#res_row + 1] = func(tmp_row[j]or 0)
		end
		res_mat[#res_mat + 1] = p.mmath.row(unpack(res_row))
	end
	return p.mmath.matrix(unpack(res_mat)) 
end

-- ================ Matrix Operations ================
function p._swap_row(mat, i, j)
	local m = p.mmath.cols(mat)
    --mw.log("Swapped rows ", i, " and ", j);
 
    for k=1,m do
        local temp = mat[i][k];
        mat[i][k] = mat[j][k];
        mat[j][k] = temp;
    end
end

function p._applyGaussElimination(matrix)
	local N = p.mmath.rows(matrix)
	local M = p.mmath.cols(matrix)
	local mat = {}
	for i = 1,N do
		local res_row = {}
		for j = 1,M do
			res_row[#res_row + 1] = (matrix[i]or{})[j]or 0
		end
		mat[#mat + 1] = p.mmath.row(unpack(res_row))
	end
	mat = p.mmath.matrix(unpack(mat))
	
    for k=1,N do
        -- Initialize maximum value and index for pivot
        local i_max = k;
        local v_max = mat[i_max][k];
 
        -- find greater amplitude for pivot if any
        for i = k+1,N do
            if cmath.abs(mat[i][k]) > cmath.abs(v_max) then v_max, i_max = mat[i][k], i end
        end
	
        --[[ if a principal diagonal element  is zero,
         * it denotes that matrix is singular, and
         * will lead to a division-by-zero later. ]]
        --if mat[k][i_max] == 0 then return k end -- Matrix is singular
 
        -- Swap the greatest value row with current row
        if i_max ~= k then p._swap_row(mat, k, i_max) end
 
        for i=k+1,N do
            --[[ factor f to set current row kth element to 0,
             * and subsequently remaining kth column to 0 ]]
	            local f = mat[i][k]/mat[k][k];
	 
	            --[[ subtract fth multiple of corresponding kth
	               row element]]
	            for j=k+1,N do 
	            	if mat[k][k] ~= 0 and mat[k][k] ~= cmath.zero then
	            		mat[i][j] = mat[i][j] - mat[k][j] * f 
	            	end
	            end
	 
	            -- filling lower triangular matrix with zeros
	            mat[i][k] = cmath.zero;

        end
 
        --mw.log(mat);        --for matrix state
    end
    --mw.log(mat);            --for matrix state
	for i = 1,N do 
		for j = 1,M do 
			local check_cNaN = tostring(mat[i][j])
			if check_cNaN == "nan" or check_cNaN == "-nan" or check_cNaN == "inf" or check_cNaN == "-inf" then mat[i][j] = cmath.zero end 
		end
	end
    return mat;
end

function p._diagonal(...)
	local diag_items = {...}
	local size = #diag_items
	local res_mat = {}
	for i = 1,size do
		local res_row = {}
		for j = 1,size do
			res_row[#res_row + 1] = (i==j)and (diag_items[i]or 0) or 0
		end
		res_mat[#res_mat + 1] = p.mmath.row(unpack(res_row))
	end
	return p.mmath.matrix(unpack(res_mat))
end

function p._determinant(matrix, n)
	local det = 0;
	local submatrix={}
	if n == 0 then return 1 
	elseif n == 1 then return matrix[1][1] 
	elseif n == 2 then return ((matrix[1][1] * matrix[2][2]) - (matrix[2][1] * matrix[1][2]));
	elseif n == 3 then
		local _x = ((matrix[2][2] * matrix[3][3]) - (matrix[3][2] * matrix[2][3]));
		local _y = ((matrix[2][1] * matrix[3][3]) - (matrix[3][1] * matrix[2][3]));
		local _z = ((matrix[2][1] * matrix[3][2]) - (matrix[3][1] * matrix[2][2]));
		return ((matrix[1][1] * _x) - (matrix[1][2] * _y) + (matrix[1][3] * _z));
	else 
		for x = 1,n do
			local subi = 1;
			for i = 2,n do
				local subj = 1;
				local res_row = {}
				for j = 1,n do
					--if (j == x) continue;
					if j ~= x then
						--submatrix[subi][subj] = matrix[i][j];
						res_row[subj] = (matrix[i]or{})[j]or 0
						subj=subj+1
					end
				end
				submatrix[subi] = p.mmath.row(unpack(res_row))
				subi=subi+1
			end
			submatrix = p.mmath.matrix(unpack(submatrix)) 
			det = det + (math.pow(-1, (x-1) ) * matrix[1][x] * p._determinant( submatrix, n - 1 ));
		end
   end
   return det;
end

return p