source: git/Singular/LIB/matrix.lib @ 49f94f

spielwiese
Last change on this file since 49f94f was 49f94f, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: typos reported by gorzel git-svn-id: file:///usr/local/Singular/svn/trunk@11088 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 38.6 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: matrix.lib,v 1.43 2008-10-01 15:29:23 Singular 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 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
35 headStand(A);          A[n-i+1,m-j+1]:=A[i,j]
36 symmetricBasis(n,k);   ring, with a basis of a k-th symmetric power of an n-dim v.space
37 symmetricPower(A,k);   module, k-th symmetric power of a module/matrix A
38 exteriorBasis(n,k);    ring, with a basis of a k-th exterior power of an n-dim v.space
39 exteriorPower(A,k);    module, 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.
533SEE ALSO:  colred
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.
623SEE ALSO: rowred
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);
672   print(addcol(A,1,xy,2));
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);
691   print(addrow(A,1,xy,3));
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.
795SEE ALSO:  gauss_row
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.
880SEE ALSO:  gauss_col
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         positife 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//////////////////////////////////////////////////////////////////////////////
1046proc headStand(matrix M)
1047"USAGE:   headStand(M);  M matrix
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);
1072   print(headStand(M));
1073}
1074//////////////////////////////////////////////////////////////////////////////
1075
1076// Symmetric/Exterior powers, thanks to Oleksandr Iena for his persistence ;-)
1077
1078proc symmetricBasis(int n, int k)
1079"USAGE:    symmetricBasis(n, k); n int, k int
1080RETURN:   ring: the basis of the k-th symmetric power of a n-dim vector space
1081NOTE:     The output consists of a polynomial ring in n variables, together with the ideal called I.
1082          The ideal I is the basis of the k-th symmetric power of a n-dim vector space (ordered lexicographically).
1083          The char. of the returned ring is the same as of current basis ring or zero.
1084SEE ALSO: exteriorBasis
1085KEYWORDS: symmetric basis
1086EXAMPLE:  example symmetricBasis; shows an example"
1087{
1088  int p = 0;
1089  if( nameof(basering) != "basering" )
1090  {
1091    p = char(basering);
1092  }
1093
1094  ring T = (p), (e(1..n)), dp;
1095  ideal I = simplify( NF(maxideal(k), std(0)), 1 + 2 + 8 );
1096  int N = ncols(I);
1097  I = sort(I)[1]; // lex
1098  I = I[N..1];
1099
1100  export I;
1101  return(T);
1102}
1103example
1104{ "EXAMPLE:"; echo = 2;
1105
1106def r = symmetricBasis(2, 3); setring r; r;
1107I; // basis of 3rd sym. power of a 2-dim v.s.
1108kill r;
1109
1110ring r = (0),(a, b, c, d), dp; r;
1111def g = symmetricBasis(3, 2); setring g; g; I;
1112
1113kill g, r;
1114
1115ring r = (32003),(a, b, c, d), dp; r;
1116def g = symmetricBasis(4, 2); setring g; g; I;
1117}
1118
1119//////////////////////////////////////////////////////////////////////////////
1120
1121proc exteriorBasis(int n, int k)
1122"USAGE:    exteriorBasis(n, k); n int, k int
1123RETURN:   ring: the basis of the k-th exterior power of a n-dim vector space
1124NOTE:     The output consists of a polynomial ring in n variables, together with the ideal called I.
1125          The ideal I is the basis of the k-th exterior power of a n-dim vector space (ordered lexicographically).
1126          The char. of the returned ring is the same as of current basis ring or zero.
1127
1128SEE ALSO: symmetricBasis
1129KEYWORDS: exterior basis
1130EXAMPLE:  example exteriorBasis; shows an example"
1131{
1132  int p = 0;
1133  if( nameof(basering) != "basering" )
1134  {
1135    p = char(basering);
1136  }
1137   
1138  ring S = (p), (e(1..n)), dp;
1139  def T = SuperCommutative(); setring T;
1140  ideal I = simplify( NF(maxideal(k), std(0)), 1 + 2 + 8 );
1141  int N = ncols(I);
1142  I = sort(I)[1];
1143  I = I[N..1 ];
1144
1145  export I;
1146  return(T);
1147}
1148example
1149{ "EXAMPLE:"; echo = 2;
1150
1151def r = exteriorBasis(2, 3); setring r; r;
1152"Basis: ", I; // basis of 3rd sym. power of a 2-dim v.s.
1153kill r;
1154
1155ring r = (0),(a, b, c, d), dp; r;
1156def g = exteriorBasis(3, 2); setring g; g; "Basis: ", I;
1157
1158kill g, r;
1159
1160ring r = (32003),(a, b, c, d), dp; r;
1161def g = exteriorBasis(4, 2); setring g; g; "Basis: ", I;
1162}
1163
1164//////////////////////////////////////////////////////////////////////////////
1165
1166proc symmetricPower(module A, int k)
1167"USAGE:    symmetricPower(A, k); A module, k int
1168RETURN:   module: the k-th symmetric power of A
1169NOTE:     the chosen bases (ordered lexicographically) and
1170          temporary data  may be shown if printlevel is big enough
1171SEE ALSO: exteriorPower
1172KEYWORDS: symmetric power
1173EXAMPLE:  example symmetricPower; shows an example"
1174{
1175  int p = printlevel - voice + 2;
1176  def save = basering;
1177
1178  int M = nrows(A);
1179  int N = ncols(A);
1180
1181  int i, j;
1182
1183  /////////////////////////////////////////////////////////////////
1184  def Tn = symmetricBasis(N, k); setring Tn;
1185  ideal Ink = I;
1186  dbprint(p-3, "Temporary Source Ring", basering);
1187  dbprint(p, "S^k(Source Basis)", Ink);
1188
1189  /////////////////////////////////////////////////////////////////
1190  def Tm = symmetricBasis(M, k); setring Tm;
1191  ideal Imk = I;
1192  dbprint(p-3, "Temporary Image Ring", basering);
1193  dbprint(p, "S^k(Image Basis)", Imk);
1194
1195  /////////////////////////////////////////////////////////////////
1196  def Rm = save + Tm;
1197  setring Rm;
1198  dbprint(p-2, "Temporary Working Ring", Rm);
1199
1200  module A = imap(save, A);
1201
1202  ideal B; poly t;
1203
1204  for( i = N; i > 0; i-- )
1205  {
1206    t = 0;
1207    for( j = M; j > 0; j-- )
1208    {
1209      t = t + A[i][j] * var( nvars(save) + j); // tensor product!!!
1210    }
1211
1212    B[i] = t;
1213  }
1214
1215  dbprint(p-1, "Matrix of simgle images", B);
1216
1217  map TMap = Tn, B; ideal C = TMap(Ink); // apply S^k A to Source basis vectors... (Ink)
1218  C = NF(C, std(0));
1219
1220  dbprint(p-1, "Image Matrix: ", C);
1221
1222  // and write it in Image basis (Imk)
1223
1224  ideal Imk = imap(Tm, Imk);
1225
1226  module D; poly lm; vector tt;
1227
1228  for( i = ncols(C); i > 0; i-- )
1229  {
1230    t = C[i];
1231    tt = 0;
1232
1233    while( t != 0 )
1234    {
1235      lm = leadmonom(t);
1236      //    lm;
1237      for( j = ncols(Imk); j > 0; j-- )
1238      {
1239        if( lm / Imk[j] != 0 )
1240        {
1241          tt = tt + (lead(t) / Imk[j]) * gen(j);
1242          break;
1243        }
1244      }
1245      t = t - lead(t);
1246    }
1247
1248    D[i] = tt;
1249//    tt;
1250  }
1251
1252  // D; print(D);
1253
1254
1255  setring save;
1256
1257  // basering;
1258
1259  module AA = imap(Rm, D);
1260
1261//  nrows(AA) = Nk; // ????
1262//  ncols(AA) = Mk;
1263
1264//  Nk, Mk;
1265
1266//  AA[Mk] = AA[Mk] * 1;
1267
1268  return(AA);
1269}
1270example
1271{ "EXAMPLE:"; echo = 2;
1272int save = printlevel; printlevel = 1;
1273
1274ring r = (0),(a, b, c, d), dp; r;
1275
1276module A = a*gen(1), b * gen(1), c*gen(1), d * gen(1);
1277
1278print(symmetricPower(A, 2));
1279print(symmetricPower(A, 3));
1280
1281module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2);
1282
1283print(symmetricPower(B, 2));
1284print(symmetricPower(B, 3));
1285
1286kill r;
1287ring r = (0),(a, b, c, d), dp;def g = SuperCommutative();setring g;
1288
1289module A = a*gen(1), b * gen(1), c*gen(1), d * gen(1);
1290
1291print(symmetricPower(A, 2));
1292print(symmetricPower(A, 3));
1293
1294module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2);
1295
1296print(symmetricPower(B, 2));
1297print(symmetricPower(B, 3));
1298
1299kill r;
1300
1301printlevel = save;
1302}
1303
1304//////////////////////////////////////////////////////////////////////////////
1305
1306proc exteriorPower(module A, int k)
1307"USAGE:    exteriorPower(A, k); A module, k int
1308RETURN:   module: the k-th exterior power of A
1309NOTE:     the chosen bases and temporary data
1310          may be shown if printlevel is big enough.
1311          Last rows may be invisible if zero.
1312SEE ALSO: symmetricPower
1313KEYWORDS: exterior power
1314EXAMPLE:  example exteriorPower; shows an example"
1315{
1316  int p = printlevel - voice + 2;
1317  def save = basering;
1318
1319  int M = nrows(A);
1320  int N = ncols(A);
1321
1322  int i, j;
1323
1324  /////////////////////////////////////////////////////////////////
1325  def Tn = exteriorBasis(N, k); setring Tn;
1326  ideal Ink = I;
1327  dbprint(p-3, "Temporary Source Ring", basering);
1328  dbprint(p, "E^k(Source Basis)", Ink);
1329
1330  /////////////////////////////////////////////////////////////////
1331  def Tm = exteriorBasis(M, k); setring Tm;
1332  ideal Imk = I;
1333  dbprint(p-3, "Temporary Image Ring", basering);
1334  dbprint(p, "E^k(Image Basis)", Imk);
1335
1336
1337  def Rm = save + Tm;
1338  setring Rm;
1339  dbprint(p-2, "Temporary Working Ring", Rm);
1340
1341  module A = imap(save, A);
1342
1343  ideal B; poly t;
1344
1345  for( i = N; i > 0; i-- )
1346  {
1347    t = 0;
1348    for( j = M; j > 0; j-- )
1349    {
1350      t = t + A[i][j] * var( nvars(save) + j); // tensor product!!!
1351    }
1352
1353    B[i] = t;
1354  }
1355
1356  dbprint(p-1, "Matrix of simgle images", B);
1357
1358  map TMap = Tn, B; ideal C = TMap(Ink); // apply S^k A to Source basis vectors... (Ink)
1359  C = NF(C, std(0));
1360
1361  dbprint(p-1, "Image Matrix: ", C);
1362
1363  // and write it in Image basis (Imk)
1364
1365  ideal Imk = imap(Tm, Imk);
1366
1367  module D; poly lm; vector tt;
1368
1369  for( i = ncols(C); i > 0; i-- )
1370  {
1371    t = C[i];
1372    tt = 0;
1373
1374    while( t != 0 )
1375    {
1376      lm = leadmonom(t);
1377      //    lm;
1378      for( j = ncols(Imk); j > 0; j-- )
1379      {
1380        if( lm / Imk[j] != 0 )
1381        {
1382          tt = tt + (lead(t) / Imk[j]) * gen(j);
1383          break;
1384        }
1385      }
1386      t = t - lead(t);
1387    }
1388
1389    D[i] = tt;
1390//    tt;
1391  }
1392
1393  // D; print(D);
1394
1395
1396  setring save;
1397
1398  // basering;
1399
1400  module AA = imap(Rm, D);
1401
1402//  nrows(AA) = Nk; // ????
1403//  ncols(AA) = Mk;
1404
1405//  Nk, Mk;
1406
1407//  AA[Mk] = AA[Mk] * 1;
1408
1409  return(AA);
1410}
1411example
1412{ "EXAMPLE:"; echo = 2;
1413  int save = printlevel; printlevel = 1;
1414
1415  ring r = (0),(a, b, c, d), dp; r;
1416
1417  module A = a*gen(1), b * gen(1), c*gen(1), d * gen(1);
1418
1419  print(exteriorPower(A, 2));
1420  print(exteriorPower(A, 3));
1421
1422  module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2);
1423
1424  print(exteriorPower(B, 2));
1425  print(exteriorPower(B, 3));
1426
1427  kill r;
1428  ring r = (0),(a, b, c, d), dp;def g = SuperCommutative();setring g;
1429
1430  module A = a*gen(1), b * gen(1), c*gen(1), d * gen(1);
1431
1432  print(exteriorPower(A, 2));
1433  print(exteriorPower(A, 3));
1434
1435  module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2);
1436
1437  print(exteriorPower(B, 2));
1438  print(exteriorPower(B, 3));
1439
1440  kill r;
1441
1442  printlevel = save;
1443}
1444
1445//////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.