--[[
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.<function> 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.<functions> 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 //--
--///////////////--