davidm/lua-matrix
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

import version 0.2.9, 2008-08-26

Browse Source 
- decoupled symbol class from matrix class:
	-   matrix.replace has new semantics, applying a function to each element.
	      For old behavior: mtx = mtx:replace(matrix.symbol.makereplacer(...))
	-   replaced mtx:gsub(a,b) with mtx = mtx:replace(symbol.gsub,a,b)
	-   replaced matrix.tosymbol(mtx) with mtx = mtx:replace(symbol)
	- eliminated dependency on complex:
	-   replaced matrix.tocomplex(mtx) with mtx = mtx:replace(complex)
	-   replaced matrix.conjugate(mtx) with mtx = mtx:replace(complex.conjugate)
	-   mulnum and divnum no longer can take num of type string
	-   complex table no longer returned on module load
	- renamed remcomplex to elementstostrings and changed it to return new
	    matrix rather than doing in-place modification
	- fixed matrix.numround (numround variable mispelled).
	    Reported by Goeff Richards.
This commit is contained in:
David Manura 2010-09-22 21:41:53 -04:00
parent 67f724300f
commit fe6837f6dc
14 changed files with 448 additions and 401 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
LICENSE.txt Normal file
View File

@ -0,0 +1,35 @@
LuaMatrix License
-----------
LuaMatrix ( http://luamatrix.luaforge.net/ ) is licensed under the
same terms as Lua (MIT license) reproduced below. This means that
LuaMatrix is free software and can be used for both academic and
commercial purposes at absolutely no cost.
For details and rationale, see http://www.lua.org/license.html .
===============================================================================
Copyright (C) 2007-2010 Michael Lutz.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
===============================================================================
(end of COPYRIGHT)

22
README.txt Normal file
View File

@ -0,0 +1,22 @@
== Description ==
LuaMatrix - Matrices and matrix operations implemented in pure Lua.
This supports operations on matrices and vectors whose elements are
real, complex, or symbolic. Implemented entirely in Lua as tables.
Includes a complex number data type too.
For details on use, see the comments in matrix.lua and complex.lua,
as well as the test suite.
To install, copy matrix.lua and complex.lua into your LUA_PATH. The
modules can alternately be installed via LuaRocks ("luarocks install
luamatrix").
== Project Page ==
http://lua-users.org/wiki/LuaMatrix
== License ==
See the LICENSE.txt file for licensing details.

5
complex_changelog.txt → doc/complex_changelog.txt
View File

@ -1,5 +1,8 @@
complex changelog
complex changelog: v 0.3.1: 2007-11-12
[ David Manura ]
rename mulconjugate to norm2
v 0.3.0: 2007-08-26 v 0.3.0: 2007-08-26
- fixed print function - fixed print function

0
fit_changelog.txt → doc/fit_changelog.txt
View File

62
matrix_changelog.txt → doc/matrix_changelog.txt
View File

@ -1,49 +1,21 @@
matrix function list: matrix changelog
matrix.add v 0.2.9: 2008-08-26
matrix.columns [ David Manura ]
matrix.concath - decoupled symbol class from matrix class:
matrix.concatv - matrix.replace has new semantics, applying a function to each element.
matrix.conjugate For old behavior: mtx = mtx:replace(matrix.symbol.makereplacer(...))
matrix.copy - replaced mtx:gsub(a,b) with mtx = mtx:replace(symbol.gsub,a,b)
matrix.cross - replaced matrix.tosymbol(mtx) with mtx = mtx:replace(symbol)
matrix.det - eliminated dependency on complex:
matrix.div - replaced matrix.tocomplex(mtx) with mtx = mtx:replace(complex)
matrix.divnum - replaced matrix.conjugate(mtx) with mtx = mtx:replace(complex.conjugate)
matrix.dogauss - mulnum and divnum no longer can take num of type string
matrix.getelement - complex table no longer returned on module load
matrix.gsub - renamed remcomplex to elementstostrings and changed it to return new
matrix.invert matrix rather than doing in-place modification
matrix.ipairs - Fixed matrix.numround (numround variable mispelled).
matrix.latex Reported by Goeff Richards.
matrix.len
matrix.mul
matrix.mulnum
matrix:new
matrix.normf
matrix.normmax
matrix.pow
matrix.print
matrix.random
matrix.remcomplex
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.tocomplex
matrix.tostring
matrix.tosymbol
matrix.transpose
matrix.type
v 0.2.8: 2007-08-26 v 0.2.8: 2007-08-26
[ Michael Lutz ] [ Michael Lutz ]

16
complex.lua → lua/complex.lua
View File

@ -1,4 +1,4 @@
-- complex 0.3.0 -- complex 0.3.1
-- Lua 5.1 -- Lua 5.1
-- 'complex' provides common tasks with complex numbers -- 'complex' provides common tasks with complex numbers
@ -15,8 +15,7 @@
-- the access is faster than in a hash table -- the access is faster than in a hash table
-- the metatable is just a add on, when it comes to speed, one is faster using a direct function call -- the metatable is just a add on, when it comes to speed, one is faster using a direct function call
-- http://luaforge.net/projects/LuaMatrix -- http://luamatrix.luaforge.net
-- http://lua-users.org/wiki/ComplexNumbers
-- Licensed under the same terms as Lua itself. -- Licensed under the same terms as Lua itself.
@ -191,9 +190,10 @@ function complex.polardeg( cx )
return math.sqrt( cx[1]^2 + cx[2]^2 ), math.atan2( cx[2], cx[1] ) / math.pi * 180 return math.sqrt( cx[1]^2 + cx[2]^2 ), math.atan2( cx[2], cx[1] ) / math.pi * 180
end end
-- complex.mulconjugate( cx ) -- complex.norm2( cx )
-- multiply with conjugate, function returning a number -- multiply with conjugate, function returning a scalar number
function complex.mulconjugate( cx ) -- norm2(x + i*y) returns x^2 + y^2
function complex.norm2( cx )
return cx[1]^2 + cx[2]^2 return cx[1]^2 + cx[2]^2
end end
@ -330,6 +330,10 @@ function complex.round( cx,idp )
math.floor( cx[2] * mult + 0.5 ) / mult }, complex_meta ) math.floor( cx[2] * mult + 0.5 ) / mult }, complex_meta )
end end
--// variables
complex.zero = complex.new(0, 0)
complex.one = complex.new(1, 0)
--// metatable functions --// metatable functions
complex_meta.__add = function( cx1,cx2 ) complex_meta.__add = function( cx1,cx2 )

530
matrix.lua → lua/matrix.lua
View File

