From 67f724300f3db59612a9cb49a4c695b576c43aae Mon Sep 17 00:00:00 2001
From: David Manura <dm.git@math2.org>
Date: Wed, 22 Sep 2010 21:00:52 -0400
Subject: [PATCH] convert from dos to unix-style line endings

---
 complex.lua           |  768 ++++++-------
 complex_changelog.txt |  108 +-
 fit.lua               |  234 ++--
 fit_changelog.txt     |   12 +-
 matrix.lua            | 2532 ++++++++++++++++++++---------------------
 matrix_changelog.txt  |  318 +++---
 test_complex.lua      |  368 +++---
 test_fit.lua          |   42 +-
 test_matrix.lua       |  514 ++++-----
 9 files changed, 2448 insertions(+), 2448 deletions(-)

diff --git a/complex.lua b/complex.lua
index e86a24b..c6dd552 100644
--- a/complex.lua
+++ b/complex.lua
@@ -1,385 +1,385 @@
--- complex 0.3.0
--- Lua 5.1
-
--- 'complex' provides common tasks with complex numbers
-
--- function complex.to( arg ); complex( arg )
--- returns a complex number on success, nil on failure
--- arg := number or { number,number } or ( "(-)<number>" and/or "(+/-)<number>i" )
---		e.g. 5; {2,3}; "2", "2+i", "-2i", "2^2*3+1/3i"
---		note: 'i' is always in the numerator, spaces are not allowed
-
--- a complex number is defined as carthesic complex number
--- complex number := { real_part, imaginary_part }
--- this gives fast access to both parts of the number for calculation
--- 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
-
--- http://luaforge.net/projects/LuaMatrix
--- http://lua-users.org/wiki/ComplexNumbers
-
--- Licensed under the same terms as Lua itself.
-
---/////////////--
---// complex //--
---/////////////--
-
--- link to complex table
-local complex = {}
-
--- link to complex metatable
-local complex_meta = {}
-
--- complex.to( arg )
--- return a complex number on success
--- return nil on failure
-local _retone = function() return 1 end
-local _retminusone = function() return -1 end
-function complex.to( num )
-	-- check for table type
-	if type( num ) == "table" then
-		-- check for a complex number
-		if getmetatable( num ) == complex_meta then
-			return num
-		end
-		local real,imag = tonumber( num[1] ),tonumber( num[2] )
-		if real and imag then
-			return setmetatable( { real,imag }, complex_meta )
-		end
-		return
-	end
-	-- check for number
-	local isnum = tonumber( num )
-	if isnum then
-		return setmetatable( { isnum,0 }, complex_meta )
-	end
-	if type( num ) == "string" then
-		-- check for real and complex
-		-- number chars [%-%+%*%^%d%./Ee]
-		local real,sign,imag = string.match( num, "^([%-%+%*%^%d%./Ee]*%d)([%+%-])([%-%+%*%^%d%./Ee]*)i$" )
-		if real then
-			if string.lower(string.sub(real,1,1)) == "e"
-			or string.lower(string.sub(imag,1,1)) == "e" then
-				return
-			end
-			if imag == "" then
-				if sign == "+" then
-					imag = _retone
-				else
-					imag = _retminusone
-				end
-			elseif sign == "+" then
-				imag = loadstring("return tonumber("..imag..")")
-			else
-				imag = loadstring("return tonumber("..sign..imag..")")
-			end
-			real = loadstring("return tonumber("..real..")")
-			if real and imag then
-				return setmetatable( { real(),imag() }, complex_meta )
-			end
-			return
-		end
-		-- check for complex
-		local imag = string.match( num,"^([%-%+%*%^%d%./Ee]*)i$" )
-		if imag then
-			if imag == "" then
-				return setmetatable( { 0,1 }, complex_meta )
-			elseif imag == "-" then
-				return setmetatable( { 0,-1 }, complex_meta )
-			end
-			if string.lower(string.sub(imag,1,1)) ~= "e" then
-				imag = loadstring("return tonumber("..imag..")")
-				if imag then
-					return setmetatable( { 0,imag() }, complex_meta )
-				end
-			end
-			return
-		end
-		-- should be real
-		local real = string.match( num,"^(%-*[%d%.][%-%+%*%^%d%./Ee]*)$" )
-		if real then
-			real = loadstring( "return tonumber("..real..")" )
-			if real then
-				return setmetatable( { real(),0 }, complex_meta )
-			end
-		end
-	end
-end
-
--- complex( arg )
--- same as complex.to( arg )
--- set __call behaviour of complex
-setmetatable( complex, { __call = function( _,num ) return complex.to( num ) end } )
-
--- complex.new( real, complex )
--- fast function to get a complex number, not invoking any checks
-function complex.new( ... )
-	return setmetatable( { ... }, complex_meta )
-end
-
--- complex.type( arg )
--- is argument of type complex
-function complex.type( arg )
-	if getmetatable( arg ) == complex_meta then
-		return "complex"
-	end
-end
-
--- complex.convpolar( r, phi )
--- convert polar coordinates ( r*e^(i*phi) ) to carthesic complex number
--- r (radius) is a number
--- phi (angle) must be in radians; e.g. [0 - 2pi]
-function complex.convpolar( radius, phi )
-	return setmetatable( { radius * math.cos( phi ), radius * math.sin( phi ) }, complex_meta )
-end
-
--- complex.convpolardeg( r, phi )
--- convert polar coordinates ( r*e^(i*phi) ) to carthesic complex number
--- r (radius) is a number
--- phi must be in degrees; e.g. [0° - 360°]
-function complex.convpolardeg( radius, phi )
-	phi = phi/180 * math.pi
-	return setmetatable( { radius * math.cos( phi ), radius * math.sin( phi ) }, complex_meta )
-end
-
---// complex number functions only
-
--- complex.tostring( cx [, formatstr] )
--- to string or real number
--- takes a complex number and returns its string value or real number value
-function complex.tostring( cx,formatstr )
-	local real,imag = cx[1],cx[2]
-	if formatstr then
-		if imag == 0 then
-			return string.format( formatstr, real )
-		elseif real == 0 then
-			return string.format( formatstr, imag ).."i"
-		elseif imag > 0 then
-			return string.format( formatstr, real ).."+"..string.format( formatstr, imag ).."i"
-		end
-		return string.format( formatstr, real )..string.format( formatstr, imag ).."i"
-	end
-	if imag == 0 then
-		return real
-	elseif real == 0 then
-		return ((imag==1 and "") or (imag==-1 and "-") or imag).."i"
-	elseif imag > 0 then
-		return real.."+"..(imag==1 and "" or imag).."i"
-	end
-	return real..(imag==-1 and "-" or imag).."i"
-end
-
--- complex.print( cx [, formatstr] )
--- print a complex number
-function complex.print( ... )
-	print( complex.tostring( ... ) )
-end
-
--- complex.polar( cx )
--- from complex number to polar coordinates
--- output in radians; [-pi,+pi]
--- returns r (radius), phi (angle)
-function complex.polar( cx )
-	return math.sqrt( cx[1]^2 + cx[2]^2 ), math.atan2( cx[2], cx[1] )
-end
-
--- complex.polardeg( cx )
--- from complex number to polar coordinates
--- output in degrees; [-180°,180°]
--- returns r (radius), phi (angle)
-function complex.polardeg( cx )
-	return math.sqrt( cx[1]^2 + cx[2]^2 ), math.atan2( cx[2], cx[1] ) / math.pi * 180
-end
-
--- complex.mulconjugate( cx )
--- multiply with conjugate, function returning a number
-function complex.mulconjugate( cx )
-	return cx[1]^2 + cx[2]^2
-end
-
--- complex.abs( cx )
--- get the absolute value of a complex number
-function complex.abs( cx )
-	return math.sqrt( cx[1]^2 + cx[2]^2 )
-end
-
--- complex.get( cx )
--- returns real_part, imaginary_part
-function complex.get( cx )
-	return cx[1],cx[2]
-end
-
--- complex.set( cx, real, imag )
--- sets real_part = real and imaginary_part = imag
-function complex.set( cx,real,imag )
-	cx[1],cx[2] = real,imag
-end
-
--- complex.is( cx, real, imag )
--- returns true if, real_part = real and imaginary_part = imag
--- else returns false
-function complex.is( cx,real,imag )
-	if cx[1] == real and cx[2] == imag then
-		return true
-	end
-	return false
-end
-
---// functions returning a new complex number
-
--- complex.copy( cx )
--- copy complex number
-function complex.copy( cx )
-	return setmetatable( { cx[1],cx[2] }, complex_meta )
-end
-
--- complex.add( cx1, cx2 )
--- add two numbers; cx1 + cx2
-function complex.add( cx1,cx2 )
-	return setmetatable( { cx1[1]+cx2[1], cx1[2]+cx2[2] }, complex_meta )
-end
-
--- complex.sub( cx1, cx2 )
--- subtract two numbers; cx1 - cx2
-function complex.sub( cx1,cx2 )
-	return setmetatable( { cx1[1]-cx2[1], cx1[2]-cx2[2] }, complex_meta )
-end
-
--- complex.mul( cx1, cx2 )
--- multiply two numbers; cx1 * cx2
-function complex.mul( cx1,cx2 )
-	return setmetatable( { cx1[1]*cx2[1] - cx1[2]*cx2[2],cx1[1]*cx2[2] + cx1[2]*cx2[1] }, complex_meta )
-end
-
--- complex.mulnum( cx, num )
--- multiply complex with number; cx1 * num
-function complex.mulnum( cx,num )
-	return setmetatable( { cx[1]*num,cx[2]*num }, complex_meta )
-end
-
--- complex.div( cx1, cx2 )
--- divide 2 numbers; cx1 / cx2
-function complex.div( cx1,cx2 )
-	-- get complex value
-	local val = cx2[1]^2 + cx2[2]^2
-	-- multiply cx1 with conjugate complex of cx2 and divide through val
-	return setmetatable( { (cx1[1]*cx2[1]+cx1[2]*cx2[2])/val,(cx1[2]*cx2[1]-cx1[1]*cx2[2])/val }, complex_meta )
-end
-
--- complex.divnum( cx, num )
--- divide through a number
-function complex.divnum( cx,num )
-	return setmetatable( { cx[1]/num,cx[2]/num }, complex_meta )
-end
-
--- complex.pow( cx, num )
--- get the power of a complex number
-function complex.pow( cx,num )
-	if math.floor( num ) == num then
-		if num < 0 then
-			local val = cx[1]^2 + cx[2]^2
-			cx = { cx[1]/val,-cx[2]/val }
-			num = -num
-		end
-		local real,imag = cx[1],cx[2]
-		for i = 2,num do
-			real,imag = real*cx[1] - imag*cx[2],real*cx[2] + imag*cx[1]
-		end
-		return setmetatable( { real,imag }, complex_meta )
-	end
-	-- we calculate the polar complex number now
-	-- since then we have the versatility to calc any potenz of the complex number
-	-- then we convert it back to a carthesic complex number, we loose precision here
-	local length,phi = math.sqrt( cx[1]^2 + cx[2]^2 )^num, math.atan2( cx[2], cx[1] )*num
-	return setmetatable( { length * math.cos( phi ), length * math.sin( phi ) }, complex_meta )
-end
-
--- complex.sqrt( cx )
--- get the first squareroot of a complex number, more accurate than cx^.5
-function complex.sqrt( cx )
-	local len = math.sqrt( cx[1]^2+cx[2]^2 )
-	local sign = (cx[2]<0 and -1) or 1
-	return setmetatable( { math.sqrt((cx[1]+len)/2), sign*math.sqrt((len-cx[1])/2) }, complex_meta )
-end
-
--- complex.ln( cx )
--- natural logarithm of cx
-function complex.ln( cx )
-	return setmetatable( { math.log(math.sqrt( cx[1]^2 + cx[2]^2 )),
-		math.atan2( cx[2], cx[1] ) }, complex_meta )
-end
-
--- complex.exp( cx )
--- exponent of cx (e^cx)
-function complex.exp( cx )
-	local expreal = math.exp(cx[1])
-	return setmetatable( { expreal*math.cos(cx[2]), expreal*math.sin(cx[2]) }, complex_meta )
-end
-
--- complex.conjugate( cx )
--- get conjugate complex of number
-function complex.conjugate( cx )
-	return setmetatable( { cx[1], -cx[2] }, complex_meta )
-end
-
--- complex.round( cx [,idp] )
--- round complex numbers, by default to 0 decimal points
-function complex.round( cx,idp )
-	local mult = 10^( idp or 0 )
-	return setmetatable( { math.floor( cx[1] * mult + 0.5 ) / mult,
-		math.floor( cx[2] * mult + 0.5 ) / mult }, complex_meta )
-end
-
---// metatable functions
-
-complex_meta.__add = function( cx1,cx2 )
-	local cx1,cx2 = complex.to( cx1 ),complex.to( cx2 )
-	return complex.add( cx1,cx2 )
-end
-complex_meta.__sub = function( cx1,cx2 )
-	local cx1,cx2 = complex.to( cx1 ),complex.to( cx2 )
-	return complex.sub( cx1,cx2 )
-end
-complex_meta.__mul = function( cx1,cx2 )
-	local cx1,cx2 = complex.to( cx1 ),complex.to( cx2 )
-	return complex.mul( cx1,cx2 )
-end
-complex_meta.__div = function( cx1,cx2 )
-	local cx1,cx2 = complex.to( cx1 ),complex.to( cx2 )
-	return complex.div( cx1,cx2 )
-end
-complex_meta.__pow = function( cx,num )
-	if num == "*" then
-		return complex.conjugate( cx )
-	end
-	return complex.pow( cx,num )
-end
-complex_meta.__unm = function( cx )
-	return setmetatable( { -cx[1], -cx[2] }, complex_meta )
-end
-complex_meta.__eq = function( cx1,cx2 )
-	if cx1[1] == cx2[1] and cx1[2] == cx2[2] then
-		return true
-	end
-	return false
-end
-complex_meta.__tostring = function( cx )
-	return tostring( complex.tostring( cx ) )
-end
-complex_meta.__concat = function( cx,cx2 )
-	return tostring(cx)..tostring(cx2)
-end
--- cx( cx, formatstr )
-complex_meta.__call = function( ... )
-	print( complex.tostring( ... ) )
-end
-complex_meta.__index = {}
-for k,v in pairs( complex ) do
-	complex_meta.__index[k] = v
-end
-
-return complex
-
---///////////////--
---// chillcode //--
+-- complex 0.3.0
+-- Lua 5.1
+
+-- 'complex' provides common tasks with complex numbers
+
+-- function complex.to( arg ); complex( arg )
+-- returns a complex number on success, nil on failure
+-- arg := number or { number,number } or ( "(-)<number>" and/or "(+/-)<number>i" )
+--		e.g. 5; {2,3}; "2", "2+i", "-2i", "2^2*3+1/3i"
+--		note: 'i' is always in the numerator, spaces are not allowed
+
+-- a complex number is defined as carthesic complex number
+-- complex number := { real_part, imaginary_part }
+-- this gives fast access to both parts of the number for calculation
+-- 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
+
+-- http://luaforge.net/projects/LuaMatrix
+-- http://lua-users.org/wiki/ComplexNumbers
+
+-- Licensed under the same terms as Lua itself.
+
+--/////////////--
+--// complex //--
+--/////////////--
+
+-- link to complex table
+local complex = {}
+
+-- link to complex metatable
+local complex_meta = {}
+
+-- complex.to( arg )
+-- return a complex number on success
+-- return nil on failure
+local _retone = function() return 1 end
+local _retminusone = function() return -1 end
+function complex.to( num )
+	-- check for table type
+	if type( num ) == "table" then
+		-- check for a complex number
+		if getmetatable( num ) == complex_meta then
+			return num
+		end
+		local real,imag = tonumber( num[1] ),tonumber( num[2] )
+		if real and imag then
+			return setmetatable( { real,imag }, complex_meta )
+		end
+		return
+	end
+	-- check for number
+	local isnum = tonumber( num )
+	if isnum then
+		return setmetatable( { isnum,0 }, complex_meta )
+	end
+	if type( num ) == "string" then
+		-- check for real and complex
+		-- number chars [%-%+%*%^%d%./Ee]
+		local real,sign,imag = string.match( num, "^([%-%+%*%^%d%./Ee]*%d)([%+%-])([%-%+%*%^%d%./Ee]*)i$" )
+		if real then
+			if string.lower(string.sub(real,1,1)) == "e"
+			or string.lower(string.sub(imag,1,1)) == "e" then
+				return
+			end
+			if imag == "" then
+				if sign == "+" then
+					imag = _retone
+				else
+					imag = _retminusone
+				end
+			elseif sign == "+" then
+				imag = loadstring("return tonumber("..imag..")")
+			else
+				imag = loadstring("return tonumber("..sign..imag..")")
+			end
+			real = loadstring("return tonumber("..real..")")
+			if real and imag then
+				return setmetatable( { real(),imag() }, complex_meta )
+			end
+			return
+		end
+		-- check for complex
+		local imag = string.match( num,"^([%-%+%*%^%d%./Ee]*)i$" )
+		if imag then
+			if imag == "" then
+				return setmetatable( { 0,1 }, complex_meta )
+			elseif imag == "-" then
+				return setmetatable( { 0,-1 }, complex_meta )
+			end
+			if string.lower(string.sub(imag,1,1)) ~= "e" then
+				imag = loadstring("return tonumber("..imag..")")
+				if imag then
+					return setmetatable( { 0,imag() }, complex_meta )
+				end
+			end
+			return
+		end
+		-- should be real
+		local real = string.match( num,"^(%-*[%d%.][%-%+%*%^%d%./Ee]*)$" )
+		if real then
+			real = loadstring( "return tonumber("..real..")" )
+			if real then
+				return setmetatable( { real(),0 }, complex_meta )
+			end
+		end
+	end
+end
+
+-- complex( arg )
+-- same as complex.to( arg )
+-- set __call behaviour of complex
+setmetatable( complex, { __call = function( _,num ) return complex.to( num ) end } )
+
+-- complex.new( real, complex )
+-- fast function to get a complex number, not invoking any checks
+function complex.new( ... )
+	return setmetatable( { ... }, complex_meta )
+end
+
+-- complex.type( arg )
+-- is argument of type complex
+function complex.type( arg )
+	if getmetatable( arg ) == complex_meta then
+		return "complex"
+	end
+end
+
+-- complex.convpolar( r, phi )
+-- convert polar coordinates ( r*e^(i*phi) ) to carthesic complex number
+-- r (radius) is a number
+-- phi (angle) must be in radians; e.g. [0 - 2pi]
+function complex.convpolar( radius, phi )
+	return setmetatable( { radius * math.cos( phi ), radius * math.sin( phi ) }, complex_meta )
+end
+
+-- complex.convpolardeg( r, phi )
+-- convert polar coordinates ( r*e^(i*phi) ) to carthesic complex number
+-- r (radius) is a number
+-- phi must be in degrees; e.g. [0° - 360°]
+function complex.convpolardeg( radius, phi )
+	phi = phi/180 * math.pi
+	return setmetatable( { radius * math.cos( phi ), radius * math.sin( phi ) }, complex_meta )
+end
+
+--// complex number functions only
+
+-- complex.tostring( cx [, formatstr] )
+-- to string or real number
+-- takes a complex number and returns its string value or real number value
+function complex.tostring( cx,formatstr )
+	local real,imag = cx[1],cx[2]
+	if formatstr then
+		if imag == 0 then
+			return string.format( formatstr, real )
+		elseif real == 0 then
+			return string.format( formatstr, imag ).."i"
+		elseif imag > 0 then
+			return string.format( formatstr, real ).."+"..string.format( formatstr, imag ).."i"
+		end
+		return string.format( formatstr, real )..string.format( formatstr, imag ).."i"
+	end
+	if imag == 0 then
+		return real
+	elseif real == 0 then
+		return ((imag==1 and "") or (imag==-1 and "-") or imag).."i"
+	elseif imag > 0 then
+		return real.."+"..(imag==1 and "" or imag).."i"
+	end
+	return real..(imag==-1 and "-" or imag).."i"
+end
+
+-- complex.print( cx [, formatstr] )
+-- print a complex number
+function complex.print( ... )
+	print( complex.tostring( ... ) )
+end
+
+-- complex.polar( cx )
+-- from complex number to polar coordinates
+-- output in radians; [-pi,+pi]
+-- returns r (radius), phi (angle)
+function complex.polar( cx )
+	return math.sqrt( cx[1]^2 + cx[2]^2 ), math.atan2( cx[2], cx[1] )
+end
+
+-- complex.polardeg( cx )
+-- from complex number to polar coordinates
+-- output in degrees; [-180°,180°]
+-- returns r (radius), phi (angle)
+function complex.polardeg( cx )
+	return math.sqrt( cx[1]^2 + cx[2]^2 ), math.atan2( cx[2], cx[1] ) / math.pi * 180
+end
+
+-- complex.mulconjugate( cx )
+-- multiply with conjugate, function returning a number
+function complex.mulconjugate( cx )
+	return cx[1]^2 + cx[2]^2
+end
+
+-- complex.abs( cx )
+-- get the absolute value of a complex number
+function complex.abs( cx )
+	return math.sqrt( cx[1]^2 + cx[2]^2 )
+end
+
+-- complex.get( cx )
+-- returns real_part, imaginary_part
+function complex.get( cx )
+	return cx[1],cx[2]
+end
+
+-- complex.set( cx, real, imag )
+-- sets real_part = real and imaginary_part = imag
+function complex.set( cx,real,imag )
+	cx[1],cx[2] = real,imag
+end
+
+-- complex.is( cx, real, imag )
+-- returns true if, real_part = real and imaginary_part = imag
+-- else returns false
+function complex.is( cx,real,imag )
+	if cx[1] == real and cx[2] == imag then
+		return true
+	end
+	return false
+end
+
+--// functions returning a new complex number
+
+-- complex.copy( cx )
+-- copy complex number
+function complex.copy( cx )
+	return setmetatable( { cx[1],cx[2] }, complex_meta )
+end
+
+-- complex.add( cx1, cx2 )
+-- add two numbers; cx1 + cx2
+function complex.add( cx1,cx2 )
+	return setmetatable( { cx1[1]+cx2[1], cx1[2]+cx2[2] }, complex_meta )
+end
+
+-- complex.sub( cx1, cx2 )
+-- subtract two numbers; cx1 - cx2
+function complex.sub( cx1,cx2 )
+	return setmetatable( { cx1[1]-cx2[1], cx1[2]-cx2[2] }, complex_meta )
+end
+
+-- complex.mul( cx1, cx2 )
+-- multiply two numbers; cx1 * cx2
+function complex.mul( cx1,cx2 )
+	return setmetatable( { cx1[1]*cx2[1] - cx1[2]*cx2[2],cx1[1]*cx2[2] + cx1[2]*cx2[1] }, complex_meta )
+end
+
+-- complex.mulnum( cx, num )
+-- multiply complex with number; cx1 * num
+function complex.mulnum( cx,num )
+	return setmetatable( { cx[1]*num,cx[2]*num }, complex_meta )
+end
+
+-- complex.div( cx1, cx2 )
+-- divide 2 numbers; cx1 / cx2
+function complex.div( cx1,cx2 )
+	-- get complex value
+	local val = cx2[1]^2 + cx2[2]^2
+	-- multiply cx1 with conjugate complex of cx2 and divide through val
+	return setmetatable( { (cx1[1]*cx2[1]+cx1[2]*cx2[2])/val,(cx1[2]*cx2[1]-cx1[1]*cx2[2])/val }, complex_meta )
+end
+
+-- complex.divnum( cx, num )
+-- divide through a number
+function complex.divnum( cx,num )
+	return setmetatable( { cx[1]/num,cx[2]/num }, complex_meta )
+end
+
+-- complex.pow( cx, num )
+-- get the power of a complex number
+function complex.pow( cx,num )
+	if math.floor( num ) == num then
+		if num < 0 then
+			local val = cx[1]^2 + cx[2]^2
+			cx = { cx[1]/val,-cx[2]/val }
+			num = -num
+		end
+		local real,imag = cx[1],cx[2]
+		for i = 2,num do
+			real,imag = real*cx[1] - imag*cx[2],real*cx[2] + imag*cx[1]
+		end
+		return setmetatable( { real,imag }, complex_meta )
+	end
+	-- we calculate the polar complex number now
+	-- since then we have the versatility to calc any potenz of the complex number
+	-- then we convert it back to a carthesic complex number, we loose precision here
+	local length,phi = math.sqrt( cx[1]^2 + cx[2]^2 )^num, math.atan2( cx[2], cx[1] )*num
+	return setmetatable( { length * math.cos( phi ), length * math.sin( phi ) }, complex_meta )
+end
+
+-- complex.sqrt( cx )
+-- get the first squareroot of a complex number, more accurate than cx^.5
+function complex.sqrt( cx )
+	local len = math.sqrt( cx[1]^2+cx[2]^2 )
+	local sign = (cx[2]<0 and -1) or 1
+	return setmetatable( { math.sqrt((cx[1]+len)/2), sign*math.sqrt((len-cx[1])/2) }, complex_meta )
+end
+
+-- complex.ln( cx )
+-- natural logarithm of cx
+function complex.ln( cx )
+	return setmetatable( { math.log(math.sqrt( cx[1]^2 + cx[2]^2 )),
+		math.atan2( cx[2], cx[1] ) }, complex_meta )
+end
+
+-- complex.exp( cx )
+-- exponent of cx (e^cx)
+function complex.exp( cx )
+	local expreal = math.exp(cx[1])
+	return setmetatable( { expreal*math.cos(cx[2]), expreal*math.sin(cx[2]) }, complex_meta )
+end
+
+-- complex.conjugate( cx )
+-- get conjugate complex of number
+function complex.conjugate( cx )
+	return setmetatable( { cx[1], -cx[2] }, complex_meta )
+end
+
+-- complex.round( cx [,idp] )
+-- round complex numbers, by default to 0 decimal points
+function complex.round( cx,idp )
+	local mult = 10^( idp or 0 )
+	return setmetatable( { math.floor( cx[1] * mult + 0.5 ) / mult,
+		math.floor( cx[2] * mult + 0.5 ) / mult }, complex_meta )
+end
+
+--// metatable functions
+
+complex_meta.__add = function( cx1,cx2 )
+	local cx1,cx2 = complex.to( cx1 ),complex.to( cx2 )
+	return complex.add( cx1,cx2 )
+end
+complex_meta.__sub = function( cx1,cx2 )
+	local cx1,cx2 = complex.to( cx1 ),complex.to( cx2 )
+	return complex.sub( cx1,cx2 )
+end
+complex_meta.__mul = function( cx1,cx2 )
+	local cx1,cx2 = complex.to( cx1 ),complex.to( cx2 )
+	return complex.mul( cx1,cx2 )
+end
+complex_meta.__div = function( cx1,cx2 )
+	local cx1,cx2 = complex.to( cx1 ),complex.to( cx2 )
+	return complex.div( cx1,cx2 )
+end
+complex_meta.__pow = function( cx,num )
+	if num == "*" then
+		return complex.conjugate( cx )
+	end
+	return complex.pow( cx,num )
+end
+complex_meta.__unm = function( cx )
+	return setmetatable( { -cx[1], -cx[2] }, complex_meta )
+end
+complex_meta.__eq = function( cx1,cx2 )
+	if cx1[1] == cx2[1] and cx1[2] == cx2[2] then
+		return true
+	end
+	return false
+end
+complex_meta.__tostring = function( cx )
+	return tostring( complex.tostring( cx ) )
+end
+complex_meta.__concat = function( cx,cx2 )
+	return tostring(cx)..tostring(cx2)
+end
+-- cx( cx, formatstr )
+complex_meta.__call = function( ... )
+	print( complex.tostring( ... ) )
+end
+complex_meta.__index = {}
+for k,v in pairs( complex ) do
+	complex_meta.__index[k] = v
+end
+
+return complex
+
+--///////////////--
+--// chillcode //--
 --///////////////--
