Tue Feb 10 12:28:14 2004 run on Linux
if lisp !*rounded then rounded_was_on := t
else rounded_was_on := nil;
mat1 := mat((1,2,3,4,5),(2,3,4,5,6),(3,4,5,6,7),(4,5,6,7,8),(5,6,7,8,9));
[1 2 3 4 5]
[ ]
[2 3 4 5 6]
[ ]
mat1 := [3 4 5 6 7]
[ ]
[4 5 6 7 8]
[ ]
[5 6 7 8 9]
mat2 := mat((1,1,1,1),(2,2,2,2),(3,3,3,3),(4,4,4,4));
[1 1 1 1]
[ ]
[2 2 2 2]
mat2 := [ ]
[3 3 3 3]
[ ]
[4 4 4 4]
mat3 := mat((x),(x),(x),(x));
[x]
[ ]
[x]
mat3 := [ ]
[x]
[ ]
[x]
mat4 := mat((3,3),(4,4),(5,5),(6,6));
[3 3]
[ ]
[4 4]
mat4 := [ ]
[5 5]
[ ]
[6 6]
mat5 := mat((1,2,1,1),(1,2,3,1),(4,5,1,2),(3,4,5,6));
[1 2 1 1]
[ ]
[1 2 3 1]
mat5 := [ ]
[4 5 1 2]
[ ]
[3 4 5 6]
mat6 := mat((i+1,i+2,i+3),(4,5,2),(1,i,0));
[i + 1 i + 2 i + 3]
[ ]
mat6 := [ 4 5 2 ]
[ ]
[ 1 i 0 ]
mat7 := mat((1,1,0),(1,3,1),(0,1,1));
[1 1 0]
[ ]
mat7 := [1 3 1]
[ ]
[0 1 1]
mat8 := mat((1,3),(-4,3));
[1 3]
mat8 := [ ]
[-4 3]
mat9 := mat((1,2,3,4),(9,8,7,6));
[1 2 3 4]
mat9 := [ ]
[9 8 7 6]
poly := x^7+x^5+4*x^4+5*x^3+12;
7 5 4 3
poly := x + x + 4*x + 5*x + 12
poly1 := x^2+x*y^3+x*y*z^3+y*x+2+y*3;
2 3 3
poly1 := x + x*y + x*y*z + x*y + 3*y + 2
on errcont;
% Basis matrix manipulations.
add_columns(mat1,1,2,5*y);
[1 5*y + 2 3 4 5]
[ ]
[2 10*y + 3 4 5 6]
[ ]
[3 15*y + 4 5 6 7]
[ ]
[4 5*(4*y + 1) 6 7 8]
[ ]
[5 25*y + 6 7 8 9]
add_rows(mat1,1,2,x);
[ 1 2 3 4 5 ]
[ ]
[x + 2 2*x + 3 3*x + 4 4*x + 5 5*x + 6]
[ ]
[ 3 4 5 6 7 ]
[ ]
[ 4 5 6 7 8 ]
[ ]
[ 5 6 7 8 9 ]
add_to_columns(mat1,3,1000);
[1 2 1003 4 5]
[ ]
[2 3 1004 5 6]
[ ]
[3 4 1005 6 7]
[ ]
[4 5 1006 7 8]
[ ]
[5 6 1007 8 9]
add_to_columns(mat1,{1,2,3},y);
[y + 1 y + 2 y + 3 4 5]
[ ]
[y + 2 y + 3 y + 4 5 6]
[ ]
[y + 3 y + 4 y + 5 6 7]
[ ]
[y + 4 y + 5 y + 6 7 8]
[ ]
[y + 5 y + 6 y + 7 8 9]
add_to_rows(mat1,2,1000);
[ 1 2 3 4 5 ]
[ ]
[1002 1003 1004 1005 1006]
[ ]
[ 3 4 5 6 7 ]
[ ]
[ 4 5 6 7 8 ]
[ ]
[ 5 6 7 8 9 ]
add_to_rows(mat1,{1,2,3},x);
[x + 1 x + 2 x + 3 x + 4 x + 5]
[ ]
[x + 2 x + 3 x + 4 x + 5 x + 6]
[ ]
[x + 3 x + 4 x + 5 x + 6 x + 7]
[ ]
[ 4 5 6 7 8 ]
[ ]
[ 5 6 7 8 9 ]
augment_columns(mat1,2);
[2]
[ ]
[3]
[ ]
[4]
[ ]
[5]
[ ]
[6]
augment_columns(mat1,{1,2,5});
[1 2 5]
[ ]
[2 3 6]
[ ]
[3 4 7]
[ ]
[4 5 8]
[ ]
[5 6 9]
stack_rows(mat1,3);
[3 4 5 6 7]
stack_rows(mat1,{1,3,5});
[1 2 3 4 5]
[ ]
[3 4 5 6 7]
[ ]
[5 6 7 8 9]
char_poly(mat1,x);
3 2
x *(x - 25*x - 50)
column_dim(mat2);
4
row_dim(mat1);
5
copy_into(mat7,mat1,2,3);
[1 2 3 4 5]
[ ]
[2 3 1 1 0]
[ ]
[3 4 1 3 1]
[ ]
[4 5 0 1 1]
[ ]
[5 6 7 8 9]
copy_into(mat7,mat1,5,5);
***** Error in copy_into: the matrix
[1 1 0]
[ ]
[1 3 1]
[ ]
[0 1 1]
does not fit into
[1 2 3 4 5]
[ ]
[2 3 4 5 6]
[ ]
[3 4 5 6 7]
[ ]
[4 5 6 7 8]
[ ]
[5 6 7 8 9]
at position 5,5.
diagonal(3);
[3]
% diagonal can take both a list of arguments or just the arguments.
diagonal({mat2,mat6});
[1 1 1 1 0 0 0 ]
[ ]
[2 2 2 2 0 0 0 ]
[ ]
[3 3 3 3 0 0 0 ]
[ ]
[4 4 4 4 0 0 0 ]
[ ]
[0 0 0 0 i + 1 i + 2 i + 3]
[ ]
[0 0 0 0 4 5 2 ]
[ ]
[0 0 0 0 1 i 0 ]
diagonal(mat1,mat2,mat5);
[1 2 3 4 5 0 0 0 0 0 0 0 0]
[ ]
[2 3 4 5 6 0 0 0 0 0 0 0 0]
[ ]
[3 4 5 6 7 0 0 0 0 0 0 0 0]
[ ]
[4 5 6 7 8 0 0 0 0 0 0 0 0]
[ ]
[5 6 7 8 9 0 0 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 1 1 1 1 0 0 0 0]
[ ]
[0 0 0 0 0 2 2 2 2 0 0 0 0]
[ ]
[0 0 0 0 0 3 3 3 3 0 0 0 0]
[ ]
[0 0 0 0 0 4 4 4 4 0 0 0 0]
[ ]
[0 0 0 0 0 0 0 0 0 1 2 1 1]
[ ]
[0 0 0 0 0 0 0 0 0 1 2 3 1]
[ ]
[0 0 0 0 0 0 0 0 0 4 5 1 2]
[ ]
[0 0 0 0 0 0 0 0 0 3 4 5 6]
extend(mat1,3,2,x);
[1 2 3 4 5 x x]
[ ]
[2 3 4 5 6 x x]
[ ]
[3 4 5 6 7 x x]
[ ]
[4 5 6 7 8 x x]
[ ]
[5 6 7 8 9 x x]
[ ]
[x x x x x x x]
[ ]
[x x x x x x x]
[ ]
[x x x x x x x]
find_companion(mat5,x);
2
x - 2*x - 2
get_columns(mat1,1);
{
[1]
[ ]
[2]
[ ]
[3]
[ ]
[4]
[ ]
[5]
}
get_columns(mat1,{1,2});
{
[1]
[ ]
[2]
[ ]
[3]
[ ]
[4]
[ ]
[5]
,
[2]
[ ]
[3]
[ ]
[4]
[ ]
[5]
[ ]
[6]
}
get_rows(mat1,3);
{
[3 4 5 6 7]
}
get_rows(mat1,{1,3});
{
[1 2 3 4 5]
,
[3 4 5 6 7]
}
hermitian_tp(mat6);
[ - i + 1 4 1 ]
[ ]
[ - i + 2 5 - i]
[ ]
[ - i + 3 2 0 ]
% matrix_augment and matrix_stack can take both a list of arguments
% or just the arguments.
matrix_augment({mat1,mat2});
***** Error in matrix_augment:
***** all input matrices must have the same row dimension.
matrix_augment(mat4,mat2,mat4);
[3 3 1 1 1 1 3 3]
[ ]
[4 4 2 2 2 2 4 4]
[ ]
[5 5 3 3 3 3 5 5]
[ ]
[6 6 4 4 4 4 6 6]
matrix_stack(mat1,mat2);
***** Error in matrix_stack:
***** all input matrices must have the same column dimension.
matrix_stack({mat6,mat((z,z,z)),mat7});
[i + 1 i + 2 i + 3]
[ ]
[ 4 5 2 ]
[ ]
[ 1 i 0 ]
[ ]
[ z z z ]
[ ]
[ 1 1 0 ]
[ ]
[ 1 3 1 ]
[ ]
[ 0 1 1 ]
minor(mat1,2,3);
[1 2 4 5]
[ ]
[3 4 6 7]
[ ]
[4 5 7 8]
[ ]
[5 6 8 9]
mult_columns(mat1,3,y);
[1 2 3*y 4 5]
[ ]
[2 3 4*y 5 6]
[ ]
[3 4 5*y 6 7]
[ ]
[4 5 6*y 7 8]
[ ]
[5 6 7*y 8 9]
mult_columns(mat1,{2,3,4},100);
[1 200 300 400 5]
[ ]
[2 300 400 500 6]
[ ]
[3 400 500 600 7]
[ ]
[4 500 600 700 8]
[ ]
[5 600 700 800 9]
mult_rows(mat1,2,x);
[ 1 2 3 4 5 ]
[ ]
[2*x 3*x 4*x 5*x 6*x]
[ ]
[ 3 4 5 6 7 ]
[ ]
[ 4 5 6 7 8 ]
[ ]
[ 5 6 7 8 9 ]
mult_rows(mat1,{1,3,5},10);
[10 20 30 40 50]
[ ]
[2 3 4 5 6 ]
[ ]
[30 40 50 60 70]
[ ]
[4 5 6 7 8 ]
[ ]
[50 60 70 80 90]
pivot(mat1,3,3);
[ - 4 - 2 2 4 ]
[------ ------ 0 --- --- ]
[ 5 5 5 5 ]
[ ]
[ - 2 - 1 1 2 ]
[------ ------ 0 --- --- ]
[ 5 5 5 5 ]
[ ]
[ 3 4 5 6 7 ]
[ ]
[ 2 1 - 1 - 2 ]
[ --- --- 0 ------ ------]
[ 5 5 5 5 ]
[ ]
[ 4 2 - 2 - 4 ]
[ --- --- 0 ------ ------]
[ 5 5 5 5 ]
rows_pivot(mat1,3,3,{1,5});
[ - 4 - 2 2 4 ]
[------ ------ 0 --- --- ]
[ 5 5 5 5 ]
[ ]
[ 2 3 4 5 6 ]
[ ]
[ 3 4 5 6 7 ]
[ ]
[ 4 5 6 7 8 ]
[ ]
[ 4 2 - 2 - 4 ]
[ --- --- 0 ------ ------]
[ 5 5 5 5 ]
remove_columns(mat1,3);
[1 2 4 5]
[ ]
[2 3 5 6]
[ ]
[3 4 6 7]
[ ]
[4 5 7 8]
[ ]
[5 6 8 9]
remove_columns(mat1,{2,3,4});
[1 5]
[ ]
[2 6]
[ ]
[3 7]
[ ]
[4 8]
[ ]
[5 9]
remove_rows(mat1,2);
[1 2 3 4 5]
[ ]
[3 4 5 6 7]
[ ]
[4 5 6 7 8]
[ ]
[5 6 7 8 9]
remove_rows(mat1,{1,3});
[2 3 4 5 6]
[ ]
[4 5 6 7 8]
[ ]
[5 6 7 8 9]
remove_rows(mat1,{1,2,3,4,5});
***** Warning in remove_rows:
all the rows have been removed. Returning nil.
swap_columns(mat1,2,4);
[1 4 3 2 5]
[ ]
[2 5 4 3 6]
[ ]
[3 6 5 4 7]
[ ]
[4 7 6 5 8]
[ ]
[5 8 7 6 9]
swap_rows(mat1,1,2);
[2 3 4 5 6]
[ ]
[1 2 3 4 5]
[ ]
[3 4 5 6 7]
[ ]
[4 5 6 7 8]
[ ]
[5 6 7 8 9]
swap_entries(mat1,{1,1},{5,5});
[9 2 3 4 5]
[ ]
[2 3 4 5 6]
[ ]
[3 4 5 6 7]
[ ]
[4 5 6 7 8]
[ ]
[5 6 7 8 1]
% Constructors - functions that create matrices.
band_matrix(x,5);
[x 0 0 0 0]
[ ]
[0 x 0 0 0]
[ ]
[0 0 x 0 0]
[ ]
[0 0 0 x 0]
[ ]
[0 0 0 0 x]
band_matrix({x,y,z},6);
[y z 0 0 0 0]
[ ]
[x y z 0 0 0]
[ ]
[0 x y z 0 0]
[ ]
[0 0 x y z 0]
[ ]
[0 0 0 x y z]
[ ]
[0 0 0 0 x y]
block_matrix(1,2,{mat1,mat2});
***** Error in block_matrix: row dimensions of
***** matrices into block_matrix are not compatible
block_matrix(2,3,{mat2,mat3,mat2,mat3,mat2,mat2});
[1 1 1 1 x 1 1 1 1]
[ ]
[2 2 2 2 x 2 2 2 2]
[ ]
[3 3 3 3 x 3 3 3 3]
[ ]
[4 4 4 4 x 4 4 4 4]
[ ]
[x 1 1 1 1 1 1 1 1]
[ ]
[x 2 2 2 2 2 2 2 2]
[ ]
[x 3 3 3 3 3 3 3 3]
[ ]
[x 4 4 4 4 4 4 4 4]
char_matrix(mat1,x);
[x - 1 -2 -3 -4 -5 ]
[ ]
[ -2 x - 3 -4 -5 -6 ]
[ ]
[ -3 -4 x - 5 -6 -7 ]
[ ]
[ -4 -5 -6 x - 7 -8 ]
[ ]
[ -5 -6 -7 -8 x - 9]
cfmat := coeff_matrix({x+y+4*z=10,y+x-z=20,x+y+4});
cfmat := {
[4 1 1]
[ ]
[-1 1 1]
[ ]
[0 1 1]
,
[z]
[ ]
[y]
[ ]
[x]
,
[10]
[ ]
[20]
[ ]
[-4]
}
first cfmat * second cfmat;
[x + y + 4*z]
[ ]
[ x + y - z ]
[ ]
[ x + y ]
third cfmat;
[10]
[ ]
[20]
[ ]
[-4]
companion(poly,x);
[0 0 0 0 0 0 -12]
[ ]
[1 0 0 0 0 0 0 ]
[ ]
[0 1 0 0 0 0 0 ]
[ ]
[0 0 1 0 0 0 -5 ]
[ ]
[0 0 0 1 0 0 -4 ]
[ ]
[0 0 0 0 1 0 -1 ]
[ ]
[0 0 0 0 0 1 0 ]
hessian(poly1,{w,x,y,z});
[0 0 0 0 ]
[ ]
[ 2 3 2 ]
[0 2 3*y + z + 1 3*y*z ]
[ ]
[ 2 3 2 ]
[0 3*y + z + 1 6*x*y 3*x*z ]
[ ]
[ 2 2 ]
[0 3*y*z 3*x*z 6*x*y*z]
hilbert(4,1);
[ 1 1 1 ]
[ 1 --- --- ---]
[ 2 3 4 ]
[ ]
[ 1 1 1 1 ]
[--- --- --- ---]
[ 2 3 4 5 ]
[ ]
[ 1 1 1 1 ]
[--- --- --- ---]
[ 3 4 5 6 ]
[ ]
[ 1 1 1 1 ]
[--- --- --- ---]
[ 4 5 6 7 ]
hilbert(3,y+x);
[ - 1 - 1 - 1 ]
[----------- ----------- -----------]
[ x + y - 2 x + y - 3 x + y - 4 ]
[ ]
[ - 1 - 1 - 1 ]
[----------- ----------- -----------]
[ x + y - 3 x + y - 4 x + y - 5 ]
[ ]
[ - 1 - 1 - 1 ]
[----------- ----------- -----------]
[ x + y - 4 x + y - 5 x + y - 6 ]
jacobian({x^4,x*y^2,x*y*z^3},{w,x,y,z});
[ 3 ]
[0 4*x 0 0 ]
[ ]
[ 2 ]
[0 y 2*x*y 0 ]
[ ]
[ 3 3 2]
[0 y*z x*z 3*x*y*z ]
jordan_block(x,5);
[x 1 0 0 0]
[ ]
[0 x 1 0 0]
[ ]
[0 0 x 1 0]
[ ]
[0 0 0 x 1]
[ ]
[0 0 0 0 x]
make_identity(11);
[1 0 0 0 0 0 0 0 0 0 0]
[ ]
[0 1 0 0 0 0 0 0 0 0 0]
[ ]
[0 0 1 0 0 0 0 0 0 0 0]
[ ]
[0 0 0 1 0 0 0 0 0 0 0]
[ ]
[0 0 0 0 1 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 1 0 0 0 0 0]
[ ]
[0 0 0 0 0 0 1 0 0 0 0]
[ ]
[0 0 0 0 0 0 0 1 0 0 0]
[ ]
[0 0 0 0 0 0 0 0 1 0 0]
[ ]
[0 0 0 0 0 0 0 0 0 1 0]
[ ]
[0 0 0 0 0 0 0 0 0 0 1]
on rounded;
% makes things a bit easier to read.
random_matrix(3,3,100);
[ - 8.11911717343 - 75.7167729277 30.62058083 ]
[ ]
[ - 50.0325962624 47.1655452861 35.8674263384 ]
[ ]
[ - 49.3715543826 - 97.5563670864 - 18.8861862756]
on not_negative;
random_matrix(3,3,100);
[43.8999853223 33.7140980286 33.75065406 ]
[ ]
[49.7333355117 98.9642944905 58.5331568816]
[ ]
[39.9146060895 67.7954727837 24.8684367642]
on only_integer;
random_matrix(3,3,100);
[16 77 49]
[ ]
[28 84 51]
[ ]
[84 56 57]
on symmetric;
random_matrix(3,3,100);
[89 74 91]
[ ]
[74 95 41]
[ ]
[91 41 87]
off symmetric;
on upper_matrix;
random_matrix(3,3,100);
[41 3 8 ]
[ ]
[0 31 80]
[ ]
[0 0 12]
off upper_matrix;
on lower_matrix;
random_matrix(3,3,100);
[69 0 0 ]
[ ]
[34 87 0 ]
[ ]
[78 72 13]
off lower_matrix;
on imaginary;
off not_negative;
random_matrix(3,3,100);
[ - 95*i - 72 - 57*i + 59 52*i + 46]
[ ]
[ - 40*i - 54 70*i 39*i + 28]
[ ]
[ - 40*i + 45 28*i - 81 9*i + 74 ]
off rounded;
% toeplitz and vandermonde can take both a list of arguments or just
% the arguments.
toeplitz({1,2,3,4,5});
[1 2 3 4 5]
[ ]
[2 1 2 3 4]
[ ]
[3 2 1 2 3]
[ ]
[4 3 2 1 2]
[ ]
[5 4 3 2 1]
toeplitz(x,y,z);
[x y z]
[ ]
[y x y]
[ ]
[z y x]
vandermonde({1,2,3,4,5});
[1 1 1 1 1 ]
[ ]
[1 2 4 8 16 ]
[ ]
[1 3 9 27 81 ]
[ ]
[1 4 16 64 256]
[ ]
[1 5 25 125 625]
vandermonde(x,y,z);
[ 2]
[1 x x ]
[ ]
[ 2]
[1 y y ]
[ ]
[ 2]
[1 z z ]
% kronecker_product
a1 := mat((1,2),(3,4),(5,6));
[1 2]
[ ]
a1 := [3 4]
[ ]
[5 6]
a2 := mat((1,x,1),(2,2,2),(3,3,3));
[1 x 1]
[ ]
a2 := [2 2 2]
[ ]
[3 3 3]
kronecker_product(a1,a2);
[1 x 1 2 2*x 2 ]
[ ]
[2 2 2 4 4 4 ]
[ ]
[3 3 3 6 6 6 ]
[ ]
[3 3*x 3 4 4*x 4 ]
[ ]
[6 6 6 8 8 8 ]
[ ]
[9 9 9 12 12 12]
[ ]
[5 5*x 5 6 6*x 6 ]
[ ]
[10 10 10 12 12 12]
[ ]
[15 15 15 18 18 18]
clear a1,a2;
% High level algorithms.
on rounded;
% makes output easier to read.
ch := cholesky(mat7);
ch := {
[1 0 0 ]
[ ]
[1 1.41421356237 0 ]
[ ]
[0 0.707106781187 0.707106781187]
,
[1 1 0 ]
[ ]
[0 1.41421356237 0.707106781187]
[ ]
[0 0 0.707106781187]
}
tp first ch - second ch;
[0 0 0]
[ ]
[0 0 0]
[ ]
[0 0 0]
tmp := first ch * second ch;
[1 1 0]
[ ]
tmp := [1 3.0 1]
[ ]
[0 1 1]
tmp - mat7;
[0 0 0]
[ ]
[0 0 0]
[ ]
[0 0 0]
off rounded;
gram_schmidt({1,0,0},{1,1,0},{1,1,1});
{{1,0,0},{0,1,0},{0,0,1}}
gram_schmidt({1,2},{3,4});
1 2 2*sqrt(5) - sqrt(5)
{{---------,---------},{-----------,------------}}
sqrt(5) sqrt(5) 5 5
on rounded;
% again, makes large quotients a bit more readable.
% The algorithm used for lu_decom sometimes swaps the rows of the input
% matrix so that (given matrix A, lu_decom(A) = {L,U,vec}), we find L*U
% does not equal A but a row equivalent of it. The call convert(A,vec)
% will return this row equivalent (ie: L*U = convert(A,vec)).
lu := lu_decom(mat5);
lu := {
[4 0 0 0 ]
[ ]
[1 0.75 0 0 ]
[ ]
[1 0.75 2.0 0 ]
[ ]
[3 0.25 4.0 4.33333333333]
,
[1 1.25 0.25 0.5 ]
[ ]
[0 1 1 0.666666666667]
[ ]
[0 0 1 0 ]
[ ]
[0 0 0 1 ]
,
[3,3,3,4]}
mat5;
[1 2 1 1]
[ ]
[1 2 3 1]
[ ]
[4 5 1 2]
[ ]
[3 4 5 6]
tmp := first lu * second lu;
[4 5.0 1 2.0]
[ ]
[1 2.0 1 1 ]
tmp := [ ]
[1 2.0 3.0 1 ]
[ ]
[3 4.0 5.0 6.0]
tmp1 := convert(mat5,third lu);
[4 5 1 2]
[ ]
[1 2 1 1]
tmp1 := [ ]
[1 2 3 1]
[ ]
[3 4 5 6]
tmp - tmp1;
[0 0 0 0]
[ ]
[0 0 0 0]
[ ]
[0 0 0 0]
[ ]
[0 0 0 0]
% and the complex case...
lu1 := lu_decom(mat6);
lu1 := {
[ 1 0 0 ]
[ ]
[ 4 - 4*i + 5 0 ]
[ ]
[i + 1 3 0.414634146341*i + 2.26829268293]
,
[1 i 0 ]
[ ]
[0 1 0.19512195122*i + 0.243902439024]
[ ]
[0 0 1 ]
,
[3,2,3]}
mat6;
[i + 1 i + 2 i + 3]
[ ]
[ 4 5 2 ]
[ ]
[ 1 i 0 ]
tmp := first lu1 * second lu1;
[ 1 i 0 ]
[ ]
tmp := [ 4 5 2.0 ]
[ ]
[i + 1 i + 2 i + 3.0]
tmp1 := convert(mat6,third lu1);
[ 1 i 0 ]
[ ]
tmp1 := [ 4 5 2 ]
[ ]
[i + 1 i + 2 i + 3]
tmp - tmp1;
[0 0 0]
[ ]
[0 0 0]
[ ]
[0 0 0]
mat9inv := pseudo_inverse(mat9);
[ - 0.199999999996 0.100000000013 ]
[ ]
[ - 0.0499999999988 0.0500000000037 ]
mat9inv := [ ]
[ 0.0999999999982 - 5.57819415659e-12]
[ ]
[ 0.249999999995 - 0.0500000000148 ]
mat9 * mat9inv;
[ 0.999999999982 - 0.0000000000557818791158]
[ ]
[5.54223333893e-12 1.00000000002 ]
simplex(min,2*x1+14*x2+36*x3,{-2*x1+x2+4*x3>=5,-x1-2*x2-3*x3<=2});
{45.0,{x1=0,x2=0,x3=1.25}}
simplex(max,10000 x1 + 1000 x2 + 100 x3 + 10 x4 + x5,{ x1 <= 1, 20 x1 +
x2 <= 100, 200 x1 + 20 x2 + x3 <= 10000, 2000 x1 + 200 x2 + 20 x3 + x4
<= 1000000, 20000 x1 + 2000 x2 + 200 x3 + 20 x4 + x5 <= 100000000});
{100000000,{x1=0,x2=0,x3=0,x4=0,x5=100000000.0}}
simplex(max, 5 x1 + 4 x2 + 3 x3,
{ 2 x1 + 3 x2 + x3 <= 5,
4 x1 + x2 + 2 x3 <= 11,
3 x1 + 4 x2 + 2 x3 <= 8 });
{13.0,{x1=2.0,x2=0,x3=1.0}}
simplex(min,3 x1 + 5 x2,{ x1 + 2 x2 >= 2, 22 x1 + x2 >= 3});
{5.04651162791,{x1=0.093023255813953,x2=0.95348837209302}}
simplex(max,10x+5y+5.5z,{5x+3z<=200,0.2x+0.1y+0.5z<=12,0.1x+0.2y+0.3z<=9,
30x+10y+50z<=1500});
{525.0,{x=40.0,y=25.0,z=0}}
%example of extra variables (>=0) being added.
simplex(min,x-y,{x>=-3});
*** Warning: variable y not defined in input. Has been defined as >=0.
***** Error in simplex: The problem is unbounded.
% unfeasible as simplex algorithm implies all x>=0.
simplex(min,x,{x<=-100});
***** Error in simplex: Problem has no feasible solution.
% three error examples.
simplex(maxx,x,{x>=5});
***** Error in simplex(first argument): must be either max or min.
simplex(max,x,x>=5);
***** Error in simplex(third argument}: must be a list.
simplex(max,x,{x<=y});
***** Error in simplex: The problem is unbounded.
simplex(max, 346 X11 + 346 X12 + 248 X21 + 248 X22 + 399 X31 + 399 X32 +
200 Y11 + 200 Y12 + 75 Y21 + 75 Y22 + 2.35 Z1 + 3.5 Z2,
{
4 X11 + 4 X12 + 2 X21 + 2 X22 + X31 + X32 + 250 Y11 + 250 Y12 + 125 Y21 +
125 Y22 <= 25000,
X11 + X12 + X21 + X22 + X31 + X32 + 2 Y11 + 2 Y12 + Y21 + Y22 <= 300,
20 X11 + 15 X12 + 30 Y11 + 20 Y21 + Z1 <= 1500,
40 X12 + 35 X22 + 50 X32 + 15 Y12 + 10 Y22 + Z2 = 5000,
X31 = 0,
Y11 + Y12 <= 50,
Y21 + Y22 <= 100
});
{99250.0,
{y21=0,
y22=0,
x31=0,
x11=75.0,
z1=0,
x21=225.0,
z2=5000.0,
x32=0,
x22=0,
x12=0,
y12=0,
y11=0}}
% from Marc van Dongen. Finding the first feasible solution for the
% solution of systems of linear diophantine inequalities.
simplex(max,0,{
3*X259+4*X261+3*X262+2*X263+X269+2*X270+3*X271+4*X272+5*X273+X229=2,
7*X259+11*X261+8*X262+5*X263+3*X269+6*X270+9*X271+12*X272+15*X273+X229=4,
2*X259+5*X261+4*X262+3*X263+3*X268+4*X269+5*X270+6*X271+7*X272+8*X273=1,
X262+2*X263+5*X268+4*X269+3*X270+2*X271+X272+2*X229=1,
X259+X262+2*X263+4*X268+3*X269+2*X270+X271-X273+3*X229=2,
X259+2*X261+2*X262+2*X263+3*X268+3*X269+3*X270+3*X271+3*X272+3*X273+X229=1,
X259+X261+X262+X263+X268+X269+X270+X271+X272+X273+X229=1});
{0,
{x229=0.5,
x259=0.5,
x261=0,
x262=0,
x263=0,
x268=0,
x269=0,
x270=0,
x271=0,
x272=0,
x273=0}}
svd_ans := svd(mat8);
svd_ans := {
[ 0.289784137735 0.957092029805]
[ ]
[ - 0.957092029805 0.289784137735]
,
[5.1491628629 0 ]
[ ]
[ 0 2.9130948854]
,
[ - 0.687215403194 0.726453707825 ]
[ ]
[ - 0.726453707825 - 0.687215403194]
}
tmp := tp first svd_ans * second svd_ans * third svd_ans;
[ 0.99999998509 2.9999999859 ]
tmp := [ ]
[ - 4.00000004924 2.99999995342]
tmp - mat8;
[ - 0.0000000149096008872 - 0.0000000141042799662]
[ ]
[ - 0.000000049243064737 - 0.000000046583274127 ]
mat9inv := pseudo_inverse(mat9);
[ - 0.199999999996 0.100000000013 ]
[ ]
[ - 0.0499999999988 0.0500000000037 ]
mat9inv := [ ]
[ 0.0999999999982 - 5.57819415659e-12]
[ ]
[ 0.249999999995 - 0.0500000000148 ]
mat9 * mat9inv;
[ 0.999999999982 - 0.0000000000557818791158]
[ ]
[5.54223333893e-12 1.00000000002 ]
% triang_adjoint(in_mat) calculates the
% triangularizing adjoint of in_mat
triang_adjoint(mat1);
[1 0 0 0 0]
[ ]
[-2 1 0 0 0]
[ ]
[-1 2 -1 0 0]
[ ]
[0 0 0 0 0]
[ ]
[0 0 0 0 0]
triang_adjoint(mat2);
[1 0 0 0]
[ ]
[-2 1 0 0]
[ ]
[0 0 0 0]
[ ]
[0 0 0 0]
triang_adjoint(mat5);
[1 0 0 0]
[ ]
[-1 1 0 0]
[ ]
[-3 3 0 0]
[ ]
[10 -12 -4 6]
triang_adjoint(mat6);
[ 1 0 0 ]
[ ]
[ -4 i + 1 0 ]
[ ]
[4*i - 5 3 i - 3]
triang_adjoint(mat7);
[1 0 0]
[ ]
[-1 1 0]
[ ]
[1 - 1 2]
triang_adjoint(mat8);
[1 0]
[ ]
[4 1]
triang_adjoint(mat9);
***** Error in triang_adjoint: input matrix should be square.
% testing triang_adjoint with random matrices
% the range of the integers is in one case from
% -1000 to 1000. in the other case it is from
% -1 to 1 so that the deteminant of the i-th
% submatrix equals very often to zero.
% random matrix contains arbitrary real values
off only_integer;
tmp:=random_matrix(5,5,1000);
tmp := mat(( - 558.996086656*i + 1.71931812953,76.9987188735*i + 1.19004104683,
- 739.283479439*i - 1.32534106204,742.101952123*i + 1.35926854848,
680.515777254*i + 1.56403177895),
( - 689.196170962*i + 1.49289170118,
- 232.584493916*i - 1.38227180395,280.109305836*i + 1.38865247861,
298.151479065*i - 1.19035182389, - 602.312143386*i - 1.82183796879),
( - 739.195910955*i - 1.45944960213,859.293884826*i + 1.70488070379,
359.856032683*i - 1.2966991869, - 105.409833087*i - 1.9360631701,
234.350529301*i - 1.15598520849),
(155.629059348*i + 1.09264385739, - 16.1559469072*i - 1.9425176505,
725.11578405*i - 1.05760723025,783.020942195*i - 1.28625265346,
- 544.129360355*i + 1.74790906085),
(373.562370318*i - 1.95218354686, - 722.109349973*i + 1.56309793677,
- 746.664959169*i - 1.9915755693,186.154794517*i - 1.09842189916,
435.90998986*i - 1.46175649496))
triang_adjoint tmp;
mat((1,0,0,0,0),
(689.196170962*i - 1.49289170118, - 558.996086656*i + 1.71931812953,0,0,0),
( - 1253.37955588*i + 7.64148089854e+5, - 1516.42713845*i - 4.23429448803e+5
,1078.01877642*i - 1.830851973e+5,0,0),
102791325687.0*i + 7.3752778526e+8
(------------------------------------,
i - 169.834887206
- 3.66748178757e+10*i - 6.62162769101e+6
-------------------------------------------,
i - 169.834887206
9.85957444629e+7*i - 1.01033337718e+6,
- 7.49414742893e+8*i - 2.25311577415e+6,0),
- 547052849318.0*i + 4.06181988045e+13
(-----------------------------------------,
i - 112.974983172
- 141265342333.0*i + 4.13350560163e+12
-----------------------------------------,
i - 112.974983172
845804392649.0*i - 9.62488227345e+13
--------------------------------------,
i - 112.974983172
876106032577.0*i - 2.66464796763e+13
--------------------------------------,
i - 112.974983172
1.47617976407e+12*i - 1.66771384004e+14
-----------------------------------------))
i - 169.834887206
tmp:=random_matrix(1,1,1000);
tmp := [ - 463.860434427*i + 1.35500571348]
triang_adjoint tmp;
[1]
% random matrix contains complex real values
on imaginary;
tmp:=random_matrix(5,5,1000);
tmp := mat((107.345792105*i - 1.98704739339,188.868545358*i + 1.22417796742,
- 630.485915434*i + 1.32195292724,
- 542.510039297*i - 1.94318764036,359.860945563*i - 1.69174206177),
( - 469.501213476*i - 1.17375946319, - 62.2197820375*i - 1.4051578261
, - 98.6604380996*i + 1.64691610034,
- 216.296595937*i + 1.56809020199,797.19877393*i - 1.31894550853),
(2.07054207792*i + 1.3891068942,393.038868455*i - 1.60894498437,
- 215.390393738*i - 1.00068640594,
- 195.674928032*i + 1.22123114986,211.921323796*i - 1.42499533273),
( - 750.357435524*i - 1.19871674827,
- 792.333836712*i - 1.63151974094, - 494.87049225*i + 1.99554801527
,638.173945822*i + 1.23793954612,111.418959978*i - 1.96273029328),
( - 255.359922267*i + 1.99035939892,
- 575.376389757*i - 1.03533681609,463.961589382*i - 1.86476410547,
83.8856338571*i + 1.10369785887, - 129.597812786*i - 1.4917934624))
triang_adjoint tmp;
mat((1,0,0,0,0),
(469.501213476*i + 1.17375946319,107.345792105*i - 1.98704739339,0,0,0),
(383.407897912*i + 1.84407237435e+5,1218.59364331*i + 41798.5118562,
769.235159465*i - 81990.7504399,0,0),
- 1.411092405e+10*i - 1.91497165215e+8
(-----------------------------------------,
i - 106.587367245
- 2.06157034475e+10*i + 1.09218575639e+8
-------------------------------------------,
i - 106.587367245
- 2.4008888901e+8*i + 13175.2571592,
- 1.02728261373e+8*i + 9.22309484944e+5,0),
- 203213290519.0*i - 3.07405185302e+12
(-----------------------------------------,
i - 184.764270765
311149245317.0*i + 2.05618234856e+13
--------------------------------------,
i - 184.764270765
212889617996.0*i - 4.13210409411e+13
--------------------------------------,
i - 184.764270765
- 7.79955805661e+10*i - 5.10418442965e+12
--------------------------------------------,
i - 184.764270765
7.62835257557e+10*i - 1.40944700076e+13
-----------------------------------------))
i - 106.587367245
tmp:=random_matrix(1,1,1000);
tmp := [276.278111177*i + 1.74724262616]
triang_adjoint tmp;
[1]
off imaginary;
% random matrix contains rounded real values
on rounded;
tmp:=random_matrix(5,5,1000);
tmp := mat(( - 293.224093687, - 99.023221037, - 819.400851656,796.020234589,
593.862087611),
( - 137.84203019,354.3234619, - 852.314261681, - 217.485901759,
256.139775139),
(324.37828726, - 56.5718498235, - 118.293003834,108.279501424,
23.2385400299),
( - 976.556156754,684.207160793,146.328625386,502.848132905,
312.766816689),
(211.783458501,166.556239469,175.715904944,251.57997022,280.123720131
))
triang_adjoint tmp;
mat((1,0,0,0,0),
(137.84203019, - 293.224093687,0,0,0),
( - 1.07136859076e+5, - 48709.2122316, - 1.17545737812e+5,0,0),
(1.27980020917e+8, - 1.64635380167e+8,4.76863677307e+8,1.43208428244e+8,0),
(5.82963241185e+10,3.9383738062e+10, - 437637051137.0, - 111757830528.0,
261327212376.0))
tmp:=random_matrix(1,1,1000);
tmp := [406.584701921]
triang_adjoint tmp;
[1]
off rounded;
% random matrix contains only integer values
on only_integer;
tmp:=random_matrix(7,7,1000);
[969 210 8 244 -887 -39 -916]
[ ]
[-774 296 -475 -694 -909 560 89 ]
[ ]
[-390 -559 -551 -567 241 -306 -655]
[ ]
tmp := [-478 809 181 -987 -144 929 -886]
[ ]
[188 267 -778 660 374 590 30 ]
[ ]
[ 73 971 -946 -43 -215 386 -365]
[ ]
[-792 -852 558 -797 343 219 110 ]
triang_adjoint tmp;
mat((1,0,0,0,0,0,0),
(774,969,0,0,0,0,0),
(548106,459771,449364,0,0,0,0),
(-108937808,399369604,-497500435,-461605941,0,0,0),
(-386678984240,-1001551613816,454549593485,637690866447,433944480084,0,0),
(-604165739229705,-320961967400919,-165015285307395,-1008712187269380,
-1670995725485274,1433408878792557,0),
(-181830640557070260,295390292387079435,816541226477288004,
850494616785589377,458176543109779557,-1709784109660828152,
-1475366833406131953))
tmp:=random_matrix(7,7,1);
[0 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 0 0]
[ ]
tmp := [0 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 0 0]
triang_adjoint tmp;
[1 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 0 0]
[ ]
[0 0 0 0 0 0 0]
% random matrix contains only complex integer
% values
on imaginary;
tmp:=random_matrix(5,5,1000);
tmp := mat((12*(38*i + 83),3*(153*i - 305),2*(141*i + 427), - 553*i + 617,
3*(83*i + 157)),
(164*i - 635, - 133*i + 991, - 373*i + 210,965*i - 608,2*(358*i - 55)
),
( - 230*i + 227,32*i + 339,2*(485*i - 219),707*i - 767, - 985*i - 51)
,
( - 609*i + 647, - 441*i + 187,930*i - 349,250*i - 211,274*i - 451),
( - 374*i - 135,793*i + 592, - 81*i - 1,89*i + 26,3*( - 40*i + 201)))
triang_adjoint tmp;
mat((1,0,0,0,0),
( - 164*i + 635,12*(38*i + 83),0,0,0),
(293397*i - 414880,9*(14243*i - 47243),3*(253651*i + 180645),0,0),
- 253324472288717*i + 71265413812547
(---------------------------------------,
253651*i + 180645
2*( - 220885726602145*i - 1441709355714)
------------------------------------------, - 1436348339*i + 1393250309,
253651*i + 180645
511458435*i - 1454012933,0),
13983048003979950612955437881*i - 71498490838832832842693585028
(-----------------------------------------------------------------,
65634686423804933*i - 9174596297286164
89295323223054915316808489269*i - 37624299403809895760446255007
-----------------------------------------------------------------,
65634686423804933*i - 9174596297286164
2*( - 71881165390656818494884812727*i - 25318671134083617432051412624)
------------------------------------------------------------------------,
65634686423804933*i - 9174596297286164
134577377248105484011524135103*i + 3495516738012600790097438251
-----------------------------------------------------------------,
65634686423804933*i - 9174596297286164
6*(65634686423804933*i - 9174596297286164)
--------------------------------------------))
253651*i + 180645
tmp:=random_matrix(5,5,2);
[i - 1 i i 0 - (i + 1)]
[ ]
[ 0 i -1 - i + 1 i + 1 ]
[ ]
tmp := [ -1 0 0 - i + 1 -1 ]
[ ]
[ -1 - i - i - i i + 1 ]
[ ]
[i - 1 0 i + 1 -1 0 ]
triang_adjoint tmp;
[ 1 0 0 0 0 ]
[ ]
[ 0 i - 1 0 0 0 ]
[ ]
[ - (i + 1) i + 1 ]
[------------ ------- - (i + 1) 0 0 ]
[ i - 1 i - 1 ]
[ ]
[ - (i + 1) 2*(2*i + 1) - 2*i ]
[------------ 0 ------------- -------- 0 ]
[ i i - 1 i - 1 ]
[ ]
[ 2*(3*i - 4) 2*(i + 2) 5*(3*i + 1) - 7*i + 1 2*(i + 2) ]
[------------- ----------- ------------- ------------ -----------]
[ 4*i + 3 i - 1 4*i + 3 4*i + 3 i - 1 ]
% Predicates.
matrixp(mat1);
t
matrixp(poly);
squarep(mat2);
t
squarep(mat3);
symmetricp(mat1);
t
symmetricp(mat3);
if not rounded_was_on then off rounded;
END;
Time for test: 200 ms