Changeset 2b9b32 in git for Singular/LIB/control.lib

Timestamp:

Mar 16, 2005, 11:16:53 PM (19 years ago)

Author:

Viktor Levandovskyy <levandov@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

42d5aa6a6beb385d8854055d86c22cf5fd0d3670

Parents:

eccbb153088fe98a0348384fdd90e9bf6414e860

Message:

*levandov: Evas code for genericity and several changes by Markus


git-svn-id: file:///usr/local/Singular/svn/trunk@7784 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• Singular/LIB/control.lib (modified) (10 diffs)

Singular/LIB/control.lib

@@ -1 @@
-version="$Id: control.lib,v 1.24 2005-03-16 16:15:14 levandov Exp $";
+version="$Id: control.lib,v 1.25 2005-03-16 22:16:53 levandov Exp $";
 category="Applications";
 info="
@@ -577 +577 @@
 "
 {
+  if( nrows(R) != rank(transpose(R)) )
+  {
+    return ("control2 cannot be applied, since R does not have full row rank");
+  }  
   intvec v=option(get);
   module R_std=std(R);
@@ -607 +611 @@
 //------------------------------------------------------------------------
 proc rank(module M)
-{
-  module M_red=bareiss(M)[1];
-  int NCols_red=ncols(M_red);
+"USAGE: proc rank(M), M a matrix/module
+PURPOSE: compute the column rank of M as of matrix
+RETURN: int
+NOTE: this procedure uses bareiss-algorithm which might not terminate
+"
+{
+  module M_red  = bareiss(M)[1];
+  int NCols_red = ncols(M_red);
   return (NCols_red);
 }
+example
+{"EXAMPLE: ";echo = 2;
+  //deRham complex
+  ring r=0,(D(1..3)),dp;
+  module R;
+  R=[0,-D(3),D(2)],
+    [D(3),0,-D(1)],
+    [-D(2),D(1),0];
+  R=transpose(R);
+  rank(R); 
+};
+  
 //------------------------------------------------------------------------
 static proc autonom_output( int i, int NVars, module RC, int R_rank ) 
 //static proc autonom_output( int i, int NVars, module RC ) 
 //static proc autonom_output( int i, int NVars ) 
-"USAGE:  proc autonom_output(i, NVars)
+"USAGE:  proc autonom_output(i, NVars, RC, R_rank)
            i:  integer, number of first nonzero Ext or 
                just number of variables in a base ring + 1 in case that all the Exts are zero
            NVars:  integer, number of variables in a base ring  
-RETURN:  list with all the autonomy properties of the system which is to be returned in 'autonom' procedure
+           RC: module, kernel-representation of controllable part of the system
+           R_rank: integer, rank of the representation matrix
+PURPOSE: compute all the autonomy properties of the system which is to be returned in 'autonom' procedure
+RETURN:  list 
 NOTE:  this procedure is used in 'autonom' procedure
 "
@@ -705 +729 @@
   module RT = transpose(R);
   module RC;  //for computation of controllable part if if exists
-  d = dim_Our( std(RT) );  //this is the dimension of the system 
-  int i=NVars-d;  //First non-zero Ext
-  if(d==0)
+  d     = dim_Our( std(RT) );  //this is the dimension of the system 
+  int i = NVars-d;  //First non-zero Ext
+  if( d==0 )
     {
       RC=LeftKernel(RightKernel(R));
@@ -728 +752 @@
 proc autonom(module R)
 "USAGE:  autonom(R);   R a module (R is a matrix of the system of equations which is to be investigated)
-RETURN:  list of all the properties concerning autonomy of the system (behavior) represented by the  matrix R
-EXAMPLE:  example autonom; shows an example
+PURPOSE: find all the properties concerning autonomy of the system (behavior) represented by the  matrix R
+RETURN:  list
+EXAMPLE: example autonom; shows an example
 "
 {
@@ -908 +933 @@
 };
 //----------------------------------------------------------
-
 proc genericity(matrix M)
 "USAGE:  genericity(M), M is a module/matrix
@@ -923 +947 @@
     return("-");
   }
-  int plevel = printlevel-voice+2;
-  // M is a matrix over a ring with params and vars;
-  ideal I = ideal(M); // a list of entries
-  I = simplify(I,2); // delete 0's
-  // decompose every coeff in every poly
-  int i;
-  int s = size(I);
-  ideal NM;
-  poly p;
-  number num;
-  int  cl=1;
-  intvec ZeroVec; ZeroVec[nvars(basering)] = 0;
-  intvec W;
-  ideal Numero, Denomiro;
-  int cNu=0; int cDe=0;
-  for (i=1; i<=s; i++)
-  {
-    // remove contents and add them as polys
-    p   = I[i];
-    W   = leadexp(p);
-    if (W == ZeroVec) // i.e. just a coef
-    {
-      num    = denominator(leadcoef(p)); // from poly.lib
-      NM[cl] = numerator(leadcoef(p));
-      dbprint(p,"numerator:");
-      dbprint(p, string(NM[cl]));
-      cNu++; Numero[cNu]= NM[cl];
-      cl++;
-      NM[cl] = num; // denominator
-      dbprint(p,"denominator:");
-      dbprint(p, string(NM[cl]));
-      cDe++; Denomiro[cDe]= NM[cl];
-      cl++;
-      p = p - lead(p); // for the next cycle
-    }      
-    if ( p!= 0)
-    {
-      num = content(p);
-      p   = p/num;
-      NM[cl] = denominator(num);
-      dbprint(p,"content denominator:");
-      dbprint(p, string(NM[cl]));
-      cNu++; Numero[cNu]= NM[cl];
-      cl++;
-      NM[cl] = numerator(num);
-      dbprint(p,"content numerator:");
-      dbprint(p, string(NM[cl]));
-      cDe++; Denomiro[cDe]= NM[cl];
-      cl++;
-    }
-    // it seems that the next elements will not have real influence
-    while( p != 0)
-    {
-      NM[cl] = leadcoef(p); // should be all integer, i.e. non-rational
-      dbprint(p,"coef:");
-      dbprint(p, string(NM[cl]));
-      cl++;
-      p = p - lead(p);
-    }
-  }
-  NM = simplify(NM,4); // delete identical
-  string newvars = parstr(basering);
-  def save = basering;
-  string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;";
-  execute(NewRing);
-  // get params as variables
-  // create a list of non-monomials
-  ideal @L;
-  ideal F;
-  ideal NM = imap(save,NM);
-  NM = simplify(NM,8); //delete multiples
-  poly p,q;
-  cl = 1;
-  int j, cf;
-  for(i=1; i<=size(NM);i++)
-  {
-    p = NM[i] - lead(NM[i]);
-    if (p!=0)
-    {
-      //      L[cl] = p;
-      F  = factorize(NM[i],1); //non-constant factors only
-      cf = 1;
-      // factorize every polynomial
-      // throw away every monomial from factorization (also constants from above ring)
-      for (j=1; j<=size(F);j++)
-      {
-        q = F[j]-lead(F[j]);
-        if (q!=0)
-        {
-          @L[cl] = F[j];
-          cl++;
-        }
-      }
-    }
-  }
-  // return the result [in string-format]
-  @L = simplify(@L,2+4+8); // skip zeroes, doubled and entries, diff. by a constant
-  list SL;
-  for (j=1; j<=size(@L);j++)
-  {
-    SL[j] = string(@L[j]);
-  }
-  setring save;
-  return(SL);
+  list RT = evas_genericity(M);
+  return(RT);
 }
 example
@@ -1049 +971 @@
 }
 //---------------------------------------------------------------
