# source:git/Singular/LIB/matrix.lib@3c4dcc

spielwiese
Last change on this file since 3c4dcc was 3c4dcc, checked in by Hans Schönemann <hannes@…>, 18 years ago
*hannes: DOS->UNIX and white space cleanup git-svn-id: file:///usr/local/Singular/svn/trunk@8073 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 28.1 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2version="\$Id: matrix.lib,v 1.29 2005-05-06 14:38:45 hannes Exp \$";
3category="Linear Algebra";
4info="
5LIBRARY:  matrix.lib    Elementary Matrix Operations
6
7PROCEDURES:
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 tensor(A,B);           matrix, tensor product of matrices A nd B
21 unitmat(n);            unit square matrix of size n
22 gauss_col(A);          transform a matrix into col-reduced Gauss normal form
23 gauss_row(A);          transform a matrix into row-reduced Gauss normal form
24 addcol(A,c1,p,c2);     add p*(c1-th col) to c2-th column of matrix A, p poly
25 addrow(A,r1,p,r2);     add p*(r1-th row) to r2-th row of matrix A, p poly
26 multcol(A,c,p);        multiply c-th column of A with poly p
27 multrow(A,r,p);        multiply r-th row of A with poly p
28 permcol(A,i,j);        permute i-th and j-th columns
29 permrow(A,i,j);        permute i-th and j-th rows
30 rowred(A[,any]);       reduction of matrix A with elementary row-operations
31 colred(A[,any]);       reduction of matrix A with elementary col-operations
32 rm_unitrow(A);         remove unit rows and associated columns of A
33 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      if ( size(m) == 0)
57      { intmat B; }
58      else
59      { intmat B[nrows(A)][size(m)]; }
60      int i,j;
61      for( i=1; i<=ncols(A); i=i+1 )
62      {
63         if( m[i]!=[0] )
64         {
65            j=j+1;
66            B[1..nrows(A),j]=A[1..nrows(A),i];
67         }
68      }
69      if( defined(C) ) { return(intvec(B)); }
70      return(B);
71    }
72   return(simplify(A,2));
73}
74example
75{ "EXAMPLE:"; echo = 2;
76   ring r=0,(x,y,z),ds;
77   matrix A[3][4]=1,0,3,0,x,0,z,0,x2,0,z2,0;
78   print(A);
79   print(compress(A));
80   module m=module(A); show(m);
81   show(compress(m));
82   intmat B[3][4]=1,0,3,0,4,0,5,0,6,0,7,0;
83   compress(B);
84   intvec C=0,0,1,2,0,3;
85   compress(C);
86}
87///////////////////////////////////////////////////////////////////////////////
88proc concat (list #)
89"USAGE:   concat(A1,A2,..); A1,A2,... matrices
90RETURN:  matrix, concatenation of A1,A2,.... Number of rows of result matrix
91         is max(nrows(A1),nrows(A2),...)
92EXAMPLE: example concat; shows an example
93"
94{
95   int i;
96   module B=module(#[1]);
97   for( i=2; i<=size(#); i=i+1 ) { B=B,module(#[i]); }
98   return(matrix(B));
99}
100example
101{ "EXAMPLE:"; echo = 2;
102   ring r=0,(x,y,z),ds;
103   matrix A[3][3]=1,2,3,x,y,z,x2,y2,z2;
104   matrix B[2][2]=1,0,2,0; matrix C[1][4]=4,5,x,y;
105   print(A);
106   print(B);
107   print(C);
108   print(concat(A,B,C));
109}
110///////////////////////////////////////////////////////////////////////////////
111
112proc diag (list #)
113"USAGE:   diag(p,n); p poly, n integer
114         diag(A);   A matrix
115RETURN:  diag(p,n): diagonal matrix, p times unit matrix of size n.
116@*       diag(A)  : n*m x n*m diagonal matrix with entries all the entries of
117                    the nxm matrix A, taken from the 1st row, 2nd row etc of A
118EXAMPLE: example diag; shows an example
119"
120{
121   if( size(#)==2 ) { return(matrix(#[1]*freemodule(#[2]))); }
122   if( size(#)==1 )
123   {
124      int i; ideal id=#[1];
125      int n=ncols(id); matrix A[n][n];
126      for( i=1; i<=n; i=i+1 ) { A[i,i]=id[i]; }
127   }
128   return(A);
129}
130example
131{ "EXAMPLE:"; echo = 2;
132   ring r = 0,(x,y,z),ds;
133   print(diag(xy,4));
134   matrix A[3][2] = 1,2,3,4,5,6;
135   print(A);
136   print(diag(A));
137}
138///////////////////////////////////////////////////////////////////////////////
139
140proc dsum (list #)
141"USAGE:   dsum(A1,A2,..); A1,A2,... matrices
142RETURN:  matrix, direct sum of A1,A2,...
143EXAMPLE: example dsum; shows an example
144"
145{
146   int i,N,a;
147   list L;
148   for( i=1; i<=size(#); i=i+1 ) { N=N+nrows(#[i]); }
149   for( i=1; i<=size(#); i=i+1 )
150   {
151      matrix B[N][ncols(#[i])];
152      B[a+1..a+nrows(#[i]),1..ncols(#[i])]=#[i];
153      a=a+nrows(#[i]);
154      L[i]=B; kill B;
155   }
156   return(concat(L));
157}
158example
159{ "EXAMPLE:"; echo = 2;
160   ring r = 0,(x,y,z),ds;
161   matrix A[3][3] = 1,2,3,4,5,6,7,8,9;
162   matrix B[2][2] = 1,x,y,z;
163   print(A);
164   print(B);
165   print(dsum(A,B));
166}
167///////////////////////////////////////////////////////////////////////////////
168
169proc flatten (matrix A)
170"USAGE:   flatten(A); A matrix
171RETURN:  ideal, generated by all entries from A
172EXAMPLE: example flatten; shows an example
173"
174{
175   return(ideal(A));
176}
177example
178{ "EXAMPLE:"; echo = 2;
179   ring r = 0,(x,y,z),ds;
180   matrix A[2][3] = 1,2,x,y,z,7;
181   print(A);
182   flatten(A);
183}
184///////////////////////////////////////////////////////////////////////////////
185
186proc genericmat (int n,int m,list #)
187"USAGE:   genericmat(n,m[,id]);  n,m=integers, id=ideal
188RETURN:  nxm matrix, with entries from id.
189NOTE:    if id has less than nxm elements, the matrix is filled with 0's,
190         (default: id=maxideal(1)).
191         genericmat(n,m); creates the generic nxm matrix
192EXAMPLE: example genericmat; shows an example
193"
194{
195   if( size(#)==0 ) { ideal id=maxideal(1); }
196   if( size(#)==1 ) { ideal id=#[1]; }
197   if( size(#)>=2 ) { "// give 3 arguments, 3-rd argument must be an ideal"; }
198   matrix B[n][m]=id;
199   return(B);
200}
201example
202{ "EXAMPLE:"; echo = 2;
203   ring R = 0,x(1..16),lp;
204   print(genericmat(3,3));      // the generic 3x3 matrix
205   ring R1 = 0,(a,b,c,d),dp;
206   matrix A = genericmat(3,4,maxideal(1)^3);
207   print(A);
208   int n,m = 3,2;
209   ideal i = ideal(randommat(1,n*m,maxideal(1),9));
210   print(genericmat(n,m,i));    // matrix of generic linear forms
211}
212///////////////////////////////////////////////////////////////////////////////
213
214proc is_complex (list c)
215"USAGE:   is_complex(c); c = list of size-compatible modules or matrices
216RETURN:  1 if c[i]*c[i+1]=0 for all i, 0 if not, hence checking whether the
217         list of matrices forms a complex.
218NOTE:    Ideals are treated internally as 1-line matrices.
219         If printlevel > 0, the position where c is not a complex is shown.
220EXAMPLE: example is_complex; shows an example
221"
222{
223   int i;
224   module @test;
225   for( i=1; i<=size(c)-1; i=i+1 )
226   {
227      c[i]=matrix(c[i]); c[i+1]=matrix(c[i+1]);
228      @test=c[i]*c[i+1];
229      if (size(@test)!=0)
230      {
231        dbprint(printlevel-voice+2,"// not a complex at position " +string(i));
232         return(0);
233      }
234   }
235   return(1);
236}
237example
238{ "EXAMPLE:";   echo = 2;
239   ring r  = 32003,(x,y,z),ds;
240   ideal i = x4+y5+z6,xyz,yx2+xz2+zy7;
241   list L  = nres(i,0);
242   is_complex(L);
243   L[4]    = matrix(i);
244   is_complex(L);
245}
246///////////////////////////////////////////////////////////////////////////////
247
248proc outer (matrix A, matrix B)
249"USAGE:   outer(A,B); A,B matrices
250RETURN:  matrix, outer (tensor) product of A and B
251EXAMPLE: example outer; shows an example
252"
253{
254   int i,j; list L;
255   int triv = nrows(B)*ncols(B);
256   if( triv==1 )
257   {
258     return(B[1,1]*A);
259   }
260   else
261   {
262     int N = nrows(A)*nrows(B);
263     matrix C[N][ncols(B)];
264     for( i=1; i<=ncols(A); i=i+1 )
265     {
266       for( j=1; j<=nrows(A); j=j+1 )
267       {
268          C[(j-1)*nrows(B)+1..j*nrows(B),1..ncols(B)]=A[j,i]*B;
269       }
270       L[i]=C;
271     }
272     return(concat(L));
273   }
274}
275example
276{ "EXAMPLE:"; echo = 2;
277   ring r=32003,(x,y,z),ds;
278   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
279   matrix B[2][2]=x,y,0,z;
280   print(A);
281   print(B);
282   print(outer(A,B));
283}
284///////////////////////////////////////////////////////////////////////////////
285
286proc power ( A, int n)
287"USAGE:   power(A,n);  A a square-matrix of type intmat or matrix, n=integer
288RETURN:  intmat resp. matrix, the n-th power of A
289NOTE:    for A=intmat and big n the result may be wrong because of int overflow
290EXAMPLE: example power; shows an example
291"
292{
293//---------------------------- type checking ----------------------------------
294   if( typeof(A)!="matrix" and typeof(A)!="intmat" )
295   {
296      "// no matrix or intmat!";
297      return (A);
298   }
299   if( ncols(A) != nrows(A) )
300   {
301      "// not a suare matrix!";
302      return();
303   }
304//---------------------------- trivial cases ----------------------------------
305   int ii;
306   if( n <= 0 )
307   {
308      if( typeof(A)=="matrix" )
309      {
310         return (unitmat(nrows(A)));
311      }
312      if( typeof(A)=="intmat" )
313      {
314         intmat B[nrows(A)][nrows(A)];
315         for( ii=1; ii<=nrows(A); ii++ )
316         {
317            B[ii,ii] = 1;
318         }
319         return (B);
320      }
321   }
322   if( n == 1 ) { return (A); }
323//---------------------------- sub procedure ----------------------------------
324   proc matpow (A, int n)
325   {
326      def B = A*A;
327      int ii= 2;
328      int jj= 4;
329      while( jj <= n )
330      {
331         B=B*B;
332         ii=jj;
333         jj=2*jj;
334      }
335      return(B,n-ii);
336   }
337//----------------------------- main program ----------------------------------
338   list L = matpow(A,n);
339   def B  = L[1];
340   ii     = L[2];
341   while( ii>=2 )
342   {
343      L = matpow(A,ii);
344      B = B*L[1];
345      ii= L[2];
346   }
347   if( ii == 0) { return(B); }
348   if( ii == 1) { return(A*B); }
349}
350example
351{ "EXAMPLE:"; echo = 2;
352   intmat A[3][3]=1,2,3,4,5,6,7,8,9;
353   print(power(A,3));"";
354   ring r=0,(x,y,z),dp;
355   matrix B[3][3]=0,x,y,z,0,0,y,z,0;
356   print(power(B,3));"";
357}
358///////////////////////////////////////////////////////////////////////////////
359
360proc skewmat (int n, list #)
361"USAGE:   skewmat(n[,id]);  n integer, id ideal
362RETURN:  skew-symmetric nxn matrix, with entries from id
363         (default: id=maxideal(1))
364         skewmat(n); creates the generic skew-symmetric matrix
365NOTE:    if id has less than n*(n-1)/2 elements, the matrix is
366         filled with 0's,
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   ring R1 = 0,(a,b,c),dp;
389   matrix A=skewmat(4,maxideal(1)^2);
390   print(A);
391   int n=3;
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
401         and columns 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   ring R1 = 0,(a,b,c),dp;
447   matrix A=symmat(4,maxideal(1)^3);
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   if (ncols(A)==0)
463   {
464     int q=nrows(A)*nrows(B);
465     matrix D[q][0];
466     return(D);
467   }
468
469   int i,j;
470   matrix C,D;
471   for( i=1; i<=nrows(A); i++ )
472   {
473     C = A[i,1]*B;
474     for( j=2; j<=ncols(A); j++ )
475     {
476       C = concat(C,A[i,j]*B);
477     }
478     D = concat(D,transpose(C));
479   }
480   D = transpose(D);
481   return(submat(D,2..nrows(D),1..ncols(D)));
482}
483example
484{ "EXAMPLE:"; echo = 2;
485   ring r=32003,(x,y,z),(c,ds);
486   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
487   matrix B[2][2]=x,y,0,z;
488   print(A);
489   print(B);
490   print(tensor(A,B));
491}
492///////////////////////////////////////////////////////////////////////////////
493
494proc unitmat (int n)
495"USAGE:   unitmat(n);  n integer >= 0
496RETURN:  nxn unit matrix
497NOTE:    needs a basering, diagonal entries are numbers (=1) in the basering
498EXAMPLE: example unitmat; shows an example
499"
500{
501   return(matrix(freemodule(n)));
502}
503example
504{ "EXAMPLE:"; echo = 2;
505   ring r=32003,(x,y,z),lp;
506   print(xyz*unitmat(4));
507   print(unitmat(5));
508}
509///////////////////////////////////////////////////////////////////////////////
510
511proc gauss_col (matrix A, list #)
512"USAGE:   gauss_col(A[,e]); A a matrix, e any type
513RETURN:  - a matrix B, if called with one argument; B is the complete column-
514           reduced upper-triangular normal form of A if A is constant,
515           (resp. as far as this is possible if A is a polynomial matrix;
516           no division by polynomials).
517@*       - a list L of two matrices, if called with two arguments;
518           L satisfies L[1] = A * L[2] with L[1] the column-reduced form of A
519           and L[2] the transformation matrix.
520NOTE:    * The procedure just applies interred to A with ordering (C,dp).
521           The transformation matrix is obtained by applying 'lift'.
522           This should be faster than the procedure colred.
523@*       * It should only be used with exact coefficient field (there is no
524           pivoting and rounding error treatment).
525@*       * Parameters are allowed. Hence, if the entries of A are parameters,
526           B is the column-reduced form of A over the rational function field.
528EXAMPLE: example gauss_col; shows an example
529"
530{
531   def R=basering; int u;
532   string mp = string(minpoly);
533   int n = nrows(A);
534   int m = ncols(A);
535   module M = A;
536   intvec v = option(get);
537//------------------------ change ordering if necessary ----------------------
538   if( ordstr(R) != ("C,dp("+string(nvars(R))+")") )
539   {
540     def @R=changeord("C,dp",R); setring @R; u=1;
541     execute("minpoly="+mp+";");
542     matrix A = imap(R,A);
543     module M = A;
544   }
545//------------------------------ start computation ---------------------------
546   option(redSB);
547   M = simplify(interred(M),1);
548   if(size(#) != 0)
549   {
550      module N = lift(A,M);
551   }
552//--------------- reset ring and options and return --------------------------
553   if ( u==1 )
554   {
555      setring R;
556      M=imap(@R,M);
557      if (size(#) != 0)
558      {
559         module N = imap(@R,N);
560      }
561      kill @R;
562   }
563   option(set,v);
564   // M = sort(M,size(M)..1)[1];
565   A = matrix(M,n,m);
566   if (size(#) != 0)
567   {
568     list L= A,matrix(N,m,m);
569     return(L);
570   }
571   return(matrix(M,n,m));
572}
573example
574{ "EXAMPLE:"; echo = 2;
575   ring r=(0,a,b),(A,B,C),dp;
576   matrix m[8][6]=
577   0,    2*C, 0,    0,  0,   0,
578   0,    -4*C,a*A,  0,  0,   0,
579   b*B,  -A,  0,    0,  0,   0,
580   -A,   B,   0,    0,  0,   0,
581   -4*C, 0,   B,    2,  0,   0,
582   2*A,  B,   0,    0,  0,   0,
583   0,    3*B, 0,    0,  2b,  0,
584   0,    AB,  0,    2*A,A,   2a;"";
585   list L=gauss_col(m,1);
586   print(L[1]);
587   print(L[2]);
588
589   ring S=0,x,(c,dp);
590   matrix A[5][4] =
591    3, 1, 1, 1,
592   13, 8, 6,-7,
593   14,10, 6,-7,
594    7, 4, 3,-3,
595    2, 1, 0, 3;
596   print(gauss_col(A));
597}
598///////////////////////////////////////////////////////////////////////////////
599
600proc gauss_row (matrix A, list #)
601"USAGE:  gauss_row(A [,e]); A matrix, e any type
602RETURN: - a matrix B, if called with one argument; B is the complete row-
603          reduced lower-triangular normal form of A if A is constant,
604          (resp. as far as this is possible if A is a polynomial matrix;
605          no division by polynomials).
606@*      - a list L of two matrices, if called with two arguments;
607          L satisfies transpose(L[2])*A=transpose(L[1])
608          with L[1] the row-reduced form of A
609          and L[2] the transformation matrix.
610NOTE:   * This procedure just applies gauss_col to the transposed matrix.
611          The transformation matrix is obtained by applying lift.
612          This should be faster than the procedure rowred.
613@*      * It should only be used with exact coefficient field (there is no
614          pivoting and rounding error treatment).
615@*      * Parameters are allowed. Hence, if the entries of A are parameters,
616          B is the row-reduced form of A over the rational function field.
618EXAMPLE: example gauss_row; shows an example
619"
620{
621   if(size(#) > 0)
622   {
623     list L = gauss_col(transpose(A),1);
624     return(L);
625   }
626   A = gauss_col(transpose(A));
627   return(transpose(A));
628}
629example
630{ "EXAMPLE:"; echo = 2;
631   ring r=(0,a,b),(A,B,C),dp;
632   matrix m[6][8]=
633   0, 0,  b*B, -A,-4C,2A,0, 0,
634   2C,-4C,-A,B, 0,  B, 3B,AB,
635   0,a*A,  0, 0, B,  0, 0, 0,
636   0, 0,  0, 0, 2,  0, 0, 2A,
637   0, 0,  0, 0, 0,  0, 2b, A,
638   0, 0,  0, 0, 0,  0, 0, 2a;"";
639   print(gauss_row(m));"";
640   ring S=0,x,dp;
641   matrix A[4][5] =  3, 1,1,-1,2,
642                    13, 8,6,-7,1,
643                    14,10,6,-7,1,
644                     7, 4,3,-3,3;
645   list L=gauss_row(A,1);
646   print(L[1]);
647   print(L[2]);
648}
649///////////////////////////////////////////////////////////////////////////////
650
651proc addcol (matrix A, int c1, poly p, int c2)
652"USAGE:   addcol(A,c1,p,c2);  A matrix, p poly, c1, c2 positive integers
653RETURN:  matrix,  A being modified by adding p times column c1 to column c2
654EXAMPLE: example addcol; shows an example
655"
656{
657   A[1..nrows(A),c2]=A[1..nrows(A),c2]+p*A[1..nrows(A),c1];
658   return(A);
659}
660example
661{ "EXAMPLE:"; echo = 2;
662   ring r=32003,(x,y,z),lp;
663   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
664   print(A);
666}
667///////////////////////////////////////////////////////////////////////////////
668
669proc addrow (matrix A, int r1, poly p, int r2)
670"USAGE:   addcol(A,r1,p,r2);  A matrix, p poly, r1, r2 positive integers
671RETURN:  matrix,  A being modified by adding p times row r1 to row r2
672EXAMPLE: example addrow; shows an example
673"
674{
675   A[r2,1..ncols(A)]=A[r2,1..ncols(A)]+p*A[r1,1..ncols(A)];
676   return(A);
677}
678example
679{ "EXAMPLE:"; echo = 2;
680   ring r=32003,(x,y,z),lp;
681   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
682   print(A);
684}
685///////////////////////////////////////////////////////////////////////////////
686
687proc multcol (matrix A, int c, poly p)
688"USAGE:   addcol(A,c,p);  A matrix, p poly, c positive integer
689RETURN:  matrix,  A being modified by multiplying column c with p
690EXAMPLE: example multcol; shows an example
691"
692{
693   A[1..nrows(A),c]=p*A[1..nrows(A),c];
694   return(A);
695}
696example
697{ "EXAMPLE:"; echo = 2;
698   ring r=32003,(x,y,z),lp;
699   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
700   print(A);
701   print(multcol(A,2,xy));
702}
703///////////////////////////////////////////////////////////////////////////////
704
705proc multrow (matrix A, int r, poly p)
706"USAGE:   multrow(A,r,p);  A matrix, p poly, r positive integer
707RETURN:  matrix,  A being modified by multiplying row r with p
708EXAMPLE: example multrow; shows an example
709"
710{
711   A[r,1..ncols(A)]=p*A[r,1..ncols(A)];
712   return(A);
713}
714example
715{ "EXAMPLE:"; echo = 2;
716   ring r=32003,(x,y,z),lp;
717   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
718   print(A);
719   print(multrow(A,2,xy));
720}
721///////////////////////////////////////////////////////////////////////////////
722
723proc permcol (matrix A, int c1, int c2)
724"USAGE:   permcol(A,c1,c2);  A matrix, c1,c2 positive integers
725RETURN:  matrix,  A being modified by permuting column c1 and c2
726EXAMPLE: example permcol; shows an example
727"
728{
729   matrix B=A;
730   B[1..nrows(B),c1]=A[1..nrows(A),c2];
731   B[1..nrows(B),c2]=A[1..nrows(A),c1];
732   return(B);
733}
734example
735{ "EXAMPLE:"; echo = 2;
736   ring r=32003,(x,y,z),lp;
737   matrix A[3][3]=1,x,3,4,y,6,7,z,9;
738   print(A);
739   print(permcol(A,2,3));
740}
741///////////////////////////////////////////////////////////////////////////////
742
743proc permrow (matrix A, int r1, int r2)
744"USAGE:   permrow(A,r1,r2);  A matrix, r1,r2 positive integers
745RETURN:  matrix,  A being modified by permuting row r1 and r2
746EXAMPLE: example permrow; shows an example
747"
748{
749   matrix B=A;
750   B[r1,1..ncols(B)]=A[r2,1..ncols(A)];
751   B[r2,1..ncols(B)]=A[r1,1..ncols(A)];
752   return(B);
753}
754example
755{ "EXAMPLE:"; echo = 2;
756   ring r=32003,(x,y,z),lp;
757   matrix A[3][3]=1,2,3,x,y,z,7,8,9;
758   print(A);
759   print(permrow(A,2,1));
760}
761///////////////////////////////////////////////////////////////////////////////
762
763proc rowred (matrix A,list #)
764"USAGE:   rowred(A[,e]);  A matrix, e any type
765RETURN:  - a matrix B, being the row reduced form of A, if rowred is called
766           with one argument.
767           (as far as this is possible over the polynomial ring; no division
768           by polynomials)
769@*       - a list L of two matrices, such that L[1] = L[2] * A with L[1]
770           the row-reduced form of A and L[2] the transformation matrix
771           (if rowred is called with two arguments).
772NOTE:    * This procedure is designed for teaching purposes mainly.
773@*       * The straight forward Gaussian algorithm is implemented in the
774           library (no standard basis computation).
775           The transformation matrix is obtained by concatenating a unit
776           matrix to A. proc gauss_row should be faster.
777@*       * It should only be used with exact coefficient field (there is no
778           pivoting) over the polynomial ring (ordering lp or dp).
779@*       * Parameters are allowed. Hence, if the entries of A are parameters
780           the computation takes place over the field of rational functions.
782EXAMPLE: example rowred; shows an example
783"
784{
785   int m,n=nrows(A),ncols(A);
786   int i,j,k,l,rk;
787   poly p;
788   matrix d[m][n];
789   for (i=1;i<=m;i=i+1)
790   {  for (j=1;j<=n;j=j+1)
791      {  p = A[i,j];
792         if (ord(p)==0)
793         {  if (deg(p)==0) { d[i,j]=p; }
794         }
795      }
796   }
797   matrix b = A;
798   if (size(#) != 0) { b = concat(b,unitmat(m)); }
799      for (l=1;l<=n;l=l+1)
800   {
801      k  = findfirst(ideal(d[l]),rk+1);
802      if (k)
803      {  rk = rk+1;
804         b  = permrow(b,rk,k);
805         p  = b[rk,l];         p = 1/p;
806         b  = multrow(b,rk,p);
807         for (i=1;i<=m;i=i+1)
808         {
809            if (rk-i) { b = addrow(b,rk,-b[i,l],i);}
810         }
811         d = 0;
812         for (i=rk+1;i<=m;i=i+1)
813         {  for (j=l+1;j<=n;j=j+1)
814            {  p = b[i,j];
815               if (ord(p)==0)
816               {  if (deg(p)==0) { d[i,j]=p; }
817               }
818            }
819         }
820
821      }
822   }
823   d = submat(b,1..m,1..n);
824   if (size(#))
825   {
826      list L=d,submat(b,1..m,n+1..n+m);
827      return(L);
828   }
829   return(d);
830}
831example
832{ "EXAMPLE:"; echo = 2;
833   ring r=(0,a,b),(A,B,C),dp;
834   matrix m[6][8]=
835   0, 0,  b*B, -A,-4C,2A,0, 0,
836   2C,-4C,-A,B, 0,  B, 3B,AB,
837   0,a*A,  0, 0, B,  0, 0, 0,
838   0, 0,  0, 0, 2,  0, 0, 2A,
839   0, 0,  0, 0, 0,  0, 2b, A,
840   0, 0,  0, 0, 0,  0, 0, 2a;"";
841   print(rowred(m));"";
842   list L=rowred(m,1);
843   print(L[1]);
844   print(L[2]);
845}
846///////////////////////////////////////////////////////////////////////////////
847
848proc colred (matrix A,list #)
849"USAGE:   colred(A[,e]);  A matrix, e any type
850RETURN:  - a matrix B, being the column reduced form of A, if colred is
851           called with one argument.
852           (as far as this is possible over the polynomial ring;
853           no division by polynomials)
854@*       - a list L of two matrices, such that L[1] = A * L[2] with L[1]
855           the column-reduced form of A and L[2] the transformation matrix
856           (if colred is called with two arguments).
857NOTE:    * This procedure is designed for teaching purposes mainly.
858@*       * It applies rowred to the transposed matrix.
859           proc gauss_col should be faster.
860@*       * It should only be used with exact coefficient field (there is no
861           pivoting) over the polynomial ring (ordering lp or dp).
862@*       * Parameters are allowed. Hence, if the entries of A are parameters
863           the computation takes place over the field of rational functions.
865EXAMPLE: example colred; shows an example
866"
867{
868   A = transpose(A);
869   if (size(#))
870   { list L = rowred(A,1); return(transpose(L[1]),transpose(L[2]));}
871   else
872   { return(transpose(rowred(A)));}
873}
874example
875{ "EXAMPLE:"; echo = 2;
876   ring r=(0,a,b),(A,B,C),dp;
877   matrix m[8][6]=
878   0,    2*C, 0,    0,  0,   0,
879   0,    -4*C,a*A,  0,  0,   0,
880   b*B,  -A,  0,    0,  0,   0,
881   -A,   B,   0,    0,  0,   0,
882   -4*C, 0,   B,    2,  0,   0,
883   2*A,  B,   0,    0,  0,   0,
884   0,    3*B, 0,    0,  2b,  0,
885   0,    AB,  0,    2*A,A,   2a;"";
886   print(colred(m));"";
887   list L=colred(m,1);
888   print(L[1]);
889   print(L[2]);
890}
891//////////////////////////////////////////////////////////////////////////////
892
893static proc findfirst (ideal i,int t)
894{
895   int n,k;
896   for (n=t;n<=ncols(i);n=n+1)
897   {
898      if (i[n]!=0) { k=n;break;}
899   }
900   return(k);
901}
902//////////////////////////////////////////////////////////////////////////////
903
904proc rm_unitcol(matrix A)
905"USAGE:   rm_unitcol(A); A matrix (being row-reduced)
906RETURN:  matrix, obtained from A by deleting unit columns (having just one 1
907         and else 0 as entries) and associated rows
908EXAMPLE: example rm_unitcol; shows an example
909"
910{
911 int l,j;
912 intvec v;
913 for (j=1;j<=ncols(A);j=j+1)
914 {
915    if (gen(l+1)==module(A)[j]) {l=l+1;}
916    else { v=v,j;}
917 }
918 if (size(v)>1)
919    {  v = v[2..size(v)];
920       return(submat(A,l+1..nrows(A),v));
921    }
922 else
923    { return(0);}
924}
925example
926{ "EXAMPLE:"; echo = 2;
927   ring r=0,(A,B,C),dp;
928   matrix m[6][8]=
929   0,  0,    A,   0, 1,0,  0,0,
930   0,  0,  -C2,   0, 0,0,  1,0,
931   0,  0,    0,1/2B, 0,0,  0,1,
932   0,  0,    B,  -A, 0,2A, 0,0,
933   2C,-4C,  -A,   B, 0,B,  0,0,
934   0,  A,    0,   0, 0,0,  0,0;
935   print(rm_unitcol(m));
936}
937//////////////////////////////////////////////////////////////////////////////
938
939proc rm_unitrow (matrix A)
940"USAGE:   rm_unitrow(A); A matrix (being col-reduced)
941RETURN:  matrix, obtained from A by deleting unit rows (having just one 1
942         and else 0 as entries) and associated columns
943EXAMPLE: example rm_unitrow; shows an example
944"
945{
946 int l,j;
947 intvec v;
948 module M = transpose(A);
949 for (j=1; j <= nrows(A); j=j+1)
950 {
951    if (gen(l+1) == M[j]) { l=l+1; }
952    else { v=v,j; }
953 }
954 if (size(v) > 1)
955    {  v = v[2..size(v)];
956       return(submat(A,v,l+1..ncols(A)));
957    }
958 else
959    { return(0);}
960}
961example
962{ "EXAMPLE:"; echo = 2;
963   ring r=0,(A,B,C),dp;
964   matrix m[8][6]=
965   0,0,  0,   0, 2C, 0,
966   0,0,  0,   0, -4C,A,
967   A,-C2,0,   B, -A, 0,
968   0,0,  1/2B,-A,B,  0,
969   1,0,  0,   0, 0,  0,
970   0,0,  0,   2A,B,  0,
971   0,1,  0,   0, 0,  0,
972   0,0,  1,   0, 0,  0;
973   print(rm_unitrow(m));
974}
975//////////////////////////////////////////////////////////////////////////////
977{
978  int i,j;
979  int n=nrows(M);
980  int m=ncols(M);
981  matrix B[n][m];
982  for(i=1;i<=n;i++)
983  {
984     for(j=1;j<=m;j++)
985     {
986        B[n-i+1,m-j+1]=M[i,j];
987     }
988  }
989  return(B);
990}
991//////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.