# source:git/Singular/LIB/matrix.lib@341696

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