diff options
author | Norbert Preining <norbert@preining.info> | 2023-02-10 03:01:41 +0000 |
---|---|---|
committer | Norbert Preining <norbert@preining.info> | 2023-02-10 03:01:41 +0000 |
commit | e10a52ebe817fe7ff0235c066d83f51f39d3e023 (patch) | |
tree | 9148206ad017e5947e391d1eb4f3858feb493648 /macros/luatex/latex/lualinalg/lualinalg.sty | |
parent | ffdef3cddca5a0b9133dd928291b000452777ad8 (diff) |
CTAN sync 202302100301
Diffstat (limited to 'macros/luatex/latex/lualinalg/lualinalg.sty')
-rw-r--r-- | macros/luatex/latex/lualinalg/lualinalg.sty | 1685 |
1 files changed, 1685 insertions, 0 deletions
diff --git a/macros/luatex/latex/lualinalg/lualinalg.sty b/macros/luatex/latex/lualinalg/lualinalg.sty new file mode 100644 index 0000000000..53fe7fd688 --- /dev/null +++ b/macros/luatex/latex/lualinalg/lualinalg.sty @@ -0,0 +1,1685 @@ + +% The lualinalg package +% Authors: Chetan Shirore and Ajit Kumar +% version 1.0 +% 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.0] +\RequirePackage{xkeyval} +\RequirePackage{amsmath} +\RequirePackage{luamaths} +\RequirePackage{luacode} +\begin{luacode*} +local complex = require "luacomplex" + +-- matrices part +matrices = {} + +matrix = {} --module + +local matrix_meta = {} + +function matrix.new(matrix, rows, columns) + 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 +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.0) +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 + +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 + +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 = {} + +vector = {} --module + +local vector_meta = {} + +function vector.new(vector, rows, columns, n) + if columns ~= "e" 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 +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.0) +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 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.0) + 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 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 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.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 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 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 + +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.__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\matrixRREFERR[1]{% + \directlua{% + tex.sprint(tostring(matrix.rref0E(matrices['#1']))) + }% +} + +\newcommand\matrixRREFE[1]{% + \directlua{% + tex.sprint(tostring(matrix.rref0E(matrices['#1']))) + }% +} + +% ========= 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\complexRound[2]{% + \directlua{% + tex.sprint(tostring(complex.round(complex('#1'),#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 |