1 | ////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="version matrix.lib 4.1.2.0 Feb_2019 "; // $Id$ |
---|
3 | category="Linear Algebra"; |
---|
4 | info=" |
---|
5 | LIBRARY: matrix.lib Elementary Matrix Operations |
---|
6 | |
---|
7 | PROCEDURES: |
---|
8 | compress(A); matrix, zero columns from A deleted |
---|
9 | concat(A1,A2,..); matrix, concatenation of matrices A1,A2,... |
---|
10 | diag(p,n); matrix, nxn diagonal matrix with entries poly p |
---|
11 | dsum(A1,A2,..); matrix, direct sum of matrices A1,A2,... |
---|
12 | flatten(A); ideal, generated by entries of matrix A |
---|
13 | genericmat(n,m[,id]); generic nxm matrix [entries from id] |
---|
14 | is_complex(c); 1 if list c is a complex, 0 if not |
---|
15 | outer(A,B); matrix, outer product of matrices A and B |
---|
16 | power(A,n); matrix/intmat, n-th power of matrix/intmat A |
---|
17 | skewmat(n[,id]); generic skew-symmetric nxn matrix [entries from id] |
---|
18 | submat(A,r,c); submatrix of A with rows/cols specified by intvec r/c |
---|
19 | symmat(n[,id]); generic symmetric nxn matrix [entries from id] |
---|
20 | unitmat(n); unit square matrix of size n |
---|
21 | gauss_col(A); transform a matrix into col-reduced Gauss normal form |
---|
22 | gauss_row(A); transform a matrix into row-reduced Gauss normal form |
---|
23 | addcol(A,c1,p,c2); add p*(c1-th col) to c2-th column of matrix A, p poly |
---|
24 | addrow(A,r1,p,r2); add p*(r1-th row) to r2-th row of matrix A, p poly |
---|
25 | multcol(A,c,p); multiply c-th column of A with poly p |
---|
26 | multrow(A,r,p); multiply r-th row of A with poly p |
---|
27 | permcol(A,i,j); permute i-th and j-th columns |
---|
28 | permrow(A,i,j); permute i-th and j-th rows |
---|
29 | rowred(A[,any]); reduction of matrix A with elementary row-operations |
---|
30 | colred(A[,any]); reduction of matrix A with elementary col-operations |
---|
31 | linear_relations(E); find linear relations between homogeneous vectors |
---|
32 | rm_unitrow(A); remove unit rows and associated columns of A |
---|
33 | rm_unitcol(A); remove unit columns and associated rows of A |
---|
34 | headStand(A); A[n-i+1,m-j+1]:=A[i,j] |
---|
35 | symmetricBasis(n,k[,s]); basis of k-th symmetric power of n-dim v.space |
---|
36 | exteriorBasis(n,k[,s]); basis of k-th exterior power of n-dim v.space |
---|
37 | symmetricPower(A,k); k-th symmetric power of a module/matrix A |
---|
38 | exteriorPower(A,k); k-th exterior power of a module/matrix A |
---|
39 | (parameters in square brackets [] are optional) |
---|
40 | "; |
---|
41 | |
---|
42 | LIB "inout.lib"; |
---|
43 | LIB "ring.lib"; |
---|
44 | LIB "random.lib"; |
---|
45 | LIB "general.lib"; // for sort |
---|
46 | LIB "nctools.lib"; // for superCommutative |
---|
47 | |
---|
48 | /////////////////////////////////////////////////////////////////////////////// |
---|
49 | |
---|
50 | proc compress (def A) |
---|
51 | "USAGE: compress(A); A matrix/ideal/module/intmat/intvec |
---|
52 | RETURN: same type, zero columns/generators from A deleted |
---|
53 | (if A=intvec, zero elements are deleted) |
---|
54 | EXAMPLE: example compress; shows an example |
---|
55 | " |
---|
56 | { |
---|
57 | if( typeof(A)=="matrix" ) { return(matrix(simplify(A,2))); } |
---|
58 | if( typeof(A)=="intmat" or typeof(A)=="intvec" ) |
---|
59 | { |
---|
60 | ring r=0,x,lp; |
---|
61 | if( typeof(A)=="intvec" ) { intmat C=transpose(A); kill A; intmat A=C; } |
---|
62 | module m = matrix(A); |
---|
63 | if ( size(m) == 0) |
---|
64 | { intmat B; } |
---|
65 | else |
---|
66 | { intmat B[nrows(A)][size(m)]; } |
---|
67 | int i,j; |
---|
68 | for( i=1; i<=ncols(A); i++ ) |
---|
69 | { |
---|
70 | if( m[i]!=[0] ) |
---|
71 | { |
---|
72 | j++; |
---|
73 | B[1..nrows(A),j]=A[1..nrows(A),i]; |
---|
74 | } |
---|
75 | } |
---|
76 | if( defined(C) ) { return(intvec(B)); } |
---|
77 | return(B); |
---|
78 | } |
---|
79 | return(simplify(A,2)); |
---|
80 | } |
---|
81 | example |
---|
82 | { "EXAMPLE:"; echo = 2; |
---|
83 | ring r=0,(x,y,z),ds; |
---|
84 | matrix A[3][4]=1,0,3,0,x,0,z,0,x2,0,z2,0; |
---|
85 | print(A); |
---|
86 | print(compress(A)); |
---|
87 | module m=module(A); show(m); |
---|
88 | show(compress(m)); |
---|
89 | intmat B[3][4]=1,0,3,0,4,0,5,0,6,0,7,0; |
---|
90 | compress(B); |
---|
91 | intvec C=0,0,1,2,0,3; |
---|
92 | compress(C); |
---|
93 | } |
---|
94 | /////////////////////////////////////////////////////////////////////////////// |
---|
95 | proc concat (list #) |
---|
96 | "USAGE: concat(A1,A2,..); A1,A2,... matrices |
---|
97 | RETURN: matrix, concatenation of A1,A2,.... Number of rows of result matrix |
---|
98 | is max(nrows(A1),nrows(A2),...) |
---|
99 | EXAMPLE: example concat; shows an example |
---|
100 | " |
---|
101 | { |
---|
102 | int i; |
---|
103 | for( i=size(#);i>0; i-- ) { #[i]=module(#[i]); } |
---|
104 | module B=#[1..size(#)]; |
---|
105 | return(matrix(B)); |
---|
106 | } |
---|
107 | example |
---|
108 | { "EXAMPLE:"; echo = 2; |
---|
109 | ring r=0,(x,y,z),ds; |
---|
110 | matrix A[3][3]=1,2,3,x,y,z,x2,y2,z2; |
---|
111 | matrix B[2][2]=1,0,2,0; matrix C[1][4]=4,5,x,y; |
---|
112 | print(A); |
---|
113 | print(B); |
---|
114 | print(C); |
---|
115 | print(concat(A,B,C)); |
---|
116 | } |
---|
117 | /////////////////////////////////////////////////////////////////////////////// |
---|
118 | |
---|
119 | proc diag (list #) |
---|
120 | "USAGE: diag(p,n); p poly, n integer |
---|
121 | diag(A); A matrix |
---|
122 | RETURN: diag(p,n): diagonal matrix, p times unit matrix of size n. |
---|
123 | @* diag(A) : n*m x n*m diagonal matrix with entries all the entries of |
---|
124 | the nxm matrix A, taken from the 1st row, 2nd row etc of A |
---|
125 | EXAMPLE: example diag; shows an example |
---|
126 | " |
---|
127 | { |
---|
128 | if( size(#)==2 ) { return(matrix(#[1]*freemodule(#[2]))); } |
---|
129 | if( size(#)==1 ) |
---|
130 | { |
---|
131 | int i; ideal id=#[1]; |
---|
132 | int n=ncols(id); matrix A[n][n]; |
---|
133 | for( i=1; i<=n; i++ ) { A[i,i]=id[i]; } |
---|
134 | } |
---|
135 | return(A); |
---|
136 | } |
---|
137 | example |
---|
138 | { "EXAMPLE:"; echo = 2; |
---|
139 | ring r = 0,(x,y,z),ds; |
---|
140 | print(diag(xy,4)); |
---|
141 | matrix A[3][2] = 1,2,3,4,5,6; |
---|
142 | print(A); |
---|
143 | print(diag(A)); |
---|
144 | } |
---|
145 | /////////////////////////////////////////////////////////////////////////////// |
---|
146 | |
---|
147 | proc dsum (list #) |
---|
148 | "USAGE: dsum(A1,A2,..); A1,A2,... matrices |
---|
149 | RETURN: matrix, direct sum of A1,A2,... |
---|
150 | EXAMPLE: example dsum; shows an example |
---|
151 | " |
---|
152 | { |
---|
153 | int i,N,a; |
---|
154 | list L; |
---|
155 | for( i=1; i<=size(#); i++ ) { N=N+nrows(#[i]); } |
---|
156 | for( i=1; i<=size(#); i++ ) |
---|
157 | { |
---|
158 | matrix B[N][ncols(#[i])]; |
---|
159 | B[a+1..a+nrows(#[i]),1..ncols(#[i])]=#[i]; |
---|
160 | a=a+nrows(#[i]); |
---|
161 | L[i]=B; kill B; |
---|
162 | } |
---|
163 | return(concat(L)); |
---|
164 | } |
---|
165 | example |
---|
166 | { "EXAMPLE:"; echo = 2; |
---|
167 | ring r = 0,(x,y,z),ds; |
---|
168 | matrix A[3][3] = 1,2,3,4,5,6,7,8,9; |
---|
169 | matrix B[2][2] = 1,x,y,z; |
---|
170 | print(A); |
---|
171 | print(B); |
---|
172 | print(dsum(A,B)); |
---|
173 | } |
---|
174 | /////////////////////////////////////////////////////////////////////////////// |
---|
175 | |
---|
176 | proc flatten (def A) |
---|
177 | "USAGE: flatten(A); A matrix/smatrix |
---|
178 | RETURN: ideal, generated by all entries from A resp. all colums of A appended |
---|
179 | EXAMPLE: example flatten; shows an example |
---|
180 | " |
---|
181 | { |
---|
182 | if (typeof(A)=="smatrix") {return(system("flatten",A)); } |
---|
183 | return(ideal(A)); |
---|
184 | } |
---|
185 | example |
---|
186 | { "EXAMPLE:"; echo = 2; |
---|
187 | ring r = 0,(x,y,z),ds; |
---|
188 | matrix A[2][3] = 1,2,x,y,z,7; |
---|
189 | print(A); |
---|
190 | flatten(A); |
---|
191 | flatten(smatrix(A)); |
---|
192 | } |
---|
193 | /////////////////////////////////////////////////////////////////////////////// |
---|
194 | |
---|
195 | proc genericmat (int n,int m,list #) |
---|
196 | "USAGE: genericmat(n,m[,id]); n,m=integers, id=ideal |
---|
197 | RETURN: nxm matrix, with entries from id. |
---|
198 | NOTE: if id has less than nxm elements, the matrix is filled with 0's, |
---|
199 | (default: id=maxideal(1)). |
---|
200 | genericmat(n,m); creates the generic nxm matrix |
---|
201 | EXAMPLE: example genericmat; shows an example |
---|
202 | " |
---|
203 | { |
---|
204 | if( size(#)==0 ) { ideal id=maxideal(1); } |
---|
205 | if( size(#)==1 ) { ideal id=#[1]; } |
---|
206 | if( size(#)>=2 ) { "// give 3 arguments, 3-rd argument must be an ideal"; } |
---|
207 | matrix B[n][m]=id; |
---|
208 | return(B); |
---|
209 | } |
---|
210 | example |
---|
211 | { "EXAMPLE:"; echo = 2; |
---|
212 | ring R = 0,x(1..16),lp; |
---|
213 | print(genericmat(3,3)); // the generic 3x3 matrix |
---|
214 | ring R1 = 0,(a,b,c,d),dp; |
---|
215 | matrix A = genericmat(3,4,maxideal(1)^3); |
---|
216 | print(A); |
---|
217 | int n,m = 3,2; |
---|
218 | ideal i = ideal(randommat(1,n*m,maxideal(1),9)); |
---|
219 | print(genericmat(n,m,i)); // matrix of generic linear forms |
---|
220 | } |
---|
221 | /////////////////////////////////////////////////////////////////////////////// |
---|
222 | |
---|
223 | proc is_complex (list c) |
---|
224 | "USAGE: is_complex(c); c = list of size-compatible modules or matrices |
---|
225 | RETURN: 1 if c[i]*c[i+1]=0 for all i, 0 if not, hence checking whether the |
---|
226 | list of matrices forms a complex. |
---|
227 | NOTE: Ideals are treated internally as 1-line matrices. |
---|
228 | If printlevel > 0, the position where c is not a complex is shown. |
---|
229 | EXAMPLE: example is_complex; shows an example |
---|
230 | " |
---|
231 | { |
---|
232 | int i; |
---|
233 | module @test; |
---|
234 | for( i=1; i<=size(c)-1; i++ ) |
---|
235 | { |
---|
236 | c[i]=matrix(c[i]); c[i+1]=matrix(c[i+1]); |
---|
237 | @test=c[i]*c[i+1]; |
---|
238 | if (size(@test)!=0) |
---|
239 | { |
---|
240 | dbprint(printlevel-voice+2,"// not a complex at position " +string(i)); |
---|
241 | return(0); |
---|
242 | } |
---|
243 | } |
---|
244 | return(1); |
---|
245 | } |
---|
246 | example |
---|
247 | { "EXAMPLE:"; echo = 2; |
---|
248 | ring r = 32003,(x,y,z),ds; |
---|
249 | ideal i = x4+y5+z6,xyz,yx2+xz2+zy7; |
---|
250 | list L = nres(i,0); |
---|
251 | is_complex(L); |
---|
252 | L[4] = matrix(i); |
---|
253 | is_complex(L); |
---|
254 | } |
---|
255 | /////////////////////////////////////////////////////////////////////////////// |
---|
256 | |
---|
257 | proc outer (matrix A, matrix B) |
---|
258 | "USAGE: outer(A,B); A,B matrices |
---|
259 | RETURN: matrix, outer (tensor) product of A and B |
---|
260 | EXAMPLE: example outer; shows an example |
---|
261 | " |
---|
262 | { |
---|
263 | int i,j; list L; |
---|
264 | int triv = nrows(B)*ncols(B); |
---|
265 | if( triv==1 ) |
---|
266 | { |
---|
267 | return(B[1,1]*A); |
---|
268 | } |
---|
269 | else |
---|
270 | { |
---|
271 | int N = nrows(A)*nrows(B); |
---|
272 | matrix C[N][ncols(B)]; |
---|
273 | for( i=ncols(A);i>0; i-- ) |
---|
274 | { |
---|
275 | for( j=1; j<=nrows(A); j++ ) |
---|
276 | { |
---|
277 | C[(j-1)*nrows(B)+1..j*nrows(B),1..ncols(B)]=A[j,i]*B; |
---|
278 | } |
---|
279 | L[i]=C; |
---|
280 | } |
---|
281 | return(concat(L)); |
---|
282 | } |
---|
283 | } |
---|
284 | example |
---|
285 | { "EXAMPLE:"; echo = 2; |
---|
286 | ring r=32003,(x,y,z),ds; |
---|
287 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
288 | matrix B[2][2]=x,y,0,z; |
---|
289 | print(A); |
---|
290 | print(B); |
---|
291 | print(outer(A,B)); |
---|
292 | } |
---|
293 | /////////////////////////////////////////////////////////////////////////////// |
---|
294 | |
---|
295 | proc power (def A, int n) |
---|
296 | "USAGE: power(A,n); A a square-matrix of type intmat or matrix, n=integer>=0 |
---|
297 | RETURN: intmat resp. matrix, the n-th power of A |
---|
298 | NOTE: for A=intmat and big n the result may be wrong because of int overflow |
---|
299 | EXAMPLE: example power; shows an example |
---|
300 | " |
---|
301 | { |
---|
302 | //---------------------------- type checking ---------------------------------- |
---|
303 | if( typeof(A)!="matrix" and typeof(A)!="intmat" ) |
---|
304 | { |
---|
305 | ERROR("no matrix or intmat!"); |
---|
306 | } |
---|
307 | if( ncols(A) != nrows(A) ) |
---|
308 | { |
---|
309 | ERROR("not a square matrix!"); |
---|
310 | } |
---|
311 | //---------------------------- trivial cases ---------------------------------- |
---|
312 | int ii; |
---|
313 | if( n <= 0 ) |
---|
314 | { |
---|
315 | if( typeof(A)=="matrix" ) |
---|
316 | { |
---|
317 | return (unitmat(nrows(A))); |
---|
318 | } |
---|
319 | if( typeof(A)=="intmat" ) |
---|
320 | { |
---|
321 | intmat B[nrows(A)][nrows(A)]; |
---|
322 | for( ii=1; ii<=nrows(A); ii++ ) |
---|
323 | { |
---|
324 | B[ii,ii] = 1; |
---|
325 | } |
---|
326 | return (B); |
---|
327 | } |
---|
328 | } |
---|
329 | if( n == 1 ) { return (A); } |
---|
330 | //---------------------------- sub procedure ---------------------------------- |
---|
331 | proc matpow (def A, int n) |
---|
332 | { |
---|
333 | def B = A*A; |
---|
334 | int ii= 2; |
---|
335 | int jj= 4; |
---|
336 | while( jj <= n ) |
---|
337 | { |
---|
338 | B=B*B; |
---|
339 | ii=jj; |
---|
340 | jj=2*jj; |
---|
341 | } |
---|
342 | return(B,n-ii); |
---|
343 | } |
---|
344 | //----------------------------- main program ---------------------------------- |
---|
345 | list L = matpow(A,n); |
---|
346 | def B = L[1]; |
---|
347 | ii = L[2]; |
---|
348 | while( ii>=2 ) |
---|
349 | { |
---|
350 | L = matpow(A,ii); |
---|
351 | B = B*L[1]; |
---|
352 | ii= L[2]; |
---|
353 | } |
---|
354 | if( ii == 0) { return(B); } |
---|
355 | if( ii == 1) { return(A*B); } |
---|
356 | } |
---|
357 | example |
---|
358 | { "EXAMPLE:"; echo = 2; |
---|
359 | intmat A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
360 | print(power(A,3));""; |
---|
361 | ring r=0,(x,y,z),dp; |
---|
362 | matrix B[3][3]=0,x,y,z,0,0,y,z,0; |
---|
363 | print(power(B,3));""; |
---|
364 | } |
---|
365 | /////////////////////////////////////////////////////////////////////////////// |
---|
366 | |
---|
367 | proc skewmat (int n, list #) |
---|
368 | "USAGE: skewmat(n[,id]); n integer, id ideal |
---|
369 | RETURN: skew-symmetric nxn matrix, with entries from id |
---|
370 | (default: id=maxideal(1)) |
---|
371 | skewmat(n); creates the generic skew-symmetric matrix |
---|
372 | NOTE: if id has less than n*(n-1)/2 elements, the matrix is |
---|
373 | filled with 0's, |
---|
374 | EXAMPLE: example skewmat; shows an example |
---|
375 | " |
---|
376 | { |
---|
377 | matrix B[n][n]; |
---|
378 | if( size(#)==0 ) { ideal id=maxideal(1); } |
---|
379 | else { ideal id=#[1]; } |
---|
380 | id = id,B[1..n,1..n]; |
---|
381 | int i,j; |
---|
382 | for( i=0; i<=n-2; i++ ) |
---|
383 | { |
---|
384 | B[i+1,i+2..n]=id[j+1..j+n-i-1]; |
---|
385 | j=j+n-i-1; |
---|
386 | } |
---|
387 | matrix A=transpose(B); |
---|
388 | B=B-A; |
---|
389 | return(B); |
---|
390 | } |
---|
391 | example |
---|
392 | { "EXAMPLE:"; echo = 2; |
---|
393 | ring R=0,x(1..5),lp; |
---|
394 | print(skewmat(4)); // the generic skew-symmetric matrix |
---|
395 | ring R1 = 0,(a,b,c),dp; |
---|
396 | matrix A=skewmat(4,maxideal(1)^2); |
---|
397 | print(A); |
---|
398 | int n=3; |
---|
399 | ideal i = ideal(randommat(1,n*(n-1) div 2,maxideal(1),9)); |
---|
400 | print(skewmat(n,i)); // skew matrix of generic linear forms |
---|
401 | kill R1; |
---|
402 | } |
---|
403 | /////////////////////////////////////////////////////////////////////////////// |
---|
404 | |
---|
405 | proc submat (matrix A, intvec r, intvec c) |
---|
406 | "USAGE: submat(A,r,c); A=matrix, r,c=intvec |
---|
407 | RETURN: matrix, submatrix of A with rows specified by intvec r |
---|
408 | and columns specified by intvec c. |
---|
409 | EXAMPLE: example submat; shows an example |
---|
410 | " |
---|
411 | { |
---|
412 | matrix B[size(r)][size(c)]=A[r,c]; |
---|
413 | return(B); |
---|
414 | } |
---|
415 | example |
---|
416 | { "EXAMPLE:"; echo = 2; |
---|
417 | ring R=32003,(x,y,z),lp; |
---|
418 | matrix A[4][4]=x,y,z,0,1,2,3,4,5,6,7,8,9,x2,y2,z2; |
---|
419 | print(A); |
---|
420 | intvec v=1,3,4; |
---|
421 | matrix B=submat(A,v,1..3); |
---|
422 | print(B); |
---|
423 | } |
---|
424 | /////////////////////////////////////////////////////////////////////////////// |
---|
425 | |
---|
426 | proc symmat (int n, list #) |
---|
427 | "USAGE: symmat(n[,id]); n integer, id ideal |
---|
428 | RETURN: symmetric nxn matrix, with entries from id (default: id=maxideal(1)) |
---|
429 | NOTE: if id has less than n*(n+1)/2 elements, the matrix is filled with 0's, |
---|
430 | symmat(n); creates the generic symmetric matrix |
---|
431 | EXAMPLE: example symmat; shows an example |
---|
432 | " |
---|
433 | { |
---|
434 | matrix B[n][n]; |
---|
435 | if( size(#)==0 ) { ideal id=maxideal(1); } |
---|
436 | else { ideal id=#[1]; } |
---|
437 | id = id,B[1..n,1..n]; |
---|
438 | int i,j; |
---|
439 | for( i=0; i<=n-1; i++ ) |
---|
440 | { |
---|
441 | B[i+1,i+1..n]=id[j+1..j+n-i]; |
---|
442 | j=j+n-i; |
---|
443 | } |
---|
444 | matrix A=transpose(B); |
---|
445 | for( i=1; i<=n; i++ ) { A[i,i]=0; } |
---|
446 | B=A+B; |
---|
447 | return(B); |
---|
448 | } |
---|
449 | example |
---|
450 | { "EXAMPLE:"; echo = 2; |
---|
451 | ring R=0,x(1..10),lp; |
---|
452 | print(symmat(4)); // the generic symmetric matrix |
---|
453 | ring R1 = 0,(a,b,c),dp; |
---|
454 | matrix A=symmat(4,maxideal(1)^3); |
---|
455 | print(A); |
---|
456 | int n=3; |
---|
457 | ideal i = ideal(randommat(1,n*(n+1) div 2,maxideal(1),9)); |
---|
458 | print(symmat(n,i)); // symmetric matrix of generic linear forms |
---|
459 | kill R1; |
---|
460 | } |
---|
461 | /////////////////////////////////////////////////////////////////////////////// |
---|
462 | |
---|
463 | proc unitmat (int n) |
---|
464 | "USAGE: unitmat(n); n integer >= 0 |
---|
465 | RETURN: nxn unit matrix |
---|
466 | NOTE: needs a basering, diagonal entries are numbers (=1) in the basering |
---|
467 | EXAMPLE: example unitmat; shows an example |
---|
468 | " |
---|
469 | { |
---|
470 | return(matrix(freemodule(n))); |
---|
471 | } |
---|
472 | example |
---|
473 | { "EXAMPLE:"; echo = 2; |
---|
474 | ring r=32003,(x,y,z),lp; |
---|
475 | print(xyz*unitmat(4)); |
---|
476 | print(unitmat(5)); |
---|
477 | } |
---|
478 | /////////////////////////////////////////////////////////////////////////////// |
---|
479 | |
---|
480 | proc gauss_col (matrix A, list #) |
---|
481 | "USAGE: gauss_col(A[,e]); A a matrix, e any type |
---|
482 | RETURN: - a matrix B, if called with one argument; B is the complete column- |
---|
483 | reduced upper-triangular normal form of A if A is constant, |
---|
484 | (resp. as far as this is possible if A is a polynomial matrix; |
---|
485 | no division by polynomials). |
---|
486 | @* - a list L of two matrices, if called with two arguments; |
---|
487 | L satisfies L[1] = A * L[2] with L[1] the column-reduced form of A |
---|
488 | and L[2] the transformation matrix. |
---|
489 | NOTE: * The procedure just applies interred to A with ordering (C,dp). |
---|
490 | The transformation matrix is obtained by applying 'lift'. |
---|
491 | This should be faster than the procedure colred. |
---|
492 | @* * It should only be used with exact coefficient field (there is no |
---|
493 | pivoting and rounding error treatment). |
---|
494 | @* * Parameters are allowed. Hence, if the entries of A are parameters, |
---|
495 | B is the column-reduced form of A over the rational function field. |
---|
496 | SEE ALSO: colred |
---|
497 | EXAMPLE: example gauss_col; shows an example |
---|
498 | " |
---|
499 | { |
---|
500 | def R=basering; int u; |
---|
501 | string mp = string(minpoly); |
---|
502 | int n = nrows(A); |
---|
503 | int m = ncols(A); |
---|
504 | module M = A; |
---|
505 | intvec v = option(get); |
---|
506 | //------------------------ change ordering if necessary ---------------------- |
---|
507 | if( ordstr(R) != ("C,dp("+string(nvars(R))+")") ) |
---|
508 | { |
---|
509 | def @R=changeord(list(list("C",0:1),list("dp",1:nvars(R))),R); |
---|
510 | setring @R; u=1; |
---|
511 | if (mp!="0") { execute("minpoly="+mp+";");} |
---|
512 | matrix A = imap(R,A); |
---|
513 | module M = A; |
---|
514 | } |
---|
515 | //------------------------------ start computation --------------------------- |
---|
516 | option(redSB); |
---|
517 | M = simplify(interred(M),1); |
---|
518 | if(size(#) != 0) |
---|
519 | { |
---|
520 | module N = lift(A,M); |
---|
521 | } |
---|
522 | //--------------- reset ring and options and return -------------------------- |
---|
523 | if ( u==1 ) |
---|
524 | { |
---|
525 | setring R; |
---|
526 | M=imap(@R,M); |
---|
527 | if (size(#) != 0) |
---|
528 | { |
---|
529 | module N = imap(@R,N); |
---|
530 | } |
---|
531 | kill @R; |
---|
532 | } |
---|
533 | option(set,v); |
---|
534 | // M = sort(M,size(M)..1)[1]; |
---|
535 | A = matrix(M,n,m); |
---|
536 | if (size(#) != 0) |
---|
537 | { |
---|
538 | list L= A,matrix(N,m,m); |
---|
539 | return(L); |
---|
540 | } |
---|
541 | return(matrix(M,n,m)); |
---|
542 | } |
---|
543 | example |
---|
544 | { "EXAMPLE:"; echo = 2; |
---|
545 | ring r=(0,a,b),(A,B,C),dp; |
---|
546 | matrix m[8][6]= |
---|
547 | 0, 2*C, 0, 0, 0, 0, |
---|
548 | 0, -4*C,a*A, 0, 0, 0, |
---|
549 | b*B, -A, 0, 0, 0, 0, |
---|
550 | -A, B, 0, 0, 0, 0, |
---|
551 | -4*C, 0, B, 2, 0, 0, |
---|
552 | 2*A, B, 0, 0, 0, 0, |
---|
553 | 0, 3*B, 0, 0, 2b, 0, |
---|
554 | 0, AB, 0, 2*A,A, 2a;""; |
---|
555 | list L=gauss_col(m,1); |
---|
556 | print(L[1]); |
---|
557 | print(L[2]); |
---|
558 | |
---|
559 | ring S=0,x,(c,dp); |
---|
560 | matrix A[5][4] = |
---|
561 | 3, 1, 1, 1, |
---|
562 | 13, 8, 6,-7, |
---|
563 | 14,10, 6,-7, |
---|
564 | 7, 4, 3,-3, |
---|
565 | 2, 1, 0, 3; |
---|
566 | print(gauss_col(A)); |
---|
567 | } |
---|
568 | /////////////////////////////////////////////////////////////////////////////// |
---|
569 | |
---|
570 | proc gauss_row (matrix A, list #) |
---|
571 | "USAGE: gauss_row(A [,e]); A matrix, e any type |
---|
572 | RETURN: - a matrix B, if called with one argument; B is the complete row- |
---|
573 | reduced lower-triangular normal form of A if A is constant, |
---|
574 | (resp. as far as this is possible if A is a polynomial matrix; |
---|
575 | no division by polynomials). |
---|
576 | @* - a list L of two matrices, if called with two arguments; |
---|
577 | L satisfies transpose(L[2])*A=transpose(L[1]) |
---|
578 | with L[1] the row-reduced form of A |
---|
579 | and L[2] the transformation matrix. |
---|
580 | NOTE: * This procedure just applies gauss_col to the transposed matrix. |
---|
581 | The transformation matrix is obtained by applying lift. |
---|
582 | This should be faster than the procedure rowred. |
---|
583 | @* * It should only be used with exact coefficient field (there is no |
---|
584 | pivoting and rounding error treatment). |
---|
585 | @* * Parameters are allowed. Hence, if the entries of A are parameters, |
---|
586 | B is the row-reduced form of A over the rational function field. |
---|
587 | SEE ALSO: rowred |
---|
588 | EXAMPLE: example gauss_row; shows an example |
---|
589 | " |
---|
590 | { |
---|
591 | if(size(#) > 0) |
---|
592 | { |
---|
593 | list L = gauss_col(transpose(A),1); |
---|
594 | return(L); |
---|
595 | } |
---|
596 | A = gauss_col(transpose(A)); |
---|
597 | return(transpose(A)); |
---|
598 | } |
---|
599 | example |
---|
600 | { "EXAMPLE:"; echo = 2; |
---|
601 | ring r=(0,a,b),(A,B,C),dp; |
---|
602 | matrix m[6][8]= |
---|
603 | 0, 0, b*B, -A,-4C,2A,0, 0, |
---|
604 | 2C,-4C,-A,B, 0, B, 3B,AB, |
---|
605 | 0,a*A, 0, 0, B, 0, 0, 0, |
---|
606 | 0, 0, 0, 0, 2, 0, 0, 2A, |
---|
607 | 0, 0, 0, 0, 0, 0, 2b, A, |
---|
608 | 0, 0, 0, 0, 0, 0, 0, 2a;""; |
---|
609 | print(gauss_row(m));""; |
---|
610 | ring S=0,x,dp; |
---|
611 | matrix A[4][5] = 3, 1,1,-1,2, |
---|
612 | 13, 8,6,-7,1, |
---|
613 | 14,10,6,-7,1, |
---|
614 | 7, 4,3,-3,3; |
---|
615 | list L=gauss_row(A,1); |
---|
616 | print(L[1]); |
---|
617 | print(L[2]); |
---|
618 | } |
---|
619 | /////////////////////////////////////////////////////////////////////////////// |
---|
620 | |
---|
621 | proc addcol (matrix A, int c1, poly p, int c2) |
---|
622 | "USAGE: addcol(A,c1,p,c2); A matrix, p poly, c1, c2 positive integers |
---|
623 | RETURN: matrix, A being modified by adding p times column c1 to column c2 |
---|
624 | EXAMPLE: example addcol; shows an example |
---|
625 | " |
---|
626 | { |
---|
627 | int k=nrows(A); |
---|
628 | A[1..k,c2]=A[1..k,c2]+p*A[1..k,c1]; |
---|
629 | return(A); |
---|
630 | } |
---|
631 | example |
---|
632 | { "EXAMPLE:"; echo = 2; |
---|
633 | ring r=32003,(x,y,z),lp; |
---|
634 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
635 | print(A); |
---|
636 | print(addcol(A,1,xy,2)); |
---|
637 | } |
---|
638 | /////////////////////////////////////////////////////////////////////////////// |
---|
639 | |
---|
640 | proc addrow (matrix A, int r1, poly p, int r2) |
---|
641 | "USAGE: addrow(A,r1,p,r2); A matrix, p poly, r1, r2 positive integers |
---|
642 | RETURN: matrix, A being modified by adding p times row r1 to row r2 |
---|
643 | EXAMPLE: example addrow; shows an example |
---|
644 | " |
---|
645 | { |
---|
646 | int k=ncols(A); |
---|
647 | A[r2,1..k]=A[r2,1..k]+p*A[r1,1..k]; |
---|
648 | return(A); |
---|
649 | } |
---|
650 | example |
---|
651 | { "EXAMPLE:"; echo = 2; |
---|
652 | ring r=32003,(x,y,z),lp; |
---|
653 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
654 | print(A); |
---|
655 | print(addrow(A,1,xy,3)); |
---|
656 | } |
---|
657 | /////////////////////////////////////////////////////////////////////////////// |
---|
658 | |
---|
659 | proc multcol (matrix A, int c, poly p) |
---|
660 | "USAGE: multcol(A,c,p); A matrix, p poly, c positive integer |
---|
661 | RETURN: matrix, A being modified by multiplying column c by p |
---|
662 | EXAMPLE: example multcol; shows an example |
---|
663 | " |
---|
664 | { |
---|
665 | int k=nrows(A); |
---|
666 | A[1..k,c]=p*A[1..k,c]; |
---|
667 | return(A); |
---|
668 | } |
---|
669 | example |
---|
670 | { "EXAMPLE:"; echo = 2; |
---|
671 | ring r=32003,(x,y,z),lp; |
---|
672 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
673 | print(A); |
---|
674 | print(multcol(A,2,xy)); |
---|
675 | } |
---|
676 | /////////////////////////////////////////////////////////////////////////////// |
---|
677 | |
---|
678 | proc multrow (matrix A, int r, poly p) |
---|
679 | "USAGE: multrow(A,r,p); A matrix, p poly, r positive integer |
---|
680 | RETURN: matrix, A being modified by multiplying row r by p |
---|
681 | EXAMPLE: example multrow; shows an example |
---|
682 | " |
---|
683 | { |
---|
684 | int k=ncols(A); |
---|
685 | A[r,1..k]=p*A[r,1..k]; |
---|
686 | return(A); |
---|
687 | } |
---|
688 | example |
---|
689 | { "EXAMPLE:"; echo = 2; |
---|
690 | ring r=32003,(x,y,z),lp; |
---|
691 | matrix A[3][3]=1,2,3,4,5,6,7,8,9; |
---|
692 | print(A); |
---|
693 | print(multrow(A,2,xy)); |
---|
694 | } |
---|
695 | /////////////////////////////////////////////////////////////////////////////// |
---|
696 | |
---|
697 | proc permcol (matrix A, int c1, int c2) |
---|
698 | "USAGE: permcol(A,c1,c2); A matrix, c1,c2 positive integers |
---|
699 | RETURN: matrix, A being modified by permuting columns c1 and c2 |
---|
700 | EXAMPLE: example permcol; shows an example |
---|
701 | " |
---|
702 | { |
---|
703 | matrix B=A; |
---|
704 | int k=nrows(B); |
---|
705 | B[1..k,c1]=A[1..k,c2]; |
---|
706 | B[1..k,c2]=A[1..k,c1]; |
---|
707 | return(B); |
---|
708 | } |
---|
709 | example |
---|
710 | { "EXAMPLE:"; echo = 2; |
---|
711 | ring r=32003,(x,y,z),lp; |
---|
712 | matrix A[3][3]=1,x,3,4,y,6,7,z,9; |
---|
713 | print(A); |
---|
714 | print(permcol(A,2,3)); |
---|
715 | } |
---|
716 | /////////////////////////////////////////////////////////////////////////////// |
---|
717 | |
---|
718 | proc permrow (matrix A, int r1, int r2) |
---|
719 | "USAGE: permrow(A,r1,r2); A matrix, r1,r2 positive integers |
---|
720 | RETURN: matrix, A being modified by permuting rows r1 and r2 |
---|
721 | EXAMPLE: example permrow; shows an example |
---|
722 | " |
---|
723 | { |
---|
724 | matrix B=A; |
---|
725 | int k=ncols(B); |
---|
726 | B[r1,1..k]=A[r2,1..k]; |
---|
727 | B[r2,1..k]=A[r1,1..k]; |
---|
728 | return(B); |
---|
729 | } |
---|
730 | example |
---|
731 | { "EXAMPLE:"; echo = 2; |
---|
732 | ring r=32003,(x,y,z),lp; |
---|
733 | matrix A[3][3]=1,2,3,x,y,z,7,8,9; |
---|
734 | print(A); |
---|
735 | print(permrow(A,2,1)); |
---|
736 | } |
---|
737 | /////////////////////////////////////////////////////////////////////////////// |
---|
738 | |
---|
739 | proc rowred (matrix A,list #) |
---|
740 | "USAGE: rowred(A[,e]); A matrix, e any type |
---|
741 | RETURN: - a matrix B, being the row reduced form of A, if rowred is called |
---|
742 | with one argument. |
---|
743 | (as far as this is possible over the polynomial ring; no division |
---|
744 | by polynomials) |
---|
745 | @* - a list L of two matrices, such that L[1] = L[2] * A with L[1] |
---|
746 | the row-reduced form of A and L[2] the transformation matrix |
---|
747 | (if rowred is called with two arguments). |
---|
748 | ASSUME: The entries of A are in the base field. It is not verified whether |
---|
749 | this assumption holds. |
---|
750 | NOTE: * This procedure is designed for teaching purposes mainly. |
---|
751 | @* * The straight forward Gaussian algorithm is implemented in the |
---|
752 | library (no standard basis computation). |
---|
753 | The transformation matrix is obtained by concatenating a unit |
---|
754 | matrix to A. proc gauss_row should be faster. |
---|
755 | @* * It should only be used with exact coefficient field (there is no |
---|
756 | pivoting) over the polynomial ring (ordering lp or dp). |
---|
757 | @* * Parameters are allowed. Hence, if the entries of A are parameters |
---|
758 | the computation takes place over the field of rational functions. |
---|
759 | SEE ALSO: gauss_row |
---|
760 | EXAMPLE: example rowred; shows an example |
---|
761 | " |
---|
762 | { |
---|
763 | int m,n=nrows(A),ncols(A); |
---|
764 | int i,j,k,l,rk; |
---|
765 | poly p; |
---|
766 | matrix d[m][n]; |
---|
767 | for (i=1;i<=m;i++) |
---|
768 | { for (j=1;j<=n;j++) |
---|
769 | { p = A[i,j]; |
---|
770 | if (ord(p)==0) |
---|
771 | { if (deg(p)==0) { d[i,j]=p; } |
---|
772 | } |
---|
773 | } |
---|
774 | } |
---|
775 | matrix b = A; |
---|
776 | if (size(#) != 0) { b = concat(b,unitmat(m)); } |
---|
777 | for (l=1;l<=n;l=l+1) |
---|
778 | { |
---|
779 | k = findfirst(ideal(d[l]),rk+1); |
---|
780 | if (k) |
---|
781 | { rk = rk+1; |
---|
782 | b = permrow(b,rk,k); |
---|
783 | p = b[rk,l]; p = 1/p; |
---|
784 | b = multrow(b,rk,p); |
---|
785 | for (i=1;i<=m;i++) |
---|
786 | { |
---|
787 | if (rk-i) { b = addrow(b,rk,-b[i,l],i);} |
---|
788 | } |
---|
789 | d = 0; |
---|
790 | for (i=rk+1;i<=m;i++) |
---|
791 | { for (j=l+1;j<=n;j++) |
---|
792 | { p = b[i,j]; |
---|
793 | if (ord(p)==0) |
---|
794 | { if (deg(p)==0) { d[i,j]=p; } |
---|
795 | } |
---|
796 | } |
---|
797 | } |
---|
798 | |
---|
799 | } |
---|
800 | } |
---|
801 | d = submat(b,1..m,1..n); |
---|
802 | if (size(#)) |
---|
803 | { |
---|
804 | list L=d,submat(b,1..m,n+1..n+m); |
---|
805 | return(L); |
---|
806 | } |
---|
807 | return(d); |
---|
808 | } |
---|
809 | example |
---|
810 | { "EXAMPLE:"; echo = 2; |
---|
811 | ring r=(0,a,b),(A,B,C),dp; |
---|
812 | matrix m[6][8]= |
---|
813 | 0, 0, b*B, -A,-4C,2A,0, 0, |
---|
814 | 2C,-4C,-A,B, 0, B, 3B,AB, |
---|
815 | 0,a*A, 0, 0, B, 0, 0, 0, |
---|
816 | 0, 0, 0, 0, 2, 0, 0, 2A, |
---|
817 | 0, 0, 0, 0, 0, 0, 2b, A, |
---|
818 | 0, 0, 0, 0, 0, 0, 0, 2a;""; |
---|
819 | print(rowred(m));""; |
---|
820 | list L=rowred(m,1); |
---|
821 | print(L[1]); |
---|
822 | print(L[2]); |
---|
823 | } |
---|
824 | /////////////////////////////////////////////////////////////////////////////// |
---|
825 | |
---|
826 | proc colred (matrix A,list #) |
---|
827 | "USAGE: colred(A[,e]); A matrix, e any type |
---|
828 | RETURN: - a matrix B, being the column reduced form of A, if colred is |
---|
829 | called with one argument. |
---|
830 | (as far as this is possible over the polynomial ring; |
---|
831 | no division by polynomials) |
---|
832 | @* - a list L of two matrices, such that L[1] = A * L[2] with L[1] |
---|
833 | the column-reduced form of A and L[2] the transformation matrix |
---|
834 | (if colred is called with two arguments). |
---|
835 | ASSUME: The entries of A are in the base field. It is not verified whether |
---|
836 | this assumption holds. |
---|
837 | NOTE: * This procedure is designed for teaching purposes mainly. |
---|
838 | @* * It applies rowred to the transposed matrix. |
---|
839 | proc gauss_col should be faster. |
---|
840 | @* * It should only be used with exact coefficient field (there is no |
---|
841 | pivoting) over the polynomial ring (ordering lp or dp). |
---|
842 | @* * Parameters are allowed. Hence, if the entries of A are parameters |
---|
843 | the computation takes place over the field of rational functions. |
---|
844 | SEE ALSO: gauss_col |
---|
845 | EXAMPLE: example colred; shows an example |
---|
846 | " |
---|
847 | { |
---|
848 | A = transpose(A); |
---|
849 | if (size(#)) |
---|
850 | { list L = rowred(A,1); return(transpose(L[1]),transpose(L[2]));} |
---|
851 | else |
---|
852 | { return(transpose(rowred(A)));} |
---|
853 | } |
---|
854 | example |
---|
855 | { "EXAMPLE:"; echo = 2; |
---|
856 | ring r=(0,a,b),(A,B,C),dp; |
---|
857 | matrix m[8][6]= |
---|
858 | 0, 2*C, 0, 0, 0, 0, |
---|
859 | 0, -4*C,a*A, 0, 0, 0, |
---|
860 | b*B, -A, 0, 0, 0, 0, |
---|
861 | -A, B, 0, 0, 0, 0, |
---|
862 | -4*C, 0, B, 2, 0, 0, |
---|
863 | 2*A, B, 0, 0, 0, 0, |
---|
864 | 0, 3*B, 0, 0, 2b, 0, |
---|
865 | 0, AB, 0, 2*A,A, 2a;""; |
---|
866 | print(colred(m));""; |
---|
867 | list L=colred(m,1); |
---|
868 | print(L[1]); |
---|
869 | print(L[2]); |
---|
870 | } |
---|
871 | ////////////////////////////////////////////////////////////////////////////// |
---|
872 | |
---|
873 | proc linear_relations(module M) |
---|
874 | "USAGE: linear_relations(M); |
---|
875 | M: a module |
---|
876 | ASSUME: All non-zero entries of M are homogeneous polynomials of the same |
---|
877 | positive degree. The base field must be an exact field (not real |
---|
878 | or complex). |
---|
879 | It is not checked whether these assumptions hold. |
---|
880 | RETURN: a maximal module R such that M*R is formed by zero vectors. |
---|
881 | EXAMPLE: example linear_relations; shows an example. |
---|
882 | " |
---|
883 | { int n = ncols(M); |
---|
884 | def BaseR = basering; |
---|
885 | def br = changeord(list(list("dp",1:nvars(basering)))); |
---|
886 | setring br; |
---|
887 | module M = imap(BaseR,M); |
---|
888 | ideal vars = maxideal(1); |
---|
889 | ring tmpR = 0, ('y(1..n)), dp; |
---|
890 | def newR = br + tmpR; |
---|
891 | setring newR; |
---|
892 | module M = imap(br,M); |
---|
893 | ideal vars = imap(br,vars); |
---|
894 | attrib(vars,"isSB",1); |
---|
895 | for (int i = 1; i<=n; i++) { |
---|
896 | M[i] = M[i] + 'y(i)*gen(1); |
---|
897 | } |
---|
898 | M = interred(M); |
---|
899 | module redM = NF(M,vars); |
---|
900 | module REL; |
---|
901 | int sizeREL; |
---|
902 | int j; |
---|
903 | for (i=1; i<=n; i++) { |
---|
904 | if (M[i][1]==redM[i][1]) { //-- relation found! |
---|
905 | sizeREL++; |
---|
906 | REL[sizeREL]=0; |
---|
907 | for (j=1; j<=n; j++) { |
---|
908 | REL[sizeREL] = REL[sizeREL] + (M[i][1]/'y(j))*gen(j); |
---|
909 | } |
---|
910 | } |
---|
911 | } |
---|
912 | setring BaseR; |
---|
913 | return(minbase(imap(newR,REL))); |
---|
914 | } |
---|
915 | example |
---|
916 | { "EXAMPLE:"; echo = 2; |
---|
917 | ring r = (3,w), (a,b,c,d),dp; |
---|
918 | minpoly = w2-w-1; |
---|
919 | module M = [a2,b2],[wab,w2c2+2b2],[(w-2)*a2+wab,wb2+w2c2]; |
---|
920 | module REL = linear_relations(M); |
---|
921 | pmat(REL); |
---|
922 | pmat(matrix(M)*matrix(REL)); |
---|
923 | } |
---|
924 | |
---|
925 | ////////////////////////////////////////////////////////////////////////////// |
---|
926 | |
---|
927 | static proc findfirst (ideal i,int t) |
---|
928 | { |
---|
929 | int n,k; |
---|
930 | for (n=t;n<=ncols(i);n=n+1) |
---|
931 | { |
---|
932 | if (i[n]!=0) { k=n;break;} |
---|
933 | } |
---|
934 | return(k); |
---|
935 | } |
---|
936 | ////////////////////////////////////////////////////////////////////////////// |
---|
937 | |
---|
938 | proc rm_unitcol(matrix A) |
---|
939 | "USAGE: rm_unitcol(A); A matrix (being row-reduced) |
---|
940 | RETURN: matrix, obtained from A by deleting unit columns (having just one 1 |
---|
941 | and else 0 as entries) and associated rows |
---|
942 | EXAMPLE: example rm_unitcol; shows an example |
---|
943 | " |
---|
944 | { |
---|
945 | int l,j; |
---|
946 | intvec v; |
---|
947 | for (j=1;j<=ncols(A);j++) |
---|
948 | { |
---|
949 | if (gen(l+1)==module(A)[j]) {l=l+1;} |
---|
950 | else { v=v,j;} |
---|
951 | } |
---|
952 | if (size(v)>1) |
---|
953 | { v = v[2..size(v)]; |
---|
954 | return(submat(A,l+1..nrows(A),v)); |
---|
955 | } |
---|
956 | else |
---|
957 | { return(0);} |
---|
958 | } |
---|
959 | example |
---|
960 | { "EXAMPLE:"; echo = 2; |
---|
961 | ring r=0,(A,B,C),dp; |
---|
962 | matrix m[6][8]= |
---|
963 | 0, 0, A, 0, 1,0, 0,0, |
---|
964 | 0, 0, -C2, 0, 0,0, 1,0, |
---|
965 | 0, 0, 0,1/2B, 0,0, 0,1, |
---|
966 | 0, 0, B, -A, 0,2A, 0,0, |
---|
967 | 2C,-4C, -A, B, 0,B, 0,0, |
---|
968 | 0, A, 0, 0, 0,0, 0,0; |
---|
969 | print(rm_unitcol(m)); |
---|
970 | } |
---|
971 | ////////////////////////////////////////////////////////////////////////////// |
---|
972 | |
---|
973 | proc rm_unitrow (matrix A) |
---|
974 | "USAGE: rm_unitrow(A); A matrix (being col-reduced) |
---|
975 | RETURN: matrix, obtained from A by deleting unit rows (having just one 1 |
---|
976 | and else 0 as entries) and associated columns |
---|
977 | EXAMPLE: example rm_unitrow; shows an example |
---|
978 | " |
---|
979 | { |
---|
980 | int l,j; |
---|
981 | intvec v; |
---|
982 | module M = transpose(A); |
---|
983 | for (j=1; j <= nrows(A); j++) |
---|
984 | { |
---|
985 | if (gen(l+1) == M[j]) { l=l+1; } |
---|
986 | else { v=v,j; } |
---|
987 | } |
---|
988 | if (size(v) > 1) |
---|
989 | { v = v[2..size(v)]; |
---|
990 | return(submat(A,v,l+1..ncols(A))); |
---|
991 | } |
---|
992 | else |
---|
993 | { return(0);} |
---|
994 | } |
---|
995 | example |
---|
996 | { "EXAMPLE:"; echo = 2; |
---|
997 | ring r=0,(A,B,C),dp; |
---|
998 | matrix m[8][6]= |
---|
999 | 0,0, 0, 0, 2C, 0, |
---|
1000 | 0,0, 0, 0, -4C,A, |
---|
1001 | A,-C2,0, B, -A, 0, |
---|
1002 | 0,0, 1/2B,-A,B, 0, |
---|
1003 | 1,0, 0, 0, 0, 0, |
---|
1004 | 0,0, 0, 2A,B, 0, |
---|
1005 | 0,1, 0, 0, 0, 0, |
---|
1006 | 0,0, 1, 0, 0, 0; |
---|
1007 | print(rm_unitrow(m)); |
---|
1008 | } |
---|
1009 | ////////////////////////////////////////////////////////////////////////////// |
---|
1010 | proc headStand(matrix M) |
---|
1011 | "USAGE: headStand(M); M matrix |
---|
1012 | RETURN: matrix B such that B[i][j]=M[n-i+1,m-j+1], n=nrows(M), m=ncols(M) |
---|
1013 | EXAMPLE: example headStand; shows an example |
---|
1014 | " |
---|
1015 | { |
---|
1016 | int i,j; |
---|
1017 | int n=nrows(M); |
---|
1018 | int m=ncols(M); |
---|
1019 | matrix B[n][m]; |
---|
1020 | for(i=1;i<=n;i++) |
---|
1021 | { |
---|
1022 | for(j=1;j<=m;j++) |
---|
1023 | { |
---|
1024 | B[n-i+1,m-j+1]=M[i,j]; |
---|
1025 | } |
---|
1026 | } |
---|
1027 | return(B); |
---|
1028 | } |
---|
1029 | example |
---|
1030 | { "EXAMPLE:"; echo = 2; |
---|
1031 | ring r=0,(A,B,C),dp; |
---|
1032 | matrix M[2][3]= |
---|
1033 | 0,A, B, |
---|
1034 | A2, B2, C; |
---|
1035 | print(M); |
---|
1036 | print(headStand(M)); |
---|
1037 | } |
---|
1038 | ////////////////////////////////////////////////////////////////////////////// |
---|
1039 | |
---|
1040 | // Symmetric/Exterior powers thanks to Oleksandr Iena for his persistence ;-) |
---|
1041 | |
---|
1042 | proc symmetricBasis(int n, int k, list #) |
---|
1043 | "USAGE: symmetricBasis(n, k[,s]); n int, k int, s string |
---|
1044 | RETURN: ring, poynomial ring containing the ideal \"symBasis\", |
---|
1045 | being a basis of the k-th symmetric power of an n-dim vector space. |
---|
1046 | NOTE: The output polynomial ring has characteristics 0 and n variables |
---|
1047 | named \"S(i)\", where the base variable name S is either given by the |
---|
1048 | optional string argument(which must not contain brackets) or equal to |
---|
1049 | "e" by default. |
---|
1050 | SEE ALSO: exteriorBasis |
---|
1051 | KEYWORDS: symmetric basis |
---|
1052 | EXAMPLE: example symmetricBasis; shows an example" |
---|
1053 | { |
---|
1054 | //------------------------ handle optional base variable name--------------- |
---|
1055 | string S = "e"; |
---|
1056 | if( size(#) > 0 ) |
---|
1057 | { |
---|
1058 | if( typeof(#[1]) != "string" ) |
---|
1059 | { |
---|
1060 | ERROR("Wrong optional argument: must be a string"); |
---|
1061 | } |
---|
1062 | S = #[1]; |
---|
1063 | if( (find(S, "(") + find(S, ")")) > 0 ) |
---|
1064 | { |
---|
1065 | ERROR("Wrong optional argument: must be a string without brackets"); |
---|
1066 | } |
---|
1067 | } |
---|
1068 | |
---|
1069 | //------------------------- create ring container for symmetric power basis- |
---|
1070 | list l2; |
---|
1071 | for (int ii = 1; ii <= n; ii++) |
---|
1072 | { |
---|
1073 | l2[ii] = S+"("+string(ii)+")"; |
---|
1074 | } |
---|
1075 | ring @@@SYM_POWER_RING_NAME = create_ring(0, l2, "dp"); |
---|
1076 | //------------------------- choose symmetric basis ------------------------- |
---|
1077 | ideal symBasis = maxideal(k); |
---|
1078 | |
---|
1079 | //------------------------- export and return ------------------------- |
---|
1080 | export symBasis; |
---|
1081 | return(basering); |
---|
1082 | } |
---|
1083 | example |
---|
1084 | { "EXAMPLE:"; echo = 2; |
---|
1085 | |
---|
1086 | // basis of the 3-rd symmetricPower of a 4-dim vector space: |
---|
1087 | def R = symmetricBasis(4, 3, "@e"); setring R; |
---|
1088 | R; // container ring: |
---|
1089 | symBasis; // symmetric basis: |
---|
1090 | } |
---|
1091 | |
---|
1092 | ////////////////////////////////////////////////////////////////////////////// |
---|
1093 | |
---|
1094 | proc exteriorBasis(int n, int k, list #) |
---|
1095 | "USAGE: exteriorBasis(n, k[,s]); n int, k int, s string |
---|
1096 | RETURN: qring, an exterior algebra containing the ideal \"extBasis\", |
---|
1097 | being a basis of the k-th exterior power of an n-dim vector space. |
---|
1098 | NOTE: The output polynomial ring has characteristics 0 and n variables |
---|
1099 | named \"S(i)\", where the base variable name S is either given by the |
---|
1100 | optional string argument(which must not contain brackets) or equal to |
---|
1101 | "e" by default. |
---|
1102 | SEE ALSO: symmetricBasis |
---|
1103 | KEYWORDS: exterior basis |
---|
1104 | EXAMPLE: example exteriorBasis; shows an example" |
---|
1105 | { |
---|
1106 | //------------------------ handle optional base variable name--------------- |
---|
1107 | string S = "e"; |
---|
1108 | if( size(#) > 0 ) |
---|
1109 | { |
---|
1110 | if( typeof(#[1]) != "string" ) |
---|
1111 | { |
---|
1112 | ERROR("Wrong optional argument: must be a string"); |
---|
1113 | } |
---|
1114 | S = #[1]; |
---|
1115 | if( (find(S, "(") + find(S, ")")) > 0 ) |
---|
1116 | { |
---|
1117 | ERROR("Wrong optional argument: must be a string without brackets"); |
---|
1118 | } |
---|
1119 | } |
---|
1120 | |
---|
1121 | //------------------------- create ring container for symmetric power basis- |
---|
1122 | list l2; |
---|
1123 | for (int ii = 1; ii <= n; ii++) |
---|
1124 | { |
---|
1125 | l2[ii] = S+"("+string(ii)+")"; |
---|
1126 | } |
---|
1127 | ring @@@EXT_POWER_RING_NAME = create_ring(0, l2, "dp"); |
---|
1128 | |
---|
1129 | //------------------------- choose exterior basis ------------------------- |
---|
1130 | def T = superCommutative(); setring T; |
---|
1131 | ideal extBasis = simplify( NF(maxideal(k), std(0)), 1 + 2 + 8 ); |
---|
1132 | |
---|
1133 | //------------------------- export and return ------------------------- |
---|
1134 | export extBasis; |
---|
1135 | return(basering); |
---|
1136 | } |
---|
1137 | example |
---|
1138 | { "EXAMPLE:"; echo = 2; |
---|
1139 | // basis of the 3-rd symmetricPower of a 4-dim vector space: |
---|
1140 | def r = exteriorBasis(4, 3, "@e"); setring r; |
---|
1141 | r; // container ring: |
---|
1142 | extBasis; // exterior basis: |
---|
1143 | } |
---|
1144 | |
---|
1145 | ////////////////////////////////////////////////////////////////////////////// |
---|
1146 | |
---|
1147 | static proc chooseSafeVarName(string prefix, string suffix) |
---|
1148 | "USAGE: give appropreate prefix for variable names |
---|
1149 | RETURN: safe variable name (repeated prefix + suffix) |
---|
1150 | " |
---|
1151 | { |
---|
1152 | string V = varstr(basering); |
---|
1153 | string S = suffix; |
---|
1154 | while( find(V, S) > 0 ) |
---|
1155 | { |
---|
1156 | S = prefix + S; |
---|
1157 | } |
---|
1158 | return(S); |
---|
1159 | } |
---|
1160 | |
---|
1161 | ////////////////////////////////////////////////////////////////////////////// |
---|
1162 | |
---|
1163 | static proc mapPower(int p, module A, int k, def Tn, def Tm) |
---|
1164 | "USAGE: by both symmetric- and exterior-Power" |
---|
1165 | NOTE: everything over the basering! |
---|
1166 | module A (matrix of the map), int k (power) |
---|
1167 | rings Tn is source- and Tm is image-ring with bases |
---|
1168 | resp. Ink and Imk. |
---|
1169 | M = max dim of Image, N - dim. of source |
---|
1170 | SEE ALSO: symmetricPower, exteriorPower" |
---|
1171 | { |
---|
1172 | def save = basering; |
---|
1173 | |
---|
1174 | int n = nvars(save); |
---|
1175 | int M = nrows(A); |
---|
1176 | int N = ncols(A); |
---|
1177 | |
---|
1178 | int i, j; |
---|
1179 | |
---|
1180 | //------------------------- compute matrix of single images ------------------ |
---|
1181 | def Rm = save + Tm; setring Rm; |
---|
1182 | dbprint(p-2, "Temporary Working Ring", Rm); |
---|
1183 | |
---|
1184 | module A = imap(save, A); |
---|
1185 | |
---|
1186 | ideal B; poly t; |
---|
1187 | |
---|
1188 | for( i = N; i > 0; i-- ) |
---|
1189 | { |
---|
1190 | t = 0; |
---|
1191 | for( j = M; j > 0; j-- ) |
---|
1192 | { |
---|
1193 | t = t + A[i][j] * var(n + j); |
---|
1194 | } |
---|
1195 | |
---|
1196 | B[i] = t; |
---|
1197 | } |
---|
1198 | |
---|
1199 | dbprint(p-1, "Matrix of single images", B); |
---|
1200 | |
---|
1201 | //------------------------- compute image --------------------- |
---|
1202 | // apply S^k(A): Tn -> Rm to Source basis vectors Ink: |
---|
1203 | map TMap = Tn, B; |
---|
1204 | |
---|
1205 | ideal C = NF(TMap(Ink), std(0)); |
---|
1206 | dbprint(p-1, "Image Matrix: ", C); |
---|
1207 | |
---|
1208 | |
---|
1209 | //------------------------- write it in Image basis --------------------- |
---|
1210 | ideal Imk = imap(Tm, Imk); |
---|
1211 | |
---|
1212 | module D; poly lm; vector tt; |
---|
1213 | |
---|
1214 | for( i = ncols(C); i > 0; i-- ) |
---|
1215 | { |
---|
1216 | t = C[i]; |
---|
1217 | tt = 0; |
---|
1218 | |
---|
1219 | while( t != 0 ) |
---|
1220 | { |
---|
1221 | lm = leadmonom(t); |
---|
1222 | // lm; |
---|
1223 | for( j = ncols(Imk); j > 0; j-- ) |
---|
1224 | { |
---|
1225 | if( lm / Imk[j] != 0 ) |
---|
1226 | { |
---|
1227 | tt = tt + (lead(t) / Imk[j]) * gen(j); |
---|
1228 | break; |
---|
1229 | } |
---|
1230 | } |
---|
1231 | t = t - lead(t); |
---|
1232 | } |
---|
1233 | |
---|
1234 | D[i] = tt; |
---|
1235 | } |
---|
1236 | |
---|
1237 | //------------------------- map it back and return --------------------- |
---|
1238 | setring save; |
---|
1239 | return( imap(Rm, D) ); |
---|
1240 | } |
---|
1241 | |
---|
1242 | |
---|
1243 | ////////////////////////////////////////////////////////////////////////////// |
---|
1244 | |
---|
1245 | proc symmetricPower(module A, int k) |
---|
1246 | "USAGE: symmetricPower(A, k); A module, k int |
---|
1247 | RETURN: module: the k-th symmetric power of A |
---|
1248 | NOTE: the chosen bases and most of intermediate data will be shown if |
---|
1249 | printlevel is big enough |
---|
1250 | SEE ALSO: exteriorPower |
---|
1251 | KEYWORDS: symmetric power |
---|
1252 | EXAMPLE: example symmetricPower; shows an example" |
---|
1253 | { |
---|
1254 | int p = printlevel - voice + 2; |
---|
1255 | |
---|
1256 | def save = basering; |
---|
1257 | |
---|
1258 | int M = nrows(A); |
---|
1259 | int N = ncols(A); |
---|
1260 | |
---|
1261 | string S = chooseSafeVarName("@", "@_e"); |
---|
1262 | |
---|
1263 | //------------------------- choose source basis ------------------------- |
---|
1264 | def Tn = symmetricBasis(N, k, S); setring Tn; |
---|
1265 | ideal Ink = symBasis; |
---|
1266 | export Ink; |
---|
1267 | dbprint(p-3, "Temporary Source Ring", basering); |
---|
1268 | dbprint(p, "S^k(Source Basis)", Ink); |
---|
1269 | |
---|
1270 | //------------------------- choose image basis ------------------------- |
---|
1271 | def Tm = symmetricBasis(M, k, S); setring Tm; |
---|
1272 | ideal Imk = symBasis; |
---|
1273 | export Imk; |
---|
1274 | dbprint(p-3, "Temporary Image Ring", basering); |
---|
1275 | dbprint(p, "S^k(Image Basis)", Imk); |
---|
1276 | |
---|
1277 | //------------------------- compute and return S^k(A) in chosen bases -- |
---|
1278 | setring save; |
---|
1279 | |
---|
1280 | return(mapPower(p, A, k, Tn, Tm)); |
---|
1281 | } |
---|
1282 | example |
---|
1283 | { "EXAMPLE:"; echo = 2; |
---|
1284 | |
---|
1285 | ring r = (0),(a, b, c, d), dp; r; |
---|
1286 | module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2); print(B); |
---|
1287 | |
---|
1288 | // symmetric power over a commutative K-algebra: |
---|
1289 | print(symmetricPower(B, 2)); |
---|
1290 | print(symmetricPower(B, 3)); |
---|
1291 | |
---|
1292 | // symmetric power over an exterior algebra: |
---|
1293 | def g = superCommutative(); setring g; g; |
---|
1294 | |
---|
1295 | module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2); print(B); |
---|
1296 | |
---|
1297 | print(symmetricPower(B, 2)); // much smaller! |
---|
1298 | print(symmetricPower(B, 3)); // zero! (over an exterior algebra!) |
---|
1299 | |
---|
1300 | } |
---|
1301 | |
---|
1302 | ////////////////////////////////////////////////////////////////////////////// |
---|
1303 | |
---|
1304 | proc exteriorPower(module A, int k) |
---|
1305 | "USAGE: exteriorPower(A, k); A module, k int |
---|
1306 | RETURN: module: the k-th exterior power of A |
---|
1307 | NOTE: the chosen bases and most of intermediate data will be shown if |
---|
1308 | printlevel is big enough. Last rows will be invisible if zero. |
---|
1309 | SEE ALSO: symmetricPower |
---|
1310 | KEYWORDS: exterior power |
---|
1311 | EXAMPLE: example exteriorPower; shows an example" |
---|
1312 | { |
---|
1313 | int p = printlevel - voice + 2; |
---|
1314 | def save = basering; |
---|
1315 | |
---|
1316 | int M = nrows(A); |
---|
1317 | int N = ncols(A); |
---|
1318 | |
---|
1319 | string S = chooseSafeVarName("@", "@_e"); |
---|
1320 | |
---|
1321 | //------------------------- choose source basis ------------------------- |
---|
1322 | def Tn = exteriorBasis(N, k, S); setring Tn; |
---|
1323 | ideal Ink = extBasis; |
---|
1324 | export Ink; |
---|
1325 | dbprint(p-3, "Temporary Source Ring", basering); |
---|
1326 | dbprint(p, "E^k(Source Basis)", Ink); |
---|
1327 | |
---|
1328 | //------------------------- choose image basis ------------------------- |
---|
1329 | def Tm = exteriorBasis(M, k, S); setring Tm; |
---|
1330 | ideal Imk = extBasis; |
---|
1331 | export Imk; |
---|
1332 | dbprint(p-3, "Temporary Image Ring", basering); |
---|
1333 | dbprint(p, "E^k(Image Basis)", Imk); |
---|
1334 | |
---|
1335 | //------------------------- compute and return E^k(A) in chosen bases -- |
---|
1336 | setring save; |
---|
1337 | return(mapPower(p, A, k, Tn, Tm)); |
---|
1338 | } |
---|
1339 | example |
---|
1340 | { "EXAMPLE:"; echo = 2; |
---|
1341 | ring r = (0),(a, b, c, d, e, f), dp; |
---|
1342 | r; "base ring:"; |
---|
1343 | |
---|
1344 | module B = a*gen(1) + c*gen(2) + e*gen(3), |
---|
1345 | b*gen(1) + d*gen(2) + f*gen(3), |
---|
1346 | e*gen(1) + f*gen(3); |
---|
1347 | |
---|
1348 | print(B); |
---|
1349 | print(exteriorPower(B, 2)); |
---|
1350 | print(exteriorPower(B, 3)); |
---|
1351 | |
---|
1352 | def g = superCommutative(); setring g; g; |
---|
1353 | |
---|
1354 | module A = a*gen(1), b * gen(1), c*gen(2), d * gen(2); |
---|
1355 | print(A); |
---|
1356 | |
---|
1357 | print(exteriorPower(A, 2)); |
---|
1358 | |
---|
1359 | module B = a*gen(1) + c*gen(2) + e*gen(3), |
---|
1360 | b*gen(1) + d*gen(2) + f*gen(3), |
---|
1361 | e*gen(1) + f*gen(3); |
---|
1362 | print(B); |
---|
1363 | |
---|
1364 | print(exteriorPower(B, 2)); |
---|
1365 | print(exteriorPower(B, 3)); |
---|
1366 | |
---|
1367 | } |
---|
1368 | |
---|
1369 | ////////////////////////////////////////////////////////////////////////////// |
---|