+proc victors_genericity(matrix M)
+"USAGE:  genericity(M), M is a module/matrix
+PURPOSE: determine parametric expressions which have been assumed to be non-zero in the process of computing the Groebner basis
+RETURN:  list (of strings) 
+NOTE: we strongly recommend to switch on the redSB and redTail options; 
+@*    the procedure is effective with the lift procedure for modules with parameters
+EXAMPLE:  example genericity; shows an example
+"
+{
+  // returns "-", if there are no parameters!
+  if (npars(basering)==0)
+  {
+    return("-");
+  }
+  int plevel = printlevel-voice+2;
+  // M is a matrix over a ring with params and vars;
+  ideal I = ideal(M); // a list of entries
+  I = simplify(I,2); // delete 0's
+  // decompose every coeff in every poly
+  int i;
+  int s = size(I);
+  ideal NM;
+  poly p;
+  number num;
+  int  cl=1;
+  intvec ZeroVec; ZeroVec[nvars(basering)] = 0;
+  intvec W;
+  ideal Numero, Denomiro;
+  int cNu=0; int cDe=0;
+  for (i=1; i<=s; i++)
+  {
+    // remove contents and add them as polys
+    p   = I[i];
+    W   = leadexp(p);
+    if (W == ZeroVec) // i.e. just a coef
+    {
+      num    = denominator(leadcoef(p)); // from poly.lib
+      NM[cl] = numerator(leadcoef(p));
+      dbprint(p,"numerator:");
+      dbprint(p, string(NM[cl]));
+      cNu++; Numero[cNu]= NM[cl];
+      cl++;
+      NM[cl] = num; // denominator
+      dbprint(p,"denominator:");
+      dbprint(p, string(NM[cl]));
+      cDe++; Denomiro[cDe]= NM[cl];
+      cl++;
+      p = p - lead(p); // for the next cycle
+    }      
+    if ( p!= 0)
+    {
+      num = content(p);
+      p   = p/num;
+      NM[cl] = denominator(num);
+      dbprint(p,"content denominator:");
+      dbprint(p, string(NM[cl]));
+      cNu++; Numero[cNu]= NM[cl];
+      cl++;
+      NM[cl] = numerator(num);
+      dbprint(p,"content numerator:");
+      dbprint(p, string(NM[cl]));
+      cDe++; Denomiro[cDe]= NM[cl];
+      cl++;
+    }
+    // it seems that the next elements will not have real influence
+    while( p != 0)
+    {
+      NM[cl] = leadcoef(p); // should be all integer, i.e. non-rational
+      dbprint(p,"coef:");
+      dbprint(p, string(NM[cl]));
+      cl++;
+      p = p - lead(p);
+    }
+  }
+  NM = simplify(NM,4); // delete identical
+  string newvars = parstr(basering);
+  def save = basering;
+  string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;";
+  execute(NewRing);
+  // get params as variables
+  // create a list of non-monomials
+  ideal @L;
+  ideal F;
+  ideal NM = imap(save,NM);
+  NM = simplify(NM,8); //delete multiples
+  poly p,q;
+  cl = 1;
+  int j, cf;
+  for(i=1; i<=size(NM);i++)
+  {
+    p = NM[i] - lead(NM[i]);
+    if (p!=0)
+    {
+      //      L[cl] = p;
+      F  = factorize(NM[i],1); //non-constant factors only
+      cf = 1;
+      // factorize every polynomial
+      // throw away every monomial from factorization (also constants from above ring)
+      for (j=1; j<=size(F);j++)
+      {
+        q = F[j]-lead(F[j]);
+        if (q!=0)
+        {
+          @L[cl] = F[j];
+          cl++;
+        }
+      }
+    }
+  }
+  // return the result [in string-format]
+  @L = simplify(@L,2+4+8); // skip zeroes, doubled and entries, diff. by a constant
+  list SL;
+  for (j=1; j<=size(@L);j++)
+  {
+    SL[j] = string(@L[j]);
+  }
+  setring save;
+  return(SL);
+}
+example
+{  // TwoPendula
+  "EXAMPLE:"; echo = 2;
+  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
+  module RR = 
+    [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2], 
+    [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2],
+    [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2];
+  module R = transpose(RR);
+  module SR = std(R);
+  matrix T = lift(R,SR);
+  genericity(T);
+  //-- The result might be different when computing reduced bases:
+  matrix T2;
+  option(redSB);
+  option(redTail);
+  module SR2 = std(R);
+  T2 =  lift(R,SR2);
+  genericity(T2);
+}
+//---------------------------------------------------------------
 proc canonize(list L)
 "USAGE:  canonize(L), L a list
