source: git/Singular/LIB/linalg.lib @ 70597c

fieker-DuValspielwiese
Last change on this file since 70597c was daa83b, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes/lossen: new ring.lib git-svn-id: file:///usr/local/Singular/svn/trunk@7920 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 44.8 KB
RevLine 
[337919]1//GMG last modified: 04/25/2000
[30f140]2//////////////////////////////////////////////////////////////////////////////
[daa83b]3version="$Id: linalg.lib,v 1.37 2005-04-28 09:22:15 Singular Exp $";
[49998f]4category="Linear Algebra";
[30f140]5info="
[275721f]6LIBRARY:  linalg.lib  Algorithmic Linear Algebra
[6188357]7AUTHORS:  Ivor Saynisch (ivs@math.tu-cottbus.de)
8@*        Mathias Schulze (mschulze@mathematik.uni-kl.de)
[30f140]9
10PROCEDURES:
[a0c62d]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
[e9124e]27 hessenberg(M);       Hessenberg form of M
[7248478]28 eigenvals(M);        eigenvalues with multiplicities of M
[2699a6]29 minipoly(M);         minimal polynomial of M
[cb40b5]30 spnf(sp);            normal form of spectrum sp
31 spprint(sp);         print spectrum sp
[275721f]32 jordan(M);           Jordan data of M
33 jordanbasis(M);      Jordan basis and weight filtration of M
[cb40b5]34 jordanmatrix(jd);    Jordan matrix with Jordan data jd
[275721f]35 jordannf(M);         Jordan normal form of M
[30f140]36";
37
38LIB "matrix.lib";
39LIB "ring.lib";
[6188357]40LIB "elim.lib";
[6d37e8]41LIB "general.lib";
[30f140]42//////////////////////////////////////////////////////////////////////////////
43// help functions
44//////////////////////////////////////////////////////////////////////////////
[6188357]45static proc const_mat(matrix A)
[30f140]46"RETURN:   1 (0) if A is (is not) a constant matrix"
47{
48  int i;
49  int n=ncols(A);
50  def BR=basering;
[daa83b]51  def @R=changeord("dp,c",BR);
52  setring @R;
[30f140]53  matrix A=fetch(BR,A);
54  for(i=1;i<=n;i=i+1){
55    if(deg(lead(A)[i])>=1){
56      //"input is not a constant matrix";
57      kill @R;
58      setring BR;
59      return(0);
60    }
61  }
62  kill @R;
63  setring BR;
[337919]64  return(1);
[30f140]65}
[8942a5]66//////////////////////////////////////////////////////////////////////////////
67static proc red(matrix A,int i,int j)
[30f140]68"USAGE:    red(A,i,j);  A = constant matrix
69          reduces column j with respect to A[i,i] and column i
70          reduces row j with respect to A[i,i] and row i
71RETURN:   matrix
72"
[337919]73{
[30f140]74  module m=module(A);
[337919]75
[30f140]76  if(A[i,i]==0){
77    m[i]=m[i]+m[j];
78    m=module(transpose(matrix(m)));
79    m[i]=m[i]+m[j];
80    m=module(transpose(matrix(m)));
81  }
[337919]82
[30f140]83  A=matrix(m);
84  m[j]=m[j]-(A[i,j]/A[i,i])*m[i];
85  m=module(transpose(matrix(m)));
86  m[j]=m[j]-(A[i,j]/A[i,i])*m[i];
87  m=module(transpose(matrix(m)));
[337919]88
[30f140]89  return(matrix(m));
90}
91
[8942a5]92//////////////////////////////////////////////////////////////////////////////
[2636865]93proc inner_product(vector v1,vector v2)
[30f140]94"RETURN:   inner product <v1,v2> "
95{
96  int k;
[337919]97  if (nrows(v2)>nrows(v1)) { k=nrows(v2); } else { k=nrows(v1); }
[30f140]98  return ((transpose(matrix(v1,k,1))*matrix(v2,k,1))[1,1]);
99}
100
101/////////////////////////////////////////////////////////////////////////////
102// user functions
103/////////////////////////////////////////////////////////////////////////////
104
[6188357]105proc inverse(matrix A, list #)
[963885]106"USAGE:    inverse(A [,opt]);  A a square matrix, opt integer
[337919]107RETURN:
108@format
[963885]109          a matrix:
110          - the inverse matrix of A, if A is invertible;
111          - the 1x1 0-matrix if A is not invertible (in the polynomial ring!).
112          There are the following options:
[337919]113          - opt=0 or not given: heuristically best option from below
[963885]114          - opt=1 : apply std to (transpose(E,A)), ordering (C,dp).
115          - opt=2 : apply interred (transpose(E,A)), ordering (C,dp).
[337919]116          - opt=3 : apply lift(A,E), ordering (C,dp).
[963885]117@end format
118NOTE:     parameters and minpoly are allowed; opt=2 is only correct for
119          matrices with entries in a field
120SEE ALSO: inverse_B, inverse_L
121EXAMPLE:  example inverse; shows an example
122"
[30f140]123{
[963885]124//--------------------------- initialization and check ------------------------
125   int ii,u,i,opt;
[337919]126   matrix invA;
[963885]127   int db = printlevel-voice+3;      //used for comments
128   def R=basering;
129   string mp = string(minpoly);
130   int n = nrows(A);
131   module M = A;
132   intvec v = option(get);           //get options to reset it later
133   if ( ncols(A)!=n )
134   {
135     ERROR("// ** no square matrix");
136   }
137//----------------------- choose heurisitically best option ------------------
138// This may change later, depending on improvements of the implemantation
139// at the monent we use if opt=0 or opt not given:
140// opt = 1 (std) for everything
141// opt = 2 (interred) for nothing, NOTE: interred is ok for constant matricea
142// opt = 3 (lift) for nothing
143// NOTE: interred is ok for constant matrices only (Beispiele am Ende der lib)
144   if(size(#) != 0) {opt = #[1];}
145   if(opt == 0)
146   {
147      if(npars(R) == 0)                       //no parameters
148      {
149         if( size(ideal(A-jet(A,0))) == 0 )   //constant matrix
150         {opt = 1;}
151         else
[337919]152         {opt = 1;}
[963885]153      }
154      else {opt = 1;}
155   }
156//------------------------- change ring if necessary -------------------------
157   if( ordstr(R) != "C,dp(nvars(R))" )
158   {
159     u=1;
[daa83b]160     def @R=changeord("C,dp",R);
161     setring @R;
[963885]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   {
[337919]181      option(redSB);
[963885]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      {
[337919]190         module D = std(transpose(B));
[963885]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 ----------------------
[337919]196      for (ii=1; ii<=n; ii++)
197      {
[963885]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
[30f140]218
219}
[2636865]220example
221{ "EXAMPLE:"; echo = 2;
[963885]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]=
[337919]229   y+1,  x+y,    y,
230   z,    z+1,    z,
[963885]231   y+z+2,x+y+z+2,y+z+1;
232  print(inverse(B));
233  print(B*inverse(B));
[30f140]234}
235
236//////////////////////////////////////////////////////////////////////////////
237proc sym_gauss(matrix A)
238"USAGE:    sym_gauss(A);  A = symmetric matrix
[2636865]239RETURN:   matrix, diagonalisation with symmetric gauss algorithm
[30f140]240EXAMPLE:  example sym_gauss; shows an example"
241{
[2636865]242  int i,j;
[30f140]243  int n=nrows(A);
244
245  if (ncols(A)!=n){
[6188357]246    "// ** input is not a square matrix";;
[30f140]247    return(A);
[337919]248  }
249
[30f140]250  if(!const_mat(A)){
[6188357]251    "// ** input is not a constant matrix";
[30f140]252    return(A);
253  }
[337919]254
255  if(deg(std(A-transpose(A))[1])!=-1){
[6188357]256    "// ** input is not a symmetric matrix";
[30f140]257    return(A);
258  }
259
[2636865]260  for(i=1; i<n; i++){
261    for(j=i+1; j<=n; j++){
[30f140]262      if(A[i,j]!=0){ A=red(A,i,j); }
263    }
264  }
[337919]265
[30f140]266  return(A);
267}
[2636865]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));
[30f140]274}
275
276//////////////////////////////////////////////////////////////////////////////
[337919]277proc orthogonalize(matrix A)
[30f140]278"USAGE:    orthogonalize(A); A = constant matrix
279RETURN:    matrix, orthogonal basis of the colum space of A
280EXAMPLE:   example orthogonalize; shows an example "
[2636865]281{
282  int i,j;
[30f140]283  int n=ncols(A);
284  poly k;
285
286  if(!const_mat(A)){
[6188357]287    "// ** input is not a constant matrix";
[30f140]288    matrix B;
289    return(B);
[337919]290  }
[30f140]291
292  module B=module(interred(A));
[337919]293
[30f140]294  for(i=1;i<=n;i=i+1) {
295    for(j=1;j<i;j=j+1) {
[2636865]296      k=inner_product(B[j],B[j]);
[6188357]297      if (k==0) { "Error: vector of length zero"; return(matrix(B)); }
[2636865]298      B[i]=B[i]-(inner_product(B[i],B[j])/k)*B[j];
[30f140]299    }
300  }
[337919]301
[30f140]302  return(matrix(B));
303}
[2636865]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));
[30f140]310}
311
312////////////////////////////////////////////////////////////////////////////
313proc diag_test(matrix A)
[337919]314"USAGE:          diag_test(A); A = const square matrix
[6188357]315RETURN:   int,  1 if A is diagonalisable, 0 if not
[963885]316               -1 no statement is possible, since A does not split.
[6188357]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).
[30f140]320EXAMPLE:  example diag_test; shows an example"
321{
[2636865]322  int i,j;
[337919]323  int n     = nrows(A);
[2636865]324  string mp = string(minpoly);
325  string cs = charstr(basering);
[30f140]326  int np=0;
327
[337919]328  if(ncols(A) != n) {
[2636865]329    "// input is not a square matrix";
[30f140]330    return(-1);
[337919]331  }
[30f140]332
333   if(!const_mat(A)){
[2636865]334    "// input is not a constant matrix";
[30f140]335    return(-1);
[337919]336  }
[30f140]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){
[337919]343        "// rings with parameters not allowed";
344        return(-1);
[30f140]345  }
346
[337919]347  //speichern des aktuellen Rings
[30f140]348  def BR=basering;
349  //setze R[t]
350  execute("ring rt=("+charstr(basering)+"),(@t,"+varstr(basering)+"),lp;");
351  execute("minpoly="+mp+";");
[337919]352  matrix A=imap(BR,A);
[30f140]353
354  intvec z;
355  intvec s;
[337919]356  poly X;         //characteristisches Polynom
[30f140]357  poly dXdt;      //Ableitung von X nach t
[337919]358  ideal g;              //ggT(X,dXdt)
359  poly b;              //Komponente der Busadjunkten-Matrix
[30f140]360  matrix E[n][n]; //Einheitsmatrix
[337919]361
[30f140]362  E=E+1;
363  A=E*@t-A;
364  X=det(A);
365
[337919]366  matrix Xfactors=matrix(factorize(X,1));        //zerfaellt die Matrtix ?
[30f140]367  int nf=ncols(Xfactors);
368
369  for(i=1;i<=nf;i++){
370    if(lead(Xfactors[1,i])>=@t^2){
[6188357]371      //" matrix does not split";
[30f140]372      setring BR;
373      return(-1);
374    }
[337919]375  }
[30f140]376
377  dXdt=diff(X,@t);
378  g=std(ideal(gcd(X,dXdt)));
379
380  //Busadjunkte
[337919]381  z=2..n;
[30f140]382  for(i=1;i<=n;i++){
[337919]383    s=2..n;
[30f140]384    for(j=1;j<=n;j++){
385      b=det(submat(A,z,s));
[337919]386
[30f140]387      if(0!=reduce(b,g)){
[337919]388        //" matrix not diagonalizable";
389        setring BR;
390        return(0);
[30f140]391      }
[337919]392
[30f140]393      s[j]=j;
394    }
395    z[i]=i;
396  }
[337919]397
398  //"Die Matrix ist diagonalisierbar";
[30f140]399  setring BR;
400  return(1);
401}
[2636865]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);
[30f140]408}
409
410//////////////////////////////////////////////////////////////////////////////
411proc busadj(matrix A)
[6188357]412"USAGE:   busadj(A);  A = square matrix (of size nxn)
413RETURN:  list L:
[337919]414@format
415         L[1] contains the (n+1) coefficients of the characteristic
[6188357]416              polynomial X of A, i.e.
[2636865]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
[6188357]419              the busadjoint bA = adjoint(E*t-A) of A, i.e.
[337919]420              bA = (Hn-1)*t^(n-1)+...+Hk*t^k+...+H0,  ( Hk=L[2][k+1] )
[6188357]421@end format
[30f140]422EXAMPLE: example busadj; shows an example"
423{
424  int k;
[2636865]425  int n    = nrows(A);
426  matrix E = unitmat(n);
[30f140]427  matrix H[n][n];
428  matrix B[n][n];
[2636865]429  list bA, X, L;
[30f140]430  poly a;
431
[337919]432  if(ncols(A) != n) {
[30f140]433    "input is not a square matrix";
434    return(L);
[337919]435  }
[30f140]436
[2636865]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);
[30f140]445  }
[2636865]446  B = A*bA[1];
447  a = -trace(B)/n;
448  X = insert(X,a);
[30f140]449
[2636865]450  L = insert(L,bA);
451  L = insert(L,X);
[30f140]452  return(L);
453}
[2636865]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;
[337919]461  matrix bA[2][2] = L[2][1]+L[2][2]*t;
[2636865]462  print(bA);               //the busadjoint of A;
[337919]463  print(bA*(t*unitmat(2)-A));
[30f140]464}
465
466//////////////////////////////////////////////////////////////////////////////
[2636865]467proc charpoly(matrix A, list #)
[6188357]468"USAGE:   charpoly(A[,v]); A square matrix, v string, name of a variable
[337919]469RETURN:  poly, the characteristic polynomial det(E*v-A)
[6188357]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.
[30f140]473EXAMPLE: example charpoly; shows an example"
[6188357]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)
[337919]481  {
[6188357]482    "// input is not a square matrix";
483    return(X);
[337919]484  }
[6188357]485  //---------------------- test for correct variable -------------------------
[337919]486  if( size(#)==0 ){
487    #[1] = varstr(z);
[6188357]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);
[337919]494  }
[6188357]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"
[30f140]531{
[2636865]532  int i,j;
533  string s,v;
[30f140]534  list L;
[2636865]535  int n     = nrows(A);
536  poly X    = 0;
537  def BR    = basering;
538  string mp = string(minpoly);
[30f140]539
[337919]540  if(ncols(A) != n){
[2636865]541    "// input is not a square matrix";
[30f140]542    return(X);
[337919]543  }
[30f140]544
[2636865]545  //test for correct variable
[337919]546  if( size(#)==0 ){
547    #[1] = varstr(nvars(BR));
[2636865]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);
[337919]555  }
556
[30f140]557  j=-1;
[2636865]558  for(i=1; i<=nvars(BR); i++){
[30f140]559    if(varstr(i)==v){j=i;}
560  }
561  if(j==-1){
[2636865]562    "// "+v+" is not a variable in the basering";
[30f140]563    return(X);
564  }
[337919]565
[30f140]566  //var can not be in A
567  s="Wp(";
[2636865]568  for( i=1; i<=nvars(BR); i++ ){
[30f140]569    if(i!=j){ s=s+"0";}
570    else{ s=s+"1";}
[2636865]571    if( i<nvars(BR)) {s=s+",";}
[30f140]572  }
573  s=s+")";
574
[daa83b]575  def @R=changeord(s);
576  setring @R;
[30f140]577  execute("minpoly="+mp+";");
[2636865]578  matrix A = imap(BR,A);
579  for(i=1; i<=n; i++){
[30f140]580    if(deg(lead(A)[i])>=1){
[6188357]581      "// matrix must not contain the variable "+v;
[30f140]582      kill @R;
583      setring BR;
584      return(X);
585    }
586  }
587
588  //get coefficients and build the char. poly
589  kill @R;
590  setring BR;
[2636865]591  L = busadj(A);
592  for(i=1; i<=n+1; i++){
[30f140]593    execute("X=X+L[1][i]*"+v+"^"+string(i-1)+";");
594  }
[337919]595
596  return(X);
[30f140]597}
[2636865]598example
599{ "EXAMPLE"; echo=2;
600  ring r=0,(x,t),dp;
601  matrix A[3][3]=1,x2,x,x2,6,4,x,4,1;
602  print(A);
[6188357]603  charpoly_B(A,"t");
[30f140]604}
605
606//////////////////////////////////////////////////////////////////////////////
[2636865]607proc adjoint(matrix A)
608"USAGE:    adjoint(A);  A = square matrix
[6188357]609RETURN:   adjoint matrix of A, i.e. Adj*A=det(A)*E
[30f140]610NOTE:     computation uses busadj(A)
[2636865]611EXAMPLE:  example adjoint; shows an example"
[30f140]612{
613  int n=nrows(A);
614  matrix Adj[n][n];
615  list L;
616
[337919]617  if(ncols(A) != n) {
[2636865]618    "// input is not a square matrix";
[30f140]619    return(Adj);
[337919]620  }
621
[2636865]622  L  = busadj(A);
623  Adj= (-1)^(n-1)*L[2][1];
[30f140]624  return(Adj);
[337919]625
[30f140]626}
[2636865]627example
628{ "EXAMPLE"; echo=2;
629  ring r=0,(t,x),lp;
630  matrix A[2][2]=1,x2,x,x2+3x;
631  print(A);
632  matrix Adj[2][2]=adjoint(A);
633  print(Adj);                    //Adj*A=det(A)*E
634  print(Adj*A);
[30f140]635}
636
637//////////////////////////////////////////////////////////////////////////////
638proc inverse_B(matrix A)
[6188357]639"USAGE:    inverse_B(A);  A = square matrix
640RETURN:   list Inv with
[337919]641          - Inv[1] = matrix I and
642          - Inv[2] = poly p
643          such that I*A = unitmat(n)*p;
[6188357]644NOTE:     p=1 if 1/det(A) is computable and  p=det(A) if not;
645          the computation uses busadj.
[963885]646SEE ALSO: inverse, inverse_L
[6188357]647EXAMPLE:  example inverse_B; shows an example"
[30f140]648{
649  int i;
650  int n=nrows(A);
651  matrix I[n][n];
652  poly factor;
653  list L;
654  list Inv;
[337919]655
656  if(ncols(A) != n) {
[30f140]657    "input is not a square matrix";
658    return(I);
[337919]659  }
660
[30f140]661  L=busadj(A);
[337919]662  I=module(-L[2][1]);        //+-Adj(A)
[30f140]663
[337919]664  if(reduce(1,std(L[1][1]))==0){
665    I=I*lift(L[1][1],1)[1][1];
[30f140]666    factor=1;
667  }
668  else{ factor=L[1][1];}     //=+-det(A) or 1
669  Inv=insert(Inv,factor);
670  Inv=insert(Inv,matrix(I));
671
672  return(Inv);
673}
[2636865]674example
675{ "EXAMPLE"; echo=2;
676  ring r=0,(x,y),lp;
[6188357]677  matrix A[3][3]=x,y,1,1,x2,y,x,6,0;
[2636865]678  print(A);
679  list Inv=inverse_B(A);
680  print(Inv[1]);
681  print(Inv[2]);
682  print(Inv[1]*A);
[30f140]683}
684
685//////////////////////////////////////////////////////////////////////////////
686proc det_B(matrix A)
[337919]687"USAGE:     det_B(A);  A any matrix
[30f140]688RETURN:    returns the determinant of A
689NOTE:      the computation uses the busadj algorithm
690EXAMPLE:   example det_B; shows an example"
691{
[337919]692  int n=nrows(A);
[30f140]693  list L;
694
695  if(ncols(A) != n){ return(0);}
696
697  L=busadj(A);
698  return((-1)^n*L[1][1]);
699}
[2636865]700example
[337919]701{ "EXAMPLE"; echo=2;
[2636865]702  ring r=0,(x),dp;
703  matrix A[10][10]=random(2,10,10)+unitmat(10)*x;
704  print(A);
[337919]705  det_B(A);
[30f140]706}
707
708//////////////////////////////////////////////////////////////////////////////
709proc inverse_L(matrix A)
710"USAGE:     inverse_L(A);  A = square matrix
[2636865]711RETURN:    list Inv representing a left inverse of A, i.e
[337919]712           - Inv[1] = matrix I and
[6188357]713           - Inv[2] = poly p
[337919]714           such that I*A = unitmat(n)*p;
[6188357]715NOTE:      p=1 if 1/det(A) is computable and p=det(A) if not;
[963885]716           the computation computes first det(A) and then uses lift
[337919]717SEE ALSO:  inverse, inverse_B
[30f140]718EXAMPLE:   example inverse_L; shows an example"
719{
720  int n=nrows(A);
721  matrix I;
722  matrix E[n][n]=unitmat(n);
723  poly factor;
724  poly d=1;
725  list Inv;
726
727  if (ncols(A)!=n){
[2636865]728    "// input is not a square matrix";
[30f140]729    return(I);
730  }
731
732  d=det(A);
733  if(d==0){
[2636865]734    "// matrix is not invertible";
[30f140]735    return(Inv);
736  }
737
738  // test if 1/det(A) exists
739  if(reduce(1,std(d))!=0){ E=E*d;}
[337919]740
[30f140]741  I=lift(A,E);
[337919]742  if(I==unitmat(n)-unitmat(n)){ //catch error in lift
[2636865]743    "// matrix is not invertible";
[30f140]744    return(Inv);
745  }
746
747  factor=d;      //=det(A) or 1
748  Inv=insert(Inv,factor);
749  Inv=insert(Inv,I);
[337919]750
[30f140]751  return(Inv);
752}
[2636865]753example
754{ "EXAMPLE"; echo=2;
755  ring r=0,(x,y),lp;
[6188357]756  matrix A[3][3]=x,y,1,1,x2,y,x,6,0;
[2636865]757  print(A);
758  list Inv=inverse_L(A);
759  print(Inv[1]);
760  print(Inv[2]);
761  print(Inv[1]*A);
[30f140]762}
763
764//////////////////////////////////////////////////////////////////////////////
765proc gaussred(matrix A)
[337919]766"USAGE:   gaussred(A);   A any constant matrix
[6188357]767RETURN:  list Z:  Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A)
768         gives a row reduced matrix S, a permutation matrix P and a
[337919]769         normalized lower triangular matrix U, with P*A=U*S
[6188357]770NOTE:    This procedure is designed for teaching purposes mainly.
[337919]771         The straight forward implementation in the interpreted library
772         is not very efficient (no standard basis computation).
[6188357]773EXAMPLE: example gaussred; shows an example"
[30f140]774{
[2636865]775  int i,j,l,k,jp,rang;
776  poly c,pivo;
[30f140]777  list Z;
[2636865]778  int n = nrows(A);
779  int m = ncols(A);
780  int mr= n; //max. rang
781  matrix P[n][n] = unitmat(n);
782  matrix U[n][n] = P;
[30f140]783
784  if(!const_mat(A)){
[2636865]785    "// input is not a constant matrix";
[30f140]786    return(Z);
787  }
788
789  if(n>m){mr=m;} //max. rang
790
791  for(i=1;i<=mr;i=i+1){
[337919]792    if((i+k)>m){break};
793
[30f140]794    //Test: Diagonalelement=0
795    if(A[i,i+k]==0){
796      jp=i;pivo=0;
797      for(j=i+1;j<=n;j=j+1){
[6d37e8]798        c=absValue(A[j,i+k]);
[337919]799        if(pivo<c){ pivo=c;jp=j;}
[30f140]800      }
[2636865]801      if(jp != i){       //Zeilentausch
[337919]802        for(j=1;j<=m;j=j+1){ //Zeilentausch in A (und U) (i-te mit jp-ter)
803          c=A[i,j];
804          A[i,j]=A[jp,j];
805          A[jp,j]=c;
806        }
807        for(j=1;j<=n;j=j+1){ //Zeilentausch in P
808          c=P[i,j];
809          P[i,j]=P[jp,j];
810          P[jp,j]=c;
811        }
[30f140]812      }
813      if(pivo==0){k++;continue;} //eine von selbst auftauchende Stufe !
814    }                          //i sollte im naechsten Lauf nicht erhoeht sein
[337919]815
816    //Eliminationsschritt
817    for(j=i+1;j<=n;j=j+1){
[30f140]818      c=A[j,i+k]/A[i,i+k];
819      for(l=i+k+1;l<=m;l=l+1){
[337919]820        A[j,l]=A[j,l]-A[i,l]*c;
[30f140]821      }
822      A[j,i+k]=0;  // nur wichtig falls k>0 ist
823      A[j,i]=c;    // bildet U
824    }
[337919]825  rang=i;
[30f140]826  }
[337919]827
[30f140]828  for(i=1;i<=mr;i=i+1){
829    for(j=i+1;j<=n;j=j+1){
830      U[j,i]=A[j,i];
831      A[j,i]=0;
832    }
833  }
[337919]834
[30f140]835  Z=insert(Z,rang);
836  Z=insert(Z,A);
837  Z=insert(Z,U);
838  Z=insert(Z,P);
[337919]839
[30f140]840  return(Z);
841}
[2636865]842example
843{ "EXAMPLE";echo=2;
844  ring r=0,(x),dp;
845  matrix A[5][4]=1,3,-1,4,2,5,-1,3,1,3,-1,4,0,4,-3,1,-3,1,-5,-2;
846  print(A);
847  list Z=gaussred(A);   //construct P,U,S s.t. P*A=U*S
848  print(Z[1]);          //P
849  print(Z[2]);          //U
850  print(Z[3]);          //S
851  print(Z[4]);          //rank
852  print(Z[1]*A);        //P*A
853  print(Z[2]*Z[3]);     //U*S
[30f140]854}
855
856//////////////////////////////////////////////////////////////////////////////
857proc gaussred_pivot(matrix A)
[2636865]858"USAGE:     gaussred_pivot(A);   A any constant matrix
[5c187b]859RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A)
860           gives n row reduced matrix S, a permutation matrix P and a
[2636865]861           normalized lower triangular matrix U, with P*A=U*S
[337919]862NOTE:      with row pivoting
[2636865]863EXAMPLE:   example gaussred_pivot; shows an example"
[30f140]864{
[2636865]865  int i,j,l,k,jp,rang;
866  poly c,pivo;
867  list Z;
[30f140]868  int n=nrows(A);
869  int m=ncols(A);
870  int mr=n; //max. rang
871  matrix P[n][n]=unitmat(n);
872  matrix U[n][n]=P;
873
874  if(!const_mat(A)){
[2636865]875    "// input is not a constant matrix";
[30f140]876    return(Z);
877  }
878
879  if(n>m){mr=m;} //max. rang
880
881  for(i=1;i<=mr;i=i+1){
[337919]882    if((i+k)>m){break};
883
[30f140]884    //Pivotisierung
[6d37e8]885    pivo=absValue(A[i,i+k]);jp=i;
[30f140]886    for(j=i+1;j<=n;j=j+1){
[6d37e8]887      c=absValue(A[j,i+k]);
[30f140]888      if(pivo<c){ pivo=c;jp=j;}
889    }
890    if(jp != i){ //Zeilentausch
891      for(j=1;j<=m;j=j+1){ //Zeilentausch in A (und U) (i-te mit jp-ter)
[337919]892        c=A[i,j];
893        A[i,j]=A[jp,j];
894        A[jp,j]=c;
[30f140]895      }
896      for(j=1;j<=n;j=j+1){ //Zeilentausch in P
[337919]897        c=P[i,j];
898        P[i,j]=P[jp,j];
899        P[jp,j]=c;
[30f140]900      }
901    }
902    if(pivo==0){k++;continue;} //eine von selbst auftauchende Stufe !
903                               //i sollte im naechsten Lauf nicht erhoeht sein
[337919]904    //Eliminationsschritt
905    for(j=i+1;j<=n;j=j+1){
[30f140]906      c=A[j,i+k]/A[i,i+k];
907      for(l=i+k+1;l<=m;l=l+1){
[337919]908        A[j,l]=A[j,l]-A[i,l]*c;
[30f140]909      }
910      A[j,i+k]=0;  // nur wichtig falls k>0 ist
911      A[j,i]=c;    // bildet U
912    }
[337919]913  rang=i;
[30f140]914  }
[337919]915
[30f140]916  for(i=1;i<=mr;i=i+1){
917    for(j=i+1;j<=n;j=j+1){
918      U[j,i]=A[j,i];
919      A[j,i]=0;
920    }
921  }
[337919]922
[30f140]923  Z=insert(Z,rang);
924  Z=insert(Z,A);
925  Z=insert(Z,U);
926  Z=insert(Z,P);
[337919]927
[30f140]928  return(Z);
929}
[2636865]930example
931{ "EXAMPLE";echo=2;
932  ring r=0,(x),dp;
933  matrix A[5][4] = 1, 3,-1,4,
934                   2, 5,-1,3,
935                   1, 3,-1,4,
936                   0, 4,-3,1,
937                  -3,1,-5,-2;
938  list Z=gaussred_pivot(A);  //construct P,U,S s.t. P*A=U*S
939  print(Z[1]);               //P
940  print(Z[2]);               //U
941  print(Z[3]);               //S
942  print(Z[4]);               //rank
943  print(Z[1]*A);             //P*A
944  print(Z[2]*Z[3]);          //U*S
[30f140]945}
946
947//////////////////////////////////////////////////////////////////////////////
[0b59f5]948proc gauss_nf(matrix A)
[2636865]949"USAGE:    gauss_nf(A); A any constant matrix
[6188357]950RETURN:   matrix; gauss normal form of A (uses gaussred)
[0b59f5]951EXAMPLE:  example gauss_nf; shows an example"
[30f140]952{
953  list Z;
954  if(!const_mat(A)){
[2636865]955    "// input is not a constant matrix";
[30f140]956    return(A);
957  }
[2636865]958  Z = gaussred(A);
[30f140]959  return(Z[3]);
960}
[2636865]961example
962{ "EXAMPLE";echo=2;
963  ring r = 0,(x),dp;
964  matrix A[4][4] = 1,4,4,7,2,5,5,4,4,1,1,3,0,2,2,7;
965  print(gauss_nf(A));
[30f140]966}
967
968//////////////////////////////////////////////////////////////////////////////
[0b59f5]969proc mat_rk(matrix A)
[2636865]970"USAGE:    mat_rk(A); A any constant matrix
[337919]971RETURN:   int, rank of A
[0b59f5]972EXAMPLE:  example mat_rk; shows an example"
[30f140]973{
974  list Z;
975  if(!const_mat(A)){
[6188357]976    "// input is not a constant matrix";
[30f140]977    return(-1);
978  }
[2636865]979  Z = gaussred(A);
[30f140]980  return(Z[4]);
981}
[2636865]982example
983{ "EXAMPLE";echo=2;
984  ring r = 0,(x),dp;
985  matrix A[4][4] = 1,4,4,7,2,5,5,4,4,1,1,3,0,2,2,7;
986  mat_rk(A);
[30f140]987}
988
989//////////////////////////////////////////////////////////////////////////////
[0b59f5]990proc U_D_O(matrix A)
[2636865]991"USAGE:     U_D_O(A);   constant invertible matrix A
[30f140]992RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=D , Z[4]=O
[337919]993           gives a permutation matrix P,
[2636865]994           a normalized lower triangular matrix U ,
[337919]995           a diagonal matrix D, and
[6188357]996           a normalized upper triangular matrix O
[2636865]997           with P*A=U*D*O
[6188357]998NOTE:      Z[1]=-1 means that A is not regular (proc uses gaussred)
[2636865]999EXAMPLE:   example U_D_O; shows an example"
[30f140]1000{
[2636865]1001  int i,j;
1002  list Z,L;
[30f140]1003  int n=nrows(A);
1004  matrix O[n][n]=unitmat(n);
1005  matrix D[n][n];
1006
1007  if (ncols(A)!=n){
[2636865]1008    "// input is not a square matrix";
[30f140]1009    return(Z);
1010  }
1011  if(!const_mat(A)){
[2636865]1012    "// input is not a constant matrix";
[30f140]1013    return(Z);
1014  }
[337919]1015
[30f140]1016  L=gaussred(A);
1017
1018  if(L[4]!=n){
[2636865]1019    "// input is not an invertible matrix";
[337919]1020    Z=insert(Z,-1);  //hint for calling procedures
[30f140]1021    return(Z);
1022  }
1023
1024  D=L[3];
1025
[2636865]1026  for(i=1; i<=n; i++){
1027    for(j=i+1; j<=n; j++){
1028      O[i,j] = D[i,j]/D[i,i];
1029      D[i,j] = 0;
[30f140]1030    }
1031  }
1032
1033  Z=insert(Z,O);
1034  Z=insert(Z,D);
1035  Z=insert(Z,L[2]);
1036  Z=insert(Z,L[1]);
1037  return(Z);
1038}
[2636865]1039example
1040{ "EXAMPLE";echo=2;
1041  ring r = 0,(x),dp;
[337919]1042  matrix A[5][5] = 10, 4,  0, -9,  8,
1043                   -3, 6, -6, -4,  9,
[2636865]1044                    0, 3, -1, -9, -8,
1045                   -4,-2, -6, -10,10,
1046                   -9, 5, -1, -6,  5;
[337919]1047  list Z = U_D_O(A);              //construct P,U,D,O s.t. P*A=U*D*O
[2636865]1048  print(Z[1]);                    //P
1049  print(Z[2]);                    //U
1050  print(Z[3]);                    //D
1051  print(Z[4]);                    //O
1052  print(Z[1]*A);                  //P*A
1053  print(Z[2]*Z[3]*Z[4]);          //U*D*O
[30f140]1054}
1055
1056//////////////////////////////////////////////////////////////////////////////
1057proc pos_def(matrix A)
[337919]1058"USAGE:     pos_def(A); A = constant, symmetric square matrix
1059RETURN:    int:
1060           1  if A is positive definit ,
1061           0  if not,
1062           -1 if unknown
[30f140]1063EXAMPLE:   example pos_def; shows an example"
1064{
1065  int j;
1066  list Z;
[2636865]1067  int n = nrows(A);
[30f140]1068  matrix H[n][n];
1069
1070  if (ncols(A)!=n){
[2636865]1071    "// input is not a square matrix";
[30f140]1072    return(0);
1073  }
1074  if(!const_mat(A)){
[2636865]1075    "// input is not a constant matrix";
[30f140]1076    return(-1);
[337919]1077  }
1078  if(deg(std(A-transpose(A))[1])!=-1){
[2636865]1079    "// input is not a hermitian (symmetric) matrix";
[30f140]1080    return(-1);
1081  }
[337919]1082
[0b59f5]1083  Z=U_D_O(A);
[30f140]1084
[2636865]1085  if(Z[1]==-1){
1086    return(0);
1087  }  //A not regular, therefore not pos. definit
[30f140]1088
1089  H=Z[1];
[337919]1090  //es fand Zeilentausch statt: also nicht positiv definit
[2636865]1091  if(deg(std(H-unitmat(n))[1])!=-1){
1092    return(0);
1093  }
[337919]1094
[30f140]1095  H=Z[3];
[337919]1096
[30f140]1097  for(j=1;j<=n;j=j+1){
[337919]1098    if(H[j,j]<=0){
[2636865]1099      return(0);
1100    } //eigenvalue<=0, not pos.definit
[30f140]1101  }
1102
1103  return(1); //positiv definit;
1104}
[2636865]1105example
1106{ "EXAMPLE"; echo=2;
1107  ring r = 0,(x),dp;
1108  matrix A[5][5] = 20,  4,  0, -9,   8,
1109                    4, 12, -6, -4,   9,
[337919]1110                    0, -6, -2, -9,  -8,
1111                   -9, -4, -9, -20, 10,
[2636865]1112                    8,  9, -8,  10, 10;
1113  pos_def(A);
1114  matrix B[3][3] =  3,  2,  0,
1115                    2, 12,  4,
1116                    0,  4,  2;
1117  pos_def(B);
[30f140]1118}
1119
1120//////////////////////////////////////////////////////////////////////////////
[2636865]1121proc linsolve(matrix A, matrix b)
1122"USAGE:     linsolve(A,b); A a constant nxm-matrix, b a constant nx1-matrix
[337919]1123RETURN:    a 1xm matrix X, solution of inhomogeneous linear system A*X = b
[6188357]1124           return the 0-matrix if system is not solvable
1125NOTE:      uses gaussred
[2636865]1126EXAMPLE:   example linsolve; shows an example"
[30f140]1127{
[2636865]1128  int i,j,k,rc,r;
[30f140]1129  poly c;
1130  list Z;
[2636865]1131  int n  = nrows(A);
1132  int m  = ncols(A);
1133  int n_b= nrows(b);
1134  matrix Ab[n][m+1];
1135  matrix X[m][1];
[337919]1136
[30f140]1137  if(ncols(b)!=1){
[2636865]1138    "// right hand side b is not a nx1 matrix";
[30f140]1139    return(X);
1140  }
1141
1142  if(!const_mat(A)){
[2636865]1143    "// input hand is not a constant matrix";
[30f140]1144    return(X);
[337919]1145  }
1146
[30f140]1147  if(n_b>n){
[2636865]1148    for(i=n; i<=n_b; i++){
[30f140]1149      if(b[i,1]!=0){
[337919]1150        "// right hand side b not in Image(A)";
1151        return X;
[30f140]1152      }
[337919]1153    }
[30f140]1154  }
[337919]1155
1156  if(n_b<n){
[30f140]1157    matrix copy[n_b][1]=b;
1158    matrix b[n][1]=0;
1159    for(i=1;i<=n_b;i=i+1){
1160      b[i,1]=copy[i,1];
1161    }
1162  }
[337919]1163
[0b59f5]1164  r=mat_rk(A);
[337919]1165
[30f140]1166  //1. b constant vector
[337919]1167  if(const_mat(b)){
[30f140]1168    //extend A with b
[2636865]1169    for(i=1; i<=n; i++){
1170      for(j=1; j<=m; j++){
1171        Ab[i,j]=A[i,j];
[30f140]1172      }
1173      Ab[i,m+1]=b[i,1];
1174    }
[337919]1175
[2636865]1176    //Gauss reduction
1177    Z  = gaussred(Ab);
1178    Ab = Z[3];  //normal form
[337919]1179    rc = Z[4];  //rank(Ab)
[2636865]1180    //print(Ab);
[30f140]1181
1182    if(r<rc){
[337919]1183        "// no solution";
1184        return(X);
[30f140]1185    }
[337919]1186    k=m;
[30f140]1187    for(i=r;i>=1;i=i-1){
[337919]1188
1189      j=1;
1190      while(Ab[i,j]==0){j=j+1;}// suche Ecke
1191
[30f140]1192      for(;k>j;k=k-1){ X[k]=0;}//springe zur Ecke
[337919]1193
[30f140]1194
1195      c=Ab[i,m+1]; //i-te Komponene von b
1196      for(j=m;j>k;j=j-1){
[337919]1197        c=c-X[j,1]*Ab[i,j];
[30f140]1198      }
1199      if(Ab[i,k]==0){
[337919]1200        X[k,1]=1; //willkuerlich
[30f140]1201      }
[337919]1202      else{
1203          X[k,1]=c/Ab[i,k];
[30f140]1204      }
1205      k=k-1;
1206      if(k==0){break;}
1207    }
[337919]1208
1209
[30f140]1210  }//endif (const b)
1211  else{  //b not constant
[2636865]1212    "// !not implemented!";
[337919]1213
[30f140]1214  }
1215
1216  return(X);
1217}
[2636865]1218example
1219{ "EXAMPLE";echo=2;
1220  ring r=0,(x),dp;
1221  matrix A[3][2] = -4,-6,
[337919]1222                    2, 3,
[2636865]1223                   -5, 7;
1224  matrix b[3][1] = 10,
1225                   -5,
1226                    2;
1227  matrix X = linsolve(A,b);
1228  print(X);
1229  print(A*X);
[30f140]1230}
[2636865]1231//////////////////////////////////////////////////////////////////////////////
[30f140]1232
[6188357]1233///////////////////////////////////////////////////////////////////////////////
[ecf3424]1234//    PROCEDURES for Jordan normal form
[6188357]1235//    AUTHOR:  Mathias Schulze, email: mschulze@mathematik.uni-kl.de
1236///////////////////////////////////////////////////////////////////////////////
1237
[e9124e]1238static proc rowcolswap(matrix M,int i,int j)
1239{
1240  if(i==j)
1241  {
1242    return(M);
1243  }
1244  poly p;
1245  for(int k=1;k<=nrows(M);k++)
1246  {
1247    p=M[i,k];
1248    M[i,k]=M[j,k];
1249    M[j,k]=p;
1250  }
1251  for(k=1;k<=ncols(M);k++)
1252  {
1253    p=M[k,i];
1254    M[k,i]=M[k,j];
1255    M[k,j]=p;
1256  }
1257  return(M);
1258}
1259//////////////////////////////////////////////////////////////////////////////
1260
1261static proc rowelim(matrix M,int i,int j,int k)
1262{
1263  if(jet(M[i,k],0)==0||jet(M[j,k],0)==0)
1264  {
1265    return(M);
1266  }
1267  number n=number(jet(M[i,k],0))/number(jet(M[j,k],0));
1268  for(int l=1;l<=ncols(M);l++)
1269  {
1270    M[i,l]=M[i,l]-n*M[j,l];
1271  }
1272  for(l=1;l<=nrows(M);l++)
1273  {
1274    M[l,j]=M[l,j]+n*M[l,i];
1275  }
1276  return(M);
1277}
1278///////////////////////////////////////////////////////////////////////////////
1279
1280static proc colelim(matrix M,int i,int j,int k)
1281{
1282  if(jet(M[k,i],0)==0||jet(M[k,j],0)==0)
1283  {
1284    return(M);
1285  }
1286  number n=number(jet(M[k,i],0))/number(jet(M[k,j],0));
1287  for(int l=1;l<=nrows(M);l++)
1288  {
1289    M[l,i]=M[l,i]-n*M[l,j];
1290  }
1291  for(l=1;l<=ncols(M);l++)
1292  {
1293    M[j,l]=M[j,l]+n*M[i,l];
1294  }
1295  return(M);
1296}
1297///////////////////////////////////////////////////////////////////////////////
1298
1299proc hessenberg(matrix M)
1300"USAGE:   hessenberg(M); matrix M
1301ASSUME:  M constant square matrix
1302RETURN:  matrix H;  Hessenberg form of M
1303EXAMPLE: example hessenberg; shows examples
1304"
1305{
1306  if(system("with","eigenval"))
1307  {
[cb40b5]1308    return(system("hessenberg",M));
[e9124e]1309  }
1310
1311  int n=ncols(M);
1312  int i,j;
1313  for(int k=1;k<n-1;k++)
1314  {
1315    j=k+1;
1316    while(j<n&&jet(M[j,k],0)==0)
1317    {
1318      j++;
1319    }
1320    if(jet(M[j,k],0)!=0)
1321    {
1322      M=rowcolswap(M,j,k+1);
1323      for(i=j+1;i<=n;i++)
1324      {
1325        M=rowelim(M,i,k+1,k);
1326      }
1327    }
1328  }
1329  return(M);
1330}
1331example
1332{ "EXAMPLE:"; echo=2;
1333  ring R=0,x,dp;
1334  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
1335  print(M);
1336  print(hessenberg(M));
1337}
1338///////////////////////////////////////////////////////////////////////////////
1339
[275721f]1340proc eigenvals(matrix M)
1341"USAGE:   eigenvals(M); matrix M
1342ASSUME:  eigenvalues of M in basefield
[ecf3424]1343RETURN:
1344@format
[275721f]1345list l;
[f91f7a6]1346  ideal l[1];
1347    number l[1][i];  i-th eigenvalue of M
[275721f]1348  intvec l[2];
[f91f7a6]1349    int l[2][i];  multiplicity of i-th eigenvalue of M
[ecf3424]1350@end format
[275721f]1351EXAMPLE: example eigenvals; shows examples
[ecf3424]1352"
1353{
[e9124e]1354  if(system("with","eigenval"))
1355  {
[cb40b5]1356    return(system("eigenvals",jet(M,0)));
[e9124e]1357  }
1358
1359  M=jet(hessenberg(M),0);
1360  int n=ncols(M);
1361  int k;
1362  ideal e;
1363  intvec m;
1364  number e0;
1365  intvec v;
1366  list l;
1367  int i,j;
1368  j=1;
1369  while(j<=n)
1370  {
1371    v=j;
1372    j++;
1373    if(j<=n)
1374    {
1375      while(j<n&&M[j,j-1]!=0)
1376      {
1377        v=v,j;
1378        j++;
1379      }
1380      if(M[j,j-1]!=0)
1381      {
1382        v=v,j;
1383        j++;
1384      }
1385    }
1386    if(size(v)==1)
1387    {
1388      k++;
1389      e[k]=M[v,v];
1390      m[k]=1;
1391    }
1392    else
1393    {
1394      l=factorize(det(submat(M,v,v)-var(1)));
1395      for(i=size(l[1]);i>=1;i--)
1396      {
1397        e0=number(jet(l[1][i]/var(1),0));
1398        if(e0!=0)
1399        {
1400          k++;
1401          e[k]=(e0*var(1)-l[1][i])/e0;
1402          m[k]=l[2][i];
1403        }
1404      }
1405    }
1406  }
[cb40b5]1407  return(spnf(list(e,m)));
[ecf3424]1408}
1409example
1410{ "EXAMPLE:"; echo=2;
1411  ring R=0,x,dp;
1412  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
1413  print(M);
[275721f]1414  eigenvals(M);
[ecf3424]1415}
1416///////////////////////////////////////////////////////////////////////////////
1417
[2699a6]1418proc minipoly(matrix M,list #)
1419"USAGE:   minpoly(M); matrix M
1420ASSUME:  eigenvalues of M in basefield
1421RETURN:
1422@format
1423list l;  minimal polynomial of M
1424  ideal l[1];
1425    number l[1][i];  i-th root of minimal polynomial of M
1426  intvec l[2];
1427    int l[2][i];  multiplicity of i-th root of minimal polynomial of M
1428@end format
1429EXAMPLE: example minipoly; shows examples
1430"
1431{
1432  if(nrows(M)==0)
1433  {
1434    ERROR("non empty expected");
1435  }
1436  if(ncols(M)!=nrows(M))
1437  {
1438    ERROR("square matrix expected");
1439  }
1440
1441  M=jet(M,0);
1442
1443  if(size(#)==0)
1444  {
1445    #=eigenvals(M);
1446  }
1447  def e0,m0=#[1..2];
1448
1449  intvec m1;
1450  matrix N0,N1;
1451  for(int i=1;i<=ncols(e0);i++)
1452  {
1453    m1[i]=1;
1454    N0=M-e0[i];
1455    N1=N0;
1456    while(size(syz(N1))<m0[i])
1457    {
1458      m1[i]=m1[i]+1;
1459      N1=N1*N0;
1460    }
1461  }
1462
1463  return(list(e0,m1));
1464}
1465example
1466{ "EXAMPLE:"; echo=2;
1467  ring R=0,x,dp;
1468  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
1469  print(M);
1470  minipoly(M);
1471}
1472///////////////////////////////////////////////////////////////////////////////
1473
[348bbc]1474proc spnf(list #)
[cb40b5]1475"USAGE:   spnf(list(a[,m])); ideal a, intvec m
1476ASSUME:  ncols(a)==size(m)
1477RETURN:  order a[i] with multiplicity m[i] lexicographically
1478EXAMPLE: example spnf; shows examples
1479"
1480{
[348bbc]1481  list sp=#;
[cb40b5]1482  ideal a=sp[1];
1483  int n=ncols(a);
1484  intvec m;
1485  list V;
1486  module v;
1487  int i,j;
1488  for(i=2;i<=size(sp);i++)
1489  {
1490    if(typeof(sp[i])=="intvec")
1491    {
1492      m=sp[i];
1493    }
1494    if(typeof(sp[i])=="module")
1495    {
1496      v=sp[i];
1497      for(j=n;j>=1;j--)
1498      {
1499        V[j]=module(v[j]);
1500      }
1501    }
1502    if(typeof(sp[i])=="list")
1503    {
1504      V=sp[i];
1505    }
1506  }
1507  if(m==0)
1508  {
1509    for(i=n;i>=1;i--)
1510    {
1511      m[i]=1;
1512    }
1513  }
1514
1515  int k;
1516  ideal a0;
1517  intvec m0;
1518  list V0;
1519  number a1;
1520  int m1;
1521  for(i=n;i>=1;i--)
1522  {
1523    if(m[i]!=0)
1524    {
1525      for(j=i-1;j>=1;j--)
1526      {
1527        if(m[j]!=0)
1528        {
1529          if(number(a[i])>number(a[j]))
1530          {
1531            a1=number(a[i]);
1532            a[i]=a[j];
1533            a[j]=a1;
1534            m1=m[i];
1535            m[i]=m[j];
1536            m[j]=m1;
1537            if(size(V)>0)
1538            {
1539              v=V[i];
1540              V[i]=V[j];
1541              V[j]=v;
1542            }
1543          }
1544          if(number(a[i])==number(a[j]))
1545          {
1546            m[i]=m[i]+m[j];
1547            m[j]=0;
1548            if(size(V)>0)
1549            {
1550              V[i]=V[i]+V[j];
1551            }
1552          }
1553        }
1554      }
1555      k++;
1556      a0[k]=a[i];
1557      m0[k]=m[i];
1558      if(size(V)>0)
1559      {
1560        V0[k]=V[i];
1561      }
1562    }
1563  }
1564
1565  if(size(V0)>0)
1566  {
1567    n=size(V0);
1568    module U=std(V0[n]);
1569    for(i=n-1;i>=1;i--)
1570    {
1571      V0[i]=simplify(reduce(V0[i],U),1);
1572      if(i>=2)
1573      {
1574        U=std(U+V0[i]);
1575      }
1576    }
1577  }
1578
1579  if(k>0)
1580  {
1581    sp=a0,m0;
1582    if(size(V0)>0)
1583    {
1584      sp[3]=V0;
1585    }
1586  }
1587  return(sp);
1588}
1589example
1590{ "EXAMPLE:"; echo=2;
1591  ring R=0,(x,y),ds;
1592  list sp=list(ideal(-1/2,-3/10,-3/10,-1/10,-1/10,0,1/10,1/10,3/10,3/10,1/2));
1593  spprint(spnf(sp));
1594}
1595///////////////////////////////////////////////////////////////////////////////
1596
1597proc spprint(list sp)
1598"USAGE:   spprint(sp); list sp
1599RETURN:  string s;  spectrum sp
1600EXAMPLE: example spprint; shows examples
1601"
1602{
1603  string s;
1604  for(int i=1;i<size(sp[2]);i++)
1605  {
1606    s=s+"("+string(sp[1][i])+","+string(sp[2][i])+"),";
1607  }
1608  s=s+"("+string(sp[1][i])+","+string(sp[2][i])+")";
1609  return(s);
1610}
1611example
1612{ "EXAMPLE:"; echo=2;
1613  ring R=0,(x,y),ds;
1614  list sp=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
1615  spprint(sp);
1616}
1617///////////////////////////////////////////////////////////////////////////////
1618
[a0c62d]1619proc jordan(matrix M,list #)
[65b27c]1620"USAGE:   jordan(M); matrix M
[275721f]1621ASSUME:  eigenvalues of M in basefield
[65b27c]1622RETURN:
[337919]1623@format
[275721f]1624list l;  Jordan data of M
1625  ideal l[1];
1626    number l[1][i];  eigenvalue of i-th Jordan block of M
1627  intvec l[2];
1628    int l[2][i];  size of i-th Jordan block of M
1629  intvec l[3];
1630    int l[3][i];  multiplicity of i-th Jordan block of M
[6188357]1631@end format
[ecf3424]1632EXAMPLE: example jordan; shows examples
[6188357]1633"
1634{
[65b27c]1635  if(nrows(M)==0)
[6188357]1636  {
[65b27c]1637    ERROR("non empty expected");
[6188357]1638  }
[65b27c]1639  if(ncols(M)!=nrows(M))
[6188357]1640  {
[65b27c]1641    ERROR("square matrix expected");
[6188357]1642  }
1643
1644  M=jet(M,0);
1645
[a0c62d]1646  if(size(#)==0)
1647  {
[275721f]1648    #=eigenvals(M);
[a0c62d]1649  }
1650  def e0,m0=#[1..2];
[6188357]1651
[65b27c]1652  int i;
[4b6c75]1653  for(i=1;i<=ncols(e0);i++)
[6188357]1654  {
[275721f]1655    if(deg(e0[i])>0)
[65b27c]1656    {
[275721f]1657
[65b27c]1658      ERROR("eigenvalues in coefficient field expected");
1659      return(list());
1660    }
[6188357]1661  }
1662
[65b27c]1663  int j,k;
[4b6c75]1664  matrix N0,N1;
[65b27c]1665  module K0;
1666  list K;
[4b6c75]1667  ideal e;
1668  intvec s,m;
[6188357]1669
[4b6c75]1670  for(i=1;i<=ncols(e0);i++)
[6188357]1671  {
[0ebbcf4]1672    N0=M-e0[i]*matrix(freemodule(ncols(M)));
[6188357]1673
[4b6c75]1674    N1=N0;
[65b27c]1675    K0=0;
1676    K=module();
[4b6c75]1677    while(size(K0)<m0[i])
[6188357]1678    {
[4b6c75]1679      K0=syz(N1);
[65b27c]1680      K=K+list(K0);
[4b6c75]1681      N1=N1*N0;
[6188357]1682    }
1683
[4b6c75]1684    for(j=2;j<size(K);j++)
[6188357]1685    {
[4b6c75]1686      if(2*size(K[j])-size(K[j-1])-size(K[j+1])>0)
[6188357]1687      {
[4b6c75]1688        k++;
1689        e[k]=e0[i];
1690        s[k]=j-1;
1691        m[k]=2*size(K[j])-size(K[j-1])-size(K[j+1]);
[6188357]1692      }
1693    }
[4b6c75]1694    if(size(K[j])-size(K[j-1])>0)
1695    {
1696      k++;
1697      e[k]=e0[i];
1698      s[k]=j-1;
1699      m[k]=size(K[j])-size(K[j-1]);
1700    }
[6188357]1701  }
1702
[4b6c75]1703  return(list(e,s,m));
[65b27c]1704}
1705example
1706{ "EXAMPLE:"; echo=2;
1707  ring R=0,x,dp;
1708  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
1709  print(M);
1710  jordan(M);
1711}
1712///////////////////////////////////////////////////////////////////////////////
1713
[a0c62d]1714proc jordanbasis(matrix M,list #)
[fa01b7]1715"USAGE:   jordanbasis(M); matrix M
[275721f]1716ASSUME:  eigenvalues of M in basefield
[4b6c75]1717RETURN:
1718@format
1719list l:
[91fc5e]1720  module l[1];  inverse(l[1])*M*l[1] in Jordan normal form
[275721f]1721  intvec l[2];
1722    int l[2][i];  weight filtration index of l[1][i]
[4b6c75]1723@end format
[ecf3424]1724EXAMPLE: example jordanbasis; shows examples
[65b27c]1725"
1726{
1727  if(nrows(M)==0)
[6188357]1728  {
[65b27c]1729    ERROR("non empty matrix expected");
[6188357]1730  }
[65b27c]1731  if(ncols(M)!=nrows(M))
[6188357]1732  {
[65b27c]1733    ERROR("square matrix expected");
[6188357]1734  }
1735
[65b27c]1736  M=jet(M,0);
1737
[a0c62d]1738  if(size(#)==0)
1739  {
[275721f]1740    #=eigenvals(M);
[a0c62d]1741  }
1742  def e,m=#[1..2];
[6188357]1743
[61549b]1744  for(int i=1;i<=ncols(e);i++)
[6188357]1745  {
[0ebbcf4]1746    if(deg(e[i])>0)
[6188357]1747    {
[65b27c]1748      ERROR("eigenvalues in coefficient field expected");
1749      return(freemodule(ncols(M)));
[6188357]1750    }
1751  }
[65b27c]1752
[61549b]1753  int j,k,l,n;
[fa01b7]1754  matrix N0,N1;
[65b27c]1755  module K0,K1;
[6188357]1756  list K;
[65b27c]1757  matrix u[ncols(M)][1];
1758  module U;
[4b6c75]1759  intvec w;
[6188357]1760
[61549b]1761  for(i=1;i<=ncols(e);i++)
[6188357]1762  {
[0ebbcf4]1763    N0=M-e[i]*matrix(freemodule(ncols(M)));
[6188357]1764
[fa01b7]1765    N1=N0;
[4b6c75]1766    K0=0;
1767    K=list();
[65b27c]1768    while(size(K0)<m[i])
[6188357]1769    {
[fa01b7]1770      K0=syz(N1);
[65b27c]1771      K=K+list(K0);
[fa01b7]1772      N1=N1*N0;
[6188357]1773    }
1774
[65b27c]1775    K1=0;
[4b6c75]1776    for(j=1;j<size(K);j++)
[6188357]1777    {
[65b27c]1778      K0=K[j];
[fa01b7]1779      K[j]=interred(reduce(K[j],std(K1+module(N0*K[j+1]))));
[65b27c]1780      K1=K0;
[6188357]1781    }
[65b27c]1782    K[j]=interred(reduce(K[j],std(K1)));
[6188357]1783
[4b6c75]1784    for(l=size(K);l>=1;l--)
[6188357]1785    {
[4b6c75]1786      for(k=size(K[l]);k>0;k--)
[6188357]1787      {
[4b6c75]1788        u=K[l][k];
1789        for(j=l;j>=1;j--)
[6188357]1790        {
[61549b]1791          U=U+module(u);
1792          n++;
1793          w[n]=2*j-l-1;
[fa01b7]1794          u=N0*u;
[6188357]1795        }
1796      }
1797    }
1798  }
[61549b]1799
[4b6c75]1800  return(list(U,w));
[6188357]1801}
1802example
1803{ "EXAMPLE:"; echo=2;
1804  ring R=0,x,dp;
1805  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
1806  print(M);
[4b6c75]1807  list l=jordanbasis(M);
1808  print(l[1]);
1809  print(l[2]);
1810  print(inverse(l[1])*M*l[1]);
[6188357]1811}
1812///////////////////////////////////////////////////////////////////////////////
1813
[cb40b5]1814proc jordanmatrix(list jd)
1815"USAGE:   jordanmatrix(list(e,s,m)); ideal e, intvec s, intvec m
[f91f7a6]1816ASSUME:  ncols(e)==size(s)==size(m)
[4b6c75]1817RETURN:
1818@format
[26a4bb]1819matrix J;  Jordan matrix with list(e,s,m)==jordan(J)
[4b6c75]1820@end format
[ecf3424]1821EXAMPLE: example jordanmatrix; shows examples
[6188357]1822"
1823{
[cb40b5]1824  ideal e=jd[1];
1825  intvec s=jd[2];
1826  intvec m=jd[3];
[4f1139]1827  if(ncols(e)!=size(s)||ncols(e)!=size(m))
[6188357]1828  {
[65b27c]1829    ERROR("arguments of equal size expected");
[6188357]1830  }
1831
[4b6c75]1832  int i,j,k,l;
1833  int n=int((transpose(matrix(s))*matrix(m))[1,1]);
[6188357]1834  matrix J[n][n];
[4b6c75]1835  for(k=1;k<=ncols(e);k++)
[6188357]1836  {
[4b6c75]1837    for(l=1;l<=m[k];l++)
[6188357]1838    {
[4b6c75]1839      j++;
1840      J[j,j]=e[k];
1841      for(i=s[k];i>=2;i--)
[6188357]1842      {
[61549b]1843        J[j+1,j]=1;
[4b6c75]1844        j++;
1845        J[j,j]=e[k];
[6188357]1846      }
1847    }
1848  }
1849
1850  return(J);
1851}
1852example
1853{ "EXAMPLE:"; echo=2;
1854  ring R=0,x,dp;
[65b27c]1855  ideal e=ideal(2,3);
[4b6c75]1856  intvec s=1,2;
1857  intvec m=1,1;
[cb40b5]1858  print(jordanmatrix(list(e,s,m)));
[6188357]1859}
1860///////////////////////////////////////////////////////////////////////////////
1861
[275721f]1862proc jordannf(matrix M,list #)
1863"USAGE:   jordannf(M); matrix M
1864ASSUME:  eigenvalues of M in basefield
1865RETURN:  matrix J; Jordan normal form of M
1866EXAMPLE: example jordannf; shows examples
[6188357]1867"
1868{
[cb40b5]1869  return(jordanmatrix(jordan(M,#)));
[6188357]1870}
1871example
1872{ "EXAMPLE:"; echo=2;
1873  ring R=0,x,dp;
1874  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
1875  print(M);
[275721f]1876  print(jordannf(M));
[6188357]1877}
[e9124e]1878
[963885]1879///////////////////////////////////////////////////////////////////////////////
1880
1881/*
1882///////////////////////////////////////////////////////////////////////////////
1883//          Auskommentierte zusaetzliche Beispiele
1884//
[6188357]1885///////////////////////////////////////////////////////////////////////////////
[963885]1886// Singular for ix86-Linux version 1-3-10  (2000121517)  Dec 15 2000 17:55:12
1887// Rechnungen auf AMD700 mit 632 MB
1888
1889 LIB "linalg.lib";
1890
18911. Sparse integer Matrizen
1892--------------------------
1893ring r1=0,(x),dp;
1894system("--random", 12345678);
1895int n = 70;
[337919]1896matrix m = sparsemat(n,n,50,100);
[963885]1897option(prot,mem);
1898
1899int t=timer;
[337919]1900matrix im = inverse(m,1)[1];
[963885]1901timer-t;
1902print(im*m);
[337919]1903//list l0 = watchdog(100,"inverse("+"m"+",3)");
[963885]1904//bricht bei 100 sec ab und gibt l0[1]: string Killed zurueck
1905
[337919]1906//inverse(m,1): std       5sec 5,5 MB
[963885]1907//inverse(m,2): interred  12sec
1908//inverse(m,2): lift      nach 180 sec 13MB abgebrochen
1909//n=60: linalgorig: 3  linalg: 5
[337919]1910//n=70: linalgorig: 6,7 linalg: 11,12
1911// aber linalgorig rechnet falsch!
[963885]1912
19132. Sparse poly Matrizen
1914-----------------------
1915ring r=(0),(a,b,c),dp;
1916system("--random", 12345678);
1917int n=6;
1918matrix m = sparsematrix(n,n,2,0,50,50,9); //matrix of polys of deg <=2
1919option(prot,mem);
1920
1921int t=timer;
[337919]1922matrix im = inverse(m);
[963885]1923timer-t;
1924print(im*m);
1925//inverse(m,1): std       0sec 1MB
1926//inverse(m,2): interred  0sec 1MB
1927//inverse(m,2): lift      nach 2000 sec 33MB abgebrochen
1928
19293. Sparse Matrizen mit Parametern
1930---------------------------------
1931//liborig rechnet hier falsch!
1932ring r=(0),(a,b),dp;
1933system("--random", 12345678);
1934int n=7;
1935matrix m = sparsematrix(n,n,1,0,40,50,9);
1936ring r1 = (0,a,b),(x),dp;
1937matrix m = imap(r,m);
1938option(prot,mem);
1939
1940int t=timer;
[337919]1941matrix im = inverse(m);
[963885]1942timer-t;
1943print(im*m);
1944//inverse(m)=inverse(m,3):15 sec inverse(m,1)=1sec inverse(m,2):>120sec
1945//Bei Parametern vergeht die Zeit beim Normieren!
1946
19473. Sparse Matrizen mit Variablen und Parametern
1948-----------------------------------------------
1949ring r=(0),(a,b),dp;
1950system("--random", 12345678);
1951int n=6;
1952matrix m = sparsematrix(n,n,1,0,35,50,9);
1953ring r1 = (0,a),(b),dp;
1954matrix m = imap(r,m);
1955option(prot,mem);
1956
1957int t=timer;
[337919]1958matrix im = inverse(m,3);
[963885]1959timer-t;
1960print(im*m);
1961//n=7: inverse(m,3):lange sec inverse(m,1)=1sec inverse(m,2):1sec
1962
19634. Ueber Polynomring invertierbare Matrizen
1964-------------------------------------------
1965LIB"random.lib"; LIB"linalg.lib";
1966system("--random", 12345678);
1967int n =3;
1968ring r= 0,(x,y,z),(C,dp);
[337919]1969matrix A=triagmatrix(n,n,1,0,0,50,2);
1970intmat B=sparsetriag(n,n,20,1);
[963885]1971matrix M = A*transpose(B);
1972M=M*transpose(M);
1973M[1,1..ncols(M)]=M[1,1..n]+xyz*M[n,1..ncols(M)];
1974print(M);
1975//M hat det=1 nach Konstruktion
1976
1977int t=timer;
1978matrix iM=inverse(M);
1979timer-t;
1980print(iM*M);                      //test
1981
1982//ACHTUNG: Interred liefert i.A. keine Inverse, Gegenbeispiel z.B.
1983//mit n=3
1984//eifacheres Gegenbeispiel:
1985matrix M =
[337919]19869yz+3y+3z+2,             9y2+6y+1,
[963885]19879xyz+3xy+3xz-9z2+2x-6z-1,9xy2+6xy-9yz+x-3y-3z
1988//det M=1, inverse(M,2); ->// ** matrix is not invertible
1989//lead(M); 9xyz*gen(2) 9xy2*gen(2) nicht teilbar!
1990
19915. charpoly:
1992-----------
[337919]1993//ring rp=(0,A,B,C),(x),dp;
1994ring r=0,(A,B,C,x),dp;
[963885]1995matrix m[12][12]=
[337919]1996AC,BC,-3BC,0,-A2+B2,-3AC+1,B2,   B2,  1,   0, -C2+1,0,
19971, 1, 2C,  0,0,     B,     -A,   -4C, 2A+1,0, 0,    0,
19980, 0, 0,   1,0,     2C+1,  -4C+1,-A,  B+1, 0, B+1,  3B,
[963885]1999AB,B2,0,   1,0,     1,     0,    1,   A,   0, 1,    B+1,
[337919]20001, 0, 1,   0,0,     1,     0,    -C2, 0,   1, 0,    1,
20010, 0, 2,   1,2A,    1,     0,    0,   0,   0, 1,    1,
20020, 1, 0,   1,1,     2,     A,    3B+1,1,   B2,1,    1,
20030, 1, 0,   1,1,     1,     1,    1,   2,   0, 0,    0,
20041, 0, 1,   0,0,     0,     1,    0,   1,   1, 0,    3,
20051, 3B,B2+1,0,0,     1,     0,    1,   0,   0, 1,    0,
20060, 0, 1,   0,0,     0,     0,    1,   0,   0, 0,    0,
20070, 1, 0,   1,1,     3,     3B+1, 0,   1,   1, 1,    0;
[963885]2008option(prot,mem);
2009
2010int t=timer;
2011poly q=charpoly(m,"x");         //1sec, charpoly_B 1sec, 16MB
2012timer-t;
2013//1sec, charpoly_B 1sec, 16MB  (gleich in r und rp)
2014
[337919]2015*/
Note: See TracBrowser for help on using the repository browser.