\ No newline at end of file
diff --git a/complex_changelog.txt b/complex_changelog.txt
index b0826c8..bd99662 100644
--- a/complex_changelog.txt
+++ b/complex_changelog.txt
@@ -1,55 +1,55 @@
-
-complex changelog:
-
-v 0.3.0: 2007-08-26
-	- fixed print function
-	- fixed call behaviour of complex number
-	- __unm now directly alters the complex number instead of calling a complex function
-	- changed cmplx to imag (imaginary), for better understand of what is the real
-		and what is the imaginary part
-	- changed convpolar, to take as first argument the radius and second the angle
-		now both have to be given but is more conform to the polar function,
-		that returns radius, phi
-	- updated test_complex.lua
-
-v 0.2.9: 2007-08-13
-	- complex.tonumber is now complex.tostring, it will still return a number
-		if it only has a real part, use tostring( cmplx ) to make sure it
-		always returns a string
-
-v 0.2.8: 2007-08-11
-	- metatable functions such as __add;__sub ... didn't return an error if it couldn't be calculated
-		now the complex.add function will return an error
-	- topolar is not polar, and topolardeg is now polardeg
-	- exp returns now e^(complexnumber)
-	- convpolar returns polar coordinates as carthesic complex number
-	- sqrt now calcs the sqrt of a complex number
-	- added changelog file and *.zip
-v 0.2.7:
-	- optimised functions with ...
-	- complex.len is now complex.abs
-v 0.2.6
-	- optimised metatable functions
-	- added __call behaviour, prints the complex number
-v 0.2.5
-	- added option formatstr in complex.tonumber
-	- added __tostring, is the same as complex.tonumber, but it always returns a string
-	- added __concat so when one does "complex: "..cx, will do the right thing
-		doesn't work with table.concat unfortunatly
-v 0.2.4
-	- changed to normal module loading via require rather than dofile
-		e.g. "local complex = require "complex", if its on the call folder
-		very good so one can decide for oneself to have complex global or local
-	- added license
-	- changed function name tocomplex to complex.to
-v 0.2.3
-	- rewritten tocomplex, clear structure, can be used as 'if tocomplex( arg ) then'
-	- added fast function complex.new( real,cmplx ), not invoking any checks
-v 0.2.2
-	- changed tocomplex, to that it only returns a complex table if the input is
-	either a table holding 2 numbers as first elemenst or a number or a string
-	representing a complex number as defined
-	else it returns nil, so one can do "if tocomplex( arg ) then"
-v 0.2.1 added cx:set( real,complex ) to set real and complex
-	real,complex = cx:get() and cx:is(real,complex)
+
+complex changelog:
+
+v 0.3.0: 2007-08-26
+	- fixed print function
+	- fixed call behaviour of complex number
+	- __unm now directly alters the complex number instead of calling a complex function
+	- changed cmplx to imag (imaginary), for better understand of what is the real
+		and what is the imaginary part
+	- changed convpolar, to take as first argument the radius and second the angle
+		now both have to be given but is more conform to the polar function,
+		that returns radius, phi
+	- updated test_complex.lua
+
+v 0.2.9: 2007-08-13
+	- complex.tonumber is now complex.tostring, it will still return a number
+		if it only has a real part, use tostring( cmplx ) to make sure it
+		always returns a string
+
+v 0.2.8: 2007-08-11
+	- metatable functions such as __add;__sub ... didn't return an error if it couldn't be calculated
+		now the complex.add function will return an error
+	- topolar is not polar, and topolardeg is now polardeg
+	- exp returns now e^(complexnumber)
+	- convpolar returns polar coordinates as carthesic complex number
+	- sqrt now calcs the sqrt of a complex number
+	- added changelog file and *.zip
+v 0.2.7:
+	- optimised functions with ...
+	- complex.len is now complex.abs
+v 0.2.6
+	- optimised metatable functions
+	- added __call behaviour, prints the complex number
+v 0.2.5
+	- added option formatstr in complex.tonumber
+	- added __tostring, is the same as complex.tonumber, but it always returns a string
+	- added __concat so when one does "complex: "..cx, will do the right thing
+		doesn't work with table.concat unfortunatly
+v 0.2.4
+	- changed to normal module loading via require rather than dofile
+		e.g. "local complex = require "complex", if its on the call folder
+		very good so one can decide for oneself to have complex global or local
+	- added license
+	- changed function name tocomplex to complex.to
+v 0.2.3
+	- rewritten tocomplex, clear structure, can be used as 'if tocomplex( arg ) then'
+	- added fast function complex.new( real,cmplx ), not invoking any checks
+v 0.2.2
+	- changed tocomplex, to that it only returns a complex table if the input is
+	either a table holding 2 numbers as first elemenst or a number or a string
+	representing a complex number as defined
+	else it returns nil, so one can do "if tocomplex( arg ) then"
+v 0.2.1 added cx:set( real,complex ) to set real and complex
+	real,complex = cx:get() and cx:is(real,complex)
 built
