Artifact c71c8eb0ac46bf7469bcceab16b748a41ae375e28c5c94a8c554c5508bd0716e:
- Executable file
r38/packages/linalg/linalg.rlg
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 42605) [annotate] [blame] [check-ins using] [more...]
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