summaryrefslogtreecommitdiff
path: root/tex/context/base/m-matrix.mkiv
diff options
context:
space:
mode:
Diffstat (limited to 'tex/context/base/m-matrix.mkiv')
-rw-r--r--tex/context/base/m-matrix.mkiv495
1 files changed, 495 insertions, 0 deletions
diff --git a/tex/context/base/m-matrix.mkiv b/tex/context/base/m-matrix.mkiv
new file mode 100644
index 000000000..ccb376e39
--- /dev/null
+++ b/tex/context/base/m-matrix.mkiv
@@ -0,0 +1,495 @@
+%D \module
+%D [ file=m-matrix,
+%D version=2014.11.04, % already a year older
+%D title=\CONTEXT\ Extra Modules,
+%D subtitle=Matrices,
+%D author={Jeong Dalyoung \& Hans Hagen},
+%D date=\currentdate,
+%D copyright={PRAGMA ADE \& \CONTEXT\ Development Team}]
+%C
+%C This module is part of the \CONTEXT\ macro||package and is
+%C therefore copyrighted by \PRAGMA. See mreadme.pdf for
+%C details.
+
+%D This code is based on a post by Dalyoung on the context list. After that
+%D we turned it into a module and improved the code a bit. Feel free to ask
+%D us for more. Once we're satisfied, a more general helper l-matrix could
+%D be made. Dalyoung does the clever bits, and Hans only cleanes up and
+%D optimizes a bit.
+
+% \registerctxluafile{l-matrix}{1.001} % not yet
+
+\startmodule[matrix]
+
+\startluacode
+
+local settings_to_hash = utilities.parsers.settings_to_hash
+local formatters = string.formatters
+local copy = table.copy
+local insert = table.insert
+local remove = table.remove
+
+local matrix = { }
+moduledata.matrix = matrix
+
+local f_matrix_slot = formatters["%s_{%s%s}"]
+
+function matrix.symbolic(sym, x, y, nx ,ny) -- symMatrix("a", "m", "n")
+ local nx = nx or 2
+ local ny = ny or nx
+ local function filled(i,y)
+ local mrow = { }
+ for j=1,nx do
+ mrow[#mrow+1] = f_matrix_slot(sym,i,j)
+ end
+ mrow[#mrow+1] = "\\cdots"
+ mrow[#mrow+1] = f_matrix_slot(sym,i,y)
+ return mrow
+ end
+ local function dummy()
+ local mrow = { }
+ for j=1,nx do
+ mrow[#mrow+1] = "\\vdots"
+ end
+ mrow[#mrow+1] = "\\ddots"
+ mrow[#mrow+1] = "\\vdots"
+ return mrow
+ end
+ --
+ local mm = { }
+ for i=1,ny do
+ mm[i] = filled(i,y)
+ end
+ mm[#mm+1] = dummy()
+ mm[#mm+1] = filled(x,y)
+ return mm
+end
+
+-- todo: define a matrix at the tex end so that we have more control
+
+local fences_p = {
+ left = "\\left(\\,",
+ right = "\\,\\right)",
+}
+
+local fences_b = {
+ left = "\\left[\\,",
+ right = "\\,\\right]",
+}
+
+function matrix.typeset(m,options)
+ local options = settings_to_hash(options or "")
+ context.startmatrix(options.determinant and fences_b or fences_p)
+ for i=1, #m do
+ local mi = m[i]
+ for j=1,#mi do
+ context.NC(mi[j])
+ end
+ context.NR()
+ end
+ context.stopmatrix()
+end
+
+-- interchange two rows (i-th, j-th)
+
+function matrix.swap(t,i,j)
+ t[i], t[j] = t[j], t[i]
+end
+
+-- replace i-th row with factor * (i-th row)
+
+function matrix.multiply(m,i,factor)
+ local mi = m[i]
+ for k=1,#mi do
+ mi[k] = factor * mi[k]
+ end
+ return m
+end
+
+-- scalar product "factor * m"
+
+function matrix.scalar(m, factor)
+ for i=1,#m do
+ local mi = m[i]
+ for j=1,#mi do
+ mi[j] = factor * mi[j]
+ end
+ end
+ return m
+end
+
+-- replace i-th row with i-th row + factor * (j-th row)
+
+function matrix.sumrow(m,i,j,factor)
+ local mi = m[i]
+ local mj = m[j]
+ for k=1,#mi do
+ mi[k] = mi[k] + factor * mj[k]
+ end
+end
+
+-- transpose of a matrix
+
+function matrix.transpose(m)
+ local t = { }
+ for j=1,#m[1] do
+ local r = { }
+ for i=1,#m do
+ r[i] = m[i][j]
+ end
+ t[j] = r
+ end
+ return t
+end
+
+-- inner product of two vectors
+
+function matrix.inner(u,v)
+ local nu = #u
+ if nu == 0 then
+ return 0
+ end
+ local nv = #v
+ if nv ~= nu then
+ return 0
+ end
+ local result = 0
+ for i=1,nu do
+ result = result + u[i] * v[i]
+ end
+ return result
+end
+
+-- product of two matrices
+
+function matrix.product(m1,m2)
+ local product = { }
+ if #m1[1] == #m2 then
+ for i=1,#m1 do
+ local m1i = m1[i]
+ local mrow = { }
+ for j=1,#m2[1] do
+ local temp = 0
+ for k=1,#m1[1] do
+ temp = temp + m1i[k] * m2[k][j]
+ end
+ mrow[j] = temp
+ end
+ product[i] = mrow
+ end
+ end
+ return product
+end
+
+local function uppertri(m,sign)
+ local temp = copy(m)
+ for i=1,#temp-1 do
+ local pivot = temp[i][i]
+ if pivot == 0 then
+ local pRow = i +1
+ while temp[pRow][i] == 0 do
+ pRow = pRow + 1
+ if pRow > #temp then -- if there is no nonzero number
+ return temp
+ end
+ end
+ temp[i], temp[pRow] = temp[pRow], temp[i]
+ if sign then
+ sign = -sign
+ end
+ end
+ local mi = temp[i]
+ for k=i+1, #temp do
+ local factor = -temp[k][i]/mi[i]
+ local mk = temp[k]
+ for l=i,#mk do
+ mk[l] = mk[l] + factor * mi[l]
+ end
+ end
+ end
+ if sign then
+ return temp, sign
+ else
+ return temp
+ end
+end
+
+matrix.uppertri = uppertri
+
+function matrix.determinant(m)
+ if #m == #m[1] then
+ local d = 1
+ local t, s = uppertri(m,1)
+ for i=1,#t do
+ d = d * t[i][i]
+ end
+ return s*d
+ else
+ return 0
+ end
+end
+
+local function rowechelon(m,r)
+ local temp = copy(m)
+ local pRow = 1
+ local pCol = 1
+ while pRow <= #temp do
+ local pivot = temp[pRow][pCol]
+ if pivot == 0 then
+ local i = pRow
+ local n = #temp
+ while temp[i][pCol] == 0 do
+ i = i + 1
+ if i > n then
+ -- no nonzero number in a column
+ pCol = pCol + 1
+ if pCol > #temp[pRow] then
+ -- there is no nonzero number in a row
+ return temp
+ end
+ i = pRow
+ end
+ end
+ temp[pRow], temp[i] = temp[i], temp[pRow]
+ end
+ local row = temp[pRow]
+ pivot = row[pCol]
+ for l=pCol,#row do
+ row[l] = row[l]/pivot
+ end
+
+ if r == 1 then
+ -- make the "reduced row echelon form"
+ local row = temp[pRow]
+ for k=1,pRow-1 do
+ local current = temp[k]
+ local factor = -current[pCol]
+ local mk = current
+ for l=pCol,#mk do
+ mk[l] = mk[l] + factor * row[l]
+ end
+ end
+ end
+ -- just make the row echelon form
+ local row = temp[pRow]
+ for k=pRow+1, #temp do
+ local current = temp[k]
+ local factor = -current[pCol]
+ local mk = current
+ for l=pCol,#mk do
+ mk[l] = mk[l] + factor * row[l]
+ end
+ end
+ pRow = pRow + 1
+ pCol = pCol + 1
+
+ if pRow > #temp or pCol > #temp[1] then
+ pRow = #temp + 1
+ end
+ end
+ return temp
+end
+
+matrix.rowechelon = rowechelon
+matrix.rowEchelon = rowechelon
+
+-- solve the linear equation m X = c
+
+local function solve(m,c)
+ local n = #m
+ if n ~= #c then
+ return copy(m)
+ end
+ local newm = copy(m)
+ local temp = copy(c)
+ for i=1,n do
+ insert(newm[i],temp[i])
+ end
+ return rowechelon(newm,1)
+end
+
+matrix.solve = solve
+
+-- find the inverse matrix of m
+
+local function inverse(m)
+ local n = #m
+ local temp = copy(m)
+ if n ~= #m[1] then
+ return temp
+ end
+ for i=1,n do
+ for j=1,n do
+ insert(temp[i],j == i and 1 or 0)
+ end
+ end
+ temp = rowechelon(temp,1)
+ for i=1,n do
+ for j=1,n do
+ remove(temp[i], 1)
+ end
+ end
+ return temp
+end
+
+matrix.inverse = inverse
+
+\stopluacode
+
+\stopmodule
+
+\unexpanded\def\ctxmodulematrix#1{\ctxlua{moduledata.matrix.#1}}
+
+\continueifinputfile{m-matrix.mkiv}
+
+\starttext
+
+\startluacode
+document.DemoMatrixA = {
+ { 0, 2, 4, -4, 1 },
+ { 0, 0, 2, 3, 4 },
+ { 2, 2, -6, 2, 4 },
+ { 2, 0, -6, 9, 7 },
+ { 2, 3, 4, 5, 6 },
+ { 6, 6, -6, 6, 6 },
+}
+
+document.DemoMatrixB = {
+ { 0, 2, 4, -4, 1 },
+ { 0, 0, 2, 3, 4 },
+ { 2, 2, -6, 2, 4 },
+ { 2, 0, -6, 9, 7 },
+ { 2, 2, -6, 2, 4 },
+ { 2, 2, -6, 2, 4 },
+}
+\stopluacode
+
+\startsubject[title={A symbolic matrix}]
+
+\ctxmodulematrix{typeset(moduledata.matrix.symbolic("a", "m", "n"))}
+\ctxmodulematrix{typeset(moduledata.matrix.symbolic("a", "m", "n", 4, 8))}
+
+\stopsubject
+
+\startsubject[title={Swap two rows (2 and 4)}]
+
+\startluacode
+moduledata.matrix.typeset(document.DemoMatrixA)
+context.blank()
+moduledata.matrix.swap(document.DemoMatrixA, 2, 4)
+context.blank()
+moduledata.matrix.typeset(document.DemoMatrixA)
+\stopluacode
+
+\stopsubject
+
+\startsubject[title={Multiply $3 \times r_2$}]
+
+\startluacode
+moduledata.matrix.typeset(document.DemoMatrixA)
+context.blank()
+moduledata.matrix.typeset(moduledata.matrix.multiply(document.DemoMatrixA, 2, 3))
+\stopluacode
+
+\stopsubject
+
+\startsubject[title={Row 2 + $3 \times r_4$}]
+
+\startluacode
+moduledata.matrix.typeset(document.DemoMatrixA)
+context.blank()
+moduledata.matrix.sumrow(document.DemoMatrixA, 2, 3, 4)
+context.blank()
+moduledata.matrix.typeset(document.DemoMatrixA)
+\stopluacode
+
+\stopsubject
+
+\startsubject[title={Transpose a matrix}]
+
+\startluacode
+moduledata.matrix.typeset(document.DemoMatrixA)
+context.blank()
+moduledata.matrix.typeset(moduledata.matrix.transpose(document.DemoMatrixA))
+\stopluacode
+
+\stopsubject
+
+\startsubject[title={The inner product of two vectors}]
+
+\startluacode
+context(moduledata.matrix.inner({ 1, 2, 3 }, { 3, 1, 2 }))
+context.blank()
+context(moduledata.matrix.inner({ 1, 2, 3 }, { 3, 1, 2, 4 }))
+\stopluacode
+
+\startsubject[title={The product of two matrices}]
+
+\startluacode
+moduledata.matrix.typeset(document.DemoMatrixA)
+context.blank()
+moduledata.matrix.typeset(moduledata.matrix.product(document.DemoMatrixA,document.DemoMatrixA))
+\stopluacode
+
+\stopsubject
+
+\startsubject[title={An Upper Triangular Matrix}]
+
+\ctxmodulematrix{typeset(moduledata.matrix.uppertri(document.DemoMatrixB))}
+
+\startsubject[title={A determinant}]
+
+\startluacode
+local m = {
+ { 1, 2, 4 },
+ { 0, 0, 2 },
+ { 2, 2, -6 },
+}
+context(moduledata.matrix.determinant(m))
+\stopluacode
+
+\stopsubject
+
+\startsubject[title={Row echelon form}]
+
+\startluacode
+local m = {
+ { 1, 3, -2, 0, 2, 0, 0 },
+ { 2, 6, -5, -2, 4, -3, -1 },
+ { 0, 0, 5, 10, 0, 15, 5 },
+ { 2, 6, 0, 8, 4, 18, 6 },
+}
+
+moduledata.matrix.typeset(m)
+moduledata.matrix.typeset(moduledata.matrix.rowechelon(m,1))
+\stopluacode
+
+\stopsubject
+
+\startsubject[title={Solving linear equation}]
+
+\startluacode
+local m = {
+ { 1, 3, -2, 0 },
+ { 2, 0, 1, 2 },
+ { 6, -5, -2, 4 },
+ { -3, -1, 5, 10 },
+}
+
+local c = { 5, 2, 6, 8 }
+
+moduledata.matrix.typeset(moduledata.matrix.solve(m,c))
+\stopluacode
+
+\stopsubject
+
+\startsubject[title={Inverse matrix}]
+
+\startcombination[2*1]
+ {\ctxlua{moduledata.matrix.typeset { { 1, 1, 1 }, { 0, 2, 3 }, { 3, 2, 1 } }}} {}
+ {\ctxlua{moduledata.matrix.typeset(moduledata.matrix.inverse { { 1, 1, 1 }, { 0, 2, 3 }, { 3, 2, 1 } })}} {}
+\stopcombination
+
+\stopsubject
+
+\stoptext
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-18 23:54:14 +0000