import version 0.2.8

master
David Manura 2010-09-22 20:20:32 -04:00
commit 5b3e7b19c5
9 changed files with 2457 additions and 0 deletions

385
complex.lua Normal file
View File

@ -0,0 +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 //--
--///////////////--

55
complex_changelog.txt Normal file
View File

@ -0,0 +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)
built

118
fit.lua Normal file
View File

@ -0,0 +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 //--
--///////////////--

7
fit_changelog.txt Normal file
View File

@ -0,0 +1,7 @@
fit changelog:
v 0.2: 2007-08-26
- some optimising, for the input
v 0.1: 2007-08-11
- built from former fit

1267
matrix.lua Normal file

File diff suppressed because it is too large Load Diff

160
matrix_changelog.txt Normal file
View File

@ -0,0 +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]
added simple operations +-..

185
test_complex.lua Normal file
View File

@ -0,0 +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"
assert( -cx == complex "-2-3i" )

22
test_fit.lua Normal file
View File

@ -0,0 +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.." )")

258
test_matrix.lua Normal file
View File

@ -0,0 +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 )
end