Changeset 0c824aa in git


Ignore:
Timestamp:
Mar 17, 2005, 7:42:48 PM (18 years ago)
Author:
Viktor Levandovskyy <levandov@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0604212ebb110535022efecad887940825b97c3f')
Children:
ce6c4df83cce5d4f8fcff0bdd37452a4649eef33
Parents:
b687a30f0e0e5c6b27bb1b68424088f6ea28d02d
Message:
*levandov: noncommutative things are thrown out, genericity and smith updated, Our_extgcd added for previous versions, minor fixes to docu


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/control.lib

    rb687a3 r0c824aa  
    1 version="$Id: control.lib,v 1.25 2005-03-16 22:16:53 levandov Exp $";
     1version="$Id: control.lib,v 1.26 2005-03-17 18:42:48 levandov Exp $";
    22category="Applications";
    33info="
    … …  
    1313
    1414PROCEDURES:
    15 control  (R);           analysis of controllability-related properties of R(using Ext modules),
    16 control2(R);            analysis of controllability-related properties of R(using dimension),
    17 autonom  (R);           analysis of autonomy-related properties of R (using Ext modules),
    18 autonom2 (R);           analysis of autonomy-related properties of R (using dimension),
    19 LeftKernel  (R);        a left kernel of R,
    20 RightKernel (R);        a right kernel of R,
    21 LeftInverse (R);        a left inverse of R,
     15control (R);           analysis of controllability-related properties of R (using Ext modules),
     16control2 (R);          analysis of controllability-related properties of R (using dimension),
     17autonom (R);           analysis of autonomy-related properties of R (using Ext modules),
     18autonom2 (R);          analysis of autonomy-related properties of R (using dimension),
     19LeftKernel (R);        a left kernel of R,
     20RightKernel (R);       a right kernel of R,
     21LeftInverse (R);       a left inverse of R,
    2222RightInverse (R);      a right inverse of R,
    23 smith (M);              a Smith form of a module M,
    24 genericity (M);         analysis of the genericity of parameters,
    25 canonize (L);           Groebnerification for modules in the output of control/autonomy procs.
     23smith (M);             a Smith form of a module M,
     24rank (M);              a column rank of M as of matrix,
     25genericity (M);        analysis of the genericity of parameters,
     26canonize (L);          Groebnerification for modules in the output of control/autonomy procs,
     27FindTorsion (R, I);    generators of the submodule of a module R, annihilated by the ideal I.
     28
    2629
    2730AUXILIARY PROCEDURES:
    … …  
    6265LIB "primdec.lib";
    6366LIB "matrix.lib";
    64 // -- noncommutative procs:
    65 LIB "ncalg.lib";
    66 LIB "nctools.lib";
    67 LIB "nchomolog.lib";
    68 LIB "gkdim.lib";
    6967
    7068//---------------------------------------------------------------
    … …  
    8179  return (v);
    8280}
     81//---------------------------------------------------------------
    8382proc declare(string NameOfRing, string Variables, list #)
    8483"USAGE: declare(NameOfRing, Variables,[Parameters, Ordering]);
    … …  
    252251PURPOSE: computes the right kernel of matrix M (a module of all elements v such that Mv=0)
    253252RETURN:  module
    254 NOTE:    in the noncommutative case RightKernel is a right module
    255253EXAMPLE: example RightKernel; shows an example
    256254"
    … …  
    297295"
    298296{
    299   // now it works also for the NC case;
    300   // two approaches are possible;
    301   // the older one is commented out
     297  // it works also for the NC case;
    302298  int NCols = ncols(M);
    303299  module Id = freemodule(NCols);
    … …  
    307303  matrix T;
    308304  // check the correctness (Id \subseteq M)
    309   // via dimension: {dim/GKdim} (M) = -1!
     305  // via dimension: dim (M) = -1!
    310306  int d = dim_Our(N);
    311307  if (d != -1)
    312308  {
    313     //    "No left inverse exists";
     309    // No left inverse exists
    314310    return(matrix(0));
    315311  }
    316   //  module N = liftstd(M, T);
    317   //  ideal TT = ideal(NF(N,Id));
    318   //  TT = simplify(TT,2);
    319   //  if (TT[1]!=0)
    320   //  {
    321   //    "No left inverse exists";
    322   //    return(matrix(0));
    323   //  }
    324312  matrix T2 = lift(N, Id);
    325   //  T2 = transpose(T2)*transpose(T);
    326313  T2 =  transpose(T2);
    327   // set options back
    328   option(set,old_opt);
     314  option(set,old_opt);  // set the options back
    329315  return(T2);
    330316}
    … …  
    359345}
    360346example
    361 { // trivial example:
    362 "EXAMPLE:";echo =2;
     347{ "EXAMPLE:";echo =2;
     348  // trivial example:
    363349  ring r = 0,(x,z),dp;
    364350  matrix M[1][2] = 1,x2+z;
    … …  
    384370    R = std(R);
    385371  }
    386   //  if (Is_NC())
    387   //  {
    388   //    d = GKdim(R);
    389   //  }
    390   //  else { d = dim(R); }
    391372  d = dim(R);
    392373  return(d);
    … …  
    395376static proc Ann_Our(module R)
    396377{
    397 //   if (Is_NC())
    398 //   {
    399 //     // TODO: use preAnn, that is an annihilator of the generators
    400 //     // what is a left ideal!
    401 //     return(-1);
    402 //   }
    403 //   else
    404 //   {
    405 //     return(Ann(R));
    406 //   }
    407378  return(Ann(R));
    408379}
    409380//-----------------------------------------------------------------------
    410381static proc Ext_Our(int i, module R, list #)
    411 // mimicking 'Ext_R' from "homolog.lib" in commutative case
    412 {
    413   //  int ISNC = Is_NC();
     382{
     383  // mimicking 'Ext_R' from homolog.lib
    414384  int ExtraArg = ( size(#)>0 );
    415   //  if (ISNC)
    416   //  {
    417   //    // ExtraArg not yet implemented and thus ignored
    418   //    return( ncExt_R(i,R) );
    419   //  }
    420   // the commutative case:
    421385  if (ExtraArg)
    422386  {
    … …  
    430394//------------------------------------------------------------------------
    431395static proc is_zero_Our
    432 //just a copy of 'is_zero' from "poly.lib" patched with GKdim
    433 {
     396{
     397  //just a copy of 'is_zero' from "poly.lib" patched with GKdim
    434398  int d = dim_Our(std(#[1]));
    435399  int a = ( d==-1 );
    … …  
    452416{
    453417  // TODO: NVars to be replaced with the global hom. dimension of basering!!!
     418  // Is not clear what to do with gl.dim of qrings
    454419  string DofS = "dimension of the system:";
    455420  string Fn   = "number of first nonzero Ext:";
    456   string Gen_mes = "The following parameter constellations might lead to a non-controllable system: ";
     421  string Gen_mes = "Parameter constellations which might lead to a non-controllable system:";
    457422
    458423  module RK = RightKernel(R);
    … …  
    539504  module R_std=std(R);
    540505  module Ext_1 = std(Ext_Our(1,R_std));
    541   //  module Ext_1 = std(Ext_Our(1,R));
    542506   
    543507  ExtIsZero=is_zero_Our(Ext_1);
    … …  
    547511    i++;
    548512    ExtIsZero = is_zero_Our( Ext_Our(i,R_std) );
    549     //    ExtIsZero = is_zero_Our( Ext_Our(i,R) );
    550513  };
    551514  matrix T=lift(R,R_std);
    … …  
    554517
    555518  return( control_output( i, NVars, R, Ext_1, l ) );
    556   //  return( control_output( i, NVars, R, Ext_1 ) );
    557519}
    558520example
    … …  
    568530  view(control(R));
    569531};
    570 //------------------------------------------------------------------------------------------
     532//--------------------------------------------------------------------------
    571533proc control2(module R)
    572534"USAGE:  control2(R); R a module (R is the matrix of the system of equations to be investigated)
    … …  
    574536RETURN: list
    575537EXAMPLE:  example control2; shows an example
    576 NOTE: same as control(R); but using dimensions
     538NOTE: same as control(R); but using dimensions, this procedure only works for full row rank matrices
    577539"
    578540{
    … …  
    623585example
    624586{"EXAMPLE: ";echo = 2;
    625   //deRham complex
     587  // de Rham complex
    626588  ring r=0,(D(1..3)),dp;
    627589  module R;
    … …  
    635597//------------------------------------------------------------------------
    636598static proc autonom_output( int i, int NVars, module RC, int R_rank )
    637 //static proc autonom_output( int i, int NVars, module RC )
    638 //static proc autonom_output( int i, int NVars )
    639599"USAGE:  proc autonom_output(i, NVars, RC, R_rank)
    640600           i:  integer, number of first nonzero Ext or
    … …  
    724684  int d;
    725685  int NVars = nvars(basering);
    726 //   R = transpose(R);
    727 //   d = dim_Our( std(R) );
    728 //   return( autonom_output(NVars-d,NVars) );
    729686  module RT = transpose(R);
    730687  module RC;  //for computation of controllable part if if exists
     688  int R_rank = ncols(R);
    731689  d     = dim_Our( std(RT) );  //this is the dimension of the system
    732690  int i = NVars-d;  //First non-zero Ext
    … …  
    734692    {
    735693      RC=LeftKernel(RightKernel(R));
    736     }
    737   return( autonom_output(i,NVars,RC) );
     694      R_rank=rank(R);
     695    }
     696  return( autonom_output(i,NVars,RC,R_rank) );
    738697}
    739698example
    … …  
    763722  int R_rank=ncols(R);
    764723  ExtIsZero=is_zero_Our(Ext_Our(0,RT));     
    765 //   R=transpose(R);
    766 //   ExtIsZero=is_zero_Our(Ext_Our(0,R));
    767724  int i=0;
    768725  while( (ExtIsZero)&&(i<=NVars) )
    … …  
    770727    i++;
    771728    ExtIsZero = is_zero_Our(Ext_Our(i,RT));
    772     //    ExtIsZero = is_zero_Our(Ext_Our(i,R));
    773729  };
    774   if(i==0)
    775     {
    776       RC=LeftKernel(RightKernel(R));
    777       R_rank=rank(R);
    778     }
     730  if (i==0)
     731  {
     732    RC=LeftKernel(RightKernel(R));
     733    R_rank=rank(R);
     734  }
    779735  return(autonom_output(i,NVars,RC,R_rank));
    780   //  return(autonom_output(i,NVars,RC));
    781   //  return(autonom_output(i,NVars));     
    782736}
    783737example
    784738{"EXAMPLE:"; echo = 2;
    785   //Cauchy
     739  // Cauchy
    786740  ring r=0,(s1,s2,s3,s4),dp;
    787741  module R= [s1,-s2],
    … …  
    834788//----------------------------------------------------------
    835789//control
    836 //
     790//----------------------------------------------------------
    837791proc ex1()
    838792{
    … …  
    849803  ring @r=0,(s1,s2,s3),dp;
    850804  module R=[0,-s3,s2],
    851             [s3,0,-s1],
    852             [-s2,s1,0];
     805           [s3,0,-s1],
     806           [-s2,s1,0];
    853807  R=transpose(R);         
    854808  export R;
    … …  
    894848};
    895849
    896 
    897850//----------------------------------------------------------
    898 
    899851proc exFlexibleRod()
    900852{
    … …  
    934886//----------------------------------------------------------
    935887proc genericity(matrix M)
    936 "USAGE:  genericity(M), M is a module/matrix
     888"USAGE:  genericity(M), M is a matrix/module
    937889PURPOSE: determine parametric expressions which have been assumed to be non-zero in the process of computing the Groebner basis
    938890RETURN:  list (of strings)
    … …  
    947899    return("-");
    948900  }
    949   list RT = evas_genericity(M);
     901  list RT = evas_genericity(M); // list of strings
     902  if ((size(RT)==1) && (RT[1] == ""))
     903  {
     904    return("-");
     905  }
    950906  return(RT);
    951907}
    … …  
    971927}
    972928//---------------------------------------------------------------
    973 proc victors_genericity(matrix M)
    974 "USAGE:  genericity(M), M is a module/matrix
    975 PURPOSE: determine parametric expressions which have been assumed to be non-zero in the process of computing the Groebner basis
    976 RETURN:  list (of strings)
    977 NOTE: we strongly recommend to switch on the redSB and redTail options;
    978 @*    the procedure is effective with the lift procedure for modules with parameters
    979 EXAMPLE:  example genericity; shows an example
    980 "
     929static proc victors_genericity(matrix M)
    981930{
    982931  // returns "-", if there are no parameters!
    … …  
    10901039  return(SL);
    10911040}
    1092 example
    1093 {  // TwoPendula
    1094   "EXAMPLE:"; echo = 2;
    1095   ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
    1096   module RR =
    1097     [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2],
    1098     [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2],
    1099     [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2];
    1100   module R = transpose(RR);
    1101   module SR = std(R);
    1102   matrix T = lift(R,SR);
    1103   genericity(T);
    1104   //-- The result might be different when computing reduced bases:
    1105   matrix T2;
    1106   option(redSB);
    1107   option(redTail);
    1108   module SR2 = std(R);
    1109   T2 =  lift(R,SR2);
    1110   genericity(T2);
    1111 }
     1041//---------------------------------------------------------------
     1042static proc evas_genericity(matrix M)
     1043{
     1044  // called from the main genericity proc
     1045  ideal I = ideal(M);
     1046  I = simplify(I,2+4);
     1047  int s = size(I);
     1048  ideal Den;
     1049  poly p;
     1050  int i;
     1051  for (i=1; i<=s; i++)
     1052  {
     1053    p = I[i];
     1054    while (p !=0)
     1055    {
     1056      Den = Den, denominator(leadcoef(p));
     1057      p   = p-lead(p);
     1058    }
     1059  }
     1060  Den = simplify(Den,2+4); 
     1061  string newvars = parstr(basering);
     1062  def save = basering;
     1063  string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;";
     1064  execute(NewRing);
     1065  ideal F;
     1066  ideal Den = imap(save,Den);
     1067  Den = simplify(Den,2);
     1068  int s1 = size(Den);
     1069  for (i=1; i<=s1; i++)
     1070  {
     1071    if (Den[i] !=1)
     1072    {
     1073      F= F, factorize(Den[i],1);
     1074    }
     1075  }
     1076  F = simplify(F, 2+4+8);
     1077  ideal @L = F;
     1078  list SL;
     1079  int c,j;
     1080  string Mono;
     1081  c = 1;
     1082  for (j=1; j<=size(@L);j++)
     1083  {
     1084    if (leadcoef(@L[j]) <0)
     1085    {
     1086      @L[j] = -1*@L[j];
     1087    }
     1088    if ( (@L[j] - lead(@L[j]))==0 ) //@L[j] is a monomial
     1089    {
     1090      Mono = Mono + string(@L[j])+ ","; // concatenation
     1091    }
     1092    else
     1093    {
     1094      c++;
     1095      SL[c] = string(@L[j]);
     1096    }
     1097  }
     1098  if (Mono!="")
     1099  {
     1100    Mono = Mono[1..size(Mono)-1]; // delete the last semicolon
     1101  }
     1102  SL[1] = Mono;
     1103  setring save;
     1104  return(SL);
     1105}
     1106
    11121107//---------------------------------------------------------------
    11131108proc canonize(list L)
    … …  
    11211116  list M = L;
    11221117  intvec v=Opt_Our();
    1123   //  option(redSB);
    1124   //  option(redTail);
    11251118  int s = size(L);
    11261119  int i;
    … …  
    11351128  }
    11361129  option(set, v); //set old values back
    1137   //  option(noredTail); // set the initial option back
    11381130  return(M);
    11391131}
    … …  
    11611153  int i,j;            // counter
    11621154
    1163   if(N==0)
    1164   {
    1165     v=1,1;
     1155  if (N==0)
     1156  {
     1157    v = 1,1;
    11661158    return(v);
    11671159  }
    11681160 
    1169   for(i=1;i<=n;i++) // hier wird ein Element ausgewaehlt(!=0) und mit dessen Grad gestartet
    1170   {
    1171     for(j=1;j<=m;j++)
     1161  for (i=1; i<=n; i++)
     1162// hier wird ein Element ausgewaehlt(!=0) und mit dessen Grad gestartet
     1163  {
     1164    for (j=1; j<=m; j++)
    11721165    {
    11731166      if( deg(N[i,j])!=-1 )
    11741167      {
    1175     d=deg(N[i,j]);
    1176     break;
     1168        d=deg(N[i,j]);
     1169        break;
    11771170      }
    11781171    }
    1179     if(d!=-1)
     1172    if (d != -1)
    11801173    {
    11811174      break;
    11821175    } 
    11831176  }
    1184   for(i =1; i<=n; i++)
    1185   {
    1186     for(j=1;j<=m;j++)
     1177  for(i=1; i<=n; i++)
     1178  {
     1179    for(j=1; j<=m; j++)
    11871180    {
    11881181      d_temp = deg(N[i,j]);
    1189       if ((d_temp < d) && (N[i,j]!=0) )
     1182      if ( (d_temp < d) && (N[i,j]!=0) )
    11901183      {
    1191     d=d_temp;
     1184        d=d_temp;
    11921185      }
    11931186    }
    … …  
    11951188  for (i=1; i<=n; i++)
    11961189  {
    1197     for (j =1; j<=m;j++)
     1190    for (j=1; j<=m;j++)
    11981191    {
    11991192      if ( (deg(N[i,j]) == d) && (N[i,j]!=0) )
    12001193      {
    1201     v = i,j;
    1202     return(v);
     1194        v = i,j;
     1195        return(v);
    12031196      }
    12041197    }
    … …  
    12111204  int i,j;
    12121205  int n = nrows(v);
    1213   for(j=1; j<=n;j++)
    1214   {
    1215     if(v[j]!=0)
     1206  for( j=1; j<=n; j++)
     1207  {
     1208    if (v[j] != 0)
    12161209    {
    12171210      i++;
    … …  
    12231216  }
    12241217  return(i);
     1218}
     1219//---------------------------------------------------------------
     1220static proc Our_extgcd(poly p, poly q)
     1221{
     1222  ideal J;   //for extgcd-computations
     1223  matrix T; //----------"------------
     1224  list L;
     1225  // the extgcd-command has a bug in versions before 2-0-7
     1226  if ( system("version")<=2006 )
     1227  {
     1228    J = p,q; // J = N[k-1,k-1],N[k,k]; //J is of type ideal
     1229    L[1] = liftstd(J,T);  //T is of type matrix
     1230    if(J[1]==p) //this is just for the case the SINGULAR swaps the
     1231    //      two elements due to ordering
     1232    { 
     1233      L[2] = T[1,1];
     1234      L[3] = T[2,1];
     1235    }
     1236    else
     1237    {
     1238      L[2] = T[2,1];
     1239      L[3] = T[1,1];
     1240    }
     1241  }
     1242  else
     1243  {
     1244    L=extgcd(p,q);
     1245    //    L=extgcd(N[k-1,k-1],N[k,k]); 
     1246    //one can use this line if extgcd-bug is fixed
     1247  }
     1248  return(L);
    12251249}
    12261250//---------------------------------------------------------------
    … …  
    12381262{
    12391263  if (nvars(basering)>1) //if more than one variable, return empty list
    1240     {
    1241       list @L;
    1242       return (@L);
    1243     }
     1264  {
     1265    list @L;
     1266    return (@L);
     1267  }
    12441268  matrix N = matrix(M);         //Typecasting
    12451269  int n = nrows(N);
    … …  
    12641288      while(k<=m )        //inner while-loop for row-operations
    12651289      {
    1266     if( (n>m) && (k < n) && (k<m))
    1267     {     
    1268       if( simplify((ideal(submat(N,k+1..n,k+1..m))),2)== 0)
    1269       {
    1270         return(N,k-1,P,Q);
    1271       }
    1272     }
     1290        if( (n>m) && (k < n) && (k<m))
     1291        {     
     1292          if( simplify((ideal(submat(N,k+1..n,k+1..m))),2)== 0)
     1293          {
     1294            return(N,k-1,P,Q);
     1295          }
     1296        }
    12731297        i,j = smdeg(submat(N,k..n,k..m)); //choose smallest degree in the remaining submatrix
    12741298        i=i+(k-1);                        //indices adjusted to the whole matrix
    … …  
    12981322        {
    12991323          N = addrow(N,k,-lambda[1,ii]*var(1)^f[1,ii],ii);
    1300       N = normalize(N);
     1324          N = normalize(N);
    13011325          P = addrow(P,k,-lambda[1,ii]*var(1)^f[1,ii],ii);
    1302     }
     1326        }
    13031327      }
    13041328      if (k>n)
    13051329      {
    1306     break;
     1330        break;
    13071331      }
    13081332      if(NoNon0Pol(transpose(N)[k])==1)
    … …  
    13271351    if( (k!=1) && (k<n) && (k<m) )
    13281352    {
    1329       L=extgcd(N[k-1,k-1],N[k,k]);  //one can use this line if extgcd-bug is fixed
    1330      
    1331       if(!(N[k-1,k-1]==L[1]))  //means that N[k-1,k-1] is not a divisor of N[k,k]
    1332         {
    1333           N=addrow(N,k-1,L[2],k);
    1334           normalize(N);
    1335           P=addrow(P,k-1,L[2],k);
    1336          
    1337           N=addcol(N,k,-L[3],k-1);
    1338           N=normalize(N);
    1339           Q=addcol(Q,k,-L[3],k-1);
    1340           k=k-2;
    1341         }
     1353      L = Our_extgcd(N[k-1,k-1],N[k,k]);     
     1354      if ( N[k-1,k-1]!=L[1] )  //means that N[k-1,k-1] is not a divisor of N[k,k]
     1355      {
     1356        N=addrow(N,k-1,L[2],k);
     1357        normalize(N);
     1358        P=addrow(P,k-1,L[2],k);
     1359       
     1360        N=addcol(N,k,-L[3],k-1);
     1361        N=normalize(N);
     1362        Q=addcol(Q,k,-L[3],k-1);
     1363        k=k-2;
     1364      }
    13421365    }
    13431366  } 
    13441367  if( (k<=n) && (k<=m) )
    1345     {
    1346       if( N[k,k]==0)
    1347       {
    1348         return(N,k-1,P,Q);
    1349       }
    1350     }
     1368  {
     1369    if( N[k,k]==0)
     1370    {
     1371      return(N,k-1,P,Q);
     1372    }
     1373  }
    13511374  return(N,k,P,Q);
    13521375}
    … …  
    13571380  option(redTail);
    13581381  ring r=0,x,dp;
    1359   // This example is trivial, but we want to see,
    1360   // that nothing is changed, when the matrix is already in Smith-Form!
     1382  // see what happens when the matrix is already in Smith-Form
    13611383  module M = [x,0,0],[0,x2,0],[0,0,x3];
    13621384  print(M);
    … …  
    14751497  }
    14761498  int nc  = ncols(S);
    1477   module Tmp = AS*freemodule(nc);
    1478   module To = modulo(Tmp,S);
    1479   To = NF(std(To+S),S);
    1480   To = simplify(To,2);
    1481   To = std(To);
     1499  module To = quotient(S,AS);
     1500  To = std(NF(To,S));
    14821501  return(To);
    14831502}
    … …  
    14901509  module RR = transpose(R);
    14911510  list L = control(RR);
    1492   view(L);
    14931511  // here, we have the annihilator:
    14941512  ideal LAnn = D1; // = L[10]
    … …  
    14971515  print(Tr); // generators of the torsion submodule
    14981516}
    1499 //---------------------------------------------------------------
    1500 proc evas_genericity(matrix M)
    1501 {
    1502 ideal I=ideal(M);
    1503 I=simplify(I,2+4);
    1504 int s = size(I);
    1505 ideal Den;
    1506 poly p;
    1507 int i;
    1508 for (i=1; i<=s; i++)
    1509  {
    1510    p=I[i];
    1511    while (p !=0)
    1512    { Den=Den,denominator(leadcoef(p));
    1513      p=p-lead(p);
    1514    }
    1515  }
    1516 Den=simplify(Den,2+4); 
    1517 string newvars = parstr(basering);
    1518 def save = basering;
    1519 string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;";
    1520 execute(NewRing);
    1521 ideal F;
    1522 ideal Den = imap(save,Den);
    1523 Den=simplify(Den,1);
    1524 int s1=size(Den);
    1525 for (i=1;i<=s1; i++)
    1526 {
    1527   if (Den[i] !=1)
    1528   {
    1529     F=F,factorize(Den[i],1);
    1530   }
    1531  }
    1532  F=simplify(F,2+4+8);
    1533  // print(F);
    1534  ideal @L = F;
    1535  list SL;
    1536  int j;
    1537  for (j=1; j<=size(@L);j++)
    1538  {
    1539    SL[j] = string(@L[j]);
    1540  }
    1541  setring save;
    1542  return(SL);
    1543 }
Note: See TracChangeset for help on using the changeset viewer.