% The	lualinalg	package.
% Authors:  Chetan  Shirore  and  Ajit  Kumar
% Version 1.9.
% Licensed under LaTeX Project Public License v1.3c or later. The complete license text is available at http://www.latex-project.org/lppl.txt.

\ProvidesPackage{lualinalg}[1.9]
\RequirePackage{xkeyval}
\RequirePackage{amsmath}
\RequirePackage{luamaths}
\RequirePackage{luacode}
\begin{luacode*}

-- matrices part
matrices = {} -- global registry to store matrices.

matrix = {} --module

local matrix_meta = {}
-- Adding functions to the matrix module.
function matrix.new(matrix, rows, columns, str)
   if type(rows) == "table" then
      for i = 1, #rows do
         if #rows[1] ~= #rows[i] then
            error("Check input matrix.")
         end
      end
      return setmetatable(rows, matrix_meta)
   end
   local mtx = {}
   if columns == "I" then
      for i = 1, rows do
         mtx[i] = {}
         for j = 1, rows do
            if i == j then
               mtx[i][j] = 1
            else
               mtx[i][j] = 0
            end
         end
      end
      return setmetatable(mtx, matrix_meta)
   end
   if str == "zero" then
      for i = 1, rows do
         mtx[i] = {}
         for j = 1, columns do
			mtx[i][j] = 0
         end
      end
      return setmetatable(mtx, matrix_meta)
   end
end

setmetatable(
matrix,
{__call = function(...)
return matrix.new(...)
end}
)

function matrix.add(m1, m2)
   local mtx = {}
   for i = 1, #m1 do
      local m3i = {}
      mtx[i] = m3i
      for j = 1, #m1[1] do
         m3i[j] = m1[i][j] + m2[i][j]
      end
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.sub(m1, m2)
   local mtx = {}
   for i = 1, #m1 do
      local m3i = {}
      mtx[i] = m3i
      for j = 1, #m1[1] do
         m3i[j] = m1[i][j] - m2[i][j]
      end
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.mulnum(m1, num)
   local mtx = {}
   -- multiply elements with number
   for i = 1, #m1 do
      mtx[i] = {}
      for j = 1, #m1[1] do
         mtx[i][j] = m1[i][j] * num
      end
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.mul(m1, m2)
   local mtx = {}
   for i = 1, #m1 do
      mtx[i] = {}
      for j = 1, #m2[1] do
         local num = m1[i][1] * m2[1][j]
         for n = 2, #m1[1] do
            num = num + m1[i][n] * m2[n][j]
         end
         mtx[i][j] = num
      end
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.swapRows(m1, p, q)
   local mtx = {}
   for i = 1, #m1 do
      mtx[i] = {}
      for j = 1, #m1[1] do
         mtx[i][j] = m1[i][j]
      end
   end
   for j = 1, #m1[1] do
      rowHold = m1[p][j]
      mtx[p][j] = m1[q][j]
      mtx[q][j] = rowHold
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.swapCols(m1, p, q)
   local mtx = {}
   for i = 1, #m1 do
      mtx[i] = {}
      for j = 1, #m1[1] do
         mtx[i][j] = m1[i][j]
      end
   end
   for j = 1, #m1 do
      rowHold = m1[j][p]
      mtx[j][p] = m1[j][q]
      mtx[j][q] = rowHold
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.mulRow(m1, p, k)
   local mtx = {}
   for i = 1, #m1 do
      mtx[i] = {}
      for j = 1, #m1[1] do
         mtx[i][j] = m1[i][j]
      end
   end
   for j = 1, #m1[1] do
      mtx[p][j] = k * m1[p][j]
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.mulAddRow(m1, k, p, q)
   if p == q then
      error("Can't operate on same row.")
   end
   local mtx = {}
   for i = 1, #m1 do
      mtx[i] = {}
      for j = 1, #m1[1] do
         mtx[i][j] = m1[i][j]
      end
   end
   for j = 1, #m1[1] do
      mtx[q][j] = k * (mtx[p][j]) + mtx[q][j]
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.mulCol(m1, p, k)
   local mtx = {}
   for i = 1, #m1 do
      mtx[i] = {}
      for j = 1, #m1[1] do
         mtx[i][j] = m1[i][j]
      end
   end
   for j = 1, #m1 do
      mtx[j][p] = k * m1[j][p]
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.mulAddCol(m1, k, p, q)
   if p == q then
      error("Can't operate on same column.")
   end
   local mtx = {}
   for i = 1, #m1 do
      mtx[i] = {}
      for j = 1, #m1[1] do
         mtx[i][j] = m1[i][j]
      end
   end
   for j = 1, #m1 do
      mtx[j][q] = k * mtx[j][p] + mtx[j][q]
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.transpose(m1)
   local mtx = {}
   for i = 1, #m1[1] do
      mtx[i] = {}
      for j = 1, #m1 do
         mtx[i][j] = m1[j][i]
      end
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.subm(m1, i1, j1, i2, j2)
   local mtx = {}
   for i = i1, i2 do
      local _i = i - i1 + 1
      mtx[_i] = {}
      for j = j1, j2 do
         local _j = j - j1 + 1
         mtx[_i][_j] = m1[i][j]
      end
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.concath(m1, m2)
   if #m1 ~= #m2 then
      error("No. of rows must be equal.")
   end
   local mtx = {}
   local offset = #m1[1]
   for i = 1, #m1 do
      mtx[i] = {}
      for j = 1, offset do
         mtx[i][j] = m1[i][j]
      end
      for j = 1, #m2[1] do
         mtx[i][j + offset] = m2[i][j]
      end
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.concatv(m1, m2)
   if #m1[1] ~= #m2[1] then
      error("No. of columns must be equal.")
   end
   local mtx = {}
   for i = 1, #m1 do
      mtx[i] = {}
      for j = 1, #m1[1] do
         mtx[i][j] = m1[i][j]
      end
   end
   local offset = #mtx
   for i = 1, #m2 do
      local _i = i + offset
      mtx[_i] = {}
      for j = 1, #m2[1] do
         mtx[_i][j] = m2[i][j]
      end
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.rows(mtx)
   return #mtx
