diff --git a/README.md b/README.md index 164e84a..a5a2146 100644 --- a/README.md +++ b/README.md @@ -8,6 +8,12 @@ If you like my contributions you may consider reading http://entuland.com/en/sup WIP MOD forum thread: https://forum.minetest.net/viewtopic.php?f=9&t=20321 +# Dependencies + +A thin wrapper around a very useful library to deal with Matrices: + +[matrix] https://github.com/entuland/lua-matrix + # Why yet another screwdriver? The default screwdriver included in minetest_game, as well as any other screwdriver mod I have found, operate differently depending on the node's direction and rotation. This means that any given click on a node may produce different results which you cannot predict at a glance. diff --git a/init.lua b/init.lua index ecbeeef..63b0ded 100644 --- a/init.lua +++ b/init.lua @@ -2,8 +2,6 @@ rhotator = {} local mod_path = minetest.get_modpath(minetest.get_current_modname()) -local matrix = dofile(mod_path .. "/lib/matrix.lua") - -- constants local POS = {} diff --git a/lib/matrix.lua b/lib/matrix.lua deleted file mode 100644 index 9e1ad67..0000000 --- a/lib/matrix.lua +++ /dev/null @@ -1,1251 +0,0 @@ ---[[ - -LUA MODULE - - matrix v$(_VERSION) - matrix functions implemented with Lua tables - -SYNOPSIS - - local matrix = require 'matrix' - m1 = matrix{{8,4,1},{6,8,3}} - m2 = matrix{{-8,1,3},{5,2,1}} - assert(m1 + m2 == matrix{{0,5,4},{11,10,4}}) - -DESCRIPTION - - With simple matrices this script is quite useful, though for more - exact calculations, one would probably use a program like Matlab instead. - Matrices of size 100x100 can still be handled very well. - The error for the determinant and the inverted matrix is around 10^-9 - with a 100x100 matrix and an element range from -100 to 100. - - Characteristics: - - - functions called via matrix. should be able to handle - any table matrix of structure t[i][j] = value - - can handle a type of complex matrix - - can handle symbolic matrices. (Symbolic matrices cannot be - used with complex matrices.) - - arithmetic functions do not change the matrix itself - but build and return a new matrix - - functions are intended to be light on checks - since one gets a Lua error on incorrect use anyways - - uses mainly Gauss-Jordan elimination - - for Lua tables optimised determinant calculation (fast) - but not invoking any checks for special types of matrices - - vectors can be set up via vec1 = matrix{{ 1,2,3 }}^'T' or matrix{1,2,3} - - vectors can be multiplied to a scalar via num = vec1^'T' * vec2 - where num will be a matrix with the result in mtx[1][1], - or use num = vec1:scalar( vec2 ), where num is a number - -API - - matrix function list: - - matrix.add - matrix.columns - matrix.concath - matrix.concatv - matrix.copy - matrix.cross - matrix.det - matrix.div - matrix.divnum - matrix.dogauss - matrix.elementstostring - matrix.getelement - matrix.gsub - matrix.invert - matrix.ipairs - matrix.latex - matrix.len - matrix.mul - matrix.mulnum - matrix:new - matrix.normf - matrix.normmax - matrix.pow - matrix.print - matrix.random - matrix.replace - matrix.root - matrix.rotl - matrix.rotr - matrix.round - matrix.rows - matrix.scalar - matrix.setelement - matrix.size - matrix.solve - matrix.sqrt - matrix.sub - matrix.subm - matrix.tostring - matrix.transpose - matrix.type - - See code and test_matrix.lua. - -DEPENDENCIES - - None (other than Lua 5.1 or 5.2). May be used with complex.lua. - -HOME PAGE - - http://luamatrix.luaforge.net - http://lua-users.org/wiki/LuaMatrix - -DOWNLOAD/INSTALL - - ./util.mk - cd tmp/* - luarocks make - -LICENSE - - Licensed under the same terms as Lua itself. - - Developers: - Michael Lutz (chillcode) - original author - David Manura http://lua-users.org/wiki/DavidManura ---]] - ---//////////// ---// matrix // ---//////////// - -local matrix = {_TYPE='module', _NAME='matrix', _VERSION='0.2.11.20120416'} - --- access to the metatable we set at the end of the file -local matrix_meta = {} - ---///////////////////////////// ---// Get 'new' matrix object // ---///////////////////////////// - ---// matrix:new ( rows [, columns [, value]] ) --- if rows is a table then sets rows as matrix --- if rows is a table of structure {1,2,3} then it sets it as a vector matrix --- if rows and columns are given and are numbers, returns a matrix with size rowsxcolumns --- if num is given then returns a matrix with given size and all values set to num --- if rows is given as number and columns is "I", will return an identity matrix of size rowsxrows -function matrix:new( rows, columns, value ) - -- check for given matrix - if type( rows ) == "table" then - -- check for vector - if type(rows[1]) ~= "table" then -- expect a vector - return setmetatable( {{rows[1]},{rows[2]},{rows[3]}},matrix_meta ) - end - return setmetatable( rows,matrix_meta ) - end - -- get matrix table - local mtx = {} - local value = value or 0 - -- build identity matrix of given rows - 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 - -- build new matrix - else - for i = 1,rows do - mtx[i] = {} - for j = 1,columns do - mtx[i][j] = value - end - end - end - -- return matrix with shared metatable - return setmetatable( mtx,matrix_meta ) -end - ---// matrix ( rows [, comlumns [, value]] ) --- set __call behaviour of matrix --- for matrix( ... ) as matrix.new( ... ) -setmetatable( matrix, { __call = function( ... ) return matrix.new( ... ) end } ) - - --- functions are designed to be light on checks --- so we get Lua errors instead on wrong input --- matrix. should handle any table of structure t[i][j] = value --- we always return a matrix with scripts metatable --- cause its faster than setmetatable( mtx, getmetatable( input matrix ) ) - ---/////////////////////////////// ---// matrix 'matrix' functions // ---/////////////////////////////// - ---// for real, complex and symbolic matrices //-- - --- note: real and complex matrices may be added, subtracted, etc. --- real and symbolic matrices may also be added, subtracted, etc. --- but one should avoid using symbolic matrices with complex ones --- since it is not clear which metatable then is used - ---// matrix.add ( m1, m2 ) --- Add two matrices; m2 may be of bigger size than m1 -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 - ---// matrix.sub ( m1 ,m2 ) --- Subtract two matrices; m2 may be of bigger size than m1 -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 - ---// matrix.mul ( m1, m2 ) --- Multiply two matrices; m1 columns must be equal to m2 rows --- e.g. #m1[1] == #m2 -function matrix.mul( m1, m2 ) - -- multiply rows with columns - 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 - ---// matrix.div ( m1, m2 ) --- Divide two matrices; m1 columns must be equal to m2 rows --- m2 must be square, to be inverted, --- if that fails returns the rank of m2 as second argument --- e.g. #m1[1] == #m2; #m2 == #m2[1] -function matrix.div( m1, m2 ) - local rank; m2,rank = matrix.invert( m2 ) - if not m2 then return m2, rank end -- singular - return matrix.mul( m1, m2 ) -end - ---// matrix.mulnum ( m1, num ) --- Multiply matrix with a number --- num may be of type 'number' or 'complex number' --- strings get converted to complex number, if that fails then to symbol -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 - ---// matrix.divnum ( m1, num ) --- Divide matrix by a number --- num may be of type 'number' or 'complex number' --- strings get converted to complex number, if that fails then to symbol -function matrix.divnum( m1, num ) - local mtx = {} - -- divide elements by number - for i = 1,#m1 do - local mtxi = {} - mtx[i] = mtxi - for j = 1,#m1[1] do - mtxi[j] = m1[i][j] / num - end - end - return setmetatable( mtx, matrix_meta ) -end - - ---// for real and complex matrices only //-- - ---// matrix.pow ( m1, num ) --- Power of matrix; mtx^(num) --- num is an integer and may be negative --- m1 has to be square --- if num is negative and inverting m1 fails --- returns the rank of matrix m1 as second argument -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 - -local function number_norm2(x) - return x * x -end - ---// matrix.det ( m1 ) --- Calculate the determinant of a matrix --- m1 needs to be square --- Can calc the det for symbolic matrices up to 3x3 too --- The function to calculate matrices bigger 3x3 --- is quite fast and for matrices of medium size ~(100x100) --- and average values quite accurate --- here we try to get the nearest element to |1|, (smallest pivot element) --- os that usually we have |mtx[i][j]/subdet| > 1 or mtx[i][j]; --- with complex matrices we use the complex.abs function to check if it is bigger or smaller -function matrix.det( m1 ) - - -- check if matrix is quadratic - 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 - - --// no symbolic matrix supported below here - 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 - - --// matrix is bigger than 3x3 - -- get determinant - -- using Gauss elimination and Laplace - -- start eliminating from below better for removals - -- get copy of matrix, set initial determinant - local mtx = matrix.copy( m1 ) - local det = 1 - -- get det up to the last element - for j = 1,#mtx[1] do - -- get smallest element so that |factor| > 1 - -- and set it as last element - local rows = #mtx - local subdet,xrow - for i = 1,rows do - -- get element - local e = mtx[i][j] - -- if no subdet has been found - if not subdet then - -- check if element it is not zero - if e ~= zero then - -- use element as new subdet - subdet,xrow = e,i - end - -- check for elements nearest to 1 or -1 - elseif e ~= zero and math.abs(norm2(e)-1) < math.abs(norm2(subdet)-1) then - subdet,xrow = e,i - end - end - -- only cary on if subdet is found - if subdet then - -- check if xrow is the last row, - -- else switch lines and multiply det by -1 - if xrow ~= rows then - mtx[rows],mtx[xrow] = mtx[xrow],mtx[rows] - det = -det - end - -- traverse all fields setting element to zero - -- we don't set to zero cause we don't use that column anymore then anyways - for i = 1,rows-1 do - -- factor is the dividor of the first element - -- if element is not already zero - if mtx[i][j] ~= zero then - local factor = mtx[i][j]/subdet - -- update all remaining fields of the matrix, with value from xrow - for n = j+1,#mtx[1] do - mtx[i][n] = mtx[i][n] - factor * mtx[rows][n] - end - end - end - -- update determinant and remove row - if math.fmod( rows,2 ) == 0 then - det = -det - end - det = det * subdet - table.remove( mtx ) - else - -- break here table det is 0 - return det * 0 - end - end - -- det ready to return - return det -end - ---// matrix.dogauss ( mtx ) --- Gauss elimination, Gauss-Jordan Method --- this function changes the matrix itself --- returns on success: true, --- returns on failure: false,'rank of matrix' - --- locals --- checking here for the element nearest but not equal to zero (smallest pivot element). --- This way the `factor` in `dogauss` will be >= 1, which --- can give better results. -local pivotOk = function( mtx,i,j,norm2 ) - -- find min value - local iMin - local normMin = math.huge - for _i = i,#mtx do - local e = mtx[_i][j] - local norm = math.abs(norm2(e)) - if norm > 0 and norm < normMin then - iMin = _i - normMin = norm - end - end - if iMin then - -- switch lines if not in position. - if iMin ~= i then - mtx[i],mtx[iMin] = mtx[iMin],mtx[i] - end - return true - end - return false -end - -local function copy(x) - return type(x) == "table" and x.copy(x) or x -end - --- note: in --// ... //-- we have a way that does no divison, --- however with big number and matrices we get problems since we do no reducing -function matrix.dogauss( mtx ) - local e = mtx[1][1] - local zero = type(e) == "table" and e.zero or 0 - local one = type(e) == "table" and e.one or 1 - local norm2 = type(e) == "table" and e.norm2 or number_norm2 - - local rows,columns = #mtx,#mtx[1] - -- stairs left -> right - for j = 1,rows do - -- check if element can be setted to one - if pivotOk( mtx,j,j,norm2 ) then - -- start parsing rows - for i = j+1,rows do - -- check if element is not already zero - if mtx[i][j] ~= zero then - -- we may add x*otherline row, to set element to zero - -- tozero - x*mtx[j][j] = 0; x = tozero/mtx[j][j] - local factor = mtx[i][j]/mtx[j][j] - --// this should not be used although it does no division, - -- yet with big matrices (since we do no reducing and other things) - -- we get too big numbers - --local factor1,factor2 = mtx[i][j],mtx[j][j] //-- - mtx[i][j] = copy(zero) - for _j = j+1,columns do - --// mtx[i][_j] = mtx[i][_j] * factor2 - factor1 * mtx[j][_j] //-- - mtx[i][_j] = mtx[i][_j] - factor * mtx[j][_j] - end - end - end - else - -- return false and the rank of the matrix - return false,j-1 - end - end - -- stairs right <- left - for j = rows,1,-1 do - -- set element to one - -- do division here - local div = mtx[j][j] - for _j = j+1,columns do - mtx[j][_j] = mtx[j][_j] / div - end - -- start parsing rows - for i = j-1,1,-1 do - -- check if element is not already zero - if mtx[i][j] ~= zero then - local factor = mtx[i][j] - for _j = j+1,columns do - mtx[i][_j] = mtx[i][_j] - factor * mtx[j][_j] - end - mtx[i][j] = copy(zero) - end - end - mtx[j][j] = copy(one) - end - return true -end - ---// matrix.invert ( m1 ) --- Get the inverted matrix or m1 --- matrix must be square and not singular --- on success: returns inverted matrix --- on failure: returns nil,'rank of matrix' -function matrix.invert( m1 ) - assert(#m1 == #m1[1], "matrix not square") - local mtx = matrix.copy( m1 ) - local ident = setmetatable( {},matrix_meta ) - local e = m1[1][1] - local zero = type(e) == "table" and e.zero or 0 - local one = type(e) == "table" and e.one or 1 - for i = 1,#m1 do - local identi = {} - ident[i] = identi - for j = 1,#m1 do - identi[j] = copy((i == j) and one or zero) - end - end - mtx = matrix.concath( mtx,ident ) - local done,rank = matrix.dogauss( mtx ) - if done then - return matrix.subm( mtx, 1,(#mtx[1]/2)+1,#mtx,#mtx[1] ) - else - return nil,rank - end -end - ---// matrix.sqrt ( m1 [,iters] ) --- calculate the square root of a matrix using "Denman Beavers square root iteration" --- condition: matrix rows == matrix columns; must have a invers matrix and a square root --- if called without additional arguments, the function finds the first nearest square root to --- input matrix, there are others but the error between them is very small --- if called with agument iters, the function will return the matrix by number of iterations --- the script returns: --- as first argument, matrix^.5 --- as second argument, matrix^-.5 --- as third argument, the average error between (matrix^.5)^2-inputmatrix --- you have to determin for yourself if the result is sufficent enough for you --- local average error -local function get_abs_avg( m1, m2 ) - local dist = 0 - local e = m1[1][1] - local abs = type(e) == "table" and e.abs or math.abs - for i=1,#m1 do - for j=1,#m1[1] do - dist = dist + abs(m1[i][j]-m2[i][j]) - end - end - -- norm by numbers of entries - return dist/(#m1*2) -end --- square root function -function matrix.sqrt( m1, iters ) - assert(#m1 == #m1[1], "matrix not square") - local iters = iters or math.huge - local y = matrix.copy( m1 ) - local z = matrix(#y, 'I') - local dist = math.huge - -- iterate, and get the average error - for n=1,iters do - local lasty,lastz = y,z - -- calc square root - -- y, z = (1/2)*(y + z^-1), (1/2)*(z + y^-1) - y, z = matrix.divnum((matrix.add(y,matrix.invert(z))),2), - matrix.divnum((matrix.add(z,matrix.invert(y))),2) - local dist1 = get_abs_avg(y,lasty) - if iters == math.huge then - if dist1 >= dist then - return lasty,lastz,get_abs_avg(matrix.mul(lasty,lasty),m1) - end - end - dist = dist1 - end - return y,z,get_abs_avg(matrix.mul(y,y),m1) -end - ---// matrix.root ( m1, root [,iters] ) --- calculate any root of a matrix --- source: http://www.dm.unipi.it/~cortona04/slides/bruno.pdf --- m1 and root have to be given;(m1 = matrix, root = number) --- conditions same as matrix.sqrt --- returns same values as matrix.sqrt -function matrix.root( m1, root, iters ) - assert(#m1 == #m1[1], "matrix not square") - local iters = iters or math.huge - local mx = matrix.copy( m1 ) - local my = matrix.mul(mx:invert(),mx:pow(root-1)) - local dist = math.huge - -- iterate, and get the average error - for n=1,iters do - local lastx,lasty = mx,my - -- calc root of matrix - --mx,my = ((p-1)*mx + my^-1)/p, - -- ((((p-1)*my + mx^-1)/p)*my^-1)^(p-2) * - -- ((p-1)*my + mx^-1)/p - mx,my = mx:mulnum(root-1):add(my:invert()):divnum(root), - my:mulnum(root-1):add(mx:invert()):divnum(root) - :mul(my:invert():pow(root-2)):mul(my:mulnum(root-1) - :add(mx:invert())):divnum(root) - local dist1 = get_abs_avg(mx,lastx) - if iters == math.huge then - if dist1 >= dist then - return lastx,lasty,get_abs_avg(matrix.pow(lastx,root),m1) - end - end - dist = dist1 - end - return mx,my,get_abs_avg(matrix.pow(mx,root),m1) -end - - ---// Norm functions //-- - ---// matrix.normf ( mtx ) --- calculates the Frobenius norm of the matrix. --- ||mtx||_F = sqrt(SUM_{i,j} |a_{i,j}|^2) --- http://en.wikipedia.org/wiki/Frobenius_norm#Frobenius_norm -function matrix.normf(mtx) - local mtype = matrix.type(mtx) - local result = 0 - for i = 1,#mtx do - for j = 1,#mtx[1] do - local e = mtx[i][j] - if mtype ~= "number" then e = e:abs() end - result = result + e^2 - end - end - local sqrt = (type(result) == "number") and math.sqrt or result.sqrt - return sqrt(result) -end - ---// matrix.normmax ( mtx ) --- calculates the max norm of the matrix. --- ||mtx||_{max} = max{|a_{i,j}|} --- Does not work with symbolic matrices --- http://en.wikipedia.org/wiki/Frobenius_norm#Max_norm -function matrix.normmax(mtx) - local abs = (matrix.type(mtx) == "number") and math.abs or mtx[1][1].abs - local result = 0 - for i = 1,#mtx do - for j = 1,#mtx[1] do - local e = abs(mtx[i][j]) - if e > result then result = e end - end - end - return result -end - - ---// only for number and complex type //-- --- Functions changing the matrix itself - ---// matrix.round ( mtx [, idp] ) --- perform round on elements -local numround = function( num,mult ) - return math.floor( num * mult + 0.5 ) / mult -end -local tround = function( t,mult ) - for i,v in ipairs(t) do - t[i] = math.floor( v * mult + 0.5 ) / mult - end - return t -end -function matrix.round( mtx, idp ) - local mult = 10^( idp or 0 ) - local fround = matrix.type( mtx ) == "number" and numround or tround - for i = 1,#mtx do - for j = 1,#mtx[1] do - mtx[i][j] = fround(mtx[i][j],mult) - end - end - return mtx -end - ---// matrix.random( mtx [,start] [, stop] [, idip] ) --- fillmatrix with random values -local numfill = function( _,start,stop,idp ) - return math.random( start,stop ) / idp -end -local tfill = function( t,start,stop,idp ) - for i in ipairs(t) do - t[i] = math.random( start,stop ) / idp - end - return t -end -function matrix.random( mtx,start,stop,idp ) - local start,stop,idp = start or -10,stop or 10,idp or 1 - local ffill = matrix.type( mtx ) == "number" and numfill or tfill - for i = 1,#mtx do - for j = 1,#mtx[1] do - mtx[i][j] = ffill( mtx[i][j], start, stop, idp ) - end - end - return mtx -end - - ---////////////////////////////// ---// Object Utility Functions // ---////////////////////////////// - ---// for all types and matrices //-- - ---// matrix.type ( mtx ) --- get type of matrix, normal/complex/symbol or tensor -function matrix.type( mtx ) - local e = mtx[1][1] - if type(e) == "table" then - if e.type then - return e:type() - end - return "tensor" - end - return "number" -end - --- local functions to copy matrix values -local num_copy = function( num ) - return num -end -local t_copy = function( t ) - local newt = setmetatable( {}, getmetatable( t ) ) - for i,v in ipairs( t ) do - newt[i] = v - end - return newt -end - ---// matrix.copy ( m1 ) --- Copy a matrix --- simple copy, one can write other functions oneself -function matrix.copy( m1 ) - local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy - local mtx = {} - for i = 1,#m1[1] do - mtx[i] = {} - for j = 1,#m1 do - mtx[i][j] = docopy( m1[i][j] ) - end - end - return setmetatable( mtx, matrix_meta ) -end - ---// matrix.transpose ( m1 ) --- Transpose a matrix --- switch rows and columns -function matrix.transpose( m1 ) - local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy - local mtx = {} - for i = 1,#m1[1] do - mtx[i] = {} - for j = 1,#m1 do - mtx[i][j] = docopy( m1[j][i] ) - end - end - return setmetatable( mtx, matrix_meta ) -end - ---// matrix.subm ( m1, i1, j1, i2, j2 ) --- Submatrix out of a matrix --- input: i1,j1,i2,j2 --- i1,j1 are the start element --- i2,j2 are the end element --- condition: i1,j1,i2,j2 are elements of the matrix -function matrix.subm( m1,i1,j1,i2,j2 ) - local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy - 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] = docopy( m1[i][j] ) - end - end - return setmetatable( mtx, matrix_meta ) -end - ---// matrix.concath( m1, m2 ) --- Concatenate two matrices, horizontal --- will return m1m2; rows have to be the same --- e.g.: #m1 == #m2 -function matrix.concath( m1,m2 ) - assert(#m1 == #m2, "matrix size mismatch") - local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy - local mtx = {} - local offset = #m1[1] - for i = 1,#m1 do - mtx[i] = {} - for j = 1,offset do - mtx[i][j] = docopy( m1[i][j] ) - end - for j = 1,#m2[1] do - mtx[i][j+offset] = docopy( m2[i][j] ) - end - end - return setmetatable( mtx, matrix_meta ) -end - ---// matrix.concatv ( m1, m2 ) --- Concatenate two matrices, vertical --- will return m1 --- m2 --- columns have to be the same; e.g.: #m1[1] == #m2[1] -function matrix.concatv( m1,m2 ) - assert(#m1[1] == #m2[1], "matrix size mismatch") - local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy - local mtx = {} - for i = 1,#m1 do - mtx[i] = {} - for j = 1,#m1[1] do - mtx[i][j] = docopy( 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] = docopy( m2[i][j] ) - end - end - return setmetatable( mtx, matrix_meta ) -end - ---// matrix.rotl ( m1 ) --- Rotate Left, 90 degrees -function matrix.rotl( m1 ) - local mtx = matrix:new( #m1[1],#m1 ) - local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy - for i = 1,#m1 do - for j = 1,#m1[1] do - mtx[#m1[1]-j+1][i] = docopy( m1[i][j] ) - end - end - return mtx -end - ---// matrix.rotr ( m1 ) --- Rotate Right, 90 degrees -function matrix.rotr( m1 ) - local mtx = matrix:new( #m1[1],#m1 ) - local docopy = matrix.type( m1 ) == "number" and num_copy or t_copy - for i = 1,#m1 do - for j = 1,#m1[1] do - mtx[j][#m1-i+1] = docopy( m1[i][j] ) - end - end - return mtx -end - -local function tensor_tostring( t,fstr ) - if not fstr then return "["..table.concat(t,",").."]" end - local tval = {} - for i,v in ipairs( t ) do - tval[i] = string.format( fstr,v ) - end - return "["..table.concat(tval,",").."]" -end -local function number_tostring( e,fstr ) - return fstr and string.format( fstr,e ) or e -end - ---// matrix.tostring ( mtx, formatstr ) --- tostring function -function matrix.tostring( mtx, formatstr ) - local ts = {} - local mtype = matrix.type( mtx ) - local e = mtx[1][1] - local tostring = mtype == "tensor" and tensor_tostring or - type(e) == "table" and e.tostring or number_tostring - for i = 1,#mtx do - local tstr = {} - for j = 1,#mtx[1] do - tstr[j] = tostring(mtx[i][j],formatstr) - end - ts[i] = table.concat(tstr, "\t") - end - return table.concat(ts, "\n") -end - ---// matrix.print ( mtx [, formatstr] ) --- print out the matrix, just calls tostring -function matrix.print( ... ) - print( matrix.tostring( ... ) ) -end - ---// matrix.latex ( mtx [, align] ) --- LaTeX output -function matrix.latex( mtx, align ) - -- align : option to align the elements - -- c = center; l = left; r = right - -- \usepackage{dcolumn}; D{.}{,}{-1}; aligns number by . replaces it with , - local align = align or "c" - local str = "$\\left( \\begin{array}{"..string.rep( align, #mtx[1] ).."}\n" - local getstr = matrix.type( mtx ) == "tensor" and tensor_tostring or number_tostring - for i = 1,#mtx do - str = str.."\t"..getstr(mtx[i][1]) - for j = 2,#mtx[1] do - str = str.." & "..getstr(mtx[i][j]) - end - -- close line - if i == #mtx then - str = str.."\n" - else - str = str.." \\\\\n" - end - end - return str.."\\end{array} \\right)$" -end - - ---// Functions not changing the matrix - ---// matrix.rows ( mtx ) --- return number of rows -function matrix.rows( mtx ) - return #mtx -end - ---// matrix.columns ( mtx ) --- return number of columns -function matrix.columns( mtx ) - return #mtx[1] -end - ---// matrix.size ( mtx ) --- get matrix size as string rows,columns -function matrix.size( mtx ) - if matrix.type( mtx ) == "tensor" then - return #mtx,#mtx[1],#mtx[1][1] - end - return #mtx,#mtx[1] -end - ---// matrix.getelement ( mtx, i, j ) --- return specific element ( row,column ) --- returns element on success and nil on failure -function matrix.getelement( mtx,i,j ) - if mtx[i] and mtx[i][j] then - return mtx[i][j] - end -end - ---// matrix.setelement( mtx, i, j, value ) --- set an element ( i, j, value ) --- returns 1 on success and nil on failure -function matrix.setelement( mtx,i,j,value ) - if matrix.getelement( mtx,i,j ) then - -- check if value type is number - mtx[i][j] = value - return 1 - end -end - ---// matrix.ipairs ( mtx ) --- iteration, same for complex -function matrix.ipairs( mtx ) - local i,j,rows,columns = 1,0,#mtx,#mtx[1] - local function iter() - j = j + 1 - if j > columns then -- return first element from next row - i,j = i + 1,1 - end - if i <= rows then - return i,j - end - end - return iter -end - ---/////////////////////////////// ---// matrix 'vector' functions // ---/////////////////////////////// - --- a vector is defined as a 3x1 matrix --- get a vector; vec = matrix{{ 1,2,3 }}^'T' - ---// matrix.scalar ( m1, m2 ) --- returns the Scalar Product of two 3x1 matrices (vectors) -function matrix.scalar( m1, m2 ) - return m1[1][1]*m2[1][1] + m1[2][1]*m2[2][1] + m1[3][1]*m2[3][1] -end - ---// matrix.cross ( m1, m2 ) --- returns the Cross Product of two 3x1 matrices (vectors) -function matrix.cross( m1, m2 ) - local mtx = {} - mtx[1] = { m1[2][1]*m2[3][1] - m1[3][1]*m2[2][1] } - mtx[2] = { m1[3][1]*m2[1][1] - m1[1][1]*m2[3][1] } - mtx[3] = { m1[1][1]*m2[2][1] - m1[2][1]*m2[1][1] } - return setmetatable( mtx, matrix_meta ) -end - ---// matrix.len ( m1 ) --- returns the Length of a 3x1 matrix (vector) -function matrix.len( m1 ) - return math.sqrt( m1[1][1]^2 + m1[2][1]^2 + m1[3][1]^2 ) -end - - ---// matrix.replace (mtx, func, ...) --- for each element e in the matrix mtx, replace it with func(mtx, ...). -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 - ---// matrix.remcomplex ( mtx ) --- set the matrix elements to strings --- IMPROVE: tostring v.s. tostringelements confusing -function matrix.elementstostrings( mtx ) - local e = mtx[1][1] - local tostring = type(e) == "table" and e.tostring or tostring - return matrix.replace(mtx, tostring) -end - ---// matrix.solve ( m1 ) --- solve; tries to solve a symbolic matrix to a number -function matrix.solve( m1 ) - assert( matrix.type( m1 ) == "symbol", "matrix not of type 'symbol'" ) - local mtx = {} - for i = 1,#m1 do - mtx[i] = {} - for j = 1,#m1[1] do - mtx[i][j] = tonumber( loadstring( "return "..m1[i][j][1] )() ) - end - end - return setmetatable( mtx, matrix_meta ) -end - ---////////////////////////-- ---// METATABLE HANDLING //-- ---////////////////////////-- - ---// MetaTable --- as we declaired on top of the page --- local/shared metatable --- matrix_meta - --- note '...' is always faster than 'arg1,arg2,...' if it can be used - --- Set add "+" behaviour -matrix_meta.__add = function( ... ) - return matrix.add( ... ) -end - --- Set subtract "-" behaviour -matrix_meta.__sub = function( ... ) - return matrix.sub( ... ) -end - --- Set multiply "*" behaviour -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 - --- Set division "/" behaviour -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 - --- Set unary minus "-" behavior -matrix_meta.__unm = function( mtx ) - return matrix.mulnum( mtx,-1 ) -end - --- Set power "^" behaviour --- if opt is any integer number will do mtx^opt --- (returning nil if answer doesn't exist) --- if opt is 'T' then it will return the transpose matrix --- only for complex: --- if opt is '*' then it returns the complex conjugate matrix - local option = { - -- only for complex - ["*"] = function( m1 ) return matrix.conjugate( m1 ) end, - -- for both - ["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 - --- Set equal "==" behaviour -matrix_meta.__eq = function( m1, m2 ) - -- check same type - if matrix.type( m1 ) ~= matrix.type( m2 ) then - return false - end - -- check same size - if #m1 ~= #m2 or #m1[1] ~= #m2[1] then - return false - end - -- check elements equal - for i = 1,#m1 do - for j = 1,#m1[1] do - if m1[i][j] ~= m2[i][j] then - return false - end - end - end - return true -end - --- Set tostring "tostring( mtx )" behaviour -matrix_meta.__tostring = function( ... ) - return matrix.tostring( ... ) -end - --- set __call "mtx( [formatstr] )" behaviour, mtx [, formatstr] -matrix_meta.__call = function( ... ) - matrix.print( ... ) -end - ---// __index handling -matrix_meta.__index = {} -for k,v in pairs( matrix ) do - matrix_meta.__index[k] = v -end - - ---///////////////////////////////// ---// symbol class implementation ---///////////////////////////////// - --- access to the symbolic metatable -local symbol_meta = {}; symbol_meta.__index = symbol_meta -local symbol = symbol_meta - -function symbol_meta.new(o) - return setmetatable({tostring(o)}, symbol_meta) -end -symbol_meta.to = symbol_meta.new - --- symbol( arg ) --- same as symbol.to( arg ) --- set __call behaviour of symbol -setmetatable( symbol_meta, { __call = function( _,s ) return symbol_meta.to( s ) end } ) - - --- Converts object to string, optionally with formatting. -function symbol_meta.tostring( e,fstr ) - return string.format( fstr,e[1] ) -end - --- Returns "symbol" if object is a symbol type, else nothing. -function symbol_meta:type() - if getmetatable(self) == symbol_meta then - return "symbol" - end -end - --- Performs string.gsub on symbol. --- for use in matrix.replace -function symbol_meta:gsub(from, to) - return symbol.to( string.gsub( self[1],from,to ) ) -end - --- creates function that replaces one letter by something else --- makereplacer( "a",4,"b",7, ... )(x) --- will replace a with 4 and b with 7 in symbol x. --- for use in matrix.replace -function symbol_meta.makereplacer( ... ) - local tosub = {} - local args = {...} - for i = 1,#args,2 do - tosub[args[i]] = args[i+1] - end - local function func( a ) return tosub[a] or a end - return function(sym) - return symbol.to( string.gsub( sym[1], "%a", func ) ) - end -end - --- applies abs function to symbol -function symbol_meta.abs(a) - return symbol.to("(" .. a[1] .. "):abs()") -end - --- applies sqrt function to symbol -function symbol_meta.sqrt(a) - return symbol.to("(" .. a[1] .. "):sqrt()") -end - -function symbol_meta.__add(a,b) - return symbol.to(a .. "+" .. b) -end - -function symbol_meta.__sub(a,b) - return symbol.to(a .. "-" .. b) -end - -function symbol_meta.__mul(a,b) - return symbol.to("(" .. a .. ")*(" .. b .. ")") -end - -function symbol_meta.__div(a,b) - return symbol.to("(" .. a .. ")/(" .. b .. ")") -end - -function symbol_meta.__pow(a,b) - return symbol.to("(" .. a .. ")^(" .. b .. ")") -end - -function symbol_meta.__eq(a,b) - return a[1] == b[1] -end - -function symbol_meta.__tostring(a) - return a[1] -end - -function symbol_meta.__concat(a,b) - return tostring(a) .. tostring(b) -end - -matrix.symbol = symbol - - --- return matrix -return matrix - ---///////////////-- ---// chillcode //-- ---///////////////-- \ No newline at end of file diff --git a/mod.conf b/mod.conf index 0327d56..aa3b0b6 100644 --- a/mod.conf +++ b/mod.conf @@ -1 +1,3 @@ name = rhotator + +depends = matrix