source: git/Singular/LIB/matrix.lib @ 7708934

spielwiese
Last change on this file since 7708934 was c2aa97, checked in by Gert-Martin Greuel <greuel@…>, 23 years ago
* GMG: Kosmetik fuer html-Hilfe, sowie proc tensor korrigiert git-svn-id: file:///usr/local/Singular/svn/trunk@4990 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 24.9 KB
Line 
1// GMG/BM, last modified: 8.10.98
2///////////////////////////////////////////////////////////////////////////////
3version="$Id: matrix.lib,v 1.16 2000-12-23 17:17:28 greuel Exp $";
4category="Linear Algebra";
5info="
6LIBRARY:  matrix.lib    Elementary Matrix Operations
7
8PROCEDURES:
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
38LIB "inout.lib";
39LIB "ring.lib";
40LIB "random.lib";
41///////////////////////////////////////////////////////////////////////////////
42
43proc compress (A)
44"USAGE:   compress(A); A matrix/ideal/module/intmat/intvec
45RETURN:  same type, zero columns/generators from A deleted
46         (if A=intvec, zero elements are deleted)
47EXAMPLE: 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}
71example
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///////////////////////////////////////////////////////////////////////////////
85proc concat (list #)
86"USAGE:   concat(A1,A2,..); A1,A2,... matrices
87RETURN:  matrix, concatenation of A1,A2,... . Number of rows of result matrix
88         is max(nrows(A1),nrows(A2),...)
89EXAMPLE: 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}
97example
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
109proc diag (list #)
110"USAGE:   diag(p,n); p poly, n integer
111         diag(A);   A matrix
112RETURN:  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
115EXAMPLE: 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}
127example
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
137proc dsum (list #)
138"USAGE:   dsum(A1,A2,..); A1,A2,... matrices
139RETURN:  matrix, direct sum of A1,A2,...
140EXAMPLE: 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}
155example
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
168proc flatten (matrix A)
169"USAGE:   flatten(A); A matrix
170RETURN:  ideal, generated by all entries from A
171EXAMPLE: example flatten; shows an example
172"
173{
174   return(ideal(A));
175}
176example
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
185proc genericmat (int n,int m,list #)
186"USAGE:   genericmat(n,m[,id]);  n,m=integers, id=ideal
187RETURN:  nxm matrix, with entries from id (default: id=maxideal(1))
188NOTE:    if id has less than nxm elements, the matrix is filled with 0's,
189         genericmat(n,m); creates the generic nxm matrix
190EXAMPLE: 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}
199example
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
213proc is_complex (list c)
214"USAGE:   is_complex(c); c = list of size-compatible modules or matrices
215RETURN:  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
217NOTE:    Ideals are treated internally as 1-line matrices
218EXAMPLE: 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}
236example
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
247proc outer (matrix A, matrix B)
248"USAGE:   outer(A,B); A,B matrices
249RETURN:  matrix, outer product of A and B
250EXAMPLE: 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}
274example
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
285proc power ( A, int n)
286"USAGE:   power(A,n);  A a square-matrix of type intmat or matrix, n=integer
287RETURN:  intmat resp. matrix, the n-th power of A
288NOTE:    for A=intmat and big n the result may be wrong because of int overflow
289EXAMPLE: 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}
349example
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
361proc skewmat (int n, list #)
362"USAGE:   skewmat(n[,id]);  n integer, id ideal
363RETURN:  skew-symmetric nxn matrix, with entries from id
364         (default: id=maxideal(1))
365NOTE:    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
367EXAMPLE: 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}
384example
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
398proc submat (matrix A, intvec r, intvec c)
399"USAGE:   submat(A,r,c);  A=matrix, r,c=intvec
400RETURN:  matrix, submatrix of A with rows specified by intvec r and columns
401         specified by intvec c
402EXAMPLE: example submat; shows an example
403"
404{
405   matrix B[size(r)][size(c)]=A[r,c];
406   return(B);
407}
408example
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
419proc symmat (int n, list #)
420"USAGE:   symmat(n[,id]);  n integer, id ideal
421RETURN:  symmetric nxn matrix, with entries from id (default: id=maxideal(1))
422NOTE:    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
424EXAMPLE: 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}
442example
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
456proc tensor (matrix A, matrix B)
457"USAGE:   tensor(A,B); A,B matrices
458RETURN:  matrix, tensor product of A and B
459EXAMPLE: 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}
476example
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
487proc unitmat (int n)
488"USAGE:   unitmat(n);  n integer >= 0
489RETURN:  nxn unit matrix
490NOTE:    needs a basering, diagonal entries are numbers (=1) in the basering
491EXAMPLE: example unitmat; shows an example
492"
493{
494   return(matrix(freemodule(n)));
495}
496example
497{ "EXAMPLE:"; echo = 2;
498   ring r=32003,(x,y,z),lp;
499   print(xyz*unitmat(4));
500   print(unitmat(5));
501}
502///////////////////////////////////////////////////////////////////////////////
503
504proc gauss_col (matrix A)
505"USAGE:   gauss_col(A); A=matrix with constant coefficients
506RETURN:  matrix = col-reduced lower-triagonal normal form of A
507NOTE:    works fine for constant matrices, use proc colred for poly matrices
508EXAMPLE: 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}
523example
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
535proc gauss_row (matrix A)
536"USAGE:   gauss_row(A); A=matrix with constant coefficients
537RETURN:  matrix = row-reduced upper-triangular normal form of A
538NOTE:    may be used to solve a system of linear equations
539         works fine for constant matrices, use proc rowred for poly matrices
540EXAMPLE: example gauss_row; shows an example
541"
542{
543   A = gauss_col(transpose(A));
544   return(transpose(A));
545}
546example
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
557proc addcol (matrix A, int c1, poly p, int c2)
558"USAGE:   addcol(A,c1,p,c2);  A matrix, p poly, c1, c2 positive integers
559RETURN:  matrix,  A being modified by adding p times column c1 to column c2
560EXAMPLE: 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}
566example
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
575proc addrow (matrix A, int r1, poly p, int r2)
576"USAGE:   addcol(A,r1,p,r2);  A matrix, p poly, r1, r2 positive integers
577RETURN:  matrix,  A being modified by adding p times row r1 to row r2
578EXAMPLE: 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}
584example
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
593proc multcol (matrix A, int c, poly p)
594"USAGE:   addcol(A,c,p);  A matrix, p poly, c positive integer
595RETURN:  matrix,  A being modified by multiplying column c with p
596EXAMPLE: example multcol; shows an example
597"
598{
599   A[1..nrows(A),c]=p*A[1..nrows(A),c];
600   return(A);
601}
602example
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
611proc multrow (matrix A, int r, poly p)
612"USAGE:   multrow(A,r,p);  A matrix, p poly, r positive integer
613RETURN:  matrix,  A being modified by multiplying row r with p
614EXAMPLE: example multrow; shows an example
615"
616{
617   A[r,1..ncols(A)]=p*A[r,1..ncols(A)];
618   return(A);
619}
620example
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
629proc permcol (matrix A, int c1, int c2)
630"USAGE:   permcol(A,c1,c2);  A matrix, c1,c2 positive integers
631RETURN:  matrix,  A being modified by permuting column c1 and c2
632EXAMPLE: 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}
640example
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
649proc permrow (matrix A, int r1, int r2)
650"USAGE:   permrow(A,r1,r2);  A matrix, r1,r2 positive integers
651RETURN:  matrix,  A being modified by permuting row r1 and r2
652EXAMPLE: 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}
660example
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
669proc rowred (matrix A,list #)
670"USAGE:   rowred(A[,any]);  A matrix
671RETURN:  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
673EXAMPLE: 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}
718example
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
740proc colred (matrix A,list #)
741"USAGE:   colred(A[,any]);  A matrix
742RETURN:  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
744EXAMPLE: 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}
753example
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
777static 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
788proc rm_unitcol(matrix A)
789"USAGE:   rm_unitcol(A); A matrix (being row reduced)
790RETURN:  matrix, obtained from A by deleting unit columns (having just one 1
791         and else 0 as entries) and associated rows
792EXAMPLE: 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}
809example
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
831proc rm_unitrow (matrix A)
832"USAGE:   rm_unitrow(A); A matrix (being colreduced)
833RETURN:  matrix, obtained from A by deleting unit rows (having just one 1
834         and else 0 as entries) and associated columns
835EXAMPLE: 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}
853example
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//////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.