@@ -1089 +1151 @@
   view(CC);
 }
-
 //---------------------------------------------------------------
 static proc smdeg(matrix N)  
@@ -1436 +1497 @@
   print(Tr); // generators of the torsion submodule
 }
+//---------------------------------------------------------------
+proc evas_genericity(matrix M)
+{
+ideal I=ideal(M);
+I=simplify(I,2+4);
+int s = size(I);
+ideal Den;
+poly p;
+int i;
+for (i=1; i<=s; i++)
+ {
+   p=I[i];
+   while (p !=0) 
+   { Den=Den,denominator(leadcoef(p)); 
+     p=p-lead(p); 
+   }
+ }
+Den=simplify(Den,2+4);  
+string newvars = parstr(basering);
+def save = basering;
+string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;";
+execute(NewRing);
+ideal F;
+ideal Den = imap(save,Den);
+Den=simplify(Den,1);
+int s1=size(Den);
+for (i=1;i<=s1; i++)
+{ 
+  if (Den[i] !=1)
+  { 
+    F=F,factorize(Den[i],1); 
+  } 
+ }
+ F=simplify(F,2+4+8); 
+ // print(F);
+ ideal @L = F;
+ list SL;
+ int j;
+ for (j=1; j<=size(@L);j++)
+ {
+   SL[j] = string(@L[j]);
+ }
+ setring save;
+ return(SL);
+}