end

function matrix.columns(mtx)
   return #mtx[1]
end

setmetatable(
matrix,
{__call = function(...)
return matrix.new(...)
end}
)

function matrix.getelement(mtx, i, j)
   if mtx[i] and mtx[i][j] then
      return mtx[i][j]
   end
end

function matrix.setelement(mtx, i, j, value)
   if matrix.getelement(mtx, i, j) then
      mtx[i][j] = value
      return value
   end
end

function matrix.invert(m1)
   if #m1 ~= #m1[1] then
      error("matrix not square")
   end
   if matrix.det(m1) == 0 then
      error("matrix not invertible")
   end
   local mtx = {}
   local idnt = matrix(#m1, "I")
   mtx = matrix.subm(matrix.rref(matrix.concath(m1, idnt)), 1, #m1 + 1, #m1, #m1 + #m1)
   return mtx
end

function matrix.trace(m1)
   if #m1 ~= #m1[1] then
      error("matrix not square")
   end
   local sum = 0

   for i = 1, #m1 do
      for j = 1, #m1[1] do
         if i == j then
            sum = sum + m1[i][j]
         end
      end
   end

   return sum
end

function matrix.normF(mtx)
   local result = 0
   for i = 1, #mtx do
      for j = 1, #mtx[1] do
         local e = mtx[i][j]
         result = result + complex.abs(complex(e)) ^ 2
      end
   end
   return complex.sqrt(complex(result))
end

function matrix.normmax(mtx)
   local result = 0
   for i = 1, #mtx do
      for j = 1, #mtx[1] do
         local e = complex.abs(complex(mtx[i][j]))
         if e > result then
            result = e
         end
      end
   end
   return result
end

function matrix.norminfty(mtx)
   local e = 0
   local result = 0
   for i = 1, #mtx do
      local e = 0
      for j = 1, #mtx[1] do
         e = e + complex.abs(complex(mtx[i][j]))
      end
      if e > result then
         result = e
      end
   end
   return result
end

function matrix.norm1(mtx)
   local e = 0
   local result = 0
   for i = 1, #mtx[1] do
      local e = 0
      for j = 1, #mtx do
         e = e + complex.abs(complex(mtx[j][i]))
      end
      if e > result then
         result = e
      end
   end
   return result
end

function matrix.conjugate(m1)
   local mtx = matrix.copy(m1)
   for i = 1, #mtx do
      for j = 1, #mtx[1] do
         mtx[i][j] = complex.conjugate(complex(mtx[i][j]))
      end
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.conjugateT(m1)
   local mtx = {}
   for i = 1, #m1[1] do
      mtx[i] = {}
      for j = 1, #m1 do
         mtx[i][j] = complex.conjugate(complex(m1[j][i]))
      end
   end
   return setmetatable(mtx, matrix_meta)
end

function copy(x)
   return type(x) == "table" and x.copy(x) or x
end

function matrix.pow(m1, num)
   assert(num == math.floor(num), "exponent not an integer")
   if num == 0 then
      return matrix:new(#m1, "I")
   end
   if num < 0 then
      local rank
      m1, rank = matrix.invert(m1)
      if not m1 then
         return m1, rank
      end -- singular
      num = -num
   end
   local mtx = matrix.copy(m1)
   for i = 2, num do
      mtx = matrix.mul(mtx, m1)
   end
   return mtx
end

function matrix.createrandom(nrow, ncol, start, stop)
   mtx = {}
   for i = 1, nrow do
      mtx[i] = {}
      for j = 1, ncol do
         mtx[i][j] = math.random(start, stop)
         mtx[i][j] = mtx[i][j] + math.min(math.random(), math.abs(mtx[i][j] - start), math.abs(stop - mtx[i][j]))
      end
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.process(m1)
   --m1=load("return "..m1)()
   return matrix.mulnum(m1, 1)
end

function matrix.det(m1)
   assert(#m1 == #m1[1], "matrix not square")

   local size = #m1

   if size == 1 then
      return m1[1][1]
   end

   if size == 2 then
      return m1[1][1] * m1[2][2] - m1[2][1] * m1[1][2]
   end

   if size == 3 then
      return (m1[1][1] * m1[2][2] * m1[3][3] + m1[1][2] * m1[2][3] * m1[3][1] + m1[1][3] * m1[2][1] * m1[3][2] -
      m1[1][3] * m1[2][2] * m1[3][1] -
      m1[1][1] * m1[2][3] * m1[3][2] -
      m1[1][2] * m1[2][1] * m1[3][3])
   end

   local e = m1[1][1]
   local zero = type(e) == "table" and e.zero or 0
   local norm2 = type(e) == "table" and e.norm2 or number_norm2

   local mtx = matrix.copy(m1)
   local det = 1

   for j = 1, #mtx[1] do
      local rows = #mtx
      local subdet, xrow
      for i = 1, rows do
         local e = mtx[i][j]
         if not subdet then
            if e ~= zero then
               subdet, xrow = e, i
            end
         elseif e ~= zero and math.abs(norm2(e) - 1) < math.abs(norm2(subdet) - 1) then
            subdet, xrow = e, i
         end
      end
      if subdet then
         if xrow ~= rows then
            mtx[rows], mtx[xrow] = mtx[xrow], mtx[rows]
            det = -det
         end

         for i = 1, rows - 1 do
            if mtx[i][j] ~= zero then
               local factor = mtx[i][j] / subdet
               for n = j + 1, #mtx[1] do
                  mtx[i][n] = mtx[i][n] - factor * mtx[rows][n]
               end
            end
         end
         if math.fmod(rows, 2) == 0 then
            det = -det
         end
         det = det * subdet
         table.remove(mtx)
      else
         return det * 0
      end
   end

   return det
end

function matrix.copy(m1)
   local mtx = {}
   for i = 1, #m1 do
      mtx[i] = {}
      for j = 1, #m1[1] do
         mtx[i][j] = m1[i][j]
      end
   end
   return setmetatable(mtx, matrix_meta)
end

norm2 = type(e) == "table" and e.norm2 or number_norm2

function number_norm2(x)
   return x * x
end

function matrix.op(exp)
   return load("return " .. exp, exp, "t", matrices)()
end

function matrix.rref(mtx)
   local mtx = matrix.copy(mtx)
   step = 1
   lead = 1
   rowCount = #mtx
   columnCount = #mtx[1]
   for r = 1, rowCount do
      if lead > columnCount then
         return mtx
      end
      i = r
      while (mtx[i][lead] == 0) do
         i = i + 1
         if (i - 1 == rowCount) then
            i = r
            if (columnCount == lead) then
               return mtx
            end
            lead = lead + 1
         end
      end

      if i ~= r then
         mtx = matrix.swapRows(mtx, i, r)
      end

      local m = mtx[r][lead]
      if (mtx[r][lead] ~= 0) then
         for u = 1, columnCount do
            mtx[r][u] = mtx[r][u] / m
         end
      end
      for i = 1, rowCount do
         local m = mtx[i][lead]
         if (i ~= r) then
            for v = 1, columnCount do
               mtx[i][v] = mtx[i][v] + ((-m) * (mtx[r][v]))
            end
         end
      end
      lead = lead + 1
   end
   return mtx
end

function matrix.rref0E(mtx, fom, dignum)
   local strng = ""
   truncate = truncate or 6
   local mtx = matrix.copy(mtx)
   step = 1
   lead = 1
   stepCnt = 0
   rowCount = #mtx
   columnCount = #mtx[1]
   for r = 1, rowCount do
      if lead > columnCount then
         return mtx
      end
      i = r
      while (mtx[i][lead] == 0) do
         i = i + 1
         if (i - 1 == rowCount) then
            i = r
            if (columnCount == lead) then
               if stepCnt == 0 then
                  stepCnt = stepCnt + 1
                  strng = strng .. "Step " .. tostring(stepCnt) ".$$" .. tostring(matrix.show(mtx, fom, dignum))
                  return strng
               end
               return strng
            end
            lead = lead + 1
         end
      end

      if i ~= r then
         mtx = matrix.swapRows(mtx, i, r)
         stepCnt = stepCnt + 1
         strng =
         strng ..
         "Step " ..
         tostring(stepCnt) ..
         ": Interchange rows " ..
         tostring(i) ..
         " and " .. tostring(r) .. ".$$" .. tostring(matrix.show(mtx, fom, dignum)) .. "$$"
      end

      local m = mtx[r][lead]
      if (mtx[r][lead] ~= 0) then
         for u = 1, columnCount do
            mtx[r][u] = mtx[r][u] / m
         end
         if m ~= 1.0 then
            if m ~= complex("1.0") then
               stepCnt = stepCnt + 1
               strng =
               strng ..
               "Step " ..
               tostring(stepCnt) ..
               ": Divide row " ..
               tostring(r) ..
               " by $" ..
               tostring(complex.round(complex(m), dignum)) ..
               "$.$$" .. tostring(matrix.show(mtx, fom, dignum)) .. "$$"
            end
         end
      end
      for i = 1, rowCount do
         local m = mtx[i][lead]
         if (i ~= r) then
            for v = 1, columnCount do
               mtx[i][v] = mtx[i][v] + ((-m) * (mtx[r][v]))
            end
            if m ~= 0 then
               if m ~= complex("0.0") then
                  stepCnt = stepCnt + 1
                  strng =
                  strng ..
                  "Step " ..
                  tostring(stepCnt) ..
                  ": Multiply row " ..
                  tostring(r) ..
                  " by $" ..
                  tostring(complex.round(complex(m), dignum)) ..
                  "$ and subtract it from row " ..
                  tostring(i) ..
                  ".$$" .. tostring(matrix.show(mtx, fom, dignum)) .. "$$"
               end
            end
         end
      end
      lead = lead + 1
   end
   return strng
end

function matrix.GaussJordan(mtx, augmt)
   local mtx = matrix.copy(mtx)
   local augmt = matrix.copy(augmt)
   step = 1
   lead = 1
   rowCount = #mtx
   columnCount = #mtx[1]
   for r = 1, rowCount do
      if lead > columnCount then
         return matrix.concath(mtx, augmt)
      end
      i = r
      while (mtx[i][lead] == 0) do
         i = i + 1
         if (i - 1 == rowCount) then
            i = r
            if (columnCount == lead) then
               return matrix.concath(mtx, augmt)
            end
            lead = lead + 1
         end
      end

      if i ~= r then
         mtx = matrix.swapRows(mtx, i, r)
         augmt = matrix.swapRows(augmt, i, r)
      end

      local m = mtx[r][lead]
      if (mtx[r][lead] ~= 0) then
         for u = 1, columnCount do
            mtx[r][u] = mtx[r][u] / m
         end
         augmt[r][1] = augmt[r][1] / m
      end
      for i = 1, rowCount do
         local m = mtx[i][lead]
         if (i ~= r) then
            for v = 1, columnCount do
               mtx[i][v] = mtx[i][v] - m * mtx[r][v]
            end
            augmt[i][1] = augmt[i][1] - m * augmt[r][1]
         end
      end
      lead = lead + 1
   end
   return matrix.concath(mtx, augmt)
end

function matrix.gauss0E(mtx, augmt, fom, dignum)
   local strng = ""
   truncate = truncate or 6
   local mtx = matrix.copy(mtx)
   local augmt = matrix.copy(augmt)
   if matrix.columns(augmt) ~= 1 then
      error("The second matrix should have only 1 column.")
   end
   step = 1
   lead = 1
   stepCnt = 0
   rowCount = #mtx
   columnCount = #mtx[1]
   for r = 1, rowCount do
      if lead > columnCount then
         return mtx
      end
      i = r
      while (mtx[i][lead] == 0) do
         i = i + 1
         if (i - 1 == rowCount) then
            i = r
            if (columnCount == lead) then
               if stepCnt == 0 then
                  stepCnt = stepCnt + 1
                  strng =
                  strng ..
                  "Step " ..
                  tostring(stepCnt) ".$$" ..
                  tostring(matrix.show(matrix.concath(mtx, augmt), fom, dignum))
                  return strng
               end
               return strng
            end
            lead = lead + 1
         end
      end

      if i ~= r then
         mtx = matrix.swapRows(mtx, i, r)
         augmt = matrix.swapRows(augmt, i, r)
         stepCnt = stepCnt + 1
         strng =
         strng ..
         "Step " ..
         tostring(stepCnt) ..
         ": Interchange rows " ..
         tostring(i) ..
         " and " .. tostring(r) .. ".$$" .. tostring(matrix.show(mtx, fom, dignum)) .. "$$"
      end

      local m = mtx[r][lead]
      if (mtx[r][lead] ~= 0) then
         for u = 1, columnCount do
            mtx[r][u] = mtx[r][u] / m
         end
         augmt[r][1] = augmt[r][1] / m
         if m ~= 1.0 then
            if m ~= complex("1.0") then
               stepCnt = stepCnt + 1
               strng =
               strng ..
               "Step " ..
               tostring(stepCnt) ..
               ": Divide row " ..
               tostring(r) ..
               " by $" ..
               tostring(complex.round(complex(m), dignum)) ..
               "$. $$" ..
               tostring(matrix.show(matrix.concath(mtx, augmt), fom, dignum)) ..
               "$$"
            end
         end
      end
      for i = 1, rowCount do
         local m = mtx[i][lead]
         if (i ~= r) then
            for v = 1, columnCount do
               mtx[i][v] = mtx[i][v] - m * mtx[r][v]
            end
            augmt[i][1] = augmt[i][1] - m * augmt[r][1]
            if m ~= 0 then
               if m ~= complex("0.0") then
                  stepCnt = stepCnt + 1
                  strng =
                  strng ..
                  "Step " ..
                  tostring(stepCnt) ..
                  ": Multiply row " ..
                  tostring(r) ..
                  " by $" ..
                  tostring(complex.round(complex(m), dignum)) ..
                  "$ and subtract it from row " ..
                  tostring(i) ..
                  ".$$" ..
                  tostring(
                  matrix.show(matrix.concath(mtx, augmt), fom, dignum)
                  ) ..
                  "$$"
               end
            end
         end
      end
      lead = lead + 1
   end
   return strng
end

function matrix.rank(m1)
   local mtx = {}
   mtx = matrix.rref(m1)
   rank = #mtx
   for i = 1, #mtx do
      if CheckEqual(mtx[i], 0) then
         rank = rank - 1
      end
   end
   return rank
end

function CheckEqual(Values, Number)
   local CheckEqual = true
   local i = 1

   while (CheckEqual and (i <= #Values)) do
      if Values[i] == Number then
         i = i + 1
      else
         CheckEqual = false
      end
   end

   return CheckEqual
end

function matrix.replace(m1, func, ...)
   local mtx = {}
   for i = 1, #m1 do
      local m1i = m1[i]
      local mtxi = {}
      for j = 1, #m1i do
         mtxi[j] = func(m1i[j], ...)
      end
      mtx[i] = mtxi
   end
   return setmetatable(mtx, matrix_meta)
end

function matrix.show(mtx, format, dig)
   mtx = matrix.process(mtx)
   local format = format or "bmatrix"
   local dig = dig or 6
   local str = "\\begin{" .. format .. "}"
   for i = 1, #mtx do
      str = str .. "\t" .. complex.round(complex(mtx[i][1]), dig)
      for j = 2, #mtx[1] do
         str = str .. " & " .. complex.round(complex(mtx[i][j]), dig)
      end
      if i == #mtx then
         str = str .. " \\\\ "
      else
         str = str .. " \\\\ "
      end
   end
   return str .. "\\end{" .. format .. "} "
end

function matrix.chqeql(m1, m2) 
   if #m1 ~= #m2 or #m1[1] ~= #m2[1] then
   return false
   end
   for i = 1, #m1 do
         for j = 1, #m1[1] do
			if not(lnumChqEql(m1[i][j],m2[i][j])) then
			return false
			end
         end
   end
   return true
end

-- Setting Meta-operations in the matrix module.
matrix_meta.__tostring = function(...)
   return matrix.show(...)
end

matrix_meta.__add = function(...)
   return matrix.add(...)
end

matrix_meta.__sub = function(...)
   return matrix.sub(...)
end

matrix_meta.__mul = function(m1, m2)
   if getmetatable(m1) ~= matrix_meta then
   return matrix.mulnum(m2, m1)
   elseif getmetatable(m2) ~= matrix_meta then
   return matrix.mulnum(m1, m2)
   end
   return matrix.mul(m1, m2)
end

matrix_meta.__div = function(m1, m2)
   if getmetatable(m1) ~= matrix_meta then
   return matrix.mulnum(matrix.invert(m2), m1)
   elseif getmetatable(m2) ~= matrix_meta then
   return matrix.divnum(m1, m2)
end
   return matrix.div(m1, m2)
end

matrix_meta.__unm = function(mtx)
   return matrix.mulnum(mtx, -1)
end

matrix_meta.__eq = function(...)
   return matrix.chqeql(...)
end

local option = {
   ["*"] = function(m1)
   return matrix.conjugate(m1)
end,
["T"] = function(m1)
   return matrix.transpose(m1)
end
}
matrix_meta.__pow = function(m1, opt)
   return option[opt] and option[opt](m1) or matrix.pow(m1, opt)
end

-- vector part

vectors = {} -- global registry to store vectors.

vector = {} --module.

local vector_meta = {}
-- Adding functions to the vector module.
function vector.new(vector, rows, columns, n)
    if columns ~= "e" and columns ~= "zero" then
        local tbl = {}
        for i = 1, #rows do
            tbl[i] = rows[i]
        end
        return setmetatable(tbl, vector_meta)
    end
    local vec = {}
    if columns == "e" then
        for i = 1, rows do
            if i == n then
                vec[i] = 1
            else
                vec[i] = 0
            end
        end
        return setmetatable(vec, vector_meta)
    end
	
	if columns == "zero" then
        for i = 1, rows do   
                vec[i] = 0
        end
        return setmetatable(vec, vector_meta)
    end
end

setmetatable(
vector,
{__call = function(...)
return vector.new(...)
end}
)

function vector.add(v1, v2)
    if #v1 ~= #v2 then
        return error("Vectors should be of same dimension.")
    end
    local vec = {}
    for i = 1, #v1 do
        vec[i] = v1[i] + v2[i]
    end
    return setmetatable(vec, vector_meta)
end

function vector.sub(v1, v2)
    if #v1 ~= #v2 then
        return error("Vectors should be of same dimension.")
    end
    local vec = {}
    for i = 1, #v1 do
        vec[i] = v1[i] - v2[i]
    end
    return setmetatable(vec, vector_meta)
end

function vector.dot(v1, v2)
    if #v1 ~= #v2 then
        return error("Vectors should be of same dimension")
    end
    local sum = 0
    for i = 1, #v1 do
        sum = sum + v1[i] * complex.conjugate(complex(v2[i]))
    end
    return sum
end

function vector.mulnum(v1, num)
    local vec = {}
    -- multiply elements with number
    for i = 1, #v1 do
        vec[i] = v1[i] * num
    end
    return setmetatable(vec, vector_meta)
end

function vector.sumnorm(v1)
    local norm = 0
    for i = 1, #v1 do
        norm = norm + complex.abs(complex(v1[i]))
    end
    return norm
end

function vector.euclidnorm(v1)
    return complex.sqrt(vector.dot(v1, v1))
end

function vector.pnorm(v1, p)
    if math.floor(p) ~= math.abs(p) or p <= 1 then
        return error("Invalid value of p")
    end
    local sum = 0
    for i = 1, #v1 do
        sum = sum + complex.abs(complex(v1[i])) ^ p
    end
    return sum ^ (1 / p)
end

function vector.supnorm(v1)
    local result = 0
    for i = 1, #v1 do
        local e = complex.abs(complex(v1[i]))
        if e > result then
            result = e
        end
    end
    return result
end

function vector.cross(v1, v2)
    if #v1 ~= 3 or #v2 ~= 3 then
        return error("Vectors should be of dimension 3")
    end
    local vec = {}
    vec[1] = v1[2] * v2[3] - v1[3] * v2[2]
    vec[2] = v1[3] * v2[1] - v1[1] * v2[3]
    vec[3] = v1[1] * v2[2] - v1[2] * v2[1]
    return setmetatable(vec, vector_meta)
end

function vector.createrandom(n, start, stop)
    start = start or 0
    stop = stop or 10
    vec = {}
    for i = 1, n do
        vec[i] = math.random(start, stop)
        vec[i] = vec[i] + math.min(math.random(), math.abs(vec[i] - start), math.abs(stop - vec[i]))
    end
    return setmetatable(vec, vector_meta)
end

function vector.getcoordinate(vec, i)
    if vec[i] then
        return vec[i]
    end
end

function vector.setcoordinate(vec, i, val)
    if vec[i] then
        vec[i] = val
        return val
    end
end

function vector.getangle(v1, v2)
    if #v1 ~= #v2 then
        return error("Vectors should be of same dimension")
    end
    local x = complex.get(vector.dot(v1, v2) / (vector.euclidnorm(v1) * vector.euclidnorm(v2)))
    return math.acos(mathround(x, 15))
end

function vector.copy(v1)
    local vec = {}
    for i = 1, #v1 do
        vec[i] = v1[i]
    end
    return setmetatable(vec, vector_meta)
end

function vector.op(exp)
    return load("return " .. exp, exp, "t", vectors)()
end

function vector.process(v1)
    --v1=load("return "..v1)()
    return vector.mulnum(v1, 1)
end

function vector.show(vec, dig)
    vec = vector.process(vec)
    local dig = dig or 6
    local str = ""
    for i = 1, #vec do
        if i == 1 then
            str = str .. complex.round(complex(vec[i]), dig)
        else
            str = str .. "," .. complex.round(complex(vec[i]), dig)
        end
    end
    return str .. ""
end

function vector.parse(vec)
    local tbl = {}
    for i = 1, #vec do
        tbl[i] = vec[i]
    end
    return "(" .. table.concat(tbl, ",") .. ")"
end

function vector.gs(inptTbl, brckt, dignum)
    local brcktR = ""
    brckt = brckt or "round"
    if brckt == "round" then
        brcktL = "("
        brcktR = ")"
    end
    if brckt == "square" then
        brcktL = "["
        brcktR = "]"
    end
    if brckt == "curly" then
        brcktL = "\\{"
        brcktR = "\\}"
    end

    local tbl = {}
    local str = ""
    k = #inptTbl

    if complex.round(vector.euclidnorm(inptTbl[1])) ~= complex("0.0") then
        tbl[1] = vector.mulnum(inptTbl[1], 1 / vector.euclidnorm(inptTbl[1]))
    else
        tbl[1] = vector.mulnum(inptTbl[1], 1)
    end
    setmetatable(tbl[1], vector_meta)
    str = str .. "$\\left" .. brcktL .. vector.show(tbl[1], dignum) .. "\\right" .. brcktR.." $"
    for i = 2, k do
        tbl[i] = inptTbl[i]
        setmetatable(tbl[i], vector_meta)
        for j = 1, i - 1 do
            setmetatable(tbl[j], vector_meta)
            tbl[i] = vector.sub(tbl[i], vector.mulnum(tbl[j], vector.dot(tbl[i], tbl[j])))
        end
        if complex.round(vector.euclidnorm(tbl[i])) ~= complex("0.0") then
            tbl[i] = vector.mulnum(tbl[i], 1 / vector.euclidnorm(tbl[i]))
        end
        tbl[i] = vector.mulnum(tbl[i], 1.0)
        str = str .. ", $\\left" .. brcktL .. vector.show(tbl[i], dignum) .. "\\right" .. brcktR.." $"
    end
    str = str
    return str
end

function vector.gsX(inptTbl, brckt, dignum)
    local brcktR = ""
    local cnt = 1
    brckt = brckt or "round"
    if brckt == "round" then
        brcktL = "\\left("
        brcktR = "\\right)$$"
    end
    if brckt == "square" then
        brcktL = "\\left["
        brcktR = "\\right]$$"
    end
    if brckt == "curly" then
        brcktL = "\\left\\{"
        brcktR = "\\right\\}$$"
    end

    local tbl = {}
    local tmpTbl = {}
    local str = ""
    k = #inptTbl
    str = str .. "\\ \\newline Take given vectors as $v_1,\\ldots, v_" .. k .. "$ in order."

    if complex.round(vector.euclidnorm(inptTbl[1])) ~= complex("0.0") then
        tbl[1] = vector.mulnum(inptTbl[1], 1.0 / vector.euclidnorm(inptTbl[1]))
    else
        tbl[1] = vector.mulnum(inptTbl[1], 1.0)
    end
    setmetatable(tbl[1], vector_meta)
    str = str .. "\\ \\newline Step " .. cnt .. ": $$ u_" .. cnt .. "=v_" .. cnt .. "="
    str = str .. brcktL .. vector.show(inptTbl[1], dignum) .. brcktR
    str = str .. " $$ e_" .. cnt .. "="
    if complex.round(vector.euclidnorm(tbl[1])) ~= complex("0.0") then
        str = str .. "\\frac{u_{" .. cnt .. "}}" .. "{||u_{" .. cnt .. "}||} ="
    end
    str = str .. brcktL .. vector.show(tbl[1], dignum) .. brcktR
    for i = 2, k do
        tbl[i] = inptTbl[i]
        setmetatable(tbl[i], vector_meta)
        for j = 1, i - 1 do
            setmetatable(tbl[j], vector_meta)
            tmpTbl[i] = vector.sub(tbl[i], vector.mulnum(tbl[j], vector.dot(tbl[i], tbl[j])))
            tbl[i] = vector.sub(tbl[i], vector.mulnum(tbl[j], vector.dot(tbl[i], tbl[j])))
        end
        if complex.round(vector.euclidnorm(tbl[i])) ~= complex("0.0") then
            tbl[i] = vector.mulnum(tbl[i], 1 / vector.euclidnorm(tbl[i]))
        else
            tbl[i] = vector.mulnum(tbl[i], 1.0)
        end

        cnt = cnt + 1
        str = str .. " Step " .. cnt
        str = str .. ": $$ u_" .. cnt .. "="
        str = str .. "v_" .. cnt .. "-\\sum_{j=1}^{" .. (cnt - 1) .. "}{{proj_{u_j}(v_" .. cnt .. ")}}="
        str = str .. brcktL .. vector.show(tmpTbl[i], dignum) .. brcktR
        str = str .. " $$ e_" .. cnt .. "="
        if vector.euclidnorm(tbl[i]) ~= complex(0.0) then
            str = str .. "\\frac{u_{" .. cnt .. "}}" .. "{||u_{" .. cnt .. "}||} ="
        end
        str = str .. brcktL .. vector.show(tbl[i], dignum) .. brcktR
    end

    return str
end

function vector.chqeql(v1, v2) 
   if #v1 ~= #v2 then
   return false
   end
	for j = 1, #v1 do
		if not(lnumChqEql(v1[j],v2[j])) then
		return false
		end
    end
   return true
end
-- Setting Meta-operations in the vector module.
vector_meta.__tostring = function(...)
    return vector.show(...)
end

vector_meta.__add = function(...)
    return vector.add(...)
end

vector_meta.__sub = function(...)
    return vector.sub(...)
end

vector_meta.__unm = function(vec)
    return vector.mulnum(vec, -1)
end

vector_meta.__eq = function(...)
    return vector.chqeql(...)
end

vector_meta.__mul = function(v1, v2)
    if getmetatable(v1) ~= vector_meta then
    return vector.mulnum(v2, v1)
    elseif getmetatable(v2) ~= vector_meta then
    return vector.mulnum(v1, v2)
    end
    return vector.dot(v1, v2)
end

function mathround(num, numDecimalPlaces)
    local mult = 10 ^ (numDecimalPlaces or 0)
    return math.floor(num * mult + 0.5) / mult
end


\end{luacode*}

% matrix latex commands

\newcommand\matrixNew[2]{%
\directlua{%
matrices['#1'] = matrix(#2)
}%
}

% ========= KEY DEFINITIONS =========
\define@key{matrixop}{type}{\def\mop@type{#1}}
\define@key{matrixop}{truncate}{\def\mop@truncate{#1}}
% ========= KEY DEFAULTS =========
\setkeys{matrixop}{type=bmatrix,truncate=6}%

\newcommand{\matrixPrint}[2][]{%
\begingroup%
\setkeys{matrixop}{#1}
\directlua{tex.sprint(matrix.show(matrices['#2'],"\mop@type",\mop@truncate))}
%
\endgroup%
}

\newcommand\matrixOp[2]{%
\directlua{%
matrices['#1'] = matrix.op('#2')
}%
}


\newcommand\matrixAdd[3]{%
\directlua{%
matrices['#1'] = matrix.add(matrices['#2'],matrices['#3'])
}%
}

\newcommand\matrixSub[3]{%
\directlua{%
matrices['#1'] = matrix.sub(matrices['#2'],matrices['#3'])
}%
}

\newcommand\matrixMulNum[3]{%
\directlua{%
matrices['#1'] = matrix.mulnum(matrices['#3'],#2)
}%
}

\newcommand\matrixMul[3]{%
\directlua{%
matrices['#1'] = matrix.mul(matrices['#2'],matrices['#3'])
}%
}

\newcommand\matrixSwapRows[4]{%
\directlua{%
matrices['#1'] = matrix.swapRows(matrices['#2'],#3,#4)
}%
}

\newcommand\matrixSwapCols[4]{%
\directlua{%
matrices['#1'] = matrix.swapCols(matrices['#2'],#3,#4)
}%
}

\newcommand\matrixMulRow[4]{%
\directlua{%
matrices['#1'] = matrix.mulRow(matrices['#2'],#3,#4)
}%
}

\newcommand\matrixMulCol[4]{%
\directlua{%
matrices['#1'] = matrix.mulCol(matrices['#2'],#3,#4)
}%
}

\newcommand\matrixMulAddRow[5]{%
\directlua{%
matrices['#1'] = matrix.mulAddRow(matrices['#2'],#4,#3,#5)
}%
}

\newcommand\matrixMulAddCol[5]{%
\directlua{%
matrices['#1'] = matrix.mulAddCol(matrices['#2'],#4,#3,#5)
}%
}

\newcommand\matrixTranspose[2]{%
\directlua{%
matrices['#1'] = matrix.transpose(matrices['#2'])
}%
}

\newcommand\matrixSubmatrix[6]{%
\directlua{%
matrices['#1'] = matrix.subm(matrices['#2'],#3,#4,#5,#6)
}%
}

\newcommand\matrixConcatH[3]{%
\directlua{%
matrices['#1'] = matrix.concath(matrices['#2'],matrices['#3'])
}%
}

\newcommand\matrixConcatV[3]{%
\directlua{%
matrices['#1'] = matrix.concatv(matrices['#2'],matrices['#3'])
}%
}

\newcommand\matrixNumRows[1]{%
\directlua{%
tex.sprint(tostring(matrix.rows(matrices['#1'])))
}%
}

\newcommand\matrixNumCols[1]{%
\directlua{%
tex.sprint(tostring(matrix.columns(matrices['#1'])))
}%
}

\newcommand\matrixGetElement[3]{%
\directlua{%
tex.sprint(tostring(matrix.getelement(matrices['#1'],#2,#3)))
}%
}

\newcommand\matrixSetElement[4]{%
\directlua{%
matrix.setelement(matrices['#1'],#2,#3,#4)
}%
}

\newcommand\matrixInvert[2]{%
\directlua{%
matrices['#1'] = matrix.invert(matrices['#2'])
}%
}

\newcommand\matrixPow[3]{%
\directlua{%
matrices['#1'] = matrix.pow(matrices['#2'],#3)
}%
}

\newcommand\matrixCreateRandom[5]{%
\directlua{%
matrices['#1'] = matrix.createrandom(#2,#3,#4,#5)
}%
}

\newcommand\matrixDet[1]{%
\directlua{%
tex.sprint(tostring(matrix.det(matrices['#1'])))
}%
}

\newcommand\matrixTrace[1]{%
\directlua{%
tex.sprint(tostring(matrix.trace(matrices['#1'])))
}%
}

\newcommand\matrixNormOne[1]{%
\directlua{%
tex.sprint(tostring(matrix.norm1(matrices['#1'])))
}%
}

\newcommand\matrixNormInfty[1]{%
\directlua{%
tex.sprint(tostring(matrix.norminfty(matrices['#1'])))
}%
}

\newcommand\matrixNormMax[1]{%
\directlua{%
tex.sprint(tostring(matrix.normmax(matrices['#1'])))
}%
}

\newcommand\matrixNormF[1]{%
\directlua{%
tex.sprint(tostring(matrix.normF(matrices['#1'])))
}%
}

\newcommand\matrixCopy[2]{%
\directlua{%
matrices['#1'] = matrix.copy(matrices['#2'])
}%
}

\newcommand\matrixRREF[2]{%
\directlua{%
matrices['#1'] = matrix.rref(matrices['#2'])
}%
}

\newcommand\matrixConjugate[2]{%
\directlua{%
matrices['#1'] = matrix.conjugate(matrices['#2'])
}%
}

\newcommand\matrixConjugateT[2]{%
\directlua{%
matrices['#1'] = matrix.conjugateT(matrices['#2'])
}%
}

\newcommand\matrixRank[1]{%
\directlua{%
tex.sprint(tostring(matrix.rank(matrices['#1'])))
}%
}

\newcommand\matrixEql[2]{%
\directlua{%
tex.sprint(tostring(matrix.chqeql(matrices['#1'],matrices['#2'])))
}%
}

% ========= KEY DEFINITIONS =========
\define@key{matrixrr}{type}{\def\moprr@type{#1}}
\define@key{matrixrr}{truncate}{\def\moprr@truncate{#1}}

% ========= KEY DEFAULTS =========
\setkeys{matrixrr}{type=bmatrix,truncate=6}%

\newcommand{\matrixRREFSteps}[2][]{%
\begingroup%
\setkeys{matrixrr}{#1}
\directlua{%
tex.sprint(matrix.rref0E(matrices['#2'],"\moprr@type",\moprr@truncate))}
%
\endgroup%
}

\newcommand\matrixGaussJordan[3]{%
\directlua{%
matrices['#1'] = matrix.GaussJordan(matrices['#2'],matrices['#3'])
}%
}

\newcommand{\matrixGaussJordanSteps}[3][]{%
\begingroup%
\setkeys{matrixrr}{#1}
\directlua{%
tex.sprint(matrix.gauss0E(matrices['#2'],matrices['#3'],"\moprr@type",\moprr@truncate))}
%
\endgroup%
}

% vector latex commands

\newcommand\vectorNew[2]{%
\directlua{%
vectors['#1'] = vector(#2)
}%
}

% ========= KEY DEFINITIONS =========
\define@key{vectorop}{truncate}{\def\vop@truncate{#1}}

% ========= KEY DEFAULTS =========
\setkeys{vectorop}{truncate=6}%

\newcommand{\vectorPrint}[2][]{%
\begingroup%
\setkeys{vectorop}{#1}
\directlua{tex.sprint(vector.show(vectors['#2'],\vop@truncate))}
%
\endgroup%
}

\newcommand\vectorParse[1]{%
\directlua{%
tex.sprint(tostring(vector.parse(vectors['#1'])))
}%
}

\newcommand\vectorOp[2]{%
\directlua{%
vectors['#1'] = vector.op('#2')
}%
}

\newcommand\vectorAdd[3]{%
\directlua{%
vectors['#1'] = vector.add(vectors['#2'],vectors['#3'])
}%
}

\newcommand\vectorSub[3]{%
\directlua{%
vectors['#1'] = vector.sub(vectors['#2'],vectors['#3'])
}%
}

\newcommand\vectorDot[2]{%
\directlua{%
tex.sprint(tostring(vector.dot(vectors['#1'],vectors['#2'])))
}%
}

\newcommand\vectorMulNum[3]{%
\directlua{%
vectors['#1'] = vector.mulnum(vectors['#2'],#3)
}%
}

\newcommand\vectorCross[3]{%
\directlua{%
vectors['#1'] = vector.cross(vectors['#2'],vectors['#3'])
}%
}

\newcommand\vectorSumNorm[1]{%
\directlua{%
tex.sprint(tostring(vector.sumnorm(vectors['#1'])))
}%
}

\newcommand\vectorEuclidNorm[1]{%
\directlua{%
tex.sprint(tostring(vector.euclidnorm(vectors['#1'])))
}%
}

\newcommand\vectorSupNorm[1]{%
\directlua{%
tex.sprint(tostring(vector.supnorm(vectors['#1'])))
}%
}

\newcommand\vectorpNorm[2]{%
\directlua{%
tex.sprint(tostring(vector.pnorm(vectors['#1'],#2)))
}%
}

\newcommand\vectorCreateRandom[4]{%
\directlua{%
vectors['#1'] = vector.createrandom(#2,#3,#4)
}%
}

\newcommand\vectorCopy[2]{%
\directlua{%
vectors['#1'] = vector.copy(vectors['#2'])
}%
}

\newcommand\vectorGetCoordinate[2]{%
\directlua{%
tex.sprint(tostring(vector.getcoordinate(vectors['#1'],#2)))
}%
}

\newcommand\vectorSetCoordinate[3]{%
\directlua{%
tex.sprint(tostring(vector.setcoordinate(vectors['#1'],#2,#3)))
}%
}

\newcommand\vectorGetAngle[2]{%
\directlua{%
tex.sprint(tostring(vector.getangle(vectors['#1'],vectors['#2'])))
}%
}

\newcommand\vectorEql[2]{%
\directlua{%
tex.sprint(tostring(vector.chqeql(vectors['#1'],vectors['#2'])))
}%
}

% ========= KEY DEFINITIONS =========
\define@key{vecrr}{brckt}{\def\voprr@brckt{#1}}
\define@key{vecrr}{truncate}{\def\voprr@truncate{#1}}

% ========= KEY DEFAULTS =========
\setkeys{vecrr}{brckt=round,truncate=6}%
\newcommand{\vectorGramSchmidt}[2][]{%
\begingroup%
\setkeys{vecrr}{#1}
\directlua{%
local tbl = #2
local outTbl={}
local sum = 0
for i=1,table.getn(tbl) do
outTbl[i] = vectors[tbl[i]]
end

tex.sprint(vector.gs(outTbl,"\voprr@brckt",\voprr@truncate))}
%
\endgroup%
}

\newcommand{\vectorGramSchmidtSteps}[2][]{%
\begingroup%
\setkeys{vecrr}{#1}
\directlua{%
local tbl = #2
local outTbl={}
local sum = 0
for i=1,table.getn(tbl) do
outTbl[i] = vectors[tbl[i]]
end
tex.sprint(vector.gsX(outTbl,"\voprr@brckt",\voprr@truncate))}
%
\endgroup%
}

\endinput