1 | // GMG/BM, last modified: 8.10.98 |
---|
2 | /////////////////////////////////////////////////////////////////////////////// |
---|
3 | version="$Id: matrix.lib,v 1.16 2000-12-23 17:17:28 greuel Exp $"; |
---|
4 | category="Linear Algebra"; |
---|
5 | info=" |
---|
6 | LIBRARY: matrix.lib Elementary Matrix Operations |
---|
7 | |
---|
8 | PROCEDURES: |
---|
9 | compress(A); matrix, zero columns from A deleted |
---|
10 | concat(A1,A2,..); matrix, concatenation of matrices A1,A2,... |
---|
11 | diag(p,n); matrix, nxn diagonal matrix with entries poly p |
---|
12 | dsum(A1,A2,..); matrix, direct sum of matrices A1,A2,... |
---|
13 | flatten(A); ideal, generated by entries of matrix A |
---|
14 | genericmat(n,m[,id]); generic nxm matrix [entries from id] |
---|
15 | is_complex(c); 1 if list c is a complex, 0 if not |
---|
16 | outer(A,B); matrix, outer product of matrices A and B |
---|
17 | power(A,n); matrix/intmat, n-th power of matrix/intmat A |
---|
18 | skewmat(n[,id]); generic skew-symmetric nxn matrix [entries from id] |
---|
19 | submat(A,r,c); submatrix of A with rows/cols specified by intvec r/c |
---|
20 | symmat(n[,id]); generic symmetric nxn matrix [entries from id] |
---|
21 | tensor(A,B); matrix, tensor product of matrices A nd B |
---|
22 | unitmat(n); unit square matrix of size n |
---|
23 | gauss_col(A); transform constant matrix A into col-reduced nf |
---|
24 | gauss_row(A); transform constant matrix A into row-reduced nf |
---|
25 | addcol(A,c1,p,c2); add p*(c1-th col) to c2-th column of matrix A, p poly |
---|
26 | addrow(A,r1,p,r2); add p*(r1-th row) to r2-th row of matrix A, p poly |
---|
27 | multcol(A,c,p); multiply c-th column of A with poly p |
---|
28 | multrow(A,r,p); multiply r-th row of A with poly p |
---|
29 | permcol(A,i,j); permute i-th and j-th columns |
---|
30 | permrow(A,i,j); permute i-th and j-th rows |
---|
31 | rowred(A[,any]); reduction of matrix A with elementary row-operations |
---|
32 | colred(A[,any]); reduction of matrix A with elementary col-operations |
---|
33 | rm_unitrow(A); remove unit rows and associated columns of A |
---|
34 | rm_unitcol(A); remove unit columns and associated rows of A |
---|
35 | (parameters in square brackets [] are optional) |
---|
36 | "; |
---|
37 | |
---|
38 | LIB "inout.lib"; |
---|
39 | LIB "ring.lib"; |
---|
40 | LIB "random.lib"; |
---|
41 | /////////////////////////////////////////////////////////////////////////////// |
---|
42 | |
---|
43 | proc compress (A) |
---|
44 | "USAGE: compress(A); A matrix/ideal/module/intmat/intvec |
---|
45 | RETURN: same type, zero columns/generators from A deleted |
---|
46 | (if A=intvec, zero elements are deleted) |
---|
47 | EXAMPLE: example compress; shows an example |
---|
48 | " |
---|
49 | { |
---|
50 | if( typeof(A)=="matrix" ) { return(matrix(simplify(A,2))); } |
---|
51 | if( typeof(A)=="intmat" or typeof(A)=="intvec" ) |
---|
52 | { |
---|
53 | ring r=0,x,lp; |
---|
54 | if( typeof(A)=="intvec" ) { intmat C=transpose(A); kill A; intmat A=C; } |
---|
55 | module m = matrix(A); |
---|
56 | intmat B[nrows(A)][size(m)]; |
---|
57 | int i,j; |
---|
58 | for( i=1; i<=ncols(A); i=i+1 ) |
---|
59 | { |
---|
60 | if( m[i]!=[0] ) |
---|
61 | { |
---|
62 | j=j+1; |
---|
63 | B[1..nrows(A),j]=A[1..nrows(A),i]; |
---|
64 | } |
---|
65 | } |
---|
66 | if( defined(C) ) { return(intvec(B)); } |
---|
67 | return(B); |
---|
68 | } |
---|
69 | return(simplify(A,2)); |
---|
70 | } |
---|
71 | example |
---|
72 | { "EXAMPLE:"; echo = 2; |
---|
73 | ring r=0,(x,y,z),ds; |
---|
74 | matrix A[3][4]=1,0,3,0,x,0,z,0,x2,0,z2,0; |
---|
75 | print(A); |
---|
76 | print(compress(A)); |
---|
77 | module m=module(A); show(m); |
---|
78 | show(compress(m)); |
---|
79 | intmat B[3][4]=1,0,3,0,4,0,5,0,6,0,7,0; |
---|
80 | compress(B); |
---|
81 | intvec C=0,0,1,2,0,3; |
---|
82 | compress(C); |
---|
83 | } |
---|
84 | /////////////////////////////////////////////////////////////////////////////// |
---|
85 | proc concat (list #) |
---|
86 | "USAGE: concat(A1,A2,..); A1,A2,... matrices |
---|
87 | RETURN: matrix, concatenation of A1,A2,... . Number of rows of result matrix |
---|
88 | is max(nrows(A1),nrows(A2),...) |
---|
89 | EXAMPLE: example concat; shows an example |
---|
90 | " |
---|
91 | { |
---|
92 | int i; |
---|
93 | module B=module(#[1]); |
---|
94 | for( i=2; i<=size(#); i=i+1 ) { B=B,module(#[i]); } |
---|
95 | return(matrix(B)); |
---|
96 | } |
---|
97 | example |
---|
98 | { "EXAMPLE:"; echo = 2; |
---|
99 | ring r=0,(x,y,z),ds; |
---|
100 | matrix A[3][3]=1,2,3,x,y,z,x2,y2,z2; |
---|
101 | matrix B[2][2]=1,0,2,0; matrix C[1][4]=4,5,x,y; |
---|
102 | print(A); |
---|
103 | print(B); |
---|
104 | print(C); |
---|
105 | print(concat(A,B,C)); |
---|
106 | } |
---|
107 | /////////////////////////////////////////////////////////////////////////////// |
---|
108 | |
---|
109 | proc diag (list #) |
---|
110 | "USAGE: diag(p,n); p poly, n integer |
---|
111 | diag(A); A matrix |
---|
112 | RETURN: diag(p,n): diagonal matrix, p times unitmatrix of size n |
---|
113 | diag(A) : n*m x n*m diagonal matrix with entries all the entries of |
---|
114 | the nxm matrix A, taken from the 1st row, 2nd row etc of A |
---|
115 | EXAMPLE: example diag; shows an example |
---|
116 | " |
---|
117 | { |
---|
118 | if( size(#)==2 ) { return(matrix(#[1]*freemodule(#[2]))); } |
---|
119 | if( size(#)==1 ) |
---|
120 | { |
---|
121 | int i; ideal id=#[1]; |
---|
122 | int n=ncols(id); matrix A[n][n]; |
---|
123 | for( i=1; i<=n; i=i+1 ) { A[i,i]=id[i]; } |
---|
124 | } |
---|
125 | return(A); |
---|
126 | } |
---|
127 | example |
---|
128 | { "EXAMPLE:"; echo = 2; |
---|
129 | ring r=0,(x,y,z),ds; |
---|
130 | print(diag(xy,4)); |
---|
131 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
132 | print(A); |
---|
133 | print(diag(A)); |
---|
134 | } |
---|
135 | /////////////////////////////////////////////////////////////////////////////// |
---|
136 | |
---|
137 | proc dsum (list #) |
---|
138 | "USAGE: dsum(A1,A2,..); A1,A2,... matrices |
---|
139 | RETURN: matrix, direct sum of A1,A2,... |
---|
140 | EXAMPLE: example dsum; shows an example |
---|
141 | " |
---|
142 | { |
---|
143 | int i,N,a; |
---|
144 | list L; |
---|
145 | for( i=1; i<=size(#); i=i+1 ) { N=N+nrows(#[i]); } |
---|
146 | for( i=1; i<=size(#); i=i+1 ) |
---|
147 | { |
---|
148 | matrix B[N][ncols(#[i])]; |
---|
149 | B[a+1..a+nrows(#[i]),1..ncols(#[i])]=#[i]; |
---|
150 | a=a+nrows(#[i]); |
---|
151 | L[i]=B; kill B; |
---|
152 | } |
---|
153 | return(concat(L)); |
---|
154 | } |
---|
155 | example |
---|
156 | { "EXAMPLE:"; echo = 2; |
---|
157 | ring r=0,(x,y,z),ds; |
---|
158 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
159 | matrix B[2][2]=1,x,y,z; |
---|
160 | matrix C[1][4]=4,5,x,y; |
---|
161 | print(A); |
---|
162 | print(B); |
---|
163 | print(C); |
---|
164 | print(dsum(A,B,C)); |
---|
165 | } |
---|
166 | /////////////////////////////////////////////////////////////////////////////// |
---|
167 | |
---|
168 | proc flatten (matrix A) |
---|
169 | "USAGE: flatten(A); A matrix |
---|
170 | RETURN: ideal, generated by all entries from A |
---|
171 | EXAMPLE: example flatten; shows an example |
---|
172 | " |
---|
173 | { |
---|
174 | return(ideal(A)); |
---|
175 | } |
---|
176 | example |
---|
177 | { "EXAMPLE:"; echo = 2; |
---|
178 | ring r=0,(x,y,z),ds; |
---|
179 | matrix A[3][3]=1,2,3,x,y,z,7,8,9; |
---|
180 | print(A); |
---|
181 | flatten(A); |
---|
182 | } |
---|
183 | /////////////////////////////////////////////////////////////////////////////// |
---|
184 | |
---|
185 | proc genericmat (int n,int m,list #) |
---|
186 | "USAGE: genericmat(n,m[,id]); n,m=integers, id=ideal |
---|
187 | RETURN: nxm matrix, with entries from id (default: id=maxideal(1)) |
---|
188 | NOTE: if id has less than nxm elements, the matrix is filled with 0's, |
---|
189 | genericmat(n,m); creates the generic nxm matrix |
---|
190 | EXAMPLE: example genericmat; shows an example |
---|
191 | " |
---|
192 | { |
---|
193 | if( size(#)==0 ) { ideal id=maxideal(1); } |
---|
194 | if( size(#)==1 ) { ideal id=#[1]; } |
---|
195 | if( size(#)>=2 ) { "// give 3 arguments, 3-rd argument must be an ideal"; } |
---|
196 | matrix B[n][m]=id; |
---|
197 | return(B); |
---|
198 | } |
---|
199 | example |
---|
200 | { "EXAMPLE:"; echo = 2; |
---|
201 | ring R=0,x(1..16),lp; |
---|
202 | print(genericmat(4,4)); // the generic 4x4 matrix |
---|
203 | changevar("R1",A_Z("a",4),R); |
---|
204 | matrix A=genericmat(4,5,maxideal(1)^3); |
---|
205 | print(A); |
---|
206 | int n,m=4,3; |
---|
207 | ideal i = ideal(randommat(1,n*m,maxideal(1),9)); |
---|
208 | print(genericmat(n,m,i)); // matrix of generic linear forms |
---|
209 | kill R1; |
---|
210 | } |
---|
211 | /////////////////////////////////////////////////////////////////////////////// |
---|
212 | |
---|
213 | proc is_complex (list c) |
---|
214 | "USAGE: is_complex(c); c = list of size-compatible modules or matrices |
---|
215 | RETURN: 1 if c[i]*c[i+1]=0 for all i, 0 if not, hence checking whether the |
---|
216 | list of matrices forms a complex |
---|
217 | NOTE: Ideals are treated internally as 1-line matrices |
---|
218 | EXAMPLE: example is_complex; shows an example |
---|
219 | " |
---|
220 | { |
---|
221 | int i; |
---|
222 | module @test; |
---|
223 | for( i=1; i<=size(c)-1; i=i+1 ) |
---|
224 | { |
---|
225 | c[i]=matrix(c[i]); c[i+1]=matrix(c[i+1]); |
---|
226 | @test=c[i]*c[i+1]; |
---|
227 | if (size(@test)!=0) |
---|
228 | { |
---|
229 | if( voice==2 ) { "// argument is not a complex at position",i; } |
---|
230 | return(0); |
---|
231 | } |
---|
232 | } |
---|
233 | if( voice==2 ) { "// argument is a complex"; } |
---|
234 | return(1); |
---|
235 | } |
---|
236 | example |
---|
237 | { "EXAMPLE:"; echo = 2; |
---|
238 | ring r=32003,(x,y,z),ds; |
---|
239 | ideal i=x4+y5+z6,xyz,yx2+xz2+zy7; |
---|
240 | list L=nres(i,0); |
---|
241 | is_complex(L); |
---|
242 | L[4]=matrix(i); |
---|
243 | is_complex(L); |
---|
244 | } |
---|
245 | /////////////////////////////////////////////////////////////////////////////// |
---|
246 | |
---|
247 | proc outer (matrix A, matrix B) |
---|
248 | "USAGE: outer(A,B); A,B matrices |
---|
249 | RETURN: matrix, outer product of A and B |
---|
250 | EXAMPLE: example outer; shows an example |
---|
251 | " |
---|
252 | { |
---|
253 | int i,j; list L; |
---|
254 | int triv = nrows(B)*ncols(B); |
---|
255 | if( triv==1 ) |
---|
256 | { |
---|
257 | return(B[1,1]*A); |
---|
258 | } |
---|
259 | else |
---|
260 | { |
---|
261 | int N = nrows(A)*nrows(B); |
---|
262 | matrix C[N][ncols(B)]; |
---|
263 | for( i=1; i<=ncols(A); i=i+1 ) |
---|
264 | { |
---|
265 | for( j=1; j<=nrows(A); j=j+1 ) |
---|
266 | { |
---|
267 | C[(j-1)*nrows(B)+1..j*nrows(B),1..ncols(B)]=A[j,i]*B; |
---|
268 | } |
---|
269 | L[i]=C; |
---|
270 | } |
---|
271 | return(concat(L)); |
---|
272 | } |
---|
273 | } |
---|
274 | example |
---|
275 | { "EXAMPLE:"; echo = 2; |
---|
276 | ring r=32003,(x,y,z),ds; |
---|
277 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
278 | matrix B[2][2]=x,y,0,z; |
---|
279 | print(A); |
---|
280 | print(B); |
---|
281 | print(outer(A,B)); |
---|
282 | } |
---|
283 | /////////////////////////////////////////////////////////////////////////////// |
---|
284 | |
---|
285 | proc power ( A, int n) |
---|
286 | "USAGE: power(A,n); A a square-matrix of type intmat or matrix, n=integer |
---|
287 | RETURN: intmat resp. matrix, the n-th power of A |
---|
288 | NOTE: for A=intmat and big n the result may be wrong because of int overflow |
---|
289 | EXAMPLE: example power; shows an example |
---|
290 | " |
---|
291 | { |
---|
292 | //---------------------------- type checking ---------------------------------- |
---|
293 | if( typeof(A)!="matrix" and typeof(A)!="intmat" ) |
---|
294 | { |
---|
295 | "// no matrix or intmat!"; |
---|
296 | return (A); |
---|
297 | } |
---|
298 | if( ncols(A) != nrows(A) ) |
---|
299 | { |
---|
300 | "// not a suare matrix!"; |
---|
301 | return(); |
---|
302 | } |
---|
303 | //---------------------------- trivial cases ---------------------------------- |
---|
304 | int ii; |
---|
305 | if( n <= 0 ) |
---|
306 | { |
---|
307 | if( typeof(A)=="matrix" ) |
---|
308 | { |
---|
309 | return (unitmat(nrows(A))); |
---|
310 | } |
---|
311 | if( typeof(A)=="intmat" ) |
---|
312 | { |
---|
313 | intmat B[nrows(A)][nrows(A)]; |
---|
314 | for( ii=1; ii<=nrows(A); ii++ ) |
---|
315 | { |
---|
316 | B[ii,ii] = 1; |
---|
317 | } |
---|
318 | return (B); |
---|
319 | } |
---|
320 | } |
---|
321 | if( n == 1 ) { return (A); } |
---|
322 | //---------------------------- sub procedure ---------------------------------- |
---|
323 | proc matpow (A, int n) |
---|
324 | { |
---|
325 | def B = A*A; |
---|
326 | int ii= 2; |
---|
327 | int jj= 4; |
---|
328 | while( jj <= n ) |
---|
329 | { |
---|
330 | B=B*B; |
---|
331 | ii=jj; |
---|
332 | jj=2*jj; |
---|
333 | } |
---|
334 | return(B,n-ii); |
---|
335 | } |
---|
336 | //----------------------------- main program ---------------------------------- |
---|
337 | list L = matpow(A,n); |
---|
338 | def B = L[1]; |
---|
339 | ii = L[2]; |
---|
340 | while( ii>=2 ) |
---|
341 | { |
---|
342 | L = matpow(A,ii); |
---|
343 | B = B*L[1]; |
---|
344 | ii= L[2]; |
---|
345 | } |
---|
346 | if( ii == 0) { return(B); } |
---|
347 | if( ii == 1) { return(A*B); } |
---|
348 | } |
---|
349 | example |
---|
350 | { "EXAMPLE:"; echo = 2; |
---|
351 | intmat A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
352 | print(power(A,3));""; |
---|
353 | ring r=0,(x,y,z),dp; |
---|
354 | matrix B[4][4]=0,x,y,z,0,0,y,z,0,0,0,z,x,y,z,0; |
---|
355 | print(power(B,3));""; |
---|
356 | matrix C[3][3]=1,2,3,4,5,6,7,8,9; |
---|
357 | power(C,50); |
---|
358 | } |
---|
359 | /////////////////////////////////////////////////////////////////////////////// |
---|
360 | |
---|
361 | proc skewmat (int n, list #) |
---|
362 | "USAGE: skewmat(n[,id]); n integer, id ideal |
---|
363 | RETURN: skew-symmetric nxn matrix, with entries from id |
---|
364 | (default: id=maxideal(1)) |
---|
365 | NOTE: if id has less than n*(n-1)/2 elements, the matrix is filled with 0's, |
---|
366 | skewmat(n); creates the generic skew-symmetric matrix |
---|
367 | EXAMPLE: example skewmat; shows an example |
---|
368 | " |
---|
369 | { |
---|
370 | matrix B[n][n]; |
---|
371 | if( size(#)==0 ) { ideal id=maxideal(1); } |
---|
372 | else { ideal id=#[1]; } |
---|
373 | id = id,B[1..n,1..n]; |
---|
374 | int i,j; |
---|
375 | for( i=0; i<=n-2; i=i+1 ) |
---|
376 | { |
---|
377 | B[i+1,i+2..n]=id[j+1..j+n-i-1]; |
---|
378 | j=j+n-i-1; |
---|
379 | } |
---|
380 | matrix A=transpose(B); |
---|
381 | B=B-A; |
---|
382 | return(B); |
---|
383 | } |
---|
384 | example |
---|
385 | { "EXAMPLE:"; echo = 2; |
---|
386 | ring R=0,x(1..5),lp; |
---|
387 | print(skewmat(4)); // the generic skew-symmetric matrix |
---|
388 | changevar("R1",A_Z("a",5),R); |
---|
389 | matrix A=skewmat(6,maxideal(1)^2); |
---|
390 | print(A); |
---|
391 | int n=4; |
---|
392 | ideal i = ideal(randommat(1,n*(n-1) div 2,maxideal(1),9)); |
---|
393 | print(skewmat(n,i)); // skew matrix of generic linear forms |
---|
394 | kill R1; |
---|
395 | } |
---|
396 | /////////////////////////////////////////////////////////////////////////////// |
---|
397 | |
---|
398 | proc submat (matrix A, intvec r, intvec c) |
---|
399 | "USAGE: submat(A,r,c); A=matrix, r,c=intvec |
---|
400 | RETURN: matrix, submatrix of A with rows specified by intvec r and columns |
---|
401 | specified by intvec c |
---|
402 | EXAMPLE: example submat; shows an example |
---|
403 | " |
---|
404 | { |
---|
405 | matrix B[size(r)][size(c)]=A[r,c]; |
---|
406 | return(B); |
---|
407 | } |
---|
408 | example |
---|
409 | { "EXAMPLE:"; echo = 2; |
---|
410 | ring R=32003,(x,y,z),lp; |
---|
411 | matrix A[4][4]=x,y,z,0,1,2,3,4,5,6,7,8,9,x2,y2,z2; |
---|
412 | print(A); |
---|
413 | intvec v=1,3,4; |
---|
414 | matrix B=submat(A,v,1..3); |
---|
415 | print(B); |
---|
416 | } |
---|
417 | /////////////////////////////////////////////////////////////////////////////// |
---|
418 | |
---|
419 | proc symmat (int n, list #) |
---|
420 | "USAGE: symmat(n[,id]); n integer, id ideal |
---|
421 | RETURN: symmetric nxn matrix, with entries from id (default: id=maxideal(1)) |
---|
422 | NOTE: if id has less than n*(n+1)/2 elements, the matrix is filled with 0's, |
---|
423 | symmat(n); creates the generic symmetric matrix |
---|
424 | EXAMPLE: example symmat; shows an example |
---|
425 | " |
---|
426 | { |
---|
427 | matrix B[n][n]; |
---|
428 | if( size(#)==0 ) { ideal id=maxideal(1); } |
---|
429 | else { ideal id=#[1]; } |
---|
430 | id = id,B[1..n,1..n]; |
---|
431 | int i,j; |
---|
432 | for( i=0; i<=n-1; i=i+1 ) |
---|
433 | { |
---|
434 | B[i+1,i+1..n]=id[j+1..j+n-i]; |
---|
435 | j=j+n-i; |
---|
436 | } |
---|
437 | matrix A=transpose(B); |
---|
438 | for( i=1; i<=n; i=i+1 ) { A[i,i]=0; } |
---|
439 | B=A+B; |
---|
440 | return(B); |
---|
441 | } |
---|
442 | example |
---|
443 | { "EXAMPLE:"; echo = 2; |
---|
444 | ring R=0,x(1..10),lp; |
---|
445 | print(symmat(4)); // the generic symmetric matrix |
---|
446 | changevar("R1",A_Z("a",5),R); |
---|
447 | matrix A=symmat(5,maxideal(1)^2); |
---|
448 | print(A); |
---|
449 | int n=3; |
---|
450 | ideal i = ideal(randommat(1,n*(n+1) div 2,maxideal(1),9)); |
---|
451 | print(symmat(n,i)); // symmetric matrix of generic linear forms |
---|
452 | kill R1; |
---|
453 | } |
---|
454 | /////////////////////////////////////////////////////////////////////////////// |
---|
455 | |
---|
456 | proc tensor (matrix A, matrix B) |
---|
457 | "USAGE: tensor(A,B); A,B matrices |
---|
458 | RETURN: matrix, tensor product of A and B |
---|
459 | EXAMPLE: example tensor; shows an example |
---|
460 | " |
---|
461 | { |
---|
462 | int i,j; |
---|
463 | matrix C,D; |
---|
464 | for( i=1; i<=nrows(A); i++ ) |
---|
465 | { |
---|
466 | C = A[i,1]*B; |
---|
467 | for( j=2; j<=ncols(A); j++ ) |
---|
468 | { |
---|
469 | C = concat(C,A[i,j]*B); |
---|
470 | } |
---|
471 | D = concat(D,transpose(C)); |
---|
472 | } |
---|
473 | D = transpose(D); |
---|
474 | return(submat(D,2..nrows(D),1..ncols(D))); |
---|
475 | } |
---|
476 | example |
---|
477 | { "EXAMPLE:"; echo = 2; |
---|
478 | ring r=32003,(x,y,z),(c,ds); |
---|
479 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
480 | matrix B[2][2]=x,y,0,z; |
---|
481 | print(A); |
---|
482 | print(B); |
---|
483 | print(tensor(A,B)); |
---|
484 | } |
---|
485 | /////////////////////////////////////////////////////////////////////////////// |
---|
486 | |
---|
487 | proc unitmat (int n) |
---|
488 | "USAGE: unitmat(n); n integer >= 0 |
---|
489 | RETURN: nxn unit matrix |
---|
490 | NOTE: needs a basering, diagonal entries are numbers (=1) in the basering |
---|
491 | EXAMPLE: example unitmat; shows an example |
---|
492 | " |
---|
493 | { |
---|
494 | return(matrix(freemodule(n))); |
---|
495 | } |
---|
496 | example |
---|
497 | { "EXAMPLE:"; echo = 2; |
---|
498 | ring r=32003,(x,y,z),lp; |
---|
499 | print(xyz*unitmat(4)); |
---|
500 | print(unitmat(5)); |
---|
501 | } |
---|
502 | /////////////////////////////////////////////////////////////////////////////// |
---|
503 | |
---|
504 | proc gauss_col (matrix A) |
---|
505 | "USAGE: gauss_col(A); A=matrix with constant coefficients |
---|
506 | RETURN: matrix = col-reduced lower-triagonal normal form of A |
---|
507 | NOTE: works fine for constant matrices, use proc colred for poly matrices |
---|
508 | EXAMPLE: example gauss_col; shows an example |
---|
509 | " |
---|
510 | { |
---|
511 | def R=basering; |
---|
512 | intvec v = option(get); |
---|
513 | changeord("@R","ds,c",R); |
---|
514 | option(redSB); option(nointStrategy); |
---|
515 | matrix A = imap(R,A); |
---|
516 | A = matrix(std(A),nrows(A),ncols(A)); |
---|
517 | setring R; |
---|
518 | A=imap(@R,A); |
---|
519 | option(set,v); |
---|
520 | kill @R; |
---|
521 | return(A); |
---|
522 | } |
---|
523 | example |
---|
524 | { "EXAMPLE:"; echo = 2; |
---|
525 | ring S=0,x,dp; |
---|
526 | matrix A[5][4] = 3, 1,1,-1, |
---|
527 | 13, 8,6,-7, |
---|
528 | 14,10,6,-7, |
---|
529 | 7, 4,3,-3, |
---|
530 | 2, 1,0, 3; |
---|
531 | print(gauss_col(A)); |
---|
532 | } |
---|
533 | /////////////////////////////////////////////////////////////////////////////// |
---|
534 | |
---|
535 | proc gauss_row (matrix A) |
---|
536 | "USAGE: gauss_row(A); A=matrix with constant coefficients |
---|
537 | RETURN: matrix = row-reduced upper-triangular normal form of A |
---|
538 | NOTE: may be used to solve a system of linear equations |
---|
539 | works fine for constant matrices, use proc rowred for poly matrices |
---|
540 | EXAMPLE: example gauss_row; shows an example |
---|
541 | " |
---|
542 | { |
---|
543 | A = gauss_col(transpose(A)); |
---|
544 | return(transpose(A)); |
---|
545 | } |
---|
546 | example |
---|
547 | { "EXAMPLE:"; echo = 2; |
---|
548 | ring S=0,x,dp; |
---|
549 | matrix A[4][5] = 3, 1,1,-1,2, |
---|
550 | 13, 8,6,-7,1, |
---|
551 | 14,10,6,-7,1, |
---|
552 | 7, 4,3,-3,3; |
---|
553 | print(gauss_row(A)); |
---|
554 | } |
---|
555 | /////////////////////////////////////////////////////////////////////////////// |
---|
556 | |
---|
557 | proc addcol (matrix A, int c1, poly p, int c2) |
---|
558 | "USAGE: addcol(A,c1,p,c2); A matrix, p poly, c1, c2 positive integers |
---|
559 | RETURN: matrix, A being modified by adding p times column c1 to column c2 |
---|
560 | EXAMPLE: example addcol; shows an example |
---|
561 | " |
---|
562 | { |
---|
563 | A[1..nrows(A),c2]=A[1..nrows(A),c2]+p*A[1..nrows(A),c1]; |
---|
564 | return(A); |
---|
565 | } |
---|
566 | example |
---|
567 | { "EXAMPLE:"; echo = 2; |
---|
568 | ring r=32003,(x,y,z),lp; |
---|
569 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
570 | print(A); |
---|
571 | print(addcol(A,1,xy,2)); |
---|
572 | } |
---|
573 | /////////////////////////////////////////////////////////////////////////////// |
---|
574 | |
---|
575 | proc addrow (matrix A, int r1, poly p, int r2) |
---|
576 | "USAGE: addcol(A,r1,p,r2); A matrix, p poly, r1, r2 positive integers |
---|
577 | RETURN: matrix, A being modified by adding p times row r1 to row r2 |
---|
578 | EXAMPLE: example addrow; shows an example |
---|
579 | " |
---|
580 | { |
---|
581 | A[r2,1..ncols(A)]=A[r2,1..ncols(A)]+p*A[r1,1..ncols(A)]; |
---|
582 | return(A); |
---|
583 | } |
---|
584 | example |
---|
585 | { "EXAMPLE:"; echo = 2; |
---|
586 | ring r=32003,(x,y,z),lp; |
---|
587 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
588 | print(A); |
---|
589 | print(addrow(A,1,xy,3)); |
---|
590 | } |
---|
591 | /////////////////////////////////////////////////////////////////////////////// |
---|
592 | |
---|
593 | proc multcol (matrix A, int c, poly p) |
---|
594 | "USAGE: addcol(A,c,p); A matrix, p poly, c positive integer |
---|
595 | RETURN: matrix, A being modified by multiplying column c with p |
---|
596 | EXAMPLE: example multcol; shows an example |
---|
597 | " |
---|
598 | { |
---|
599 | A[1..nrows(A),c]=p*A[1..nrows(A),c]; |
---|
600 | return(A); |
---|
601 | } |
---|
602 | example |
---|
603 | { "EXAMPLE:"; echo = 2; |
---|
604 | ring r=32003,(x,y,z),lp; |
---|
605 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
606 | print(A); |
---|
607 | print(multcol(A,2,xy)); |
---|
608 | } |
---|
609 | /////////////////////////////////////////////////////////////////////////////// |
---|
610 | |
---|
611 | proc multrow (matrix A, int r, poly p) |
---|
612 | "USAGE: multrow(A,r,p); A matrix, p poly, r positive integer |
---|
613 | RETURN: matrix, A being modified by multiplying row r with p |
---|
614 | EXAMPLE: example multrow; shows an example |
---|
615 | " |
---|
616 | { |
---|
617 | A[r,1..ncols(A)]=p*A[r,1..ncols(A)]; |
---|
618 | return(A); |
---|
619 | } |
---|
620 | example |
---|
621 | { "EXAMPLE:"; echo = 2; |
---|
622 | ring r=32003,(x,y,z),lp; |
---|
623 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
624 | print(A); |
---|
625 | print(multrow(A,2,xy)); |
---|
626 | } |
---|
627 | /////////////////////////////////////////////////////////////////////////////// |
---|
628 | |
---|
629 | proc permcol (matrix A, int c1, int c2) |
---|
630 | "USAGE: permcol(A,c1,c2); A matrix, c1,c2 positive integers |
---|
631 | RETURN: matrix, A being modified by permuting column c1 and c2 |
---|
632 | EXAMPLE: example permcol; shows an example |
---|
633 | " |
---|
634 | { |
---|
635 | matrix B=A; |
---|
636 | B[1..nrows(B),c1]=A[1..nrows(A),c2]; |
---|
637 | B[1..nrows(B),c2]=A[1..nrows(A),c1]; |
---|
638 | return(B); |
---|
639 | } |
---|
640 | example |
---|
641 | { "EXAMPLE:"; echo = 2; |
---|
642 | ring r=32003,(x,y,z),lp; |
---|
643 | matrix A[3][3]=1,x,3,4,y,6,7,z,9; |
---|
644 | print(A); |
---|
645 | print(permcol(A,2,3)); |
---|
646 | } |
---|
647 | /////////////////////////////////////////////////////////////////////////////// |
---|
648 | |
---|
649 | proc permrow (matrix A, int r1, int r2) |
---|
650 | "USAGE: permrow(A,r1,r2); A matrix, r1,r2 positive integers |
---|
651 | RETURN: matrix, A being modified by permuting row r1 and r2 |
---|
652 | EXAMPLE: example permrow; shows an example |
---|
653 | " |
---|
654 | { |
---|
655 | matrix B=A; |
---|
656 | B[r1,1..ncols(B)]=A[r2,1..ncols(A)]; |
---|
657 | B[r2,1..ncols(B)]=A[r1,1..ncols(A)]; |
---|
658 | return(B); |
---|
659 | } |
---|
660 | example |
---|
661 | { "EXAMPLE:"; echo = 2; |
---|
662 | ring r=32003,(x,y,z),lp; |
---|
663 | matrix A[3][3]=1,2,3,x,y,z,7,8,9; |
---|
664 | print(A); |
---|
665 | print(permrow(A,2,1)); |
---|
666 | } |
---|
667 | /////////////////////////////////////////////////////////////////////////////// |
---|
668 | |
---|
669 | proc rowred (matrix A,list #) |
---|
670 | "USAGE: rowred(A[,any]); A matrix |
---|
671 | RETURN: matrix B, being the row reduced form of A, |
---|
672 | if any defined: a list of two matrices B,C, such that B = C * A |
---|
673 | EXAMPLE: example rowred; shows an example |
---|
674 | " |
---|
675 | { |
---|
676 | int m,n=nrows(A),ncols(A); |
---|
677 | int i,j,k,l,rk; |
---|
678 | poly p; |
---|
679 | matrix d[m][n]; |
---|
680 | for (i=1;i<=m;i=i+1) |
---|
681 | { for (j=1;j<=n;j=j+1) |
---|
682 | { p = A[i,j]; |
---|
683 | if (ord(p)==0) |
---|
684 | { if (deg(p)==0) { d[i,j]=p; } |
---|
685 | } |
---|
686 | } |
---|
687 | } |
---|
688 | matrix b = A; |
---|
689 | if (size(#)) { b = concat(b,unitmat(m)); } |
---|
690 | for (l=1;l<=n;l=l+1) |
---|
691 | { |
---|
692 | k = findfirst(ideal(d[l]),rk+1); |
---|
693 | if (k) |
---|
694 | { rk = rk+1; |
---|
695 | b = permrow(b,rk,k); |
---|
696 | p = b[rk,l]; p = 1/p; |
---|
697 | b = multrow(b,rk,p); |
---|
698 | for (i=1;i<=m;i=i+1) |
---|
699 | { |
---|
700 | if (rk-i) { b = addrow(b,rk,-b[i,l],i);} |
---|
701 | } |
---|
702 | d = 0; |
---|
703 | for (i=rk+1;i<=m;i=i+1) |
---|
704 | { for (j=l+1;j<=n;j=j+1) |
---|
705 | { p = b[i,j]; |
---|
706 | if (ord(p)==0) |
---|
707 | { if (deg(p)==0) { d[i,j]=p; } |
---|
708 | } |
---|
709 | } |
---|
710 | } |
---|
711 | |
---|
712 | } |
---|
713 | } |
---|
714 | d = submat(b,1..m,1..n); |
---|
715 | if (size(#)) {return(d,submat(b,1..m,n+1..n+m));} |
---|
716 | return(d); |
---|
717 | } |
---|
718 | example |
---|
719 | { "EXAMPLE:"; echo = 2; |
---|
720 | ring r=0,(A,B,C),dp; |
---|
721 | matrix m[12][14]= |
---|
722 | AC,BC,-3BC,0, B2-A2,-3AC,B2,B2,0, 0, -C2,0, 0, 0, |
---|
723 | 2C,0, 0, B, -A, -4C, 2A,0, 0, 0, 0, 0, 0, 0, |
---|
724 | 0, 2C,-4C, -A,B, 0, B, 3B,AB,B2,0, 0, 0, 0, |
---|
725 | 0, 0, A, 0, 0, B, 0, 0, 0, 0, 0, 0, 0, -C2, |
---|
726 | 0, 0, 0, 0, 0, 0, 2, 0, 2A,0, 0, 0, 0, 0, |
---|
727 | 0, 0, 0, 0, 0, 0, 0, 2, A, 3B,0, B2,0, 0, |
---|
728 | 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, |
---|
729 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 3B,B2,0, |
---|
730 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, |
---|
731 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3B,0, |
---|
732 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, |
---|
733 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1; |
---|
734 | print(rowred(m)); |
---|
735 | list L=rowred(m,1); |
---|
736 | print(L[1]);"---------";print(L[2]); |
---|
737 | } |
---|
738 | /////////////////////////////////////////////////////////////////////////////// |
---|
739 | |
---|
740 | proc colred (matrix A,list #) |
---|
741 | "USAGE: colred(A[,any]); A matrix |
---|
742 | RETURN: matrix B, being the column reduced form of A, |
---|
743 | if any defined: a list of two matrices B,C, such that B = A * C |
---|
744 | EXAMPLE: example colred; shows an example |
---|
745 | " |
---|
746 | { |
---|
747 | A = transpose(A); |
---|
748 | if (size(#)) |
---|
749 | { list L = rowred(A,1); return(transpose(L[1]),transpose(L[2]));} |
---|
750 | else |
---|
751 | { return(transpose(rowred(A)));} |
---|
752 | } |
---|
753 | example |
---|
754 | { "EXAMPLE:"; echo = 2; |
---|
755 | ring r=0,(A,B,C),dp; |
---|
756 | matrix m[14][12]= |
---|
757 | AC, 2C, 0, 0, 0, 0, 0,0, 0,0, 0,0, |
---|
758 | BC, 0, 2C, 0, 0, 0, 0,0, 0,0, 0,0, |
---|
759 | -3BC, 0, -4C,A, 0, 0, 0,0, 0,0, 0,0, |
---|
760 | 0, B, -A, 0, 0, 0, 0,0, 0,0, 0,0, |
---|
761 | -A2+B2,-A, B, 0, 0, 0, 0,0, 0,0, 0,0, |
---|
762 | -3AC, -4C,0, B, 0, 0, 0,0, 0,0, 0,0, |
---|
763 | B2, 2A, B, 0, 2, 0, 0,0, 0,0, 0,0, |
---|
764 | B2, 0, 3B, 0, 0, 2, 0,0, 0,0, 0,0, |
---|
765 | 0, 0, AB, 0, 2A,A, 2,0, 0,0, 0,0, |
---|
766 | 0, 0, B2, 0, 0, 3B,0,2, 0,0, 0,0, |
---|
767 | -C2, 0, 0, 0, 0, 0, 0,0, 1,0, 0,0, |
---|
768 | 0, 0, 0, 0, 0, B2,0,3B,0,2, 0,0, |
---|
769 | 0, 0, 0, 0, 0, 0, 0,B2,0,3B,2,0, |
---|
770 | 0, 0, 0, -C2,0, 0, 0,0, 0,0, 0,1; |
---|
771 | print(colred(m)); |
---|
772 | list L=colred(m,1); |
---|
773 | print(L[1]);"---------";print(L[2]); |
---|
774 | } |
---|
775 | ////////////////////////////////////////////////////////////////////////////// |
---|
776 | |
---|
777 | static proc findfirst (ideal i,int t) |
---|
778 | { |
---|
779 | int n,k; |
---|
780 | for (n=t;n<=ncols(i);n=n+1) |
---|
781 | { |
---|
782 | if (i[n]!=0) { k=n;break;} |
---|
783 | } |
---|
784 | return(k); |
---|
785 | } |
---|
786 | ////////////////////////////////////////////////////////////////////////////// |
---|
787 | |
---|
788 | proc rm_unitcol(matrix A) |
---|
789 | "USAGE: rm_unitcol(A); A matrix (being row reduced) |
---|
790 | RETURN: matrix, obtained from A by deleting unit columns (having just one 1 |
---|
791 | and else 0 as entries) and associated rows |
---|
792 | EXAMPLE: example rm_unitcol; shows an example |
---|
793 | " |
---|
794 | { |
---|
795 | int l,j; |
---|
796 | intvec v; |
---|
797 | for (j=1;j<=ncols(A);j=j+1) |
---|
798 | { |
---|
799 | if (gen(l+1)==module(A)[j]) {l=l+1;} |
---|
800 | else { v=v,j;} |
---|
801 | } |
---|
802 | if (size(v)>1) |
---|
803 | { v = v[2..size(v)]; |
---|
804 | return(submat(A,l+1..nrows(A),v)); |
---|
805 | } |
---|
806 | else |
---|
807 | { return(0);} |
---|
808 | } |
---|
809 | example |
---|
810 | { "EXAMPLE:"; echo = 2; |
---|
811 | ring r=0,(A,B,C),dp; |
---|
812 | matrix m[14][12]= |
---|
813 | 1/2B2, A, 1/2B, 0, 1,0,0,0,0,0,0,0, |
---|
814 | 1/2B2, 0, 3/2B, 0, 0,1,0,0,0,0,0,0, |
---|
815 | -3/4AB2, -A2,-3/4AB, 0, 0,0,1,0,0,0,0,0, |
---|
816 | -3/4B3, 0, -7/4B2, 0, 0,0,0,1,0,0,0,0, |
---|
817 | -C2, 0, 0, 0, 0,0,0,0,1,0,0,0, |
---|
818 | 7/8B4, 0, 15/8B3, 0, 0,0,0,0,0,1,0,0, |
---|
819 | -15/16B5,0, -31/16B4,0, 0,0,0,0,0,0,1,0, |
---|
820 | 0, 0, 0, -C2,0,0,0,0,0,0,0,1, |
---|
821 | -3BC, 0, -4C, A, 0,0,0,0,0,0,0,0, |
---|
822 | 0, B, -A, 0, 0,0,0,0,0,0,0,0, |
---|
823 | -A2+B2, -A, B, 0, 0,0,0,0,0,0,0,0, |
---|
824 | -3AC, -4C,0, B, 0,0,0,0,0,0,0,0, |
---|
825 | AC, 2C, 0, 0, 0,0,0,0,0,0,0,0, |
---|
826 | BC, 0, 2C, 0, 0,0,0,0,0,0,0,0; |
---|
827 | print(rm_unitcol(m)); |
---|
828 | } |
---|
829 | ////////////////////////////////////////////////////////////////////////////// |
---|
830 | |
---|
831 | proc rm_unitrow (matrix A) |
---|
832 | "USAGE: rm_unitrow(A); A matrix (being colreduced) |
---|
833 | RETURN: matrix, obtained from A by deleting unit rows (having just one 1 |
---|
834 | and else 0 as entries) and associated columns |
---|
835 | EXAMPLE: example rm_unitrow; shows an example |
---|
836 | " |
---|
837 | { |
---|
838 | int l,j; |
---|
839 | intvec v; |
---|
840 | module M = transpose(A); |
---|
841 | for (j=1; j <= nrows(A); j=j+1) |
---|
842 | { |
---|
843 | if (gen(l+1) == M[j]) { l=l+1; } |
---|
844 | else { v=v,j; } |
---|
845 | } |
---|
846 | if (size(v) > 1) |
---|
847 | { v = v[2..size(v)]; |
---|
848 | return(submat(A,v,l+1..ncols(A))); |
---|
849 | } |
---|
850 | else |
---|
851 | { return(0);} |
---|
852 | } |
---|
853 | example |
---|
854 | { "EXAMPLE:"; echo = 2; |
---|
855 | ring r=0,(A,B,C),dp; |
---|
856 | matrix m[12][14]= |
---|
857 | 1/2B2,1/2B2,-3/4AB2,-3/4B3,-C2,7/8B4, -15/16B5,0, -3BC,0, B2-A2,-3AC,AC,BC, |
---|
858 | A, 0, -A2, 0, 0, 0, 0, 0, 0, B, -A, -4C, 2C,0, |
---|
859 | 1/2B, 3/2B, -3/4AB, -7/4B2,0, 15/8B3,-31/16B4,0, -4C, -A,B, 0, 0, 2C, |
---|
860 | 0, 0, 0, 0, 0, 0, 0, -C2,A, 0, 0, B, 0, 0, |
---|
861 | 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, |
---|
862 | 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, |
---|
863 | 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, |
---|
864 | 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, |
---|
865 | 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, |
---|
866 | 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, |
---|
867 | 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, |
---|
868 | 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0; |
---|
869 | print(rm_unitrow(m)); |
---|
870 | } |
---|
871 | ////////////////////////////////////////////////////////////////////////////// |
---|