@ -1,5 +1,5 @@
--[[ --[[
matrix v 0.2.8 matrix v 0.2.9
Lua 5.1 compatible Lua 5.1 compatible
@ -30,10 +30,53 @@
where num will be a matrix with the result in mtx[1][1], where num will be a matrix with the result in mtx[1][1],
or use num = vec1:scalar( vec2 ), where num is a number or use num = vec1:scalar( vec2 ), where num is a number
Sites: Site: http://luamatrix.luaforge.net
http://luaforge.net/projects/LuaMatrix
http://lua-users.org/wiki/SimpleMatrix 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
Licensed under the same terms as Lua itself. Licensed under the same terms as Lua itself.
Developers: Developers:
@ -41,10 +84,6 @@
David Manura http://lua-users.org/wiki/DavidManura David Manura http://lua-users.org/wiki/DavidManura
]]-- ]]--
-- for speed and clearer code load the complex function table
-- in there we define the complex number
local complex = require "complex"
--//////////// --////////////
--// matrix // --// matrix //
--//////////// --////////////
@ -54,13 +93,6 @@ local matrix = {}
-- access to the metatable we set at the end of the file -- access to the metatable we set at the end of the file
local matrix_meta = {} local matrix_meta = {}
-- access to the symbolic metatable
local symbol_meta = {}; symbol_meta.__index = symbol_meta
-- set up a symbol type
local function newsymbol(o)
return setmetatable({tostring(o)}, symbol_meta)
end
--///////////////////////////// --/////////////////////////////
--// Get 'new' matrix object // --// Get 'new' matrix object //
--///////////////////////////// --/////////////////////////////
@ -124,7 +156,7 @@ setmetatable( matrix, { __call = function( ... ) return matrix.new( ... ) end }
--// matrix 'matrix' functions // --// matrix 'matrix' functions //
--/////////////////////////////// --///////////////////////////////
--// for real, complx and symbolic matrices //-- --// for real, complex and symbolic matrices //--
-- note: real and complex matrices may be added, subtracted, etc. -- note: real and complex matrices may be added, subtracted, etc.
-- real and symbolic matrices may also be added, subtracted, etc. -- real and symbolic matrices may also be added, subtracted, etc.
@ -132,33 +164,35 @@ setmetatable( matrix, { __call = function( ... ) return matrix.new( ... ) end }
-- since it is not clear which metatable then is used -- since it is not clear which metatable then is used
--// matrix.add ( m1, m2 ) --// matrix.add ( m1, m2 )
-- Add 2 matrices; m2 may be of bigger size than m1 -- Add two matrices; m2 may be of bigger size than m1
function matrix.add( m1, m2 ) function matrix.add( m1, m2 )
local mtx = {} local mtx = {}
for i = 1,#m1 do for i = 1,#m1 do
mtx[i] = {} local m3i = {}
mtx[i] = m3i
for j = 1,#m1[1] do for j = 1,#m1[1] do
mtx[i][j] = m1[i][j] + m2[i][j] m3i[j] = m1[i][j] + m2[i][j]
end end
end end
return setmetatable( mtx, matrix_meta ) return setmetatable( mtx, matrix_meta )
end end
--// matrix.sub ( m1 ,m2 ) --// matrix.sub ( m1 ,m2 )
-- Subtract 2 matrices; m2 may be of bigger size than m1 -- Subtract two matrices; m2 may be of bigger size than m1
function matrix.sub( m1, m2 ) function matrix.sub( m1, m2 )
local mtx = {} local mtx = {}
for i = 1,#m1 do for i = 1,#m1 do
mtx[i] = {} local m3i = {}
mtx[i] = m3i
for j = 1,#m1[1] do for j = 1,#m1[1] do
mtx[i][j] = m1[i][j] - m2[i][j] m3i[j] = m1[i][j] - m2[i][j]
end end
end end
return setmetatable( mtx, matrix_meta ) return setmetatable( mtx, matrix_meta )
end end
--// matrix.mul ( m1, m2 ) --// matrix.mul ( m1, m2 )
-- Multiply 2 matrices; m1 columns must be equal to m2 rows -- Multiply two matrices; m1 columns must be equal to m2 rows
-- e.g. #m1[1] == #m2 -- e.g. #m1[1] == #m2
function matrix.mul( m1, m2 ) function matrix.mul( m1, m2 )
-- multiply rows with columns -- multiply rows with columns
@ -177,7 +211,7 @@ function matrix.mul( m1, m2 )
end end
--// matrix.div ( m1, m2 ) --// matrix.div ( m1, m2 )
-- Divide 2 matrices; m1 columns must be equal to m2 rows -- Divide two matrices; m1 columns must be equal to m2 rows
-- m2 must be square, to be inverted, -- m2 must be square, to be inverted,
-- if that fails returns the rank of m2 as second argument -- if that fails returns the rank of m2 as second argument
-- e.g. #m1[1] == #m2; #m2 == #m2[1] -- e.g. #m1[1] == #m2; #m2 == #m2[1]
@ -189,12 +223,9 @@ end
--// matrix.mulnum ( m1, num ) --// matrix.mulnum ( m1, num )
-- Multiply matrix with a number -- Multiply matrix with a number
-- num may be of type 'number','complex number' or 'string' -- num may be of type 'number' or 'complex number'
-- strings get converted to complex number, if that fails then to symbol -- strings get converted to complex number, if that fails then to symbol
function matrix.mulnum( m1, num ) function matrix.mulnum( m1, num )
if type(num) == "string" then
num = complex.to(num) or newsymbol(num)
end
local mtx = {} local mtx = {}
-- multiply elements with number -- multiply elements with number
for i = 1,#m1 do for i = 1,#m1 do
@ -208,18 +239,16 @@ end
--// matrix.divnum ( m1, num ) --// matrix.divnum ( m1, num )
-- Divide matrix by a number -- Divide matrix by a number
-- num may be of type 'number','complex number' or 'string' -- num may be of type 'number' or 'complex number'
-- strings get converted to complex number, if that fails then to symbol -- strings get converted to complex number, if that fails then to symbol
function matrix.divnum( m1, num ) function matrix.divnum( m1, num )
if type(num) == "string" then
num = complex.to(num) or newsymbol(num)
end
local mtx = {} local mtx = {}
-- divide elements by number -- divide elements by number
for i = 1,#m1 do for i = 1,#m1 do
mtx[i] = {} local mtxi = {}
mtx[i] = mtxi
for j = 1,#m1[1] do for j = 1,#m1[1] do
mtx[i][j] = m1[i][j] / num mtxi[j] = m1[i][j] / num
end end
end end
return setmetatable( mtx, matrix_meta ) return setmetatable( mtx, matrix_meta )
@ -251,6 +280,10 @@ function matrix.pow( m1, num )
return mtx return mtx
end end
local function number_norm2(x)
return x * x
end
--// matrix.det ( m1 ) --// matrix.det ( m1 )
-- Calculate the determinant of a matrix -- Calculate the determinant of a matrix
-- m1 needs to be square -- m1 needs to be square
@ -261,8 +294,6 @@ end
-- here we try to get the nearest element to |1|, (smallest pivot element) -- 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]; -- 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 -- with complex matrices we use the complex.abs function to check if it is bigger or smaller
local fiszerocomplex = function( cx ) return complex.is(cx,0,0) end
local fiszeronumber = function( num ) return num == 0 end
function matrix.det( m1 ) function matrix.det( m1 )
-- check if matrix is quadratic -- check if matrix is quadratic
@ -284,16 +315,10 @@ function matrix.det( m1 )
end end
--// no symbolic matrix supported below here --// no symbolic matrix supported below here
local e = m1[1][1]
local fiszero, abs local zero = type(e) == "table" and e.zero or 0
if matrix.type( m1 ) == "complex" then local norm2 = type(e) == "table" and e.norm2 or number_norm2
fiszero = fiszerocomplex
abs = complex.mulconjugate
else
fiszero = fiszeronumber
abs = math.abs
end
--// matrix is bigger than 3x3 --// matrix is bigger than 3x3
-- get determinant -- get determinant
-- using Gauss elimination and Laplace -- using Gauss elimination and Laplace
@ -313,12 +338,12 @@ function matrix.det( m1 )
-- if no subdet has been found -- if no subdet has been found
if not subdet then if not subdet then
-- check if element it is not zero -- check if element it is not zero
if not fiszero(e) then if e ~= zero then
-- use element as new subdet -- use element as new subdet
subdet,xrow = e,i subdet,xrow = e,i
end end
-- check for elements nearest to 1 or -1 -- check for elements nearest to 1 or -1
elseif (not fiszero(e)) and math.abs(abs(e)-1) < math.abs(abs(subdet)-1) then elseif e ~= zero and math.abs(norm2(e)-1) < math.abs(norm2(subdet)-1) then
subdet,xrow = e,i subdet,xrow = e,i
end end
end end
@ -335,7 +360,7 @@ function matrix.det( m1 )
for i = 1,rows-1 do for i = 1,rows-1 do
-- factor is the dividor of the first element -- factor is the dividor of the first element
-- if element is not already zero -- if element is not already zero
if not fiszero( mtx[i][j] ) then if mtx[i][j] ~= zero then
local factor = mtx[i][j]/subdet local factor = mtx[i][j]/subdet
-- update all remaining fields of the matrix, with value from xrow -- update all remaining fields of the matrix, with value from xrow
for n = j+1,#mtx[1] do for n = j+1,#mtx[1] do
@ -366,23 +391,23 @@ end
-- locals -- locals
-- checking here for the nearest element to 1 or -1; (smallest pivot element) -- checking here for the nearest element to 1 or -1; (smallest pivot element)
-- this way the factor of the evolving number division should be > 1 or the divided number itself, -- this way the factor of the evolving number division should be > 1 or the
-- what gives better results -- divided number itself, what gives better results
local setelementtosmallest = function( mtx,i,j,fiszero,fisone,abs ) local setelementtosmallest = function( mtx,i,j,zero,one,norm2 )
-- check if element is one -- check if element is one
if fisone(mtx[i][j]) then return true end if mtx[i][j] == one then return true end
-- check for lowest value -- check for lowest value
local _ilow local _ilow
for _i = i,#mtx do for _i = i,#mtx do
local e = mtx[_i][j] local e = mtx[_i][j]
if fisone(e) then if e == one then
break break
end end
if not _ilow then if not _ilow then
if not fiszero(e) then if e ~= zero then
_ilow = _i _ilow = _i
end end
elseif (not fiszero(e)) and math.abs(abs(e)-1) < math.abs(abs(mtx[_ilow][j])-1) then elseif (e ~= zero) and math.abs(norm2(e)-1) < math.abs(norm2(mtx[_ilow][j])-1) then
_ilow = _i _ilow = _i
end end
end end
@ -395,40 +420,28 @@ local setelementtosmallest = function( mtx,i,j,fiszero,fisone,abs )
return true return true
end end
end end
local cxfiszero = function( cx ) return complex.is(cx,0,0) end
local cxfsetzero = function( mtx,i,j ) complex.set(mtx[i][j],0,0) end local function copy(x)
local cxfisone = function( cx ) return complex.abs(cx) == 1 end return type(x) == "table" and x.copy(x) or x
local cxfsetone = function( mtx,i,j ) complex.set(mtx[i][j],1,0) end end
local numfiszero = function( num ) return num == 0 end
local numfsetzero = function( mtx,i,j ) mtx[i][j] = 0 end
local numfisone = function( num ) return math.abs(num) == 1 end
local numfsetone = function( mtx,i,j ) mtx[i][j] = 1 end
-- note: in --// ... //-- we have a way that does no divison, -- note: in --// ... //-- we have a way that does no divison,
-- however with big number and matrices we get problems since we do no reducing -- however with big number and matrices we get problems since we do no reducing
function matrix.dogauss( mtx ) function matrix.dogauss( mtx )
local fiszero,fsetzero,fisone,fsetone,abs local e = mtx[1][1]
if matrix.type( mtx ) == "complex" then local zero = type(e) == "table" and e.zero or 0
fiszero = cxfiszero local one = type(e) == "table" and e.one or 1
fsetzero = cxfsetzero local norm2 = type(e) == "table" and e.norm2 or number_norm2
fisone = cxfisone
fsetone = cxfsetone
abs = complex.mulconjugate
else
fiszero = numfiszero
fsetzero = numfsetzero
fisone = numfisone
fsetone = numfsetone
abs = math.abs
end
local rows,columns = #mtx,#mtx[1] local rows,columns = #mtx,#mtx[1]
-- stairs left -> right -- stairs left -> right
for j = 1,rows do for j = 1,rows do
-- check if element can be setted to one -- check if element can be setted to one
if setelementtosmallest( mtx,j,j,fiszero,fisone,abs ) then if setelementtosmallest( mtx,j,j,zero,one,norm2 ) then
-- start parsing rows -- start parsing rows
for i = j+1,rows do for i = j+1,rows do
-- check if element is not already zero -- check if element is not already zero
if not fiszero(mtx[i][j]) then if mtx[i][j] ~= zero then
-- we may add x*otherline row, to set element to zero -- we may add x*otherline row, to set element to zero
-- tozero - x*mtx[j][j] = 0; x = tozero/mtx[j][j] -- tozero - x*mtx[j][j] = 0; x = tozero/mtx[j][j]
local factor = mtx[i][j]/mtx[j][j] local factor = mtx[i][j]/mtx[j][j]
@ -436,7 +449,7 @@ function matrix.dogauss( mtx )
-- yet with big matrices (since we do no reducing and other things) -- yet with big matrices (since we do no reducing and other things)
-- we get too big numbers -- we get too big numbers
--local factor1,factor2 = mtx[i][j],mtx[j][j] //-- --local factor1,factor2 = mtx[i][j],mtx[j][j] //--
fsetzero(mtx,i,j) mtx[i][j] = copy(zero)
for _j = j+1,columns do for _j = j+1,columns do
--// mtx[i][_j] = mtx[i][_j] * factor2 - factor1 * mtx[j][_j] //-- --// mtx[i][_j] = mtx[i][_j] * factor2 - factor1 * mtx[j][_j] //--
mtx[i][_j] = mtx[i][_j] - factor * mtx[j][_j] mtx[i][_j] = mtx[i][_j] - factor * mtx[j][_j]
@ -459,15 +472,15 @@ function matrix.dogauss( mtx )
-- start parsing rows -- start parsing rows
for i = j-1,1,-1 do for i = j-1,1,-1 do
-- check if element is not already zero -- check if element is not already zero
if not fiszero(mtx[i][j]) then if mtx[i][j] ~= zero then
local factor = mtx[i][j] local factor = mtx[i][j]
for _j = j+1,columns do for _j = j+1,columns do
mtx[i][_j] = mtx[i][_j] - factor * mtx[j][_j] mtx[i][_j] = mtx[i][_j] - factor * mtx[j][_j]
end end
fsetzero(mtx,i,j) mtx[i][j] = copy(zero)
end end
end end
fsetone(mtx,j,j) mtx[j][j] = copy(one)
end end
return true return true
end end
@ -481,27 +494,14 @@ function matrix.invert( m1 )
assert(#m1 == #m1[1], "matrix not square") assert(#m1 == #m1[1], "matrix not square")
local mtx = matrix.copy( m1 ) local mtx = matrix.copy( m1 )
local ident = setmetatable( {},matrix_meta ) local ident = setmetatable( {},matrix_meta )
if matrix.type( mtx ) == "complex" then local e = m1[1][1]
for i = 1,#m1 do local zero = type(e) == "table" and e.zero or 0
ident[i] = {} local one = type(e) == "table" and e.one or 1
for j = 1,#m1 do for i = 1,#m1 do
if i == j then local identi = {}
ident[i][j] = complex.new( 1,0 ) ident[i] = identi
else for j = 1,#m1 do
ident[i][j] = complex.new( 0,0 ) identi[j] = copy((i == j) and one or zero)
end
end
end
else
for i = 1,#m1 do
ident[i] = {}
for j = 1,#m1 do
if i == j then
ident[i][j] = 1
else
ident[i][j] = 0
end
end
end end
end end
mtx = matrix.concath( mtx,ident ) mtx = matrix.concath( mtx,ident )
@ -527,7 +527,8 @@ end
-- local average error -- local average error
local function get_abs_avg( m1, m2 ) local function get_abs_avg( m1, m2 )
local dist = 0 local dist = 0
local abs = matrix.type(m1) == "complex" and complex.abs or math.abs local e = m1[1][1]
local abs = type(e) == "table" and e.abs or math.abs
for i=1,#m1 do for i=1,#m1 do
for j=1,#m1[1] do for j=1,#m1[1] do
dist = dist + abs(m1[i][j]-m2[i][j]) dist = dist + abs(m1[i][j]-m2[i][j])
@ -650,7 +651,7 @@ local tround = function( t,mult )
end end
function matrix.round( mtx, idp ) function matrix.round( mtx, idp )
local mult = 10^( idp or 0 ) local mult = 10^( idp or 0 )
local fround = matrix.type( mtx ) == "number" and numound or tround local fround = matrix.type( mtx ) == "number" and numround or tround
for i = 1,#mtx do for i = 1,#mtx do
for j = 1,#mtx[1] do for j = 1,#mtx[1] do
mtx[i][j] = fround(mtx[i][j],mult) mtx[i][j] = fround(mtx[i][j],mult)
@ -691,12 +692,10 @@ end
--// matrix.type ( mtx ) --// matrix.type ( mtx )
-- get type of matrix, normal/complex/symbol or tensor -- get type of matrix, normal/complex/symbol or tensor
function matrix.type( mtx ) function matrix.type( mtx )
if type(mtx[1][1]) == "table" then local e = mtx[1][1]
if complex.type(mtx[1][1]) then if type(e) == "table" then
return "complex" if e.type then
end return e:type()
if getmetatable(mtx[1][1]) == symbol_meta then
return "symbol"
end end
return "tensor" return "tensor"
end end
@ -766,7 +765,7 @@ function matrix.subm( m1,i1,j1,i2,j2 )
end end
--// matrix.concath( m1, m2 ) --// matrix.concath( m1, m2 )
-- Concatenate 2 matrices, horizontal -- Concatenate two matrices, horizontal
-- will return m1m2; rows have to be the same -- will return m1m2; rows have to be the same
-- e.g.: #m1 == #m2 -- e.g.: #m1 == #m2
function matrix.concath( m1,m2 ) function matrix.concath( m1,m2 )
@ -787,7 +786,7 @@ function matrix.concath( m1,m2 )
end end
--// matrix.concatv ( m1, m2 ) --// matrix.concatv ( m1, m2 )
-- Concatenate 2 matrices, vertical -- Concatenate two matrices, vertical
-- will return m1 -- will return m1
-- m2 -- m2
-- columns have to be the same; e.g.: #m1[1] == #m2[1] -- columns have to be the same; e.g.: #m1[1] == #m2[1]
@ -838,59 +837,32 @@ function matrix.rotr( m1 )
return mtx return mtx
end end
-- local get_elemnts in string local function tensor_tostring( t,fstr )
local get_tstr = function( t ) if not fstr then return "["..table.concat(t,",").."]" end
return "["..table.concat(t,",").."]"
end
local get_str = function( e )
return tostring(e)
end
-- local get_elemnts in string and formated
local getf_tstr = function( t,fstr )
local tval = {} local tval = {}
for i,v in ipairs( t ) do for i,v in ipairs( t ) do
tval[i] = string.format( fstr,v ) tval[i] = string.format( fstr,v )
end end
return "["..table.concat(tval,",").."]" return "["..table.concat(tval,",").."]"
end end
local getf_cxstr = function( e,fstr ) local function number_tostring( e,fstr )
return complex.tostring( e,fstr ) return fstr and string.format( fstr,e ) or e
end
local getf_symstr = function( e,fstr )
return string.format( fstr,e[1] )
end
local getf_str = function( e,fstr )
return string.format( fstr,e )
end end
--// matrix.tostring ( mtx, formatstr ) --// matrix.tostring ( mtx, formatstr )
-- tostring function -- tostring function
function matrix.tostring( mtx, formatstr ) function matrix.tostring( mtx, formatstr )
local ts = {} local ts = {}
local getstr local mtype = matrix.type( mtx )
if formatstr then -- get str formatted local e = mtx[1][1]
local mtype = matrix.type( mtx ) local tostring = mtype == "tensor" and tensor_tostring or
if mtype == "tensor" then getstr = getf_tstr type(e) == "table" and e.tostring or number_tostring
elseif mtype == "complex" then getstr = getf_cxstr for i = 1,#mtx do
elseif mtype == "symbol" then getstr = getf_symstr local tstr = {}
else getstr = getf_str end for j = 1,#mtx[1] do
-- iteratr tstr[j] = tostring(mtx[i][j],formatstr)
for i = 1,#mtx do
local tstr = {}
for j = 1,#mtx[1] do
tstr[j] = getstr(mtx[i][j],formatstr)
end
ts[i] = table.concat(tstr, "\t")
end
else
getstr = matrix.type( mtx ) == "tensor" and get_tstr or get_str
for i = 1,#mtx do
local tstr = {}
for j = 1,#mtx[1] do
tstr[j] = getstr(mtx[i][j])
end
ts[i] = table.concat(tstr, "\t")
end end
ts[i] = table.concat(tstr, "\t")
end end
return table.concat(ts, "\n") return table.concat(ts, "\n")
end end
@ -909,7 +881,7 @@ function matrix.latex( mtx, align )
-- \usepackage{dcolumn}; D{.}{,}{-1}; aligns number by . replaces it with , -- \usepackage{dcolumn}; D{.}{,}{-1}; aligns number by . replaces it with ,
local align = align or "c" local align = align or "c"
local str = "$\\left( \\begin{array}{"..string.rep( align, #mtx[1] ).."}\n" local str = "$\\left( \\begin{array}{"..string.rep( align, #mtx[1] ).."}\n"
local getstr = matrix.type( mtx ) == "tensor" and get_tstr or get_str local getstr = matrix.type( mtx ) == "tensor" and tensor_tostring or number_tostring
for i = 1,#mtx do for i = 1,#mtx do
str = str.."\t"..getstr(mtx[i][1]) str = str.."\t"..getstr(mtx[i][1])
for j = 2,#mtx[1] do for j = 2,#mtx[1] do
@ -1014,96 +986,29 @@ function matrix.len( m1 )
return math.sqrt( m1[1][1]^2 + m1[2][1]^2 + m1[3][1]^2 ) return math.sqrt( m1[1][1]^2 + m1[2][1]^2 + m1[3][1]^2 )
end end
--////////////////////////////////
--// matrix 'complex' functions //
--////////////////////////////////
--// matrix.tocomplex ( mtx ) --// matrix.replace (mtx, func, ...)
-- we set now all elements to a complex number -- for each element e in the matrix mtx, replace it with func(mtx, ...).
-- also set the metatable function matrix.replace( m1, func, ... )
function matrix.tocomplex( mtx ) local mtx = {}
assert( matrix.type(mtx) == "number", "matrix not of type 'number'" ) for i = 1,#m1 do
for i = 1,#mtx do local m1i = m1[i]
for j = 1,#mtx[1] do local mtxi = {}
mtx[i][j] = complex.to( mtx[i][j] ) for j = 1,#m1i do
mtxi[j] = func( m1i[j], ... )
end end
mtx[i] = mtxi
end end
return setmetatable( mtx, matrix_meta ) return setmetatable( mtx, matrix_meta )
end end
--// matrix.remcomplex ( mtx ) --// matrix.remcomplex ( mtx )
-- set the matrix elements to a number or complex number string -- set the matrix elements to strings
function matrix.remcomplex( mtx ) -- IMPROVE: tostring v.s. tostringelements confusing
assert( matrix.type(mtx) == "complex", "matrix not of type 'complex'" ) function matrix.elementstostrings( mtx )
for i = 1,#mtx do local e = mtx[1][1]
for j = 1,#mtx[1] do local tostring = type(e) == "table" and e.tostring or tostring
mtx[i][j] = complex.tostring( mtx[i][j] ) return matrix.replace(mtx, tostring)
end
end
return setmetatable( mtx, matrix_meta )
end
--// matrix.conjugate ( m1 )
-- get the conjugate complex matrix
function matrix.conjugate( m1 )
assert( matrix.type(m1) == "complex", "matrix not of type 'complex'" )
local mtx = {}
for i = 1,#m1 do
mtx[i] = {}
for j = 1,#m1[1] do
mtx[i][j] = complex.conjugate( m1[i][j] )
end
end
return setmetatable( mtx, matrix_meta )
end
--/////////////////////////////////
--// matrix 'symbol' functions //
--/////////////////////////////////
--// matrix.tosymbol ( mtx )
-- set the matrix elements to symbolic values
function matrix.tosymbol( mtx )
assert( matrix.type( mtx ) ~= "tensor", "cannot convert type 'tensor' to 'symbol'" )
for i = 1,#mtx do
for j = 1,#mtx[1] do
mtx[i][j] = newsymbol( mtx[i][j] )
end
end
return setmetatable( mtx, matrix_meta )
end
--// matrix.gsub( m1, from, to )
-- perform gsub on all elements
function matrix.gsub( m1,from,to )
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] = newsymbol( string.gsub( m1[i][j][1],from,to ) )
end
end
return setmetatable( mtx, matrix_meta )
end
--// matrix.replace ( m1, ... )
-- replace one letter by something else
-- replace( "a",4,"b",7, ... ) will replace a with 4 and b with 7
function matrix.replace( m1,... )
assert( matrix.type( m1 ) == "symbol", "matrix not of type 'symbol'" )
local tosub,args = {},{...}
for i = 1,#args,2 do
tosub[args[i]] = args[i+1]
end
local mtx = {}
for i = 1,#m1 do
mtx[i] = {}
for j = 1,#m1[1] do
mtx[i][j] = newsymbol( string.gsub( m1[i][j][1], "%a", function( a ) return tosub[a] or a end ) )
end
end
return setmetatable( mtx, matrix_meta )
end end
--// matrix.solve ( m1 ) --// matrix.solve ( m1 )
@ -1120,46 +1025,6 @@ function matrix.solve( m1 )
return setmetatable( mtx, matrix_meta ) return setmetatable( mtx, matrix_meta )
end end
function symbol_meta.__add(a,b)
return newsymbol(a .. "+" .. b)
end
function symbol_meta.__sub(a,b)
return newsymbol(a .. "-" .. b)
end
function symbol_meta.__mul(a,b)
return newsymbol("(" .. a .. ")*(" .. b .. ")")
end
function symbol_meta.__div(a,b)
return newsymbol("(" .. a .. ")/(" .. b .. ")")
end
function symbol_meta.__pow(a,b)
return newsymbol("(" .. 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
function symbol_meta.abs(a)
return newsymbol("(" .. a[1] .. "):abs()")
end
function symbol_meta.sqrt(a)
return newsymbol("(" .. a[1] .. "):sqrt()")
end
--////////////////////////-- --////////////////////////--
--// METATABLE HANDLING //-- --// METATABLE HANDLING //--
--////////////////////////-- --////////////////////////--
@ -1232,7 +1097,7 @@ matrix_meta.__eq = function( m1, m2 )
if #m1 ~= #m2 or #m1[1] ~= #m2[1] then if #m1 ~= #m2 or #m1[1] ~= #m2[1] then
return false return false
end end
-- check normal,complex and symbolic -- check elements equal
for i = 1,#m1 do for i = 1,#m1 do
for j = 1,#m1[1] do for j = 1,#m1[1] do
if m1[i][j] ~= m2[i][j] then if m1[i][j] ~= m2[i][j] then
@ -1259,8 +1124,107 @@ for k,v in pairs( matrix ) do
matrix_meta.__index[k] = v matrix_meta.__index[k] = v
end end
-- return the matrix and complex
return matrix, complex --/////////////////////////////////
--// 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 //-- --// chillcode //--

31
rockspec Normal file
View File

@ -0,0 +1,31 @@
package = "LuaMatrix"
version = "[VERSION]"
source = {
url = "[URL]",
}
description = {
summary = "Matrices and matrix operations implemented in pure Lua.",
detailed = [[
This supports operations on matrices and vectors whose elements are
real, complex, or symbolic. Implemented entirely in Lua as tables.
Includes a complex number data type too.
]],
license = "MIT/X11",
homepage = "http://luamatrix.luaforge.net/",
maintainer = "David Manura <http://lua-users.org/wiki/DavidManura>",
}
dependencies = {
"lua >= 5.1",
}
build = {
type = "none",
install = {
lua = {
["complex"] = "lua/complex.lua",
["matrix"] = "lua/matrix.lua",
}
},
copy_directories = {"doc", "samples", "tests"},
}
-- test: tests/test.lua

3
fit.lua → samples/fit.lua
View File

@ -8,8 +8,7 @@
-- little add-on to the matrix module, to show some curve fitting -- little add-on to the matrix module, to show some curve fitting
-- http://luaforge.net/projects/LuaMatrix -- http://luamatrix.luaforge.net
-- http://lua-users.org/wiki/SimpleFit
-- Licensed under the same terms as Lua itself. -- Licensed under the same terms as Lua itself.

22
test_fit.lua
View File

@ -1,22 +0,0 @@
-- require fit
local fit = require "fit"
print( "Fit a straight line " )
-- x(i) = 2 | 3 | 4 | 5
-- y(i) = 5 | 9 | 15 | 21
-- model = y = a + b * x
-- r(i) = y(i) - ( a + b * x(i) )
local a,b = fit.linear( { 2,3, 4, 5 },
{ 5,9,15,21 } )
print( "=> y = ( "..a.." ) + ( "..b.." ) * x")
print( "Fit a parabola " )
local a, b, c = fit.parabola( { 0,1,2,4,6 },
{ 3,1,0,1,4 } )
print( "=> y = ( "..a.." ) + ( "..b.." ) * x + ( "..c.." ) * x²")
print( "Fit exponential" )
local a, b = fit.exponential( {1, 2, 3, 4, 5},
{1,3.1,5.6,9.1,12.9} )
print( "=> y = ( "..a.." ) * x^( "..b.." )")

8
tests/test.lua Normal file
View File

@ -0,0 +1,8 @@
-- test suite (run from parent directory).
package.path = './lua/?.lua;' .. package.path
dofile 'tests/test_complex.lua'
dofile 'tests/test_matrix.lua'
dofile 'tests/test_fit.lua'
print 'ALL PASSED'

10
test_complex.lua → tests/test_complex.lua
View File

@ -82,9 +82,9 @@ local r,phi = complex.polardeg( {0,-3} )
assert( r == 3 ) assert( r == 3 )
assert( phi == -90 ) assert( phi == -90 )
-- complex.mulconjugate( cx ) -- complex.norm2( cx )
cx = complex "2+3i" cx = complex "2+3i"
assert( complex.mulconjugate( cx ) == 13 ) assert( complex.norm2( cx ) == 13 )
-- complex.abs( cx ) -- complex.abs( cx )
cx = complex "3+4i" cx = complex "3+4i"
@ -168,7 +168,7 @@ assert( cx:ln():round( 4 ) == complex "1.6094+0.9273i" )
-- complex.exp( cx ) -- complex.exp( cx )
cx = complex "2+3i" cx = complex "2+3i"
assert( cx:ln():exp() == complex "2+3i" ) assert( cx.abs( cx:ln():exp() - complex "2+3i" ) < 1e-7 )
-- complex.conjugate( cx ) -- complex.conjugate( cx )
cx = complex "2+3i" cx = complex "2+3i"
@ -182,4 +182,6 @@ assert( cx+2 == complex "4+3i" )
-- __unm -- __unm
cx = complex "2+3i" cx = complex "2+3i"
assert( -cx == complex "-2-3i" ) assert( -cx == complex "-2-3i" )
print 'PASSED'

25
tests/test_fit.lua Normal file
View File

@ -0,0 +1,25 @@
-- require fit
package.path = "samples/?.lua;" .. package.path
local fit = require "fit"
-- Fit a straight line
-- x(i) = 2 | 3 | 4 | 5
-- y(i) = 5 | 9 | 15 | 21
-- model = y = a + b * x
-- r(i) = y(i) - ( a + b * x(i) )
local a,b = fit.linear( { 2,3, 4, 5 }, { 5,9,15,21 } )
assert(math.abs(a - -6.4) < 0.001)
assert(math.abs(b - 5.4) < 0.001)
-- Fit a parabola
local a, b, c = fit.parabola( { 0,1,2,4,6 }, { 3,1,0,1,4 } )
assert(math.abs(a - 2.8251599147122) < 0.001)
assert(math.abs(b - -2.0490405117271) < 0.001)
assert(math.abs(c - 0.3773987206823) < 0.001)
-- Fit exponential
local a, b = fit.exponential( {1, 2, 3, 4, 5}, {1,3.1,5.6,9.1,12.9} )
assert(math.abs(a - 1.0077958966968) < 0.001)
assert(math.abs(b - 1.5834684450364) < 0.001)
print 'PASSED'

80
test_matrix.lua → tests/test_matrix.lua
View File

@ -1,4 +1,6 @@
local matrix, complex = require "matrix" local matrix = require "matrix"
local complex = require "complex"
local symbol = matrix.symbol
local mtx, m1,m2,m3,m4,m5, ms,ms1,ms2,ms3,ms4 local mtx, m1,m2,m3,m4,m5, ms,ms1,ms2,ms3,ms4
@ -24,42 +26,42 @@ m1 = matrix{{8,4,1},{6,8,3}}
m2 = matrix{{-8,1,3},{5,2,1}} m2 = matrix{{-8,1,3},{5,2,1}}
assert(m1 + m2 == matrix{{0,5,4},{11,10,4}}) assert(m1 + m2 == matrix{{0,5,4},{11,10,4}})
-- matrix.add; complex -- matrix.add; complex
m1 = matrix{{10,"2+6i",1},{5,1,"4-2i"}}:tocomplex() m1 = matrix{{10,"2+6i",1},{5,1,"4-2i"}}:replace(complex)
m2 = matrix{{3,4,5},{"2+3i",4,"6i"}}:tocomplex() m2 = matrix{{3,4,5},{"2+3i",4,"6i"}}:replace(complex)
assert(m1 + m2 == matrix{{13,"6+6i",6},{"7+3i",5,"4+4i"}}:tocomplex()) assert(m1 + m2 == matrix{{13,"6+6i",6},{"7+3i",5,"4+4i"}}:replace(complex))
-- matrix.add; symbol -- matrix.add; symbol
m1 = matrix{{8,4,1},{6,8,3}}:tosymbol() m1 = matrix{{8,4,1},{6,8,3}}:replace(symbol)
m2 = matrix{{-8,1,3},{5,2,1}}:tosymbol() m2 = matrix{{-8,1,3},{5,2,1}}:replace(symbol)
assert(m1 + m2 == matrix{{"8+-8","4+1","1+3"},{"6+5","8+2","3+1"}}:tosymbol()) assert(m1 + m2 == matrix{{"8+-8","4+1","1+3"},{"6+5","8+2","3+1"}}:replace(symbol))
-- matrix.sub; number -- matrix.sub; number
m1 = matrix{{8,4,1},{6,8,3}} m1 = matrix{{8,4,1},{6,8,3}}
m2 = matrix{{-8,1,3},{5,2,1}} m2 = matrix{{-8,1,3},{5,2,1}}
assert(m1 - m2 == matrix{{16,3,-2},{1,6,2}}) assert(m1 - m2 == matrix{{16,3,-2},{1,6,2}})
-- matrix.sub; complex -- matrix.sub; complex
m1 = matrix{{10,"2+6i",1},{5,1,"4-2i"}}:tocomplex() m1 = matrix{{10,"2+6i",1},{5,1,"4-2i"}}:replace(complex)
m2 = matrix{{3,4,5},{"2+3i",4,"6i"}}:tocomplex() m2 = matrix{{3,4,5},{"2+3i",4,"6i"}}:replace(complex)
assert(m1 - m2 == matrix{{7,"-2+6i",-4},{"3-3i",-3,"4-8i"}}:tocomplex()) assert(m1 - m2 == matrix{{7,"-2+6i",-4},{"3-3i",-3,"4-8i"}}:replace(complex))
-- matrix.sub; symbol -- matrix.sub; symbol
m1 = matrix{{8,4,1},{6,8,3}}:tosymbol() m1 = matrix{{8,4,1},{6,8,3}}:replace(symbol)
m2 = matrix{{-8,1,3},{5,2,1}}:tosymbol() m2 = matrix{{-8,1,3},{5,2,1}}:replace(symbol)
assert(m1 - m2 == matrix{{"8--8","4-1","1-3"},{"6-5","8-2","3-1"}}:tosymbol()) assert(m1 - m2 == matrix{{"8--8","4-1","1-3"},{"6-5","8-2","3-1"}}:replace(symbol))
-- matrix.mul; number -- matrix.mul; number
m1 = matrix{{8,4,1},{6,8,3}} m1 = matrix{{8,4,1},{6,8,3}}
m2 = matrix{{3,1},{2,5},{7,4}} m2 = matrix{{3,1},{2,5},{7,4}}
assert(m1 * m2 == matrix{{39,32},{55,58}}) assert(m1 * m2 == matrix{{39,32},{55,58}})
-- matrix.mul; complex -- matrix.mul; complex
m1 = matrix{{"1+2i","3-i"},{"2-2i","1+i"}}:tocomplex() m1 = matrix{{"1+2i","3-i"},{"2-2i","1+i"}}:replace(complex)
m2 = matrix{{"i","5-i"},{2,"1-i"}}:tocomplex() m2 = matrix{{"i","5-i"},{2,"1-i"}}:replace(complex)
assert( m1*m2 == matrix{{"4-i","9+5i"},{"4+4i","10-12i"}}:tocomplex() ) assert( m1*m2 == matrix{{"4-i","9+5i"},{"4+4i","10-12i"}}:replace(complex) )
-- matrix.mul; symbol -- matrix.mul; symbol
m1 = matrix{{8,4,1},{6,8,3}}:tosymbol() m1 = matrix{{8,4,1},{6,8,3}}:replace(symbol)
m2 = matrix{{3,1},{2,5},{7,4}}:tosymbol() m2 = matrix{{3,1},{2,5},{7,4}}:replace(symbol)
assert(m1 * m2 == matrix{ assert(m1 * m2 == matrix{
{"(8)*(3)+(4)*(2)+(1)*(7)", "(8)*(1)+(4)*(5)+(1)*(4)"}, {"(8)*(3)+(4)*(2)+(1)*(7)", "(8)*(1)+(4)*(5)+(1)*(4)"},
{"(6)*(3)+(8)*(2)+(3)*(7)", "(6)*(1)+(8)*(5)+(3)*(4)"} {"(6)*(3)+(8)*(2)+(3)*(7)", "(6)*(1)+(8)*(5)+(3)*(4)"}
}:tosymbol()) }:replace(symbol))
-- matrix.div; number, same for complex, not for symbol -- matrix.div; number, same for complex, not for symbol
m1 = matrix {{1,2},{3,4}} m1 = matrix {{1,2},{3,4}}
@ -72,9 +74,9 @@ assert( m2/2 == matrix{{2,2.5},{3,3.5}} )
mtx = matrix {{3,5,1},{2,4,5},{1,2,2}} mtx = matrix {{3,5,1},{2,4,5},{1,2,2}}
assert( 2 / mtx == matrix{{4,16,-42},{-2,-10,26},{0,2,-4}} ) assert( 2 / mtx == matrix{{4,16,-42},{-2,-10,26},{0,2,-4}} )
-- matrix.mulnum; symbol -- matrix.mulnum; symbol
m1 = m1:tosymbol() m1 = m1:replace(symbol)
assert( m1*2 == matrix{{"(1)*(2)","(2)*(2)"},{"(3)*(2)","(4)*(2)"}}:tosymbol() ) assert( m1*2 == matrix{{"(1)*(2)","(2)*(2)"},{"(3)*(2)","(4)*(2)"}}:replace(symbol) )
assert( m1/2 == matrix{{"(1)/(2)","(2)/(2)"},{"(3)/(2)","(4)/(2)"}}:tosymbol() ) assert( m1/2 == matrix{{"(1)/(2)","(2)/(2)"},{"(3)/(2)","(4)/(2)"}}:replace(symbol) )
-- matrix.pow; number, same complex -- matrix.pow; number, same complex
mtx = matrix{{3,5,1},{2,4,5},{1,2,2}} mtx = matrix{{3,5,1},{2,4,5},{1,2,2}}
@ -94,8 +96,8 @@ assert(select(2, pcall(function() return mtx^-1 end))
mtx = matrix {{1,4,3,2},{2,1,-1,-1},{-3,2,2,-2},{-1,-5,-4,1}} mtx = matrix {{1,4,3,2},{2,1,-1,-1},{-3,2,2,-2},{-1,-5,-4,1}}
assert( mtx:det() == 78 ) assert( mtx:det() == 78 )
-- matrix.det; complex -- matrix.det; complex
m1 = matrix{{"1+2i","3-i"},{"2-2i","1+i"}}:tocomplex() m1 = matrix{{"1+2i","3-i"},{"2-2i","1+i"}}:replace(complex)
m2 = matrix{{"i","5-i"},{2,"1-i"}}:tocomplex() m2 = matrix{{"i","5-i"},{2,"1-i"}}:replace(complex)
m3 = m1*m2 m3 = m1*m2
-- (checked in maple) -- (checked in maple)
assert( m3:det() == complex "12-114i" ) assert( m3:det() == complex "12-114i" )
@ -104,7 +106,7 @@ mtx = {{"2+3i","1+4i","-2i",3,2},
{3,"-2i",6,"4+5i",0}, {3,"-2i",6,"4+5i",0},
{1,"1+2i",3,5,7}, {1,"1+2i",3,5,7},
{"-3+3i","3+3i",3,-8,2}} {"-3+3i","3+3i",3,-8,2}}
matrix(mtx):tocomplex() mtx = matrix(mtx):replace(complex)
-- (checked in maple) -- (checked in maple)
assert( mtx:det():round(10) == complex "5527+2687i" ) assert( mtx:det():round(10) == complex "5527+2687i" )
@ -123,7 +125,7 @@ mtx = {
{3,"-1+5i",-3}, {3,"-1+5i",-3},
{4,0,7}, {4,0,7},
} }
matrix.tocomplex( mtx ) mtx = matrix.replace( mtx, complex )
local mtxinv = mtx^-1 local mtxinv = mtx^-1
local mtxinvcomp = { local mtxinvcomp = {
{"0.13349-0.07005i","0.14335+0.03609i","0.04237+0.02547i"}, {"0.13349-0.07005i","0.14335+0.03609i","0.04237+0.02547i"},
@ -131,7 +133,7 @@ local mtxinvcomp = {
{"-0.07628+0.04003i","-0.08192-0.02062i","0.11865-0.01456i"},} {"-0.07628+0.04003i","-0.08192-0.02062i","0.11865-0.01456i"},}
mtxinvcomp = matrix( mtxinvcomp ) mtxinvcomp = matrix( mtxinvcomp )
mtxinv:round( 5 ) mtxinv:round( 5 )
mtxinv:remcomplex() mtxinv = mtxinv:elementstostrings()
assert( mtxinvcomp == mtxinv ) assert( mtxinvcomp == mtxinv )
-- matrix.sqrt; number -- matrix.sqrt; number
@ -140,7 +142,7 @@ local m2 = m1*m1
local msqrt = m2:sqrt() local msqrt = m2:sqrt()
assert((m2 - msqrt^2):normmax() < 1E-12) assert((m2 - msqrt^2):normmax() < 1E-12)
-- matrix.sqrt; complex -- matrix.sqrt; complex
local m1 = matrix{{4,"2+i",1},{1,5,"4-2i"},{1,"5+3i",2}}:tocomplex() local m1 = matrix{{4,"2+i",1},{1,5,"4-2i"},{1,"5+3i",2}}:replace(complex)
local m2 = m1*m1 local m2 = m1*m1
local msqrt = m2:sqrt() local msqrt = m2:sqrt()
assert((m2 - msqrt^2):normmax() < 1E-12) assert((m2 - msqrt^2):normmax() < 1E-12)
@ -152,7 +154,7 @@ local m2 = m1^p
local mroot = m2:root(p) local mroot = m2:root(p)
assert((m2 - mroot^p):normmax() < 1E-7) assert((m2 - mroot^p):normmax() < 1E-7)
-- matrix.root; complex -- matrix.root; complex
local m1 = matrix{{4,"2+i",1},{1,5,"4-2i"},{1,"5+3i",2}}:tocomplex() local m1 = matrix{{4,"2+i",1},{1,5,"4-2i"},{1,"5+3i",2}}:replace(complex)
local m2 = m1^p local m2 = m1^p
local mroot = m2:root(p) local mroot = m2:root(p)
assert((m2 - mroot^p):normmax() < 1E-7) assert((m2 - mroot^p):normmax() < 1E-7)
@ -160,16 +162,16 @@ assert((m2 - mroot^p):normmax() < 1E-7)
-- matrix.normf -- matrix.normf
mtx = matrix{{2,3},{-2,-3}} mtx = matrix{{2,3},{-2,-3}}
assert(mtx:normf() == math.sqrt(2^2+3^2+2^2+3^2)) assert(mtx:normf() == math.sqrt(2^2+3^2+2^2+3^2))
mtx = matrix{{'2i','3'},{'-2i','-3'}}:tocomplex() mtx = matrix{{'2i','3'},{'-2i','-3'}}:replace(complex)
assert(mtx:normf() == math.sqrt(2^2+3^2+2^2+3^2)) assert(mtx:normf() == math.sqrt(2^2+3^2+2^2+3^2))
mtx = matrix{{'a','b'},{'c','d'}}:tosymbol() mtx = matrix{{'a','b'},{'c','d'}}:replace(symbol)
assert(tostring(mtx:normf()) == "(0+((a):abs())^(2)+((b):abs())^(2)+((c):abs())^(2)+((d):abs())^(2)):sqrt()") assert(tostring(mtx:normf()) == "(0+((a):abs())^(2)+((b):abs())^(2)+((c):abs())^(2)+((d):abs())^(2)):sqrt()")
-- matrix.normmax -- matrix.normmax
-- note: symbolic matrices not supported -- note: symbolic matrices not supported
mtx = matrix{{2,3},{-2,-4}} mtx = matrix{{2,3},{-2,-4}}
assert(mtx:normmax() == 4) assert(mtx:normmax() == 4)
mtx = matrix{{'2i','3'},{'-2i','-4i'}}:tocomplex() mtx = matrix{{'2i','3'},{'-2i','-4i'}}:replace(complex)
assert(mtx:normmax() == 4) assert(mtx:normmax() == 4)
-- matrix.transpose -- matrix.transpose
@ -187,7 +189,7 @@ assert( m1:rotr() == matrix{{6,4,2},{7,5,3}} )
mtx = matrix{{4,2,-3},{3,-5,2}} mtx = matrix{{4,2,-3},{3,-5,2}}
assert(tostring(mtx) == "4\t2\t-3\n3\t-5\t2" ) assert(tostring(mtx) == "4\t2\t-3\n3\t-5\t2" )
-- matrix.tostring; complex -- matrix.tostring; complex
mtx = matrix{{4,"2+i"},{"3-4i",5}}:tocomplex() mtx = matrix{{4,"2+i"},{"3-4i",5}}:replace(complex)
assert(tostring(mtx) == "4\t2+i\n3-4i\t5" ) assert(tostring(mtx) == "4\t2+i\n3-4i\t5" )
-- matrix.tostring; tensor -- matrix.tostring; tensor
local mt = matrix {{{1,2},{3,4}},{{5,6},{7,8}}} local mt = matrix {{{1,2},{3,4}},{{5,6},{7,8}}}
@ -216,14 +218,14 @@ assert( vx:scalar( v2 ) == 0 )
assert( v2:len() == math.sqrt( 3^2+4^2 ) ) assert( v2:len() == math.sqrt( 3^2+4^2 ) )
--// test symbolic --// test symbolic
ms = matrix {{ "a",1 },{2,"b"}}:tosymbol() ms = matrix {{ "a",1 },{2,"b"}}:replace(symbol)
ms2 = matrix {{ "a",2 },{"b",3}}:tosymbol() ms2 = matrix {{ "a",2 },{"b",3}}:replace(symbol)
ms3 = ms2+ms ms3 = ms2+ms
ms3 = ms3:replace( "a",4,"b",2 ) ms3 = ms3:replace( symbol.makereplacer( "a",4,"b",2 ) )
ms3 = ms3:solve() ms3 = ms3:solve()
assert( ms3 == matrix {{8,3},{4,5}} ) assert( ms3 == matrix {{8,3},{4,5}} )
ms4 = ms2*ms ms4 = ms2*ms
ms4 = ms4:replace( "a",4,"b",2 ) ms4 = ms4:replace( symbol.makereplacer( "a",4,"b",2 ) )
ms4 = ms4:solve() ms4 = ms4:solve()
assert( ms4 == matrix {{20,8},{14,8}} ) assert( ms4 == matrix {{20,8},{14,8}} )
@ -255,4 +257,6 @@ end
table.sort( t ) table.sort( t )
for i,v in ipairs( t ) do for i,v in ipairs( t ) do
--print( "matrix."..v ) --print( "matrix."..v )
end end
print 'PASSED'