source: git/Singular/LIB/linalg.lib @ 46de44

spielwiese
Last change on this file since 46de44 was 46de44, checked in by Mathias Schulze <mschulze@…>, 22 years ago
*mschulze: ring change in eigenvals not needed git-svn-id: file:///usr/local/Singular/svn/trunk@5885 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 38.6 KB
Line 
1//GMG last modified: 04/25/2000
2//////////////////////////////////////////////////////////////////////////////
3version="$Id: linalg.lib,v 1.24 2002-02-16 13:50:19 mschulze Exp $";
4category="Linear Algebra";
5info="
6LIBRARY:  linalg.lib  Algorithmic Linear Algebra
7AUTHORS:  Ivor Saynisch (ivs@math.tu-cottbus.de)
8@*        Mathias Schulze (mschulze@mathematik.uni-kl.de)
9
10PROCEDURES:
11 inverse(A);          matrix, the inverse of A
12 inverse_B(A);        list(matrix Inv,poly p),Inv*A=p*En ( using busadj(A) )
13 inverse_L(A);        list(matrix Inv,poly p),Inv*A=p*En ( using lift )
14 sym_gauss(A);        symmetric gaussian algorithm
15 orthogonalize(A);    Gram-Schmidt orthogonalization
16 diag_test(A);        test whether A can be diagnolized
17 busadj(A);           coefficients of Adj(E*t-A) and coefficients of det(E*t-A)
18 charpoly(A,v);       characteristic polynomial of A ( using busadj(A) )
19 adjoint(A);          adjoint of A ( using busadj(A) )
20 det_B(A);            determinant of A ( using busadj(A) )
21 gaussred(A);         gaussian reduction: P*A=U*S, S a row reduced form of A
22 gaussred_pivot(A);   gaussian reduction: P*A=U*S, uses row pivoting
23 gauss_nf(A);         gaussian normal form of A
24 mat_rk(A);           rank of constant matrix A
25 U_D_O(A);            P*A=U*D*O, P,D,U,O=permutaion,diag,lower-,upper-triang
26 pos_def(A,i);        test symmetric matrix for positive definiteness
27 eigenvals(M);        eigenvalues and multiplicities of M
28 jordan(M);           Jordan data of M
29 jordanbasis(M);      Jordan basis and weight filtration of M
30 jordanmatrix(e,s,m); Jordan matrix with Jordan data (e,s,m)
31 jordannf(M);         Jordan normal form of M
32";
33
34LIB "matrix.lib";
35LIB "ring.lib";
36LIB "elim.lib";
37//////////////////////////////////////////////////////////////////////////////
38// help functions
39//////////////////////////////////////////////////////////////////////////////
40static proc abs(poly c)
41"RETURN: absolut value of c, c must be constants"
42{
43  if(c>=0){ return(c);}
44  else{ return(-c);}
45}
46
47static proc const_mat(matrix A)
48"RETURN:   1 (0) if A is (is not) a constant matrix"
49{
50  int i;
51  int n=ncols(A);
52  def BR=basering;
53  changeord("@R","dp,c",BR);
54  matrix A=fetch(BR,A);
55  for(i=1;i<=n;i=i+1){
56    if(deg(lead(A)[i])>=1){
57      //"input is not a constant matrix";
58      kill @R;
59      setring BR;
60      return(0);
61    }
62  }
63  kill @R;
64  setring BR;
65  return(1);
66}
67//////////////////////////////////////////////////////////////////////////////
68static proc red(matrix A,int i,int j)
69"USAGE:    red(A,i,j);  A = constant matrix
70          reduces column j with respect to A[i,i] and column i
71          reduces row j with respect to A[i,i] and row i
72RETURN:   matrix
73"
74{
75  module m=module(A);
76
77  if(A[i,i]==0){
78    m[i]=m[i]+m[j];
79    m=module(transpose(matrix(m)));
80    m[i]=m[i]+m[j];
81    m=module(transpose(matrix(m)));
82  }
83
84  A=matrix(m);
85  m[j]=m[j]-(A[i,j]/A[i,i])*m[i];
86  m=module(transpose(matrix(m)));
87  m[j]=m[j]-(A[i,j]/A[i,i])*m[i];
88  m=module(transpose(matrix(m)));
89
90  return(matrix(m));
91}
92
93//////////////////////////////////////////////////////////////////////////////
94proc inner_product(vector v1,vector v2)
95"RETURN:   inner product <v1,v2> "
96{
97  int k;
98  if (nrows(v2)>nrows(v1)) { k=nrows(v2); } else { k=nrows(v1); }
99  return ((transpose(matrix(v1,k,1))*matrix(v2,k,1))[1,1]);
100}
101
102/////////////////////////////////////////////////////////////////////////////
103// user functions
104/////////////////////////////////////////////////////////////////////////////
105
106proc inverse(matrix A, list #)
107"USAGE:    inverse(A [,opt]);  A a square matrix, opt integer
108RETURN:
109@format
110          a matrix:
111          - the inverse matrix of A, if A is invertible;
112          - the 1x1 0-matrix if A is not invertible (in the polynomial ring!).
113          There are the following options:
114          - opt=0 or not given: heuristically best option from below
115          - opt=1 : apply std to (transpose(E,A)), ordering (C,dp).
116          - opt=2 : apply interred (transpose(E,A)), ordering (C,dp).
117          - opt=3 : apply lift(A,E), ordering (C,dp).
118@end format
119NOTE:     parameters and minpoly are allowed; opt=2 is only correct for
120          matrices with entries in a field
121SEE ALSO: inverse_B, inverse_L
122EXAMPLE:  example inverse; shows an example
123"
124{
125//--------------------------- initialization and check ------------------------
126   int ii,u,i,opt;
127   matrix invA;
128   int db = printlevel-voice+3;      //used for comments
129   def R=basering;
130   string mp = string(minpoly);
131   int n = nrows(A);
132   module M = A;
133   intvec v = option(get);           //get options to reset it later
134   if ( ncols(A)!=n )
135   {
136     ERROR("// ** no square matrix");
137   }
138//----------------------- choose heurisitically best option ------------------
139// This may change later, depending on improvements of the implemantation
140// at the monent we use if opt=0 or opt not given:
141// opt = 1 (std) for everything
142// opt = 2 (interred) for nothing, NOTE: interred is ok for constant matricea
143// opt = 3 (lift) for nothing
144// NOTE: interred is ok for constant matrices only (Beispiele am Ende der lib)
145   if(size(#) != 0) {opt = #[1];}
146   if(opt == 0)
147   {
148      if(npars(R) == 0)                       //no parameters
149      {
150         if( size(ideal(A-jet(A,0))) == 0 )   //constant matrix
151         {opt = 1;}
152         else
153         {opt = 1;}
154      }
155      else {opt = 1;}
156   }
157//------------------------- change ring if necessary -------------------------
158   if( ordstr(R) != "C,dp(nvars(R))" )
159   {
160     u=1;
161     changeord("@R","C,dp",R);
162     module M = fetch(R,M);
163     execute("minpoly="+mp+";");
164   }
165//----------------------------- opt=3: use lift ------------------------------
166   if( opt==3 )
167   {
168      module D2;
169      D2 = lift(M,freemodule(n));
170      if (size(ideal(D2))==0)
171      {                                               //catch error in lift
172         dbprint(db,"// ** matrix is not invertible");
173         setring R;
174         if (u==1) { kill @R;}
175         return(invA);
176      }
177   }
178//-------------- opt = 1 resp. opt = 2: use std resp. interred --------------
179   if( opt==1 or opt==2 )
180   {
181      option(redSB);
182      module B = freemodule(n),M;
183      if(opt == 2)
184      {
185         module D = interred(transpose(B));
186         D = transpose(simplify(D,1));
187      }
188      if(opt == 1)
189      {
190         module D = std(transpose(B));
191         D = transpose(simplify(D,1));
192      }
193      module D2 = D[1..n];
194      module D1 = D[n+1..2*n];
195//----------------------- check if matrix is invertible ----------------------
196      for (ii=1; ii<=n; ii++)
197      {
198         if ( D1[ii] != gen(ii) )
199         {
200            dbprint(db,"// ** matrix is not invertible");
201            i = 1;
202            break;
203         }
204      }
205   }
206   option(set,v);
207//------------------ return to basering and return result ---------------------
208   if ( u==1 )
209   {
210      setring R;
211      module D2 = fetch(@R,D2);
212      if( opt==1 or opt==2 )
213      { module D1 = fetch(@R,D1);}
214      kill @R;
215   }
216   if( i==1 ) { return(invA); }     //matrix not invetible
217   else { return(matrix(D2)); }     //matrix invertible with inverse D2
218
219}
220example
221{ "EXAMPLE:"; echo = 2;
222  ring r=0,(x,y,z),lp;
223  matrix A[3][3]=
224   1,4,3,
225   1,5,7,
226   0,4,17;
227  print(inverse(A));"";
228  matrix B[3][3]=
229   y+1,  x+y,    y,
230   z,    z+1,    z,
231   y+z+2,x+y+z+2,y+z+1;
232  print(inverse(B));
233  print(B*inverse(B));
234}
235
236//////////////////////////////////////////////////////////////////////////////
237proc sym_gauss(matrix A)
238"USAGE:    sym_gauss(A);  A = symmetric matrix
239RETURN:   matrix, diagonalisation with symmetric gauss algorithm
240EXAMPLE:  example sym_gauss; shows an example"
241{
242  int i,j;
243  int n=nrows(A);
244
245  if (ncols(A)!=n){
246    "// ** input is not a square matrix";;
247    return(A);
248  }
249
250  if(!const_mat(A)){
251    "// ** input is not a constant matrix";
252    return(A);
253  }
254
255  if(deg(std(A-transpose(A))[1])!=-1){
256    "// ** input is not a symmetric matrix";
257    return(A);
258  }
259
260  for(i=1; i<n; i++){
261    for(j=i+1; j<=n; j++){
262      if(A[i,j]!=0){ A=red(A,i,j); }
263    }
264  }
265
266  return(A);
267}
268example
269{"EXAMPLE:"; echo = 2;
270  ring r=0,(x),lp;
271  matrix A[2][2]=1,4,4,15;
272  print(A);
273  print(sym_gauss(A));
274}
275
276//////////////////////////////////////////////////////////////////////////////
277proc orthogonalize(matrix A)
278"USAGE:    orthogonalize(A); A = constant matrix
279RETURN:    matrix, orthogonal basis of the colum space of A
280EXAMPLE:   example orthogonalize; shows an example "
281{
282  int i,j;
283  int n=ncols(A);
284  poly k;
285
286  if(!const_mat(A)){
287    "// ** input is not a constant matrix";
288    matrix B;
289    return(B);
290  }
291
292  module B=module(interred(A));
293
294  for(i=1;i<=n;i=i+1) {
295    for(j=1;j<i;j=j+1) {
296      k=inner_product(B[j],B[j]);
297      if (k==0) { "Error: vector of length zero"; return(matrix(B)); }
298      B[i]=B[i]-(inner_product(B[i],B[j])/k)*B[j];
299    }
300  }
301
302  return(matrix(B));
303}
304example
305{ "EXAMPLE:"; echo = 2;
306  ring r=0,(x),lp;
307  matrix A[4][4]=5,6,12,4,7,3,2,6,12,1,1,2,6,4,2,10;
308  print(A);
309  print(orthogonalize(A));
310}
311
312////////////////////////////////////////////////////////////////////////////
313proc diag_test(matrix A)
314"USAGE:          diag_test(A); A = const square matrix
315RETURN:   int,  1 if A is diagonalisable, 0 if not
316               -1 no statement is possible, since A does not split.
317NOTE:     The test works only for split matrices, i.e if eigenvalues of A
318          are in the ground field.
319          Does not work with parameters (uses factorize,gcd).
320EXAMPLE:  example diag_test; shows an example"
321{
322  int i,j;
323  int n     = nrows(A);
324  string mp = string(minpoly);
325  string cs = charstr(basering);
326  int np=0;
327
328  if(ncols(A) != n) {
329    "// input is not a square matrix";
330    return(-1);
331  }
332
333   if(!const_mat(A)){
334    "// input is not a constant matrix";
335    return(-1);
336  }
337
338  //Parameterring wegen factorize nicht erlaubt
339  for(i=1;i<size(cs);i=i+1){
340    if(cs[i]==","){np=np+1;} //Anzahl der Parameter
341  }
342  if(np>0){
343        "// rings with parameters not allowed";
344        return(-1);
345  }
346
347  //speichern des aktuellen Rings
348  def BR=basering;
349  //setze R[t]
350  execute("ring rt=("+charstr(basering)+"),(@t,"+varstr(basering)+"),lp;");
351  execute("minpoly="+mp+";");
352  matrix A=imap(BR,A);
353
354  intvec z;
355  intvec s;
356  poly X;         //characteristisches Polynom
357  poly dXdt;      //Ableitung von X nach t
358  ideal g;              //ggT(X,dXdt)
359  poly b;              //Komponente der Busadjunkten-Matrix
360  matrix E[n][n]; //Einheitsmatrix
361
362  E=E+1;
363  A=E*@t-A;
364  X=det(A);
365
366  matrix Xfactors=matrix(factorize(X,1));        //zerfaellt die Matrtix ?
367  int nf=ncols(Xfactors);
368
369  for(i=1;i<=nf;i++){
370    if(lead(Xfactors[1,i])>=@t^2){
371      //" matrix does not split";
372      setring BR;
373      return(-1);
374    }
375  }
376
377  dXdt=diff(X,@t);
378  g=std(ideal(gcd(X,dXdt)));
379
380  //Busadjunkte
381  z=2..n;
382  for(i=1;i<=n;i++){
383    s=2..n;
384    for(j=1;j<=n;j++){
385      b=det(submat(A,z,s));
386
387      if(0!=reduce(b,g)){
388        //" matrix not diagonalizable";
389        setring BR;
390        return(0);
391      }
392
393      s[j]=j;
394    }
395    z[i]=i;
396  }
397
398  //"Die Matrix ist diagonalisierbar";
399  setring BR;
400  return(1);
401}
402example
403{ "EXAMPLE:"; echo = 2;
404  ring r=0,(x),dp;
405  matrix A[4][4]=6,0,0,0,0,0,6,0,0,6,0,0,0,0,0,6;
406  print(A);
407  diag_test(A);
408}
409
410//////////////////////////////////////////////////////////////////////////////
411proc busadj(matrix A)
412"USAGE:   busadj(A);  A = square matrix (of size nxn)
413RETURN:  list L:
414@format
415         L[1] contains the (n+1) coefficients of the characteristic
416              polynomial X of A, i.e.
417              X = L[1][1]+..+L[1][k]*t^(k-1)+..+(L[1][n+1])*t^n
418         L[2] contains the n (nxn)-matrices Hk which are the coefficients of
419              the busadjoint bA = adjoint(E*t-A) of A, i.e.
420              bA = (Hn-1)*t^(n-1)+...+Hk*t^k+...+H0,  ( Hk=L[2][k+1] )
421@end format
422EXAMPLE: example busadj; shows an example"
423{
424  int k;
425  int n    = nrows(A);
426  matrix E = unitmat(n);
427  matrix H[n][n];
428  matrix B[n][n];
429  list bA, X, L;
430  poly a;
431
432  if(ncols(A) != n) {
433    "input is not a square matrix";
434    return(L);
435  }
436
437  bA   = E;
438  X[1] = 1;
439  for(k=1; k<n; k++){
440    B  = A*bA[1];              //bA[1] is the last H
441    a  = -trace(B)/k;
442    H  = B+a*E;
443    bA = insert(bA,H);
444    X  = insert(X,a);
445  }
446  B = A*bA[1];
447  a = -trace(B)/n;
448  X = insert(X,a);
449
450  L = insert(L,bA);
451  L = insert(L,X);
452  return(L);
453}
454example
455{ "EXAMPLE"; echo = 2;
456  ring r = 0,(t,x),lp;
457  matrix A[2][2] = 1,x2,x,x2+3x;
458  print(A);
459  list L = busadj(A);
460  poly X = L[1][1]+L[1][2]*t+L[1][3]*t2; X;
461  matrix bA[2][2] = L[2][1]+L[2][2]*t;
462  print(bA);               //the busadjoint of A;
463  print(bA*(t*unitmat(2)-A));
464}
465
466//////////////////////////////////////////////////////////////////////////////
467proc charpoly(matrix A, list #)
468"USAGE:   charpoly(A[,v]); A square matrix, v string, name of a variable
469RETURN:  poly, the characteristic polynomial det(E*v-A)
470         (default: v=name of last variable)
471NOTE:    A must be independent of the variable v. The computation uses det.
472         If printlevel>0, det(E*v-A) is diplayed recursively.
473EXAMPLE: example charpoly; shows an example"
474{
475  int n = nrows(A);
476  int z = nvars(basering);
477  int i,j;
478  string v;
479  poly X;
480  if(ncols(A) != n)
481  {
482    "// input is not a square matrix";
483    return(X);
484  }
485  //---------------------- test for correct variable -------------------------
486  if( size(#)==0 ){
487    #[1] = varstr(z);
488  }
489  if( typeof(#[1]) == "string") { v = #[1]; }
490  else
491  {
492    "// 2nd argument must be a name of a variable not contained in the matrix";
493    return(X);
494  }
495  j=-1;
496  for(i=1; i<=z; i++)
497  {
498    if(varstr(i)==v){j=i;}
499  }
500  if(j==-1)
501  {
502    "// "+v+" is not a variable in the basering";
503    return(X);
504  }
505  if ( size(select1(module(A),j)) != 0 )
506  {
507    "// matrix must not contain the variable "+v;
508    "// change to a ring with an extra variable, if necessary."
509    return(X);
510  }
511  //-------------------------- compute charpoly ------------------------------
512  X = det(var(j)*unitmat(n)-A);
513  if( printlevel-voice+2 >0) { showrecursive(X,var(j));}
514  return(X);
515}
516example
517{ "EXAMPLE"; echo=2;
518  ring r=0,(x,t),dp;
519  matrix A[3][3]=1,x2,x,x2,6,4,x,4,1;
520  print(A);
521  charpoly(A,"t");
522}
523
524//////////////////////////////////////////////////////////////////////////////
525proc charpoly_B(matrix A, list #)
526"USAGE:   charpoly_B(A[,v]); A square matrix, v string, name of a variable
527RETURN:  poly, the characteristic polynomial det(E*v-A)
528         (default: v=name of last variable)
529NOTE:    A must be constant in the variable v. The computation uses busadj(A).
530EXAMPLE: example charpoly_B; shows an example"
531{
532  int i,j;
533  string s,v;
534  list L;
535  int n     = nrows(A);
536  poly X    = 0;
537  def BR    = basering;
538  string mp = string(minpoly);
539
540  if(ncols(A) != n){
541    "// input is not a square matrix";
542    return(X);
543  }
544
545  //test for correct variable
546  if( size(#)==0 ){
547    #[1] = varstr(nvars(BR));
548  }
549  if( typeof(#[1]) == "string"){
550     v = #[1];
551  }
552  else{
553    "// 2nd argument must be a name of a variable not contained in the matrix";
554    return(X);
555  }
556
557  j=-1;
558  for(i=1; i<=nvars(BR); i++){
559    if(varstr(i)==v){j=i;}
560  }
561  if(j==-1){
562    "// "+v+" is not a variable in the basering";
563    return(X);
564  }
565
566  //var can not be in A
567  s="Wp(";
568  for( i=1; i<=nvars(BR); i++ ){
569    if(i!=j){ s=s+"0";}
570    else{ s=s+"1";}
571    if( i<nvars(BR)) {s=s+",";}
572  }
573  s=s+")";
574
575  changeord("@R",s);
576  execute("minpoly="+mp+";");
577  matrix A = imap(BR,A);
578  for(i=1; i<=n; i++){
579    if(deg(lead(A)[i])>=1){
580      "// matrix must not contain the variable "+v;
581      kill @R;
582      setring BR;
583      return(X);
584    }
585  }
586
587  //get coefficients and build the char. poly
588  kill @R;
589  setring BR;
590  L = busadj(A);
591  for(i=1; i<=n+1; i++){
592    execute("X=X+L[1][i]*"+v+"^"+string(i-1)+";");
593  }
594
595  return(X);
596}
597example
598{ "EXAMPLE"; echo=2;
599  ring r=0,(x,t),dp;
600  matrix A[3][3]=1,x2,x,x2,6,4,x,4,1;
601  print(A);
602  charpoly_B(A,"t");
603}
604
605//////////////////////////////////////////////////////////////////////////////
606proc adjoint(matrix A)
607"USAGE:    adjoint(A);  A = square matrix
608RETURN:   adjoint matrix of A, i.e. Adj*A=det(A)*E
609NOTE:     computation uses busadj(A)
610EXAMPLE:  example adjoint; shows an example"
611{
612  int n=nrows(A);
613  matrix Adj[n][n];
614  list L;
615
616  if(ncols(A) != n) {
617    "// input is not a square matrix";
618    return(Adj);
619  }
620
621  L  = busadj(A);
622  Adj= (-1)^(n-1)*L[2][1];
623  return(Adj);
624
625}
626example
627{ "EXAMPLE"; echo=2;
628  ring r=0,(t,x),lp;
629  matrix A[2][2]=1,x2,x,x2+3x;
630  print(A);
631  matrix Adj[2][2]=adjoint(A);
632  print(Adj);                    //Adj*A=det(A)*E
633  print(Adj*A);
634}
635
636//////////////////////////////////////////////////////////////////////////////
637proc inverse_B(matrix A)
638"USAGE:    inverse_B(A);  A = square matrix
639RETURN:   list Inv with
640          - Inv[1] = matrix I and
641          - Inv[2] = poly p
642          such that I*A = unitmat(n)*p;
643NOTE:     p=1 if 1/det(A) is computable and  p=det(A) if not;
644          the computation uses busadj.
645SEE ALSO: inverse, inverse_L
646EXAMPLE:  example inverse_B; shows an example"
647{
648  int i;
649  int n=nrows(A);
650  matrix I[n][n];
651  poly factor;
652  list L;
653  list Inv;
654
655  if(ncols(A) != n) {
656    "input is not a square matrix";
657    return(I);
658  }
659
660  L=busadj(A);
661  I=module(-L[2][1]);        //+-Adj(A)
662
663  if(reduce(1,std(L[1][1]))==0){
664    I=I*lift(L[1][1],1)[1][1];
665    factor=1;
666  }
667  else{ factor=L[1][1];}     //=+-det(A) or 1
668  Inv=insert(Inv,factor);
669  Inv=insert(Inv,matrix(I));
670
671  return(Inv);
672}
673example
674{ "EXAMPLE"; echo=2;
675  ring r=0,(x,y),lp;
676  matrix A[3][3]=x,y,1,1,x2,y,x,6,0;
677  print(A);
678  list Inv=inverse_B(A);
679  print(Inv[1]);
680  print(Inv[2]);
681  print(Inv[1]*A);
682}
683
684//////////////////////////////////////////////////////////////////////////////
685proc det_B(matrix A)
686"USAGE:     det_B(A);  A any matrix
687RETURN:    returns the determinant of A
688NOTE:      the computation uses the busadj algorithm
689EXAMPLE:   example det_B; shows an example"
690{
691  int n=nrows(A);
692  list L;
693
694  if(ncols(A) != n){ return(0);}
695
696  L=busadj(A);
697  return((-1)^n*L[1][1]);
698}
699example
700{ "EXAMPLE"; echo=2;
701  ring r=0,(x),dp;
702  matrix A[10][10]=random(2,10,10)+unitmat(10)*x;
703  print(A);
704  det_B(A);
705}
706
707//////////////////////////////////////////////////////////////////////////////
708proc inverse_L(matrix A)
709"USAGE:     inverse_L(A);  A = square matrix
710RETURN:    list Inv representing a left inverse of A, i.e
711           - Inv[1] = matrix I and
712           - Inv[2] = poly p
713           such that I*A = unitmat(n)*p;
714NOTE:      p=1 if 1/det(A) is computable and p=det(A) if not;
715           the computation computes first det(A) and then uses lift
716SEE ALSO:  inverse, inverse_B
717EXAMPLE:   example inverse_L; shows an example"
718{
719  int n=nrows(A);
720  matrix I;
721  matrix E[n][n]=unitmat(n);
722  poly factor;
723  poly d=1;
724  list Inv;
725
726  if (ncols(A)!=n){
727    "// input is not a square matrix";
728    return(I);
729  }
730
731  d=det(A);
732  if(d==0){
733    "// matrix is not invertible";
734    return(Inv);
735  }
736
737  // test if 1/det(A) exists
738  if(reduce(1,std(d))!=0){ E=E*d;}
739
740  I=lift(A,E);
741  if(I==unitmat(n)-unitmat(n)){ //catch error in lift
742    "// matrix is not invertible";
743    return(Inv);
744  }
745
746  factor=d;      //=det(A) or 1
747  Inv=insert(Inv,factor);
748  Inv=insert(Inv,I);
749
750  return(Inv);
751}
752example
753{ "EXAMPLE"; echo=2;
754  ring r=0,(x,y),lp;
755  matrix A[3][3]=x,y,1,1,x2,y,x,6,0;
756  print(A);
757  list Inv=inverse_L(A);
758  print(Inv[1]);
759  print(Inv[2]);
760  print(Inv[1]*A);
761}
762
763//////////////////////////////////////////////////////////////////////////////
764proc gaussred(matrix A)
765"USAGE:   gaussred(A);   A any constant matrix
766RETURN:  list Z:  Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A)
767         gives a row reduced matrix S, a permutation matrix P and a
768         normalized lower triangular matrix U, with P*A=U*S
769NOTE:    This procedure is designed for teaching purposes mainly.
770         The straight forward implementation in the interpreted library
771         is not very efficient (no standard basis computation).
772EXAMPLE: example gaussred; shows an example"
773{
774  int i,j,l,k,jp,rang;
775  poly c,pivo;
776  list Z;
777  int n = nrows(A);
778  int m = ncols(A);
779  int mr= n; //max. rang
780  matrix P[n][n] = unitmat(n);
781  matrix U[n][n] = P;
782
783  if(!const_mat(A)){
784    "// input is not a constant matrix";
785    return(Z);
786  }
787
788  if(n>m){mr=m;} //max. rang
789
790  for(i=1;i<=mr;i=i+1){
791    if((i+k)>m){break};
792
793    //Test: Diagonalelement=0
794    if(A[i,i+k]==0){
795      jp=i;pivo=0;
796      for(j=i+1;j<=n;j=j+1){
797        c=abs(A[j,i+k]);
798        if(pivo<c){ pivo=c;jp=j;}
799      }
800      if(jp != i){       //Zeilentausch
801        for(j=1;j<=m;j=j+1){ //Zeilentausch in A (und U) (i-te mit jp-ter)
802          c=A[i,j];
803          A[i,j]=A[jp,j];
804          A[jp,j]=c;
805        }
806        for(j=1;j<=n;j=j+1){ //Zeilentausch in P
807          c=P[i,j];
808          P[i,j]=P[jp,j];
809          P[jp,j]=c;
810        }
811      }
812      if(pivo==0){k++;continue;} //eine von selbst auftauchende Stufe !
813    }                          //i sollte im naechsten Lauf nicht erhoeht sein
814
815    //Eliminationsschritt
816    for(j=i+1;j<=n;j=j+1){
817      c=A[j,i+k]/A[i,i+k];
818      for(l=i+k+1;l<=m;l=l+1){
819        A[j,l]=A[j,l]-A[i,l]*c;
820      }
821      A[j,i+k]=0;  // nur wichtig falls k>0 ist
822      A[j,i]=c;    // bildet U
823    }
824  rang=i;
825  }
826
827  for(i=1;i<=mr;i=i+1){
828    for(j=i+1;j<=n;j=j+1){
829      U[j,i]=A[j,i];
830      A[j,i]=0;
831    }
832  }
833
834  Z=insert(Z,rang);
835  Z=insert(Z,A);
836  Z=insert(Z,U);
837  Z=insert(Z,P);
838
839  return(Z);
840}
841example
842{ "EXAMPLE";echo=2;
843  ring r=0,(x),dp;
844  matrix A[5][4]=1,3,-1,4,2,5,-1,3,1,3,-1,4,0,4,-3,1,-3,1,-5,-2;
845  print(A);
846  list Z=gaussred(A);   //construct P,U,S s.t. P*A=U*S
847  print(Z[1]);          //P
848  print(Z[2]);          //U
849  print(Z[3]);          //S
850  print(Z[4]);          //rank
851  print(Z[1]*A);        //P*A
852  print(Z[2]*Z[3]);     //U*S
853}
854
855//////////////////////////////////////////////////////////////////////////////
856proc gaussred_pivot(matrix A)
857"USAGE:     gaussred_pivot(A);   A any constant matrix
858RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A)
859           gives n row reduced matrix S, a permutation matrix P and a
860           normalized lower triangular matrix U, with P*A=U*S
861NOTE:      with row pivoting
862EXAMPLE:   example gaussred_pivot; shows an example"
863{
864  int i,j,l,k,jp,rang;
865  poly c,pivo;
866  list Z;
867  int n=nrows(A);
868  int m=ncols(A);
869  int mr=n; //max. rang
870  matrix P[n][n]=unitmat(n);
871  matrix U[n][n]=P;
872
873  if(!const_mat(A)){
874    "// input is not a constant matrix";
875    return(Z);
876  }
877
878  if(n>m){mr=m;} //max. rang
879
880  for(i=1;i<=mr;i=i+1){
881    if((i+k)>m){break};
882
883    //Pivotisierung
884    pivo=abs(A[i,i+k]);jp=i;
885    for(j=i+1;j<=n;j=j+1){
886      c=abs(A[j,i+k]);
887      if(pivo<c){ pivo=c;jp=j;}
888    }
889    if(jp != i){ //Zeilentausch
890      for(j=1;j<=m;j=j+1){ //Zeilentausch in A (und U) (i-te mit jp-ter)
891        c=A[i,j];
892        A[i,j]=A[jp,j];
893        A[jp,j]=c;
894      }
895      for(j=1;j<=n;j=j+1){ //Zeilentausch in P
896        c=P[i,j];
897        P[i,j]=P[jp,j];
898        P[jp,j]=c;
899      }
900    }
901    if(pivo==0){k++;continue;} //eine von selbst auftauchende Stufe !
902                               //i sollte im naechsten Lauf nicht erhoeht sein
903    //Eliminationsschritt
904    for(j=i+1;j<=n;j=j+1){
905      c=A[j,i+k]/A[i,i+k];
906      for(l=i+k+1;l<=m;l=l+1){
907        A[j,l]=A[j,l]-A[i,l]*c;
908      }
909      A[j,i+k]=0;  // nur wichtig falls k>0 ist
910      A[j,i]=c;    // bildet U
911    }
912  rang=i;
913  }
914
915  for(i=1;i<=mr;i=i+1){
916    for(j=i+1;j<=n;j=j+1){
917      U[j,i]=A[j,i];
918      A[j,i]=0;
919    }
920  }
921
922  Z=insert(Z,rang);
923  Z=insert(Z,A);
924  Z=insert(Z,U);
925  Z=insert(Z,P);
926
927  return(Z);
928}
929example
930{ "EXAMPLE";echo=2;
931  ring r=0,(x),dp;
932  matrix A[5][4] = 1, 3,-1,4,
933                   2, 5,-1,3,
934                   1, 3,-1,4,
935                   0, 4,-3,1,
936                  -3,1,-5,-2;
937  list Z=gaussred_pivot(A);  //construct P,U,S s.t. P*A=U*S
938  print(Z[1]);               //P
939  print(Z[2]);               //U
940  print(Z[3]);               //S
941  print(Z[4]);               //rank
942  print(Z[1]*A);             //P*A
943  print(Z[2]*Z[3]);          //U*S
944}
945
946//////////////////////////////////////////////////////////////////////////////
947proc gauss_nf(matrix A)
948"USAGE:    gauss_nf(A); A any constant matrix
949RETURN:   matrix; gauss normal form of A (uses gaussred)
950EXAMPLE:  example gauss_nf; shows an example"
951{
952  list Z;
953  if(!const_mat(A)){
954    "// input is not a constant matrix";
955    return(A);
956  }
957  Z = gaussred(A);
958  return(Z[3]);
959}
960example
961{ "EXAMPLE";echo=2;
962  ring r = 0,(x),dp;
963  matrix A[4][4] = 1,4,4,7,2,5,5,4,4,1,1,3,0,2,2,7;
964  print(gauss_nf(A));
965}
966
967//////////////////////////////////////////////////////////////////////////////
968proc mat_rk(matrix A)
969"USAGE:    mat_rk(A); A any constant matrix
970RETURN:   int, rank of A
971EXAMPLE:  example mat_rk; shows an example"
972{
973  list Z;
974  if(!const_mat(A)){
975    "// input is not a constant matrix";
976    return(-1);
977  }
978  Z = gaussred(A);
979  return(Z[4]);
980}
981example
982{ "EXAMPLE";echo=2;
983  ring r = 0,(x),dp;
984  matrix A[4][4] = 1,4,4,7,2,5,5,4,4,1,1,3,0,2,2,7;
985  mat_rk(A);
986}
987
988//////////////////////////////////////////////////////////////////////////////
989proc U_D_O(matrix A)
990"USAGE:     U_D_O(A);   constant invertible matrix A
991RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=D , Z[4]=O
992           gives a permutation matrix P,
993           a normalized lower triangular matrix U ,
994           a diagonal matrix D, and
995           a normalized upper triangular matrix O
996           with P*A=U*D*O
997NOTE:      Z[1]=-1 means that A is not regular (proc uses gaussred)
998EXAMPLE:   example U_D_O; shows an example"
999{
1000  int i,j;
1001  list Z,L;
1002  int n=nrows(A);
1003  matrix O[n][n]=unitmat(n);
1004  matrix D[n][n];
1005
1006  if (ncols(A)!=n){
1007    "// input is not a square matrix";
1008    return(Z);
1009  }
1010  if(!const_mat(A)){
1011    "// input is not a constant matrix";
1012    return(Z);
1013  }
1014
1015  L=gaussred(A);
1016
1017  if(L[4]!=n){
1018    "// input is not an invertible matrix";
1019    Z=insert(Z,-1);  //hint for calling procedures
1020    return(Z);
1021  }
1022
1023  D=L[3];
1024
1025  for(i=1; i<=n; i++){
1026    for(j=i+1; j<=n; j++){
1027      O[i,j] = D[i,j]/D[i,i];
1028      D[i,j] = 0;
1029    }
1030  }
1031
1032  Z=insert(Z,O);
1033  Z=insert(Z,D);
1034  Z=insert(Z,L[2]);
1035  Z=insert(Z,L[1]);
1036  return(Z);
1037}
1038example
1039{ "EXAMPLE";echo=2;
1040  ring r = 0,(x),dp;
1041  matrix A[5][5] = 10, 4,  0, -9,  8,
1042                   -3, 6, -6, -4,  9,
1043                    0, 3, -1, -9, -8,
1044                   -4,-2, -6, -10,10,
1045                   -9, 5, -1, -6,  5;
1046  list Z = U_D_O(A);              //construct P,U,D,O s.t. P*A=U*D*O
1047  print(Z[1]);                    //P
1048  print(Z[2]);                    //U
1049  print(Z[3]);                    //D
1050  print(Z[4]);                    //O
1051  print(Z[1]*A);                  //P*A
1052  print(Z[2]*Z[3]*Z[4]);          //U*D*O
1053}
1054
1055//////////////////////////////////////////////////////////////////////////////
1056proc pos_def(matrix A)
1057"USAGE:     pos_def(A); A = constant, symmetric square matrix
1058RETURN:    int:
1059           1  if A is positive definit ,
1060           0  if not,
1061           -1 if unknown
1062EXAMPLE:   example pos_def; shows an example"
1063{
1064  int j;
1065  list Z;
1066  int n = nrows(A);
1067  matrix H[n][n];
1068
1069  if (ncols(A)!=n){
1070    "// input is not a square matrix";
1071    return(0);
1072  }
1073  if(!const_mat(A)){
1074    "// input is not a constant matrix";
1075    return(-1);
1076  }
1077  if(deg(std(A-transpose(A))[1])!=-1){
1078    "// input is not a hermitian (symmetric) matrix";
1079    return(-1);
1080  }
1081
1082  Z=U_D_O(A);
1083
1084  if(Z[1]==-1){
1085    return(0);
1086  }  //A not regular, therefore not pos. definit
1087
1088  H=Z[1];
1089  //es fand Zeilentausch statt: also nicht positiv definit
1090  if(deg(std(H-unitmat(n))[1])!=-1){
1091    return(0);
1092  }
1093
1094  H=Z[3];
1095
1096  for(j=1;j<=n;j=j+1){
1097    if(H[j,j]<=0){
1098      return(0);
1099    } //eigenvalue<=0, not pos.definit
1100  }
1101
1102  return(1); //positiv definit;
1103}
1104example
1105{ "EXAMPLE"; echo=2;
1106  ring r = 0,(x),dp;
1107  matrix A[5][5] = 20,  4,  0, -9,   8,
1108                    4, 12, -6, -4,   9,
1109                    0, -6, -2, -9,  -8,
1110                   -9, -4, -9, -20, 10,
1111                    8,  9, -8,  10, 10;
1112  pos_def(A);
1113  matrix B[3][3] =  3,  2,  0,
1114                    2, 12,  4,
1115                    0,  4,  2;
1116  pos_def(B);
1117}
1118
1119//////////////////////////////////////////////////////////////////////////////
1120proc linsolve(matrix A, matrix b)
1121"USAGE:     linsolve(A,b); A a constant nxm-matrix, b a constant nx1-matrix
1122RETURN:    a 1xm matrix X, solution of inhomogeneous linear system A*X = b
1123           return the 0-matrix if system is not solvable
1124NOTE:      uses gaussred
1125EXAMPLE:   example linsolve; shows an example"
1126{
1127  int i,j,k,rc,r;
1128  poly c;
1129  list Z;
1130  int n  = nrows(A);
1131  int m  = ncols(A);
1132  int n_b= nrows(b);
1133  matrix Ab[n][m+1];
1134  matrix X[m][1];
1135
1136  if(ncols(b)!=1){
1137    "// right hand side b is not a nx1 matrix";
1138    return(X);
1139  }
1140
1141  if(!const_mat(A)){
1142    "// input hand is not a constant matrix";
1143    return(X);
1144  }
1145
1146  if(n_b>n){
1147    for(i=n; i<=n_b; i++){
1148      if(b[i,1]!=0){
1149        "// right hand side b not in Image(A)";
1150        return X;
1151      }
1152    }
1153  }
1154
1155  if(n_b<n){
1156    matrix copy[n_b][1]=b;
1157    matrix b[n][1]=0;
1158    for(i=1;i<=n_b;i=i+1){
1159      b[i,1]=copy[i,1];
1160    }
1161  }
1162
1163  r=mat_rk(A);
1164
1165  //1. b constant vector
1166  if(const_mat(b)){
1167    //extend A with b
1168    for(i=1; i<=n; i++){
1169      for(j=1; j<=m; j++){
1170        Ab[i,j]=A[i,j];
1171      }
1172      Ab[i,m+1]=b[i,1];
1173    }
1174
1175    //Gauss reduction
1176    Z  = gaussred(Ab);
1177    Ab = Z[3];  //normal form
1178    rc = Z[4];  //rank(Ab)
1179    //print(Ab);
1180
1181    if(r<rc){
1182        "// no solution";
1183        return(X);
1184    }
1185    k=m;
1186    for(i=r;i>=1;i=i-1){
1187
1188      j=1;
1189      while(Ab[i,j]==0){j=j+1;}// suche Ecke
1190
1191      for(;k>j;k=k-1){ X[k]=0;}//springe zur Ecke
1192
1193
1194      c=Ab[i,m+1]; //i-te Komponene von b
1195      for(j=m;j>k;j=j-1){
1196        c=c-X[j,1]*Ab[i,j];
1197      }
1198      if(Ab[i,k]==0){
1199        X[k,1]=1; //willkuerlich
1200      }
1201      else{
1202          X[k,1]=c/Ab[i,k];
1203      }
1204      k=k-1;
1205      if(k==0){break;}
1206    }
1207
1208
1209  }//endif (const b)
1210  else{  //b not constant
1211    "// !not implemented!";
1212
1213  }
1214
1215  return(X);
1216}
1217example
1218{ "EXAMPLE";echo=2;
1219  ring r=0,(x),dp;
1220  matrix A[3][2] = -4,-6,
1221                    2, 3,
1222                   -5, 7;
1223  matrix b[3][1] = 10,
1224                   -5,
1225                    2;
1226  matrix X = linsolve(A,b);
1227  print(X);
1228  print(A*X);
1229}
1230//////////////////////////////////////////////////////////////////////////////
1231
1232///////////////////////////////////////////////////////////////////////////////
1233//    PROCEDURES for Jordan normal form
1234//    AUTHOR:  Mathias Schulze, email: mschulze@mathematik.uni-kl.de
1235///////////////////////////////////////////////////////////////////////////////
1236
1237proc eigenvals(matrix M)
1238"USAGE:   eigenvals(M); matrix M
1239ASSUME:  eigenvalues of M in basefield
1240RETURN:
1241@format
1242list l;
1243  ideal l[1];
1244    number l[1][i];  i-th eigenvalue of M
1245  intvec l[2];
1246    int l[2][i];  multiplicity of i-th eigenvalue of M
1247@end format
1248EXAMPLE: example eigenvals; shows examples
1249"
1250{
1251  return(system("eigenvalues",M));
1252}
1253example
1254{ "EXAMPLE:"; echo=2;
1255  ring R=0,x,dp;
1256  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
1257  print(M);
1258  eigenvals(M);
1259}
1260///////////////////////////////////////////////////////////////////////////////
1261
1262proc jordan(matrix M,list #)
1263"USAGE:   jordan(M); matrix M
1264ASSUME:  eigenvalues of M in basefield
1265RETURN:
1266@format
1267list l;  Jordan data of M
1268  ideal l[1];
1269    number l[1][i];  eigenvalue of i-th Jordan block of M
1270  intvec l[2];
1271    int l[2][i];  size of i-th Jordan block of M
1272  intvec l[3];
1273    int l[3][i];  multiplicity of i-th Jordan block of M
1274@end format
1275EXAMPLE: example jordan; shows examples
1276"
1277{
1278  if(nrows(M)==0)
1279  {
1280    ERROR("non empty expected");
1281  }
1282  if(ncols(M)!=nrows(M))
1283  {
1284    ERROR("square matrix expected");
1285  }
1286
1287  M=jet(M,0);
1288
1289  if(size(#)==0)
1290  {
1291    #=eigenvals(M);
1292  }
1293  def e0,m0=#[1..2];
1294
1295  int i;
1296  for(i=1;i<=ncols(e0);i++)
1297  {
1298    if(deg(e0[i])>0)
1299    {
1300
1301      ERROR("eigenvalues in coefficient field expected");
1302      return(list());
1303    }
1304  }
1305
1306  int j,k;
1307  matrix N0,N1;
1308  module K0;
1309  list K;
1310  ideal e;
1311  intvec s,m;
1312
1313  for(i=1;i<=ncols(e0);i++)
1314  {
1315    N0=M-e0[i]*freemodule(ncols(M));
1316
1317    N1=N0;
1318    K0=0;
1319    K=module();
1320    while(size(K0)<m0[i])
1321    {
1322      K0=syz(N1);
1323      K=K+list(K0);
1324      N1=N1*N0;
1325    }
1326
1327    for(j=2;j<size(K);j++)
1328    {
1329      if(2*size(K[j])-size(K[j-1])-size(K[j+1])>0)
1330      {
1331        k++;
1332        e[k]=e0[i];
1333        s[k]=j-1;
1334        m[k]=2*size(K[j])-size(K[j-1])-size(K[j+1]);
1335      }
1336    }
1337    if(size(K[j])-size(K[j-1])>0)
1338    {
1339      k++;
1340      e[k]=e0[i];
1341      s[k]=j-1;
1342      m[k]=size(K[j])-size(K[j-1]);
1343    }
1344  }
1345
1346  return(list(e,s,m));
1347}
1348example
1349{ "EXAMPLE:"; echo=2;
1350  ring R=0,x,dp;
1351  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
1352  print(M);
1353  jordan(M);
1354}
1355///////////////////////////////////////////////////////////////////////////////
1356
1357proc jordanbasis(matrix M,list #)
1358"USAGE:   jordanbasis(M); matrix M
1359ASSUME:  eigenvalues of M in basefield
1360RETURN:
1361@format
1362list l:
1363  module l[1];  inverse(l[1])*M*l[1] Jordan normal form
1364  intvec l[2];
1365    int l[2][i];  weight filtration index of l[1][i]
1366@end format
1367EXAMPLE: example jordanbasis; shows examples
1368"
1369{
1370  if(nrows(M)==0)
1371  {
1372    ERROR("non empty matrix expected");
1373  }
1374  if(ncols(M)!=nrows(M))
1375  {
1376    ERROR("square matrix expected");
1377  }
1378
1379  M=jet(M,0);
1380
1381  if(size(#)==0)
1382  {
1383    #=eigenvals(M);
1384  }
1385  def e,m=#[1..2];
1386
1387  for(int i=1;i<=ncols(e);i++)
1388  {
1389    if(deg(e[i]>0))
1390    {
1391      ERROR("eigenvalues in coefficient field expected");
1392      return(freemodule(ncols(M)));
1393    }
1394  }
1395
1396  int j,k,l,n;
1397  matrix N0,N1;
1398  module K0,K1;
1399  list K;
1400  matrix u[ncols(M)][1];
1401  module U;
1402  intvec w;
1403
1404  for(i=1;i<=ncols(e);i++)
1405  {
1406    N0=M-e[i]*freemodule(ncols(M));
1407
1408    N1=N0;
1409    K0=0;
1410    K=list();
1411    while(size(K0)<m[i])
1412    {
1413      K0=syz(N1);
1414      K=K+list(K0);
1415      N1=N1*N0;
1416    }
1417
1418    K1=0;
1419    for(j=1;j<size(K);j++)
1420    {
1421      K0=K[j];
1422      K[j]=interred(reduce(K[j],std(K1+module(N0*K[j+1]))));
1423      K1=K0;
1424    }
1425    K[j]=interred(reduce(K[j],std(K1)));
1426
1427    for(l=size(K);l>=1;l--)
1428    {
1429      for(k=size(K[l]);k>0;k--)
1430      {
1431        u=K[l][k];
1432        for(j=l;j>=1;j--)
1433        {
1434          U=U+module(u);
1435          n++;
1436          w[n]=2*j-l-1;
1437          u=N0*u;
1438        }
1439      }
1440    }
1441  }
1442
1443  return(list(U,w));
1444}
1445example
1446{ "EXAMPLE:"; echo=2;
1447  ring R=0,x,dp;
1448  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
1449  print(M);
1450  list l=jordanbasis(M);
1451  print(l[1]);
1452  print(l[2]);
1453  print(inverse(l[1])*M*l[1]);
1454}
1455///////////////////////////////////////////////////////////////////////////////
1456
1457proc jordanmatrix(ideal e,intvec s,intvec m)
1458"USAGE:   jordanmatrix(e,s,m); ideal e, intvec s, intvec m
1459ASSUME:  ncols(e)==size(s)==size(m)
1460RETURN:
1461@format
1462matrix J;  list(e,s,m)==jordan(J)
1463@end format
1464EXAMPLE: example jordanmatrix; shows examples
1465"
1466{
1467  if(ncols(e)!=size(s)||size(e)!=size(m))
1468  {
1469    ERROR("arguments of equal size expected");
1470  }
1471
1472  int i,j,k,l;
1473  int n=int((transpose(matrix(s))*matrix(m))[1,1]);
1474  matrix J[n][n];
1475  for(k=1;k<=ncols(e);k++)
1476  {
1477    for(l=1;l<=m[k];l++)
1478    {
1479      j++;
1480      J[j,j]=e[k];
1481      for(i=s[k];i>=2;i--)
1482      {
1483        J[j+1,j]=1;
1484        j++;
1485        J[j,j]=e[k];
1486      }
1487    }
1488  }
1489
1490  return(J);
1491}
1492example
1493{ "EXAMPLE:"; echo=2;
1494  ring R=0,x,dp;
1495  ideal e=ideal(2,3);
1496  intvec s=1,2;
1497  intvec m=1,1;
1498  print(jordanmatrix(e,s,m));
1499}
1500///////////////////////////////////////////////////////////////////////////////
1501
1502proc jordannf(matrix M,list #)
1503"USAGE:   jordannf(M); matrix M
1504ASSUME:  eigenvalues of M in basefield
1505RETURN:  matrix J; Jordan normal form of M
1506EXAMPLE: example jordannf; shows examples
1507"
1508{
1509  list l=jordan(M,#);
1510  return(jordanmatrix(l[1],l[2]));
1511}
1512example
1513{ "EXAMPLE:"; echo=2;
1514  ring R=0,x,dp;
1515  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
1516  print(M);
1517  print(jordannf(M));
1518}
1519///////////////////////////////////////////////////////////////////////////////
1520
1521/*
1522///////////////////////////////////////////////////////////////////////////////
1523//          Auskommentierte zusaetzliche Beispiele
1524//
1525///////////////////////////////////////////////////////////////////////////////
1526// Singular for ix86-Linux version 1-3-10  (2000121517)  Dec 15 2000 17:55:12
1527// Rechnungen auf AMD700 mit 632 MB
1528
1529 LIB "linalg.lib";
1530
15311. Sparse integer Matrizen
1532--------------------------
1533ring r1=0,(x),dp;
1534system("--random", 12345678);
1535int n = 70;
1536matrix m = sparsemat(n,n,50,100);
1537option(prot,mem);
1538
1539int t=timer;
1540matrix im = inverse(m,1)[1];
1541timer-t;
1542print(im*m);
1543//list l0 = watchdog(100,"inverse("+"m"+",3)");
1544//bricht bei 100 sec ab und gibt l0[1]: string Killed zurueck
1545
1546//inverse(m,1): std       5sec 5,5 MB
1547//inverse(m,2): interred  12sec
1548//inverse(m,2): lift      nach 180 sec 13MB abgebrochen
1549//n=60: linalgorig: 3  linalg: 5
1550//n=70: linalgorig: 6,7 linalg: 11,12
1551// aber linalgorig rechnet falsch!
1552
15532. Sparse poly Matrizen
1554-----------------------
1555ring r=(0),(a,b,c),dp;
1556system("--random", 12345678);
1557int n=6;
1558matrix m = sparsematrix(n,n,2,0,50,50,9); //matrix of polys of deg <=2
1559option(prot,mem);
1560
1561int t=timer;
1562matrix im = inverse(m);
1563timer-t;
1564print(im*m);
1565//inverse(m,1): std       0sec 1MB
1566//inverse(m,2): interred  0sec 1MB
1567//inverse(m,2): lift      nach 2000 sec 33MB abgebrochen
1568
15693. Sparse Matrizen mit Parametern
1570---------------------------------
1571//liborig rechnet hier falsch!
1572ring r=(0),(a,b),dp;
1573system("--random", 12345678);
1574int n=7;
1575matrix m = sparsematrix(n,n,1,0,40,50,9);
1576ring r1 = (0,a,b),(x),dp;
1577matrix m = imap(r,m);
1578option(prot,mem);
1579
1580int t=timer;
1581matrix im = inverse(m);
1582timer-t;
1583print(im*m);
1584//inverse(m)=inverse(m,3):15 sec inverse(m,1)=1sec inverse(m,2):>120sec
1585//Bei Parametern vergeht die Zeit beim Normieren!
1586
15873. Sparse Matrizen mit Variablen und Parametern
1588-----------------------------------------------
1589ring r=(0),(a,b),dp;
1590system("--random", 12345678);
1591int n=6;
1592matrix m = sparsematrix(n,n,1,0,35,50,9);
1593ring r1 = (0,a),(b),dp;
1594matrix m = imap(r,m);
1595option(prot,mem);
1596
1597int t=timer;
1598matrix im = inverse(m,3);
1599timer-t;
1600print(im*m);
1601//n=7: inverse(m,3):lange sec inverse(m,1)=1sec inverse(m,2):1sec
1602
16034. Ueber Polynomring invertierbare Matrizen
1604-------------------------------------------
1605LIB"random.lib"; LIB"linalg.lib";
1606system("--random", 12345678);
1607int n =3;
1608ring r= 0,(x,y,z),(C,dp);
1609matrix A=triagmatrix(n,n,1,0,0,50,2);
1610intmat B=sparsetriag(n,n,20,1);
1611matrix M = A*transpose(B);
1612M=M*transpose(M);
1613M[1,1..ncols(M)]=M[1,1..n]+xyz*M[n,1..ncols(M)];
1614print(M);
1615//M hat det=1 nach Konstruktion
1616
1617int t=timer;
1618matrix iM=inverse(M);
1619timer-t;
1620print(iM*M);                      //test
1621
1622//ACHTUNG: Interred liefert i.A. keine Inverse, Gegenbeispiel z.B.
1623//mit n=3
1624//eifacheres Gegenbeispiel:
1625matrix M =
16269yz+3y+3z+2,             9y2+6y+1,
16279xyz+3xy+3xz-9z2+2x-6z-1,9xy2+6xy-9yz+x-3y-3z
1628//det M=1, inverse(M,2); ->// ** matrix is not invertible
1629//lead(M); 9xyz*gen(2) 9xy2*gen(2) nicht teilbar!
1630
16315. charpoly:
1632-----------
1633//ring rp=(0,A,B,C),(x),dp;
1634ring r=0,(A,B,C,x),dp;
1635matrix m[12][12]=
1636AC,BC,-3BC,0,-A2+B2,-3AC+1,B2,   B2,  1,   0, -C2+1,0,
16371, 1, 2C,  0,0,     B,     -A,   -4C, 2A+1,0, 0,    0,
16380, 0, 0,   1,0,     2C+1,  -4C+1,-A,  B+1, 0, B+1,  3B,
1639AB,B2,0,   1,0,     1,     0,    1,   A,   0, 1,    B+1,
16401, 0, 1,   0,0,     1,     0,    -C2, 0,   1, 0,    1,
16410, 0, 2,   1,2A,    1,     0,    0,   0,   0, 1,    1,
16420, 1, 0,   1,1,     2,     A,    3B+1,1,   B2,1,    1,
16430, 1, 0,   1,1,     1,     1,    1,   2,   0, 0,    0,
16441, 0, 1,   0,0,     0,     1,    0,   1,   1, 0,    3,
16451, 3B,B2+1,0,0,     1,     0,    1,   0,   0, 1,    0,
16460, 0, 1,   0,0,     0,     0,    1,   0,   0, 0,    0,
16470, 1, 0,   1,1,     3,     3B+1, 0,   1,   1, 1,    0;
1648option(prot,mem);
1649
1650int t=timer;
1651poly q=charpoly(m,"x");         //1sec, charpoly_B 1sec, 16MB
1652timer-t;
1653//1sec, charpoly_B 1sec, 16MB  (gleich in r und rp)
1654
1655*/
Note: See TracBrowser for help on using the repository browser.