\ No newline at end of file
diff --git a/fit.lua b/fit.lua
index 3a3ff2a..de58a38 100644
--- a/fit.lua
+++ b/fit.lua
@@ -1,118 +1,118 @@
---///////////////////--
---// Curve Fitting //--
---///////////////////--
-
--- v 0.2
-
--- Lua 5.1 compatible
-
--- little add-on to the matrix module, to show some curve fitting
-
--- http://luaforge.net/projects/LuaMatrix
--- http://lua-users.org/wiki/SimpleFit
-
--- Licensed under the same terms as Lua itself.
-
--- requires matrix module
-local matrix = require "matrix"
-
--- The Fit Table
-local fit = {}
-
--- Note all these Algos use the Gauss-Jordan Method to caculate equation systems
-
--- function to get the results
-local function getresults( mtx )
-	assert( #mtx+1 == #mtx[1], "Cannot calculate Results" )
-	mtx:dogauss()
-	-- tresults
-	local cols = #mtx[1]
-	local tres = {}
-	for i = 1,#mtx do
-		tres[i] = mtx[i][cols]
-	end
-	return unpack( tres )
-end
-
--- fit.linear ( x_values, y_values )
--- fit a straight line
--- model (  y = a + b * x  )
--- returns a, b
-function fit.linear( x_values,y_values )
-	-- x_values = { x1,x2,x3,...,xn }
-	-- y_values = { y1,y2,y3,...,yn }
-	
-	-- values for A matrix
-	local a_vals = {}
-	-- values for Y vector
-	local y_vals = {}
-
-	for i,v in ipairs( x_values ) do
-		a_vals[i] = { 1, v }
-		y_vals[i] = { y_values[i] }
-	end
-
-	-- create both Matrixes
-	local A = matrix:new( a_vals )
-	local Y = matrix:new( y_vals )
-
-	local ATA = matrix.mul( matrix.transpose(A), A )
-	local ATY = matrix.mul( matrix.transpose(A), Y )
-
-	local ATAATY = matrix.concath(ATA,ATY)
-
-	return getresults( ATAATY )
-end
-
--- fit.parabola ( x_values, y_values )
--- Fit a parabola
--- model (  y = a + b * x + c * x² )
--- returns a, b, c
-function fit.parabola( x_values,y_values )
-	-- x_values = { x1,x2,x3,...,xn }
-	-- y_values = { y1,y2,y3,...,yn }
-
-	-- values for A matrix
-	local a_vals = {}
-	-- values for Y vector
-	local y_vals = {}
-
-	for i,v in ipairs( x_values ) do
-		a_vals[i] = { 1, v, v*v }
-		y_vals[i] = { y_values[i] }
-	end
-
-	-- create both Matrixes
-	local A = matrix:new( a_vals )
-	local Y = matrix:new( y_vals )
-
-	local ATA = matrix.mul( matrix.transpose(A), A )
-	local ATY = matrix.mul( matrix.transpose(A), Y )
-
-	local ATAATY = matrix.concath(ATA,ATY)
-
-	return getresults( ATAATY )
-end
-
--- fit.exponential ( x_values, y_values )
--- Fit exponential
--- model (  y = a * x^b )
--- returns a, b
-function fit.exponential( x_values,y_values )
-	-- convert to linear problem
-	-- ln(y) = ln(a) + b * ln(x)
-	for i,v in ipairs( x_values ) do
-		x_values[i] = math.log( v )
-		y_values[i] = math.log( y_values[i] )
-	end
-
-	local a,b = fit.linear( x_values,y_values )
-
-	return math.exp(a), b
-end
-
-return fit
-
---///////////////--
---// chillcode //--
+--///////////////////--
+--// Curve Fitting //--
+--///////////////////--
+
+-- v 0.2
+
+-- Lua 5.1 compatible
+
+-- little add-on to the matrix module, to show some curve fitting
+
+-- http://luaforge.net/projects/LuaMatrix
+-- http://lua-users.org/wiki/SimpleFit
+
+-- Licensed under the same terms as Lua itself.
+
+-- requires matrix module
+local matrix = require "matrix"
+
+-- The Fit Table
+local fit = {}
+
+-- Note all these Algos use the Gauss-Jordan Method to caculate equation systems
+
+-- function to get the results
+local function getresults( mtx )
+	assert( #mtx+1 == #mtx[1], "Cannot calculate Results" )
+	mtx:dogauss()
+	-- tresults
+	local cols = #mtx[1]
+	local tres = {}
+	for i = 1,#mtx do
+		tres[i] = mtx[i][cols]
+	end
+	return unpack( tres )
+end
+
+-- fit.linear ( x_values, y_values )
+-- fit a straight line
+-- model (  y = a + b * x  )
+-- returns a, b
+function fit.linear( x_values,y_values )
+	-- x_values = { x1,x2,x3,...,xn }
+	-- y_values = { y1,y2,y3,...,yn }
+	
+	-- values for A matrix
+	local a_vals = {}
+	-- values for Y vector
+	local y_vals = {}
+
+	for i,v in ipairs( x_values ) do
+		a_vals[i] = { 1, v }
+		y_vals[i] = { y_values[i] }
+	end
+
+	-- create both Matrixes
+	local A = matrix:new( a_vals )
+	local Y = matrix:new( y_vals )
+
+	local ATA = matrix.mul( matrix.transpose(A), A )
+	local ATY = matrix.mul( matrix.transpose(A), Y )
+
+	local ATAATY = matrix.concath(ATA,ATY)
+
+	return getresults( ATAATY )
+end
+
+-- fit.parabola ( x_values, y_values )
+-- Fit a parabola
+-- model (  y = a + b * x + c * x² )
+-- returns a, b, c
+function fit.parabola( x_values,y_values )
+	-- x_values = { x1,x2,x3,...,xn }
+	-- y_values = { y1,y2,y3,...,yn }
+
+	-- values for A matrix
+	local a_vals = {}
+	-- values for Y vector
+	local y_vals = {}
+
+	for i,v in ipairs( x_values ) do
+		a_vals[i] = { 1, v, v*v }
+		y_vals[i] = { y_values[i] }
+	end
+
+	-- create both Matrixes
+	local A = matrix:new( a_vals )
+	local Y = matrix:new( y_vals )
+
+	local ATA = matrix.mul( matrix.transpose(A), A )
+	local ATY = matrix.mul( matrix.transpose(A), Y )
+
+	local ATAATY = matrix.concath(ATA,ATY)
+
+	return getresults( ATAATY )
+end
+
+-- fit.exponential ( x_values, y_values )
+-- Fit exponential
+-- model (  y = a * x^b )
+-- returns a, b
+function fit.exponential( x_values,y_values )
+	-- convert to linear problem
+	-- ln(y) = ln(a) + b * ln(x)
+	for i,v in ipairs( x_values ) do
+		x_values[i] = math.log( v )
+		y_values[i] = math.log( y_values[i] )
+	end
+
+	local a,b = fit.linear( x_values,y_values )
+
+	return math.exp(a), b
+end
+
+return fit
+
+--///////////////--
+--// chillcode //--
 --///////////////--
\ No newline at end of file
diff --git a/fit_changelog.txt b/fit_changelog.txt
index 64a8797..69b62af 100644
--- a/fit_changelog.txt
+++ b/fit_changelog.txt
@@ -1,7 +1,7 @@
-fit changelog:
-
-v 0.2: 2007-08-26
-	- some optimising, for the input
-
-v 0.1: 2007-08-11
+fit changelog:
+
+v 0.2: 2007-08-26
+	- some optimising, for the input
+
+v 0.1: 2007-08-11
 	- built from former fit
\ No newline at end of file
diff --git a/matrix.lua b/matrix.lua
index 88f2bd9..8096097 100644
--- a/matrix.lua
+++ b/matrix.lua
@@ -1,1267 +1,1267 @@
---[[
-	matrix v 0.2.8
-	
-	Lua 5.1 compatible
-	
-	'matrix' provides a good selection of matrix functions.
-
-	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 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
-		
-	Sites:
-		http://luaforge.net/projects/LuaMatrix
-		http://lua-users.org/wiki/SimpleMatrix
-		
-	Licensed under the same terms as Lua itself.
-	
-	Developers:
-		Michael Lutz (chillcode)
-		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 //
---////////////
-
-local matrix = {}
-
--- access to the metatable we set at the end of the file
-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 //
---/////////////////////////////
-
---// matrix:new ( rows [, comlumns [, 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, complx 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 2 matrices; m2 may be of bigger size than m1
-function matrix.add( m1, m2 )
-	local mtx = {}
-	for i = 1,#m1 do
-		mtx[i] = {}
-		for j = 1,#m1[1] do
-			mtx[i][j] = m1[i][j] + m2[i][j]
-		end
-	end
-	return setmetatable( mtx, matrix_meta )
-end
-
---// matrix.sub ( m1 ,m2 )
--- Subtract 2 matrices; m2 may be of bigger size than m1
-function matrix.sub( m1, m2 )
-	local mtx = {}
-	for i = 1,#m1 do
-		mtx[i] = {}
-		for j = 1,#m1[1] do
-			mtx[i][j] = m1[i][j] - m2[i][j]
-		end
-	end
-	return setmetatable( mtx, matrix_meta )
-end
-
---// matrix.mul ( m1, m2 )
--- Multiply 2 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 2 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','complex number' or 'string'
--- strings get converted to complex number, if that fails then to symbol
-function matrix.mulnum( m1, num )
-	if type(num) == "string" then
-		num = complex.to(num) or newsymbol(num)
-	end
-	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','complex number' or 'string'
--- strings get converted to complex number, if that fails then to symbol
-function matrix.divnum( m1, num )
-	if type(num) == "string" then
-		num = complex.to(num) or newsymbol(num)
-	end
-	local mtx = {}
-	-- divide elements by 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
-
-
---// 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
-
---// 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
-local fiszerocomplex = function( cx ) return complex.is(cx,0,0) end
-local fiszeronumber = function( num ) return num == 0 end
-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 fiszero, abs
-	if matrix.type( m1 ) == "complex" then
-		fiszero = fiszerocomplex
-		abs = complex.mulconjugate
-	else
-		fiszero = fiszeronumber
-		abs = math.abs
-	end
-	
-	--// 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 not fiszero(e) then
-					-- use element as new subdet
-					subdet,xrow = e,i
-				end
-			-- check for elements nearest to 1 or -1
-			elseif (not fiszero(e)) and math.abs(abs(e)-1) < math.abs(abs(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 not fiszero( mtx[i][j] ) 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 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,
--- what gives better results
-local setelementtosmallest = function( mtx,i,j,fiszero,fisone,abs )
-	-- check if element is one
-	if fisone(mtx[i][j]) then return true end
-	-- check for lowest value
-	local _ilow
-	for _i = i,#mtx do
-		local e = mtx[_i][j]
-		if fisone(e) then
-			break
-		end
-		if not _ilow then
-			if not fiszero(e) then
-				_ilow = _i
-			end
-		elseif (not fiszero(e)) and math.abs(abs(e)-1) < math.abs(abs(mtx[_ilow][j])-1) then
-			_ilow = _i
-		end
-	end
-	if _ilow then
-		-- switch lines if not input line
-		-- legal operation
-		if _ilow ~= i then
-			mtx[i],mtx[_ilow] = mtx[_ilow],mtx[i]
-		end
-		return true
-	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 cxfisone = function( cx ) return complex.abs(cx) == 1 end
-local cxfsetone = function( mtx,i,j ) complex.set(mtx[i][j],1,0) 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,
--- however with big number and matrices we get problems since we do no reducing
-function matrix.dogauss( mtx )
-	local fiszero,fsetzero,fisone,fsetone,abs
-	if matrix.type( mtx ) == "complex" then
-		fiszero = cxfiszero
-		fsetzero = cxfsetzero
-		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]
-	-- stairs left -> right
-	for j = 1,rows do
-		-- check if element can be setted to one
-		if setelementtosmallest( mtx,j,j,fiszero,fisone,abs ) then
-			-- start parsing rows
-			for i = j+1,rows do
-				-- check if element is not already zero
-				if not fiszero(mtx[i][j]) 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] //--
-					fsetzero(mtx,i,j)
-					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 not fiszero(mtx[i][j]) then
-				local factor = mtx[i][j]
-				for _j = j+1,columns do
-					mtx[i][_j] = mtx[i][_j] - factor * mtx[j][_j]
-				end
-				fsetzero(mtx,i,j)
-			end
-		end
-		fsetone(mtx,j,j)
-	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 )
-	if matrix.type( mtx ) == "complex" then
-		for i = 1,#m1 do
-			ident[i] = {}
-			for j = 1,#m1 do
-				if i == j then
-					ident[i][j] = complex.new( 1,0 )
-				else
-					ident[i][j] = complex.new( 0,0 )
-				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
-	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 abs = matrix.type(m1) == "complex" and complex.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 numound 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 )
-	if type(mtx[1][1]) == "table" then
-		if complex.type(mtx[1][1]) then
-			return "complex"
-		end
-		if getmetatable(mtx[1][1]) == symbol_meta then
-			return "symbol"
-		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 2 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 2 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 get_elemnts in string
-local get_tstr = function( t )
-	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 = {}
-	for i,v in ipairs( t ) do
-		tval[i] = string.format( fstr,v )
-	end
-	return "["..table.concat(tval,",").."]"
-end
-local getf_cxstr = function( e,fstr )
-	return complex.tostring( e,fstr )
-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
-
---// matrix.tostring ( mtx, formatstr )
--- tostring function
-function matrix.tostring( mtx, formatstr )
-	local ts = {}
-	local getstr
-	if formatstr then -- get str formatted
-		local mtype = matrix.type( mtx )
-		if mtype == "tensor" then getstr = getf_tstr
-		elseif mtype == "complex" then getstr = getf_cxstr
-		elseif mtype == "symbol" then getstr = getf_symstr
-		else getstr = getf_str end
-		-- iteratr
-		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
-	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 get_tstr or get_str
-	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 'complex' functions //
---////////////////////////////////
-
---// matrix.tocomplex ( mtx )
--- we set now all elements to a complex number
--- also set the metatable
-function matrix.tocomplex( mtx )
-	assert( matrix.type(mtx) == "number", "matrix not of type 'number'" )
-	for i = 1,#mtx do
-		for j = 1,#mtx[1] do
-			mtx[i][j] = complex.to( mtx[i][j] )
-		end
-	end
-	return setmetatable( mtx, matrix_meta )
-end
-
---// matrix.remcomplex ( mtx )
--- set the matrix elements to a number or complex number string
-function matrix.remcomplex( mtx )
-	assert( matrix.type(mtx) == "complex", "matrix not of type 'complex'" )
-	for i = 1,#mtx do
-		for j = 1,#mtx[1] do
-			mtx[i][j] = complex.tostring( mtx[i][j] )
-		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
-
---// 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
-
-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
--- 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 normal,complex and symbolic
-	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
-
--- return the matrix and complex
-return matrix, complex
-
---///////////////--
---// chillcode //--
+--[[
+	matrix v 0.2.8
+	
+	Lua 5.1 compatible
+	
+	'matrix' provides a good selection of matrix functions.
+
+	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 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
+		
+	Sites:
+		http://luaforge.net/projects/LuaMatrix
+		http://lua-users.org/wiki/SimpleMatrix
+		
+	Licensed under the same terms as Lua itself.
+	
+	Developers:
+		Michael Lutz (chillcode)
+		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 //
+--////////////
+
+local matrix = {}
+
+-- access to the metatable we set at the end of the file
+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 //
+--/////////////////////////////
+
+--// matrix:new ( rows [, comlumns [, 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, complx 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 2 matrices; m2 may be of bigger size than m1
+function matrix.add( m1, m2 )
+	local mtx = {}
+	for i = 1,#m1 do
+		mtx[i] = {}
+		for j = 1,#m1[1] do
+			mtx[i][j] = m1[i][j] + m2[i][j]
+		end
+	end
+	return setmetatable( mtx, matrix_meta )
+end
+
+--// matrix.sub ( m1 ,m2 )
+-- Subtract 2 matrices; m2 may be of bigger size than m1
+function matrix.sub( m1, m2 )
+	local mtx = {}
+	for i = 1,#m1 do
+		mtx[i] = {}
+		for j = 1,#m1[1] do
+			mtx[i][j] = m1[i][j] - m2[i][j]
+		end
+	end
+	return setmetatable( mtx, matrix_meta )
+end
+
+--// matrix.mul ( m1, m2 )
+-- Multiply 2 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 2 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','complex number' or 'string'
+-- strings get converted to complex number, if that fails then to symbol
+function matrix.mulnum( m1, num )
+	if type(num) == "string" then
+		num = complex.to(num) or newsymbol(num)
+	end
+	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','complex number' or 'string'
+-- strings get converted to complex number, if that fails then to symbol
+function matrix.divnum( m1, num )
+	if type(num) == "string" then
+		num = complex.to(num) or newsymbol(num)
+	end
+	local mtx = {}
+	-- divide elements by 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
+
+
+--// 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
+
+--// 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
+local fiszerocomplex = function( cx ) return complex.is(cx,0,0) end
+local fiszeronumber = function( num ) return num == 0 end
+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 fiszero, abs
+	if matrix.type( m1 ) == "complex" then
+		fiszero = fiszerocomplex
+		abs = complex.mulconjugate
+	else
+		fiszero = fiszeronumber
+		abs = math.abs
+	end
+	
+	--// 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 not fiszero(e) then
+					-- use element as new subdet
+					subdet,xrow = e,i
+				end
+			-- check for elements nearest to 1 or -1
+			elseif (not fiszero(e)) and math.abs(abs(e)-1) < math.abs(abs(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 not fiszero( mtx[i][j] ) 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 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,
+-- what gives better results
+local setelementtosmallest = function( mtx,i,j,fiszero,fisone,abs )
+	-- check if element is one
+	if fisone(mtx[i][j]) then return true end
+	-- check for lowest value
+	local _ilow
+	for _i = i,#mtx do
+		local e = mtx[_i][j]
+		if fisone(e) then
+			break
+		end
+		if not _ilow then
+			if not fiszero(e) then
+				_ilow = _i
+			end
+		elseif (not fiszero(e)) and math.abs(abs(e)-1) < math.abs(abs(mtx[_ilow][j])-1) then
+			_ilow = _i
+		end
+	end
+	if _ilow then
+		-- switch lines if not input line
+		-- legal operation
+		if _ilow ~= i then
+			mtx[i],mtx[_ilow] = mtx[_ilow],mtx[i]
+		end
+		return true
+	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 cxfisone = function( cx ) return complex.abs(cx) == 1 end
+local cxfsetone = function( mtx,i,j ) complex.set(mtx[i][j],1,0) 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,
+-- however with big number and matrices we get problems since we do no reducing
+function matrix.dogauss( mtx )
+	local fiszero,fsetzero,fisone,fsetone,abs
+	if matrix.type( mtx ) == "complex" then
+		fiszero = cxfiszero
+		fsetzero = cxfsetzero
+		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]
+	-- stairs left -> right
+	for j = 1,rows do
+		-- check if element can be setted to one
+		if setelementtosmallest( mtx,j,j,fiszero,fisone,abs ) then
+			-- start parsing rows
+			for i = j+1,rows do
+				-- check if element is not already zero
+				if not fiszero(mtx[i][j]) 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] //--
+					fsetzero(mtx,i,j)
+					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 not fiszero(mtx[i][j]) then
+				local factor = mtx[i][j]
+				for _j = j+1,columns do
+					mtx[i][_j] = mtx[i][_j] - factor * mtx[j][_j]
+				end
+				fsetzero(mtx,i,j)
+			end
+		end
+		fsetone(mtx,j,j)
+	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 )
+	if matrix.type( mtx ) == "complex" then
+		for i = 1,#m1 do
+			ident[i] = {}
+			for j = 1,#m1 do
+				if i == j then
+					ident[i][j] = complex.new( 1,0 )
+				else
+					ident[i][j] = complex.new( 0,0 )
+				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
+	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 abs = matrix.type(m1) == "complex" and complex.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 numound 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 )
+	if type(mtx[1][1]) == "table" then
+		if complex.type(mtx[1][1]) then
+			return "complex"
+		end
+		if getmetatable(mtx[1][1]) == symbol_meta then
+			return "symbol"
+		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 2 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 2 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 get_elemnts in string
+local get_tstr = function( t )
+	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 = {}
+	for i,v in ipairs( t ) do
+		tval[i] = string.format( fstr,v )
+	end
+	return "["..table.concat(tval,",").."]"
+end
+local getf_cxstr = function( e,fstr )
+	return complex.tostring( e,fstr )
+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
+
+--// matrix.tostring ( mtx, formatstr )
+-- tostring function
+function matrix.tostring( mtx, formatstr )
+	local ts = {}
+	local getstr
+	if formatstr then -- get str formatted
+		local mtype = matrix.type( mtx )
+		if mtype == "tensor" then getstr = getf_tstr
+		elseif mtype == "complex" then getstr = getf_cxstr
+		elseif mtype == "symbol" then getstr = getf_symstr
+		else getstr = getf_str end
+		-- iteratr
+		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
+	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 get_tstr or get_str
+	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 'complex' functions //
+--////////////////////////////////
+
+--// matrix.tocomplex ( mtx )
+-- we set now all elements to a complex number
+-- also set the metatable
+function matrix.tocomplex( mtx )
+	assert( matrix.type(mtx) == "number", "matrix not of type 'number'" )
+	for i = 1,#mtx do
+		for j = 1,#mtx[1] do
+			mtx[i][j] = complex.to( mtx[i][j] )
+		end
+	end
+	return setmetatable( mtx, matrix_meta )
+end
+
+--// matrix.remcomplex ( mtx )
+-- set the matrix elements to a number or complex number string
+function matrix.remcomplex( mtx )
+	assert( matrix.type(mtx) == "complex", "matrix not of type 'complex'" )
+	for i = 1,#mtx do
+		for j = 1,#mtx[1] do
+			mtx[i][j] = complex.tostring( mtx[i][j] )
+		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
+
+--// 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
+
+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
+-- 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 normal,complex and symbolic
+	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
+
+-- return the matrix and complex
+return matrix, complex
+
+--///////////////--
+--// chillcode //--
 --///////////////--
\ No newline at end of file
diff --git a/matrix_changelog.txt b/matrix_changelog.txt
index 24e878f..fc9ee31 100644
--- a/matrix_changelog.txt
+++ b/matrix_changelog.txt
@@ -1,160 +1,160 @@
-matrix function list:
-
-matrix.add
-matrix.columns
-matrix.concath
-matrix.concatv
-matrix.conjugate
-matrix.copy
-matrix.cross
-matrix.det
-matrix.div
-matrix.divnum
-matrix.dogauss
-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.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
-	[ Michael Lutz ]
-	- fixed rotr and rotl for rotating type 'complex','symbol' and 'tensor'
-	- mulscalar is now mulnum, divscalar is now divnum, mul/div a complex number or a string,
-		strings first get checked if they can be converted to a complex number (what returns a
-		complex matrix) and if that fails they get converted to a symbol (what returns a symbol matrix)
-	- require "matrix" will return 'matrix, complex' just for convinience
-	- matrix.size will now returns the correct size of tensors
-	- function matrix.div returns rank on failed invertion of m2
-	- function matrix.det, was adjusted to better combine complex and number matrices
-		and made sure it finds the element nearest to 1 or -1
-	- function dogauss, was updated to handle complex and number matrices in one
-	- tweaked some utility functions, should speed up number matrices at least
-	- updated matrix.latex to support all types of matrices
-	- matrix.tostring can now also format the output, and was updated to handle all types better
-	- matrix.print now just calls matrix.tostring
-	- updated test_matrix.lua
-	- added fit (curve fitting to LuaMatrix package)
-	[ David Manura ]
-	- tweaked matrix.sqrt and matrix.root function;
-		replaced "dist1 > dist or dist1 == dist" with "dist1 >= dist"
-	- tweaked get_abs_avg function
-	- added function matrix.normf ( mtx ), returns the norm abs of a matrix
-	- added function matrix.normmax ( mtx ), returns the biggest abs(element)
-	- added __pow to symbol
-	- added abs() and sqrt() to symbol
-	- fixed some global variables, that were allocated
-
-v 0.2.7: 2007-08-19
-	- added __div to metatable and the corresponding matrix functions( matrix.div(m1,m2); matrix.divscalar(m1,num) )
-	- updated square root function, now returns to the matrix the average error of the calculated to the original matrix
-	- added matrix.root function (from David Manura/http://www.dm.unipi.it/~cortona04/slides/bruno.pdf)
-		to calculate any root of a matrix, returns same values as matrix.sqrt
-	- added function matrix.rotl and matrix.rotr, for rotate left and rotate right
-
-v 0.2.6: 2007-08-12
-	- added patch#5 from DavidManura, fixes symbolic matrices handling
-	- added sqrt function to function list, thx David for the hint, that some matrices don't convergent
-	- added print function for rank3 tensors
-	- added solve function, for symbolic matrices, tries to solve a symbolic matrix so that is numeric again
-	- tocomplex converts to a complex matrix, remcomplex removes the complex metatable and returns its
-		string/number value
-
-v 0.2.5: 2007-08-11
-	- added path#4 from DavidManura
-		It contains some doc updates, fixed handling mtx^-1 for singular matrices,
-		and a few checks that can be useful.
-	- added setting up of a vector as matrix{1,2,3}
-
-v 0.2.4
-	- added patch#3 from DavidManura, fixes negative exponent
-	- removed concat, added concath (concat horizontal) and concatv (concat vertical)
-	- removed get from all commands, all commands return a matrix even those only changing
-		the input matrix itself, except getelement and setelement
-	- submatrix is now subm, and dogaussjordan is now dogauss
-	- __tostring returns matrix.tostring, matrix.tostring returns a simple string with the matrix elements
-	- __call; e.g. mtx(arg), will return matrix.print(mtx,arg)
-	- added function matrix.cross( v1,v2 ) to get the cross product of 2 matrices e.g. m1[3][1], m2[3][1]
-		to get the scalar product of these you have to do m1:scalar( m2 );
-		to create a vector do matrix{{ 1,2,3 }}^'T'
-
-v 0.2.3
-	- some optimising, in parts where complex and normal matrix are separate anyways
-	- added patch from DavidManura; matrix^0 returning the identitiy matrix, makes a lot sense
-		print and printf merged, scalar multiply fixed, and __unm and __tostring added, 'good to know'
-
-v 0.2.2
-	updated matrix function to suit new complex functions, function names
-		should become more stable now
-	matrix can now be loaded via local matrix = require "matrix"
-		or local matrix = dofile( "matrixfile" ), but require is better and faster
-		when loading multiple times
-
-v 0.2.1
-	updated matrix functions to suit new complex functions
-
-v 0.2.0
-	changed matrix functions to use the updated
-		complex number functions; makes a lot less code
-		now complex and normal matrices share the same functions
-		one can add/sub/mul complex and normal matrices, returning a complex matrix
-
-v 0.1.6
-	optimised functions a bit
-	changed matrix.getinv to matrix.invert
-	matrix:new is now matrix.new, 'matrix' should only provide functions
-		added mtx = matrix( ... ); as matrix.new( ... )
-	matrix functions only return tables with metatable of first argument matrix
-	changed complex functions to
-		matrix.c<func_name>; e.g. matirx.cadd; matrix.csub
-	added: function matrix.cconjugate, returns the conjuagte complex matrix
-
-v 0.1.5
-	written complex functions from normal matrices
-	added metatable handling for overloading operators
-
-v 0.1.4
-	added complex add on, to be able to handle
-	complex elements in the format declaired in 'complex'
-
-v 0.1.3
-	added many functions also getdet and dogaussjordan
-		only defined one way to get the determinant, so we are
-		slow here in certain cases where matrices have a triangle shape for exsample
-
-v 0.1.2
-	object returns mainly internal errors for lighter code,
-	structure should be no checks
-
-v 0.1.0
-	matrix
-	structure of matrix m[i][j]
+matrix function list:
+
+matrix.add
+matrix.columns
+matrix.concath
+matrix.concatv
+matrix.conjugate
+matrix.copy
+matrix.cross
+matrix.det
+matrix.div
+matrix.divnum
+matrix.dogauss
+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.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
+	[ Michael Lutz ]
+	- fixed rotr and rotl for rotating type 'complex','symbol' and 'tensor'
+	- mulscalar is now mulnum, divscalar is now divnum, mul/div a complex number or a string,
+		strings first get checked if they can be converted to a complex number (what returns a
+		complex matrix) and if that fails they get converted to a symbol (what returns a symbol matrix)
+	- require "matrix" will return 'matrix, complex' just for convinience
+	- matrix.size will now returns the correct size of tensors
+	- function matrix.div returns rank on failed invertion of m2
+	- function matrix.det, was adjusted to better combine complex and number matrices
+		and made sure it finds the element nearest to 1 or -1
+	- function dogauss, was updated to handle complex and number matrices in one
+	- tweaked some utility functions, should speed up number matrices at least
+	- updated matrix.latex to support all types of matrices
+	- matrix.tostring can now also format the output, and was updated to handle all types better
+	- matrix.print now just calls matrix.tostring
+	- updated test_matrix.lua
+	- added fit (curve fitting to LuaMatrix package)
+	[ David Manura ]
+	- tweaked matrix.sqrt and matrix.root function;
+		replaced "dist1 > dist or dist1 == dist" with "dist1 >= dist"
+	- tweaked get_abs_avg function
+	- added function matrix.normf ( mtx ), returns the norm abs of a matrix
+	- added function matrix.normmax ( mtx ), returns the biggest abs(element)
+	- added __pow to symbol
+	- added abs() and sqrt() to symbol
+	- fixed some global variables, that were allocated
+
+v 0.2.7: 2007-08-19
+	- added __div to metatable and the corresponding matrix functions( matrix.div(m1,m2); matrix.divscalar(m1,num) )
+	- updated square root function, now returns to the matrix the average error of the calculated to the original matrix
+	- added matrix.root function (from David Manura/http://www.dm.unipi.it/~cortona04/slides/bruno.pdf)
+		to calculate any root of a matrix, returns same values as matrix.sqrt
+	- added function matrix.rotl and matrix.rotr, for rotate left and rotate right
+
+v 0.2.6: 2007-08-12
+	- added patch#5 from DavidManura, fixes symbolic matrices handling
+	- added sqrt function to function list, thx David for the hint, that some matrices don't convergent
+	- added print function for rank3 tensors
+	- added solve function, for symbolic matrices, tries to solve a symbolic matrix so that is numeric again
+	- tocomplex converts to a complex matrix, remcomplex removes the complex metatable and returns its
+		string/number value
+
+v 0.2.5: 2007-08-11
+	- added path#4 from DavidManura
+		It contains some doc updates, fixed handling mtx^-1 for singular matrices,
+		and a few checks that can be useful.
+	- added setting up of a vector as matrix{1,2,3}
+
+v 0.2.4
+	- added patch#3 from DavidManura, fixes negative exponent
+	- removed concat, added concath (concat horizontal) and concatv (concat vertical)
+	- removed get from all commands, all commands return a matrix even those only changing
+		the input matrix itself, except getelement and setelement
+	- submatrix is now subm, and dogaussjordan is now dogauss
+	- __tostring returns matrix.tostring, matrix.tostring returns a simple string with the matrix elements
+	- __call; e.g. mtx(arg), will return matrix.print(mtx,arg)
+	- added function matrix.cross( v1,v2 ) to get the cross product of 2 matrices e.g. m1[3][1], m2[3][1]
+		to get the scalar product of these you have to do m1:scalar( m2 );
+		to create a vector do matrix{{ 1,2,3 }}^'T'
+
+v 0.2.3
+	- some optimising, in parts where complex and normal matrix are separate anyways
+	- added patch from DavidManura; matrix^0 returning the identitiy matrix, makes a lot sense
+		print and printf merged, scalar multiply fixed, and __unm and __tostring added, 'good to know'
+
+v 0.2.2
+	updated matrix function to suit new complex functions, function names
+		should become more stable now
+	matrix can now be loaded via local matrix = require "matrix"
+		or local matrix = dofile( "matrixfile" ), but require is better and faster
+		when loading multiple times
+
+v 0.2.1
+	updated matrix functions to suit new complex functions
+
+v 0.2.0
+	changed matrix functions to use the updated
+		complex number functions; makes a lot less code
+		now complex and normal matrices share the same functions
+		one can add/sub/mul complex and normal matrices, returning a complex matrix
+
+v 0.1.6
+	optimised functions a bit
+	changed matrix.getinv to matrix.invert
+	matrix:new is now matrix.new, 'matrix' should only provide functions
+		added mtx = matrix( ... ); as matrix.new( ... )
+	matrix functions only return tables with metatable of first argument matrix
+	changed complex functions to
+		matrix.c<func_name>; e.g. matirx.cadd; matrix.csub
+	added: function matrix.cconjugate, returns the conjuagte complex matrix
+
+v 0.1.5
+	written complex functions from normal matrices
+	added metatable handling for overloading operators
+
+v 0.1.4
+	added complex add on, to be able to handle
+	complex elements in the format declaired in 'complex'
+
+v 0.1.3
+	added many functions also getdet and dogaussjordan
+		only defined one way to get the determinant, so we are
+		slow here in certain cases where matrices have a triangle shape for exsample
+
+v 0.1.2
+	object returns mainly internal errors for lighter code,
+	structure should be no checks
+
+v 0.1.0
+	matrix
+	structure of matrix m[i][j]
 	added simple operations +-..
\ No newline at end of file
diff --git a/test_complex.lua b/test_complex.lua
index 63e6d9f..0b1aa5c 100644
--- a/test_complex.lua
+++ b/test_complex.lua
@@ -1,185 +1,185 @@
-local complex = require "complex"
-
-local cx,cx1,cx2,re,im
-
--- complex.to / complex call
-cx = complex { 2,3 }
-assert( tostring( cx ) == "2+3i" )
-cx = complex ( 2 )
-assert( tostring( cx ) == "2" )
-assert( cx:tostring() == 2 )
-cx = complex "2^2+3/2i"
-assert( tostring( cx ) == "4+1.5i" )
-cx = complex ".5-2E-3i"
-assert( tostring( cx ) == "0.5-0.002i" )
-cx = complex "3i"
-assert( tostring( cx ) == "3i" )
-cx = complex "2"
-assert( tostring( cx ) == "2" )
-assert( cx:tostring() == 2 )
-assert( complex "2 + 4i" == nil )
-
--- complex.new
-cx = complex.new( 2,3 )
-assert( tostring( cx ) == "2+3i" )
-
--- complex.type
-assert( complex.type( cx ) == "complex" )
-assert( complex.type( {} ) == nil )
-
--- complex.convpolar( radius, phi )
-assert( complex.convpolar( 3, 0 ):round(10) == complex "3" )
-assert( complex.convpolar( 3, math.pi/2 ):round(10) == complex "3i" )
-assert( complex.convpolar( 3, math.pi ):round(10) == complex "-3" )
-assert( complex.convpolar( 3, math.pi*3/2 ):round(10) == complex "-3i" )
-assert( complex.convpolar( 3, math.pi*2 ):round(10) == complex "3" )
-
--- complex.convpolardeg( radius, phi )
-assert( complex.convpolardeg( 3, 0 ):round(10) == complex "3" )
-assert( complex.convpolardeg( 3, 90 ):round(10) == complex "3i" )
-assert( complex.convpolardeg( 3, 180 ):round(10) == complex "-3" )
-assert( complex.convpolardeg( 3, 270 ):round(10) == complex "-3i" )
-assert( complex.convpolardeg( 3, 360 ):round(10) == complex "3" )
-
--- complex.tostring( cx,formatstr )
-cx = complex "2+3i"
-assert( complex.tostring( cx ) == "2+3i" )
-assert( complex.tostring( cx, "%.2f" ) == "2.00+3.00i" )
--- does not support a second argument
-assert( tostring( cx, "%.2f" ) == "2+3i" )
-
--- complex.polar( cx )
-local r,phi = complex.polar( {3,0} )
-assert( r == 3 )
-assert( phi == 0 )
-
-local r,phi = complex.polar( {0,3} )
-assert( r == 3 )
-assert( phi == math.pi/2 )
-
-local r,phi = complex.polar( {-3,0} )
-assert( r == 3 )
-assert( phi == math.pi )
-
-local r,phi = complex.polar( {0,-3} )
-assert( r == 3 )
-assert( phi == -math.pi/2 )
-
--- complex.polardeg( cx )
-local r,phi = complex.polardeg( {3,0} )
-assert( r == 3 )
-assert( phi == 0 )
-
-local r,phi = complex.polardeg( {0,3} )
-assert( r == 3 )
-assert( phi == 90 )
-
-local r,phi = complex.polardeg( {-3,0} )
-assert( r == 3 )
-assert( phi == 180 )
-
-local r,phi = complex.polardeg( {0,-3} )
-assert( r == 3 )
-assert( phi == -90 )
-
--- complex.mulconjugate( cx )
-cx = complex "2+3i"
-assert( complex.mulconjugate( cx ) == 13 )
-
--- complex.abs( cx )
-cx = complex "3+4i"
-assert( complex.abs( cx ) == 5 )
-
--- complex.get( cx )
-cx = complex "2+3i"
-re,im = complex.get( cx )
-assert( re == 2 )
-assert( im == 3 )
-
--- complex.set( cx, re, im )
-cx = complex "2+3i"
-complex.set( cx, 3, 2 )
-assert( cx == complex "3+2i" )
-
--- complex.is( cx, re, im )
-cx = complex "2+3i"
-assert( complex.is( cx, 2, 3 ) == true )
-assert( complex.is( cx, 3, 2 ) == false )
-
--- complex.copy( cx )
-cx = complex "2+3i"
-cx1 = complex.copy( cx )
-complex.set( cx, 1, 1 )
-assert( cx1 == complex "2+3i" )
-
--- complex.add( cx1, cx2 )
-cx1 = complex "2+3i"
-cx2 = complex "3+2i"
-assert( complex.add(cx1,cx2) == complex "5+5i" )
-
--- complex.sub( cx1, cx2 )
-cx1 = complex "2+3i"
-cx2 = complex "3+2i"
-assert( complex.sub(cx1,cx2) == complex "-1+1i" )
-
--- complex.mul( cx1, cx2 )
-cx1 = complex "2+3i"
-cx2 = complex "3+2i"
-assert( complex.mul(cx1,cx2) == complex "13i" )
-
--- complex.mulnum( cx, num )
-cx = complex "2+3i"
-assert( complex.mulnum( cx, 2 ) == complex "4+6i" )
-
--- complex.div( cx1, cx2 )
-cx1 = complex "2+3i"
-cx2 = complex "3-2i"
-assert( complex.div(cx1,cx2) == complex "i" )
-
--- complex.divnum( cx, num )
-cx = complex "2+3i"
-assert( complex.divnum( cx, 2 ) == complex "1+1.5i" )
-
--- complex.pow( cx, num )
-cx = complex "2+3i"
-assert( complex.pow( cx, 3 ) == complex "-46+9i" )
-
-cx = complex( -121 )
-cx = cx^.5
--- we have to round here due to the polar calculation of the squareroot
-cx = cx:round( 10 )
-assert( cx == complex "11i" )
-
-cx = complex"2+3i"
-assert( cx^-2 ~= 1/cx^2 )
-assert( cx^-2 == (cx^-1)^2 )
-assert( tostring( cx^-2 ) == tostring( 1/cx^2 ) )
-
--- complex.sqrt( cx )
-cx = complex( -121 )
-assert( complex.sqrt( cx ) == complex "11i" )
-cx = complex "2-3i"
-cx = cx^2
-assert( cx:sqrt() == complex "2-3i" )
-
--- complex.ln( cx )
-cx = complex "3+4i"
-assert( cx:ln():round( 4 ) == complex "1.6094+0.9273i" )
-
--- complex.exp( cx )
-cx = complex "2+3i"
-assert( cx:ln():exp() == complex "2+3i" )
-
--- complex.conjugate( cx )
-cx = complex "2+3i"
-assert( cx:conjugate() == complex "2-3i" )
-
--- metatable
-
--- __add
-cx = complex "2+3i"
-assert( cx+2 == complex "4+3i" )
-
--- __unm
-cx = complex "2+3i"
+local complex = require "complex"
+
+local cx,cx1,cx2,re,im
+
+-- complex.to / complex call
+cx = complex { 2,3 }
+assert( tostring( cx ) == "2+3i" )
+cx = complex ( 2 )
+assert( tostring( cx ) == "2" )
+assert( cx:tostring() == 2 )
+cx = complex "2^2+3/2i"
+assert( tostring( cx ) == "4+1.5i" )
+cx = complex ".5-2E-3i"
+assert( tostring( cx ) == "0.5-0.002i" )
+cx = complex "3i"
+assert( tostring( cx ) == "3i" )
+cx = complex "2"
+assert( tostring( cx ) == "2" )
+assert( cx:tostring() == 2 )
+assert( complex "2 + 4i" == nil )
+
+-- complex.new
+cx = complex.new( 2,3 )
+assert( tostring( cx ) == "2+3i" )
+
+-- complex.type
+assert( complex.type( cx ) == "complex" )
+assert( complex.type( {} ) == nil )
+
+-- complex.convpolar( radius, phi )
+assert( complex.convpolar( 3, 0 ):round(10) == complex "3" )
+assert( complex.convpolar( 3, math.pi/2 ):round(10) == complex "3i" )
+assert( complex.convpolar( 3, math.pi ):round(10) == complex "-3" )
+assert( complex.convpolar( 3, math.pi*3/2 ):round(10) == complex "-3i" )
+assert( complex.convpolar( 3, math.pi*2 ):round(10) == complex "3" )
+
+-- complex.convpolardeg( radius, phi )
+assert( complex.convpolardeg( 3, 0 ):round(10) == complex "3" )
+assert( complex.convpolardeg( 3, 90 ):round(10) == complex "3i" )
+assert( complex.convpolardeg( 3, 180 ):round(10) == complex "-3" )
+assert( complex.convpolardeg( 3, 270 ):round(10) == complex "-3i" )
+assert( complex.convpolardeg( 3, 360 ):round(10) == complex "3" )
+
+-- complex.tostring( cx,formatstr )
+cx = complex "2+3i"
+assert( complex.tostring( cx ) == "2+3i" )
+assert( complex.tostring( cx, "%.2f" ) == "2.00+3.00i" )
+-- does not support a second argument
+assert( tostring( cx, "%.2f" ) == "2+3i" )
+
+-- complex.polar( cx )
+local r,phi = complex.polar( {3,0} )
+assert( r == 3 )
+assert( phi == 0 )
+
+local r,phi = complex.polar( {0,3} )
+assert( r == 3 )
+assert( phi == math.pi/2 )
+
+local r,phi = complex.polar( {-3,0} )
+assert( r == 3 )
+assert( phi == math.pi )
+
+local r,phi = complex.polar( {0,-3} )
+assert( r == 3 )
+assert( phi == -math.pi/2 )
+
+-- complex.polardeg( cx )
+local r,phi = complex.polardeg( {3,0} )
+assert( r == 3 )
+assert( phi == 0 )
+
+local r,phi = complex.polardeg( {0,3} )
+assert( r == 3 )
+assert( phi == 90 )
+
+local r,phi = complex.polardeg( {-3,0} )
+assert( r == 3 )
+assert( phi == 180 )
+
+local r,phi = complex.polardeg( {0,-3} )
+assert( r == 3 )
+assert( phi == -90 )
+
+-- complex.mulconjugate( cx )
+cx = complex "2+3i"
+assert( complex.mulconjugate( cx ) == 13 )
+
+-- complex.abs( cx )
+cx = complex "3+4i"
+assert( complex.abs( cx ) == 5 )
+
+-- complex.get( cx )
+cx = complex "2+3i"
+re,im = complex.get( cx )
+assert( re == 2 )
+assert( im == 3 )
+
+-- complex.set( cx, re, im )
+cx = complex "2+3i"
+complex.set( cx, 3, 2 )
+assert( cx == complex "3+2i" )
+
+-- complex.is( cx, re, im )
+cx = complex "2+3i"
+assert( complex.is( cx, 2, 3 ) == true )
+assert( complex.is( cx, 3, 2 ) == false )
+
+-- complex.copy( cx )
+cx = complex "2+3i"
+cx1 = complex.copy( cx )
+complex.set( cx, 1, 1 )
+assert( cx1 == complex "2+3i" )
+
+-- complex.add( cx1, cx2 )
+cx1 = complex "2+3i"
+cx2 = complex "3+2i"
+assert( complex.add(cx1,cx2) == complex "5+5i" )
+
+-- complex.sub( cx1, cx2 )
+cx1 = complex "2+3i"
+cx2 = complex "3+2i"
+assert( complex.sub(cx1,cx2) == complex "-1+1i" )
+
+-- complex.mul( cx1, cx2 )
+cx1 = complex "2+3i"
+cx2 = complex "3+2i"
+assert( complex.mul(cx1,cx2) == complex "13i" )
+
+-- complex.mulnum( cx, num )
+cx = complex "2+3i"
+assert( complex.mulnum( cx, 2 ) == complex "4+6i" )
+
+-- complex.div( cx1, cx2 )
+cx1 = complex "2+3i"
+cx2 = complex "3-2i"
+assert( complex.div(cx1,cx2) == complex "i" )
+
+-- complex.divnum( cx, num )
+cx = complex "2+3i"
+assert( complex.divnum( cx, 2 ) == complex "1+1.5i" )
+
+-- complex.pow( cx, num )
+cx = complex "2+3i"
+assert( complex.pow( cx, 3 ) == complex "-46+9i" )
+
+cx = complex( -121 )
+cx = cx^.5
+-- we have to round here due to the polar calculation of the squareroot
+cx = cx:round( 10 )
+assert( cx == complex "11i" )
+
+cx = complex"2+3i"
+assert( cx^-2 ~= 1/cx^2 )
+assert( cx^-2 == (cx^-1)^2 )
+assert( tostring( cx^-2 ) == tostring( 1/cx^2 ) )
+
+-- complex.sqrt( cx )
+cx = complex( -121 )
+assert( complex.sqrt( cx ) == complex "11i" )
+cx = complex "2-3i"
+cx = cx^2
+assert( cx:sqrt() == complex "2-3i" )
+
+-- complex.ln( cx )
+cx = complex "3+4i"
+assert( cx:ln():round( 4 ) == complex "1.6094+0.9273i" )
+
+-- complex.exp( cx )
+cx = complex "2+3i"
+assert( cx:ln():exp() == complex "2+3i" )
+
+-- complex.conjugate( cx )
+cx = complex "2+3i"
+assert( cx:conjugate() == complex "2-3i" )
+
+-- metatable
+
+-- __add
+cx = complex "2+3i"
+assert( cx+2 == complex "4+3i" )
+
+-- __unm
+cx = complex "2+3i"
 assert( -cx == complex "-2-3i" )
\ No newline at end of file
diff --git a/test_fit.lua b/test_fit.lua
index 3ba243a..80ccb84 100644
--- a/test_fit.lua
+++ b/test_fit.lua
@@ -1,22 +1,22 @@
--- 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.." )")
+-- 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.." )")
 	
\ No newline at end of file
diff --git a/test_matrix.lua b/test_matrix.lua
index 2fd597c..1302bf2 100644
--- a/test_matrix.lua
+++ b/test_matrix.lua
@@ -1,258 +1,258 @@
-local matrix, complex = require "matrix"
-
-local mtx, m1,m2,m3,m4,m5, ms,ms1,ms2,ms3,ms4
-
--- test matrix:new/matrix call function
--- normal matrix
-mtx = matrix {{1,2},{3,4}}
-assert( tostring(mtx) == "1\t2\n3\t4" )
--- vector
-mtx = matrix {1,2,3}
-assert( tostring(mtx) == "1\n2\n3" )
--- matrix with size 2x2
-mtx = matrix (2,2)
-assert( tostring(mtx) == "0\t0\n0\t0" )
--- matrix with size 2x2 and default value 1/3
-mtx = matrix (2,2,1/3)
-assert( tostring(mtx) == "0.33333333333333\t0.33333333333333\n0.33333333333333\t0.33333333333333" )
--- identity matrix with size 2x2
-mtx = matrix (2,"I")
-assert( tostring(mtx) == "1\t0\n0\t1" )
-
--- matrix.add; number
-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}})
--- matrix.add; complex
-m1 = matrix{{10,"2+6i",1},{5,1,"4-2i"}}:tocomplex()
-m2 = matrix{{3,4,5},{"2+3i",4,"6i"}}:tocomplex()
-assert(m1 + m2 == matrix{{13,"6+6i",6},{"7+3i",5,"4+4i"}}:tocomplex())
--- matrix.add; symbol
-m1 = matrix{{8,4,1},{6,8,3}}:tosymbol()
-m2 = matrix{{-8,1,3},{5,2,1}}:tosymbol()
-assert(m1 + m2 == matrix{{"8+-8","4+1","1+3"},{"6+5","8+2","3+1"}}:tosymbol())
-
--- matrix.sub; number
-m1 = matrix{{8,4,1},{6,8,3}}
-m2 = matrix{{-8,1,3},{5,2,1}}
-assert(m1 - m2 == matrix{{16,3,-2},{1,6,2}})
--- matrix.sub; complex
-m1 = matrix{{10,"2+6i",1},{5,1,"4-2i"}}:tocomplex()
-m2 = matrix{{3,4,5},{"2+3i",4,"6i"}}:tocomplex()
-assert(m1 - m2 == matrix{{7,"-2+6i",-4},{"3-3i",-3,"4-8i"}}:tocomplex())
--- matrix.sub; symbol
-m1 = matrix{{8,4,1},{6,8,3}}:tosymbol()
-m2 = matrix{{-8,1,3},{5,2,1}}:tosymbol()
-assert(m1 - m2 == matrix{{"8--8","4-1","1-3"},{"6-5","8-2","3-1"}}:tosymbol())
-
--- matrix.mul; number
-m1 = matrix{{8,4,1},{6,8,3}}
-m2 = matrix{{3,1},{2,5},{7,4}}
-assert(m1 * m2 == matrix{{39,32},{55,58}})
--- matrix.mul; complex
-m1 = matrix{{"1+2i","3-i"},{"2-2i","1+i"}}:tocomplex()
-m2 = matrix{{"i","5-i"},{2,"1-i"}}:tocomplex()
-assert( m1*m2 == matrix{{"4-i","9+5i"},{"4+4i","10-12i"}}:tocomplex() )
--- matrix.mul; symbol
-m1 = matrix{{8,4,1},{6,8,3}}:tosymbol()
-m2 = matrix{{3,1},{2,5},{7,4}}:tosymbol()
-assert(m1 * m2 == matrix{
-  {"(8)*(3)+(4)*(2)+(1)*(7)", "(8)*(1)+(4)*(5)+(1)*(4)"},
-  {"(6)*(3)+(8)*(2)+(3)*(7)", "(6)*(1)+(8)*(5)+(3)*(4)"}
-}:tosymbol())
-
--- matrix.div; number, same for complex, not for symbol
-m1 = matrix {{1,2},{3,4}}
-m2 = matrix {{4,5},{6,7}}
-assert( m1*m2^-1 == m1/m2 )
-
--- matrix.divnum; number, same complex
-m2 = matrix {{4,5},{6,7}}
-assert( m2/2 == matrix{{2,2.5},{3,3.5}} )
-mtx = matrix {{3,5,1},{2,4,5},{1,2,2}}
-assert( 2 / mtx == matrix{{4,16,-42},{-2,-10,26},{0,2,-4}} )
--- matrix.mulnum; symbol
-m1 = m1:tosymbol()
-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)"}}:tosymbol() )
-
--- matrix.pow; number, same complex
-mtx = matrix{{3,5,1},{2,4,5},{1,2,2}}
-assert(mtx^3 == matrix{{164,308,265},{161,303,263},{76,143,124}})
-assert(mtx^1 == mtx)
-assert(mtx^0 == matrix{{1,0,0},{0,1,0},{0,0,1}} )
-assert(mtx^-1 == mtx:invert())
-assert(mtx^-3 == (mtx^-1)^3)
-mtx = matrix{{1,2,3},{1,2,3},{3,2,1}} -- singular
-assert(mtx^-1 == nil)
-local m1,rank = mtx:invert(); assert(m1 == nil and rank == 2)
-mtx = matrix{{1,2},{3,4},{5,6}} -- non-square
-assert(select(2, pcall(function() return mtx^-1 end))
-       :find("matrix not square"))
-
--- matrix.det; number
-mtx = matrix {{1,4,3,2},{2,1,-1,-1},{-3,2,2,-2},{-1,-5,-4,1}}
-assert( mtx:det() == 78 )
--- matrix.det; complex
-m1 = matrix{{"1+2i","3-i"},{"2-2i","1+i"}}:tocomplex()
-m2 = matrix{{"i","5-i"},{2,"1-i"}}:tocomplex()
-m3 = m1*m2
--- (checked in maple)
-assert( m3:det() == complex "12-114i" )
-mtx = {{"2+3i","1+4i","-2i",3,2},
-	{"2i",3,"2+3i",0,3},
-	{3,"-2i",6,"4+5i",0},
-	{1,"1+2i",3,5,7},
-	{"-3+3i","3+3i",3,-8,2}}
-matrix(mtx):tocomplex()
--- (checked in maple)
-assert( mtx:det():round(10) == complex "5527+2687i" )
-
--- matrix.invert; number
---1
-mtx = matrix{{3,5,1},{2,4,5},{1,2,2}}
-local mtxinv = matrix{{2,8,-21},{-1,-5,13},{0,1,-2}}
-assert( mtx^-1 == mtxinv )
---2
-mtx = matrix{{1,0,2},{4,1,1},{3,2,-7}} 
-local mtxinv = matrix{{-9,4,-2},{31,-13,7},{5,-2,1}}
-assert( mtx^-1 == mtxinv )
--- matrix.invert; complex
-mtx = {
-{ 3,"4-3i",1},
-{3,"-1+5i",-3},
-{4,0,7},
-}
-matrix.tocomplex( mtx )
-local mtxinv = mtx^-1
-local mtxinvcomp = {
-{"0.13349-0.07005i","0.14335+0.03609i","0.04237+0.02547i"},
-{"0.08771+0.10832i","-0.04519-0.0558i","-0.0319-0.03939i"},
-{"-0.07628+0.04003i","-0.08192-0.02062i","0.11865-0.01456i"},}
-mtxinvcomp = matrix( mtxinvcomp )
-mtxinv:round( 5 )
-mtxinv:remcomplex()
-assert( mtxinvcomp == mtxinv  )
-
--- matrix.sqrt; number
-local m1 = matrix{{4,2,1},{1,5,4},{1,5,2}}
-local m2 = m1*m1
-local msqrt = m2:sqrt()
-assert((m2 - msqrt^2):normmax() < 1E-12)
--- matrix.sqrt; complex
-local m1 = matrix{{4,"2+i",1},{1,5,"4-2i"},{1,"5+3i",2}}:tocomplex()
-local m2 = m1*m1
-local msqrt = m2:sqrt()
-assert((m2 - msqrt^2):normmax() < 1E-12)
-
--- matrix.root; number
-local p = 3
-local m1 = matrix {{4,2,3},{1,9,7},{6,5,8}}
-local m2 = m1^p
-local mroot = m2:root(p)
-assert((m2 - mroot^p):normmax() < 1E-7)
--- matrix.root; complex
-local m1 = matrix{{4,"2+i",1},{1,5,"4-2i"},{1,"5+3i",2}}:tocomplex()
-local m2 = m1^p
-local mroot = m2:root(p)
-assert((m2 - mroot^p):normmax() < 1E-7)
-
--- matrix.normf
-mtx = matrix{{2,3},{-2,-3}}
-assert(mtx:normf() == math.sqrt(2^2+3^2+2^2+3^2))
-mtx = matrix{{'2i','3'},{'-2i','-3'}}:tocomplex()
-assert(mtx:normf() == math.sqrt(2^2+3^2+2^2+3^2))
-mtx = matrix{{'a','b'},{'c','d'}}:tosymbol()
-assert(tostring(mtx:normf()) == "(0+((a):abs())^(2)+((b):abs())^(2)+((c):abs())^(2)+((d):abs())^(2)):sqrt()")
-
--- matrix.normmax
--- note: symbolic matrices not supported
-mtx = matrix{{2,3},{-2,-4}}
-assert(mtx:normmax() == 4)
-mtx = matrix{{'2i','3'},{'-2i','-4i'}}:tocomplex()
-assert(mtx:normmax() == 4)
-
--- matrix.transpose
--- test transpose; number, same for all others
-mtx = matrix{{3,5,1},{2,4,5},{1,2,2}}
-assert(mtx^'T' == matrix{{3,2,1},{5,4,2},{1,5,2}})
-
--- matrix.rotl; number, same for all others
-local m1 = matrix{{2,3},{4,5},{6,7}}
-assert( m1:rotl() == matrix{{3,5,7},{2,4,6}} )
--- matris.rotr; number, same for all others
-assert( m1:rotr() == matrix{{6,4,2},{7,5,3}} )
-
--- matrix.tostring; number
-mtx = matrix{{4,2,-3},{3,-5,2}}
-assert(tostring(mtx) == "4\t2\t-3\n3\t-5\t2" )
--- matrix.tostring; complex
-mtx = matrix{{4,"2+i"},{"3-4i",5}}:tocomplex()
-assert(tostring(mtx) == "4\t2+i\n3-4i\t5" )
--- matrix.tostring; tensor
-local mt = matrix {{{1,2},{3,4}},{{5,6},{7,8}}}
-assert( tostring(mt) == "[1,2]\t[3,4]\n[5,6]\t[7,8]" )
-local i,j,k = mt:size()
-assert( i == 2 ); assert( j == 2 );assert( k == 2 )
-
--- matrix.latex; tensor
-local mt = matrix {{{1,2},{3,4}},{{5,6},{7,8}}}
-assert( mt:latex() == "$\\left( \\begin{array}{cc}\n\t[1,2] & [3,4] \\\\\n\t[5,6] & [7,8]\n\\end{array} \\right)$" )
-
--- matrix.cross
-local e1 = matrix{ 1,0,0 }
-local e2 = matrix{ 0,1,0 }
-local e3 = e1:cross( e2 )
-assert( e3 == matrix{ 0,0,1 } )
-
--- matrix.scalar
-local v1 = matrix{ 2,3,0, }
-local v2 = matrix{ 0,3,4 }
-local vx = v1:cross( v2 )
-assert( v1:scalar( vx ) == 0 )
-assert( vx:scalar( v2 ) == 0 )
-
--- matrix.len
-assert( v2:len() == math.sqrt( 3^2+4^2 ) )
-
---// test symbolic
-ms = matrix {{ "a",1 },{2,"b"}}:tosymbol()
-ms2 = matrix {{ "a",2 },{"b",3}}:tosymbol()
-ms3 = ms2+ms
-ms3 = ms3:replace( "a",4,"b",2 )
-ms3 = ms3:solve()
-assert( ms3 == matrix {{8,3},{4,5}} )
-ms4 = ms2*ms
-ms4 = ms4:replace( "a",4,"b",2 )
-ms4 = ms4:solve()
-assert( ms4 == matrix {{20,8},{14,8}} )
-
--- __unm
-mtx = matrix{{4,2},{3,5}}
-assert(-mtx == matrix{{-4,-2},{-3,-5}})
-
--- test inverted with big table
---[[mtx = matrix:new( 100,100 )
-mtx:random( -100,100 )
---mtx:print()
-t1 = os.clock()
-identm = mtx*mtx^-1
-print("Diff:",os.clock()-t1 )
--- round to 8 decimal points
--- this depends on the numbers used, small, big, range
--- the average error in this calculation was 10^-9
-identm:round( 8 )
---identm:print()
-ident = matrix:new( 100, "I" )
-assert( identm == ident )--]]
-
-
-
-local t = {}
-for i,v in pairs( matrix ) do
-	table.insert( t, i )
-end
-table.sort( t )
-for i,v in ipairs( t ) do
-	--print( "matrix."..v )
+local matrix, complex = require "matrix"
+
+local mtx, m1,m2,m3,m4,m5, ms,ms1,ms2,ms3,ms4
+
+-- test matrix:new/matrix call function
+-- normal matrix
+mtx = matrix {{1,2},{3,4}}
+assert( tostring(mtx) == "1\t2\n3\t4" )
+-- vector
+mtx = matrix {1,2,3}
+assert( tostring(mtx) == "1\n2\n3" )
+-- matrix with size 2x2
+mtx = matrix (2,2)
+assert( tostring(mtx) == "0\t0\n0\t0" )
+-- matrix with size 2x2 and default value 1/3
+mtx = matrix (2,2,1/3)
+assert( tostring(mtx) == "0.33333333333333\t0.33333333333333\n0.33333333333333\t0.33333333333333" )
+-- identity matrix with size 2x2
+mtx = matrix (2,"I")
+assert( tostring(mtx) == "1\t0\n0\t1" )
+
+-- matrix.add; number
+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}})
+-- matrix.add; complex
+m1 = matrix{{10,"2+6i",1},{5,1,"4-2i"}}:tocomplex()
+m2 = matrix{{3,4,5},{"2+3i",4,"6i"}}:tocomplex()
+assert(m1 + m2 == matrix{{13,"6+6i",6},{"7+3i",5,"4+4i"}}:tocomplex())
+-- matrix.add; symbol
+m1 = matrix{{8,4,1},{6,8,3}}:tosymbol()
+m2 = matrix{{-8,1,3},{5,2,1}}:tosymbol()
+assert(m1 + m2 == matrix{{"8+-8","4+1","1+3"},{"6+5","8+2","3+1"}}:tosymbol())
+
+-- matrix.sub; number
+m1 = matrix{{8,4,1},{6,8,3}}
+m2 = matrix{{-8,1,3},{5,2,1}}
+assert(m1 - m2 == matrix{{16,3,-2},{1,6,2}})
+-- matrix.sub; complex
+m1 = matrix{{10,"2+6i",1},{5,1,"4-2i"}}:tocomplex()
+m2 = matrix{{3,4,5},{"2+3i",4,"6i"}}:tocomplex()
+assert(m1 - m2 == matrix{{7,"-2+6i",-4},{"3-3i",-3,"4-8i"}}:tocomplex())
+-- matrix.sub; symbol
+m1 = matrix{{8,4,1},{6,8,3}}:tosymbol()
+m2 = matrix{{-8,1,3},{5,2,1}}:tosymbol()
+assert(m1 - m2 == matrix{{"8--8","4-1","1-3"},{"6-5","8-2","3-1"}}:tosymbol())
+
+-- matrix.mul; number
+m1 = matrix{{8,4,1},{6,8,3}}
+m2 = matrix{{3,1},{2,5},{7,4}}
+assert(m1 * m2 == matrix{{39,32},{55,58}})
+-- matrix.mul; complex
+m1 = matrix{{"1+2i","3-i"},{"2-2i","1+i"}}:tocomplex()
+m2 = matrix{{"i","5-i"},{2,"1-i"}}:tocomplex()
+assert( m1*m2 == matrix{{"4-i","9+5i"},{"4+4i","10-12i"}}:tocomplex() )
+-- matrix.mul; symbol
+m1 = matrix{{8,4,1},{6,8,3}}:tosymbol()
+m2 = matrix{{3,1},{2,5},{7,4}}:tosymbol()
+assert(m1 * m2 == matrix{
+  {"(8)*(3)+(4)*(2)+(1)*(7)", "(8)*(1)+(4)*(5)+(1)*(4)"},
+  {"(6)*(3)+(8)*(2)+(3)*(7)", "(6)*(1)+(8)*(5)+(3)*(4)"}
+}:tosymbol())
+
+-- matrix.div; number, same for complex, not for symbol
+m1 = matrix {{1,2},{3,4}}
+m2 = matrix {{4,5},{6,7}}
+assert( m1*m2^-1 == m1/m2 )
+
+-- matrix.divnum; number, same complex
+m2 = matrix {{4,5},{6,7}}
+assert( m2/2 == matrix{{2,2.5},{3,3.5}} )
+mtx = matrix {{3,5,1},{2,4,5},{1,2,2}}
+assert( 2 / mtx == matrix{{4,16,-42},{-2,-10,26},{0,2,-4}} )
+-- matrix.mulnum; symbol
+m1 = m1:tosymbol()
+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)"}}:tosymbol() )
+
+-- matrix.pow; number, same complex
+mtx = matrix{{3,5,1},{2,4,5},{1,2,2}}
+assert(mtx^3 == matrix{{164,308,265},{161,303,263},{76,143,124}})
+assert(mtx^1 == mtx)
+assert(mtx^0 == matrix{{1,0,0},{0,1,0},{0,0,1}} )
+assert(mtx^-1 == mtx:invert())
+assert(mtx^-3 == (mtx^-1)^3)
+mtx = matrix{{1,2,3},{1,2,3},{3,2,1}} -- singular
+assert(mtx^-1 == nil)
+local m1,rank = mtx:invert(); assert(m1 == nil and rank == 2)
+mtx = matrix{{1,2},{3,4},{5,6}} -- non-square
+assert(select(2, pcall(function() return mtx^-1 end))
+       :find("matrix not square"))
+
+-- matrix.det; number
+mtx = matrix {{1,4,3,2},{2,1,-1,-1},{-3,2,2,-2},{-1,-5,-4,1}}
+assert( mtx:det() == 78 )
+-- matrix.det; complex
+m1 = matrix{{"1+2i","3-i"},{"2-2i","1+i"}}:tocomplex()
+m2 = matrix{{"i","5-i"},{2,"1-i"}}:tocomplex()
+m3 = m1*m2
+-- (checked in maple)
+assert( m3:det() == complex "12-114i" )
+mtx = {{"2+3i","1+4i","-2i",3,2},
+	{"2i",3,"2+3i",0,3},
+	{3,"-2i",6,"4+5i",0},
+	{1,"1+2i",3,5,7},
+	{"-3+3i","3+3i",3,-8,2}}
+matrix(mtx):tocomplex()
+-- (checked in maple)
+assert( mtx:det():round(10) == complex "5527+2687i" )
+
+-- matrix.invert; number
+--1
+mtx = matrix{{3,5,1},{2,4,5},{1,2,2}}
+local mtxinv = matrix{{2,8,-21},{-1,-5,13},{0,1,-2}}
+assert( mtx^-1 == mtxinv )
+--2
+mtx = matrix{{1,0,2},{4,1,1},{3,2,-7}} 
+local mtxinv = matrix{{-9,4,-2},{31,-13,7},{5,-2,1}}
+assert( mtx^-1 == mtxinv )
+-- matrix.invert; complex
+mtx = {
+{ 3,"4-3i",1},
+{3,"-1+5i",-3},
+{4,0,7},
+}
+matrix.tocomplex( mtx )
+local mtxinv = mtx^-1
+local mtxinvcomp = {
+{"0.13349-0.07005i","0.14335+0.03609i","0.04237+0.02547i"},
+{"0.08771+0.10832i","-0.04519-0.0558i","-0.0319-0.03939i"},
+{"-0.07628+0.04003i","-0.08192-0.02062i","0.11865-0.01456i"},}
+mtxinvcomp = matrix( mtxinvcomp )
+mtxinv:round( 5 )
+mtxinv:remcomplex()
+assert( mtxinvcomp == mtxinv  )
+
+-- matrix.sqrt; number
+local m1 = matrix{{4,2,1},{1,5,4},{1,5,2}}
+local m2 = m1*m1
+local msqrt = m2:sqrt()
+assert((m2 - msqrt^2):normmax() < 1E-12)
+-- matrix.sqrt; complex
+local m1 = matrix{{4,"2+i",1},{1,5,"4-2i"},{1,"5+3i",2}}:tocomplex()
+local m2 = m1*m1
+local msqrt = m2:sqrt()
+assert((m2 - msqrt^2):normmax() < 1E-12)
+
+-- matrix.root; number
+local p = 3
+local m1 = matrix {{4,2,3},{1,9,7},{6,5,8}}
+local m2 = m1^p
+local mroot = m2:root(p)
+assert((m2 - mroot^p):normmax() < 1E-7)
+-- matrix.root; complex
+local m1 = matrix{{4,"2+i",1},{1,5,"4-2i"},{1,"5+3i",2}}:tocomplex()
+local m2 = m1^p
+local mroot = m2:root(p)
+assert((m2 - mroot^p):normmax() < 1E-7)
+
+-- matrix.normf
+mtx = matrix{{2,3},{-2,-3}}
+assert(mtx:normf() == math.sqrt(2^2+3^2+2^2+3^2))
+mtx = matrix{{'2i','3'},{'-2i','-3'}}:tocomplex()
+assert(mtx:normf() == math.sqrt(2^2+3^2+2^2+3^2))
+mtx = matrix{{'a','b'},{'c','d'}}:tosymbol()
+assert(tostring(mtx:normf()) == "(0+((a):abs())^(2)+((b):abs())^(2)+((c):abs())^(2)+((d):abs())^(2)):sqrt()")
+
+-- matrix.normmax
+-- note: symbolic matrices not supported
+mtx = matrix{{2,3},{-2,-4}}
+assert(mtx:normmax() == 4)
+mtx = matrix{{'2i','3'},{'-2i','-4i'}}:tocomplex()
+assert(mtx:normmax() == 4)
+
+-- matrix.transpose
+-- test transpose; number, same for all others
+mtx = matrix{{3,5,1},{2,4,5},{1,2,2}}
+assert(mtx^'T' == matrix{{3,2,1},{5,4,2},{1,5,2}})
+
+-- matrix.rotl; number, same for all others
+local m1 = matrix{{2,3},{4,5},{6,7}}
+assert( m1:rotl() == matrix{{3,5,7},{2,4,6}} )
+-- matris.rotr; number, same for all others
+assert( m1:rotr() == matrix{{6,4,2},{7,5,3}} )
+
+-- matrix.tostring; number
+mtx = matrix{{4,2,-3},{3,-5,2}}
+assert(tostring(mtx) == "4\t2\t-3\n3\t-5\t2" )
+-- matrix.tostring; complex
+mtx = matrix{{4,"2+i"},{"3-4i",5}}:tocomplex()
+assert(tostring(mtx) == "4\t2+i\n3-4i\t5" )
+-- matrix.tostring; tensor
+local mt = matrix {{{1,2},{3,4}},{{5,6},{7,8}}}
+assert( tostring(mt) == "[1,2]\t[3,4]\n[5,6]\t[7,8]" )
+local i,j,k = mt:size()
+assert( i == 2 ); assert( j == 2 );assert( k == 2 )
+
+-- matrix.latex; tensor
+local mt = matrix {{{1,2},{3,4}},{{5,6},{7,8}}}
+assert( mt:latex() == "$\\left( \\begin{array}{cc}\n\t[1,2] & [3,4] \\\\\n\t[5,6] & [7,8]\n\\end{array} \\right)$" )
+
+-- matrix.cross
+local e1 = matrix{ 1,0,0 }
+local e2 = matrix{ 0,1,0 }
+local e3 = e1:cross( e2 )
+assert( e3 == matrix{ 0,0,1 } )
+
+-- matrix.scalar
+local v1 = matrix{ 2,3,0, }
+local v2 = matrix{ 0,3,4 }
+local vx = v1:cross( v2 )
+assert( v1:scalar( vx ) == 0 )
+assert( vx:scalar( v2 ) == 0 )
+
+-- matrix.len
+assert( v2:len() == math.sqrt( 3^2+4^2 ) )
+
+--// test symbolic
+ms = matrix {{ "a",1 },{2,"b"}}:tosymbol()
+ms2 = matrix {{ "a",2 },{"b",3}}:tosymbol()
+ms3 = ms2+ms
+ms3 = ms3:replace( "a",4,"b",2 )
+ms3 = ms3:solve()
+assert( ms3 == matrix {{8,3},{4,5}} )
+ms4 = ms2*ms
+ms4 = ms4:replace( "a",4,"b",2 )
+ms4 = ms4:solve()
+assert( ms4 == matrix {{20,8},{14,8}} )
+
+-- __unm
+mtx = matrix{{4,2},{3,5}}
+assert(-mtx == matrix{{-4,-2},{-3,-5}})
+
+-- test inverted with big table
+--[[mtx = matrix:new( 100,100 )
+mtx:random( -100,100 )
+--mtx:print()
+t1 = os.clock()
+identm = mtx*mtx^-1
+print("Diff:",os.clock()-t1 )
+-- round to 8 decimal points
+-- this depends on the numbers used, small, big, range
+-- the average error in this calculation was 10^-9
+identm:round( 8 )
+--identm:print()
+ident = matrix:new( 100, "I" )
+assert( identm == ident )--]]
+
+
+
+local t = {}
+for i,v in pairs( matrix ) do
+	table.insert( t, i )
+end
+table.sort( t )
+for i,v in ipairs( t ) do
+	--print( "matrix."..v )
 end
\ No newline at end of file