Changeset 1e1ec4 in git


Ignore:
Timestamp:
Jan 4, 2013, 5:54:18 PM (10 years ago)
Author:
Oleksandr Motsak <motsak@…>
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
42ea852aa2e1e683808b1ac3305dda96677af761
Parents:
8f296a6216092a84f1ebb509dbcda5fe428004f7
git-author:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2013-01-04 17:54:18+01:00
git-committer:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2013-01-15 20:41:56+01:00
Message:
Updated LIBs according to master

add: new LIBs from master
fix: updated LIBs due to minpoly/(de)numerator changes
fix: -> $Id$
fix: Fixing wrong rebase of SW on master (LIBs)
Location:
Singular/LIB
Files:
11 added
59 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/COPYING

    r8f296a r1e1ec4  
    88           Authors: see the header of the libraries or the list below
    99
    10                         Copyright (C) 1986-2011
     10                        Copyright (C) 1986-2012
    1111
    1212This directory contains all Singular libraries,
    … …  
    139139                Gert-Martin Greuel              greuel@mathematik.uni-kl.de
    140140                Anne Fruehbis-Krueger           anne@mathematik.uni-kl.de
    141 polymake.lib    Thomas Markwig                  keilen@mathematik.uni-kl.de
     141oldpolymake.lib Thomas Markwig                  keilen@mathematik.uni-kl.de
    142142presolve.lib    Gert-Martin Greuel              greuel@mathematik.uni-kl.de
    143143primdec.lib     Gerhard Pfister                 pfister@mathematik.uni-kl.de

  • Singular/LIB/aksaka.lib

    r8f296a r1e1ec4  
    325325  r=2;
    326326
    327   if(gcdN(n,r)!=1)   //falls ggt ungleich eins, Teiler von n gefunden,
     327  if(gcd(n,r)!=1)   //falls ggt ungleich eins, Teiler von n gefunden,
    328328                    // Schritt 3 des ASK-Alg.
    329329  {
    … …  
    358358      }
    359359      r=r+1;
    360       if(gcdN(n,r)!=1)   //falls ggt ungleich eins, Teiler von n gefunden,
     360      if(gcd(n,r)!=1)   //falls ggt ungleich eins, Teiler von n gefunden,
    361361                         //wird aufgrund von Schritt 3 des ASK-Alg. fuer
    362362                         //jeden Kandidaten r getestet

  • Singular/LIB/alexpoly.lib

    r8f296a r1e1ec4  
    17461746  int initial_tm=0;                  // total multipl. of the endpoint of Puiseux chain P_i-1
    17471747  int g=size(charexp);
    1748   list es=divsequence(charexp[2],charexp[1]);   // keeps the lenghts of the Puiseuxchainparts s_i,j
     1748  list es=divsequence(charexp[2],charexp[1]);   // keeps the lengths of the Puiseuxchainparts s_i,j
    17491749  intvec divseq=es[1];
    17501750  int r=es[2];

  • Singular/LIB/algebra.lib

    r8f296a r1e1ec4  
    33//new proc  nonZeroEntry(id), used to fix a bug in proc finitenessTest
    44///////////////////////////////////////////////////////////////////////////////
    5 version="$Id: algebra.lib 14661 2012-03-05 11:00:14Z motsak $";
     5version="$Id$";
    66category="Commutative Algebra";
    77info="
    … …  
    2323 mapIsFinite(R,phi,I);  query for finiteness of map phi:R --> basering/I
    2424
    25 AUXILIARY PROCEDURES:
    2625 finitenessTest(i,z);   find variables which occur as pure power in lead(i)
    2726 nonZeroEntry(id);      list describing non-zero entries of an identifier
    … …  
    816815   //a procedure from ring.lib changing the order to dp creating a new
    817816   //basering @R in order to compute the dimension d of i
    818    def @R=changeord("dp");
     817   def @R=changeord(list(list("dp",1:nvars(basering))));
    819818   setring @R;
    820819   ideal i = imap(r,i);
    … …  
    830829   // ------------------------ change of ordering ---------------------------
    831830   //Now change the order to (dp(n-d),lp) creating a new basering @S
    832    string s ="dp("+string(n-d)+"),lp";
    833    def @S=changeord(s);
     831   def @S=changeord(list(list("dp",1:(n-d)),list("lp",1:d)));
    834832   setring @S;
    835833   ideal m;

  • Singular/LIB/arcpoint.lib

    r8f296a r1e1ec4  
    128128
    129129
    130    // consider all sequences of lenght <step-1> giving rise to a
     130   // consider all sequences of length <step-1> giving rise to a
    131131   // family...
    132132

  • Singular/LIB/assprimeszerodim.lib

    r8f296a r1e1ec4  
    243243         int alg = 1;
    244244         if(n1 >= 10)
    245          { 
     245         {
    246246            int n2 = n1 + 1;
    247247            int n3 = n1;
    … …  
    258258         int alg = #[3];
    259259         if(n1 >= 10)
    260          { 
     260         {
    261261            int n2 = n1 + 1;
    262262            int n3 = n1;
    … …  
    276276      int n3 = 10;
    277277   }
    278    
     278
    279279   if(printlevel >= 10)
    280280   {
    … …  
    292292   ideal F;
    293293   poly F1;
    294    
     294
    295295   if(printlevel >= 10) { "========== Start modStd =========="; }
    296296   I = modStd(I,n1);
    … …  
    415415   int tt = timer;
    416416   int rt = rtimer;
    417    
     417
    418418   while(1)
    419419   {
    420420      tt = timer;
    421421      rt = rtimer;
    422      
     422
    423423      if(printlevel >= 9) { "size(L) = "+string(size(L)); }
    424      
     424
    425425      if(n1 > 1)
    426426      {
    … …  
    471471         }
    472472      }
    473      
     473
    474474      if(printlevel >= 9)
    475475      {
    … …  
    490490                             +" seconds"; }
    491491
    492       if(pTestPoly(F[1], ringL, alg, L))                             
     492      if(pTestPoly(F[1], ringL, alg, L))
    493493      {
    494494         F = cleardenom(F[1]);
    … …  
    513513               phi = rHelp,var(nvars(SPR));
    514514               H = phi(H);
    515                
     515
    516516               if(printlevel >= 9)
    517517               {
    … …  
    519519                  "CPU-time without test is "+string(timer - T)+" seconds.";
    520520               }
    521                
     521
    522522               T = timer;
    523523               RT = rtimer;
    … …  
    647647   setring R;
    648648   map phi = Rhelp,p;
    649    def Rlp = changeord("dp("+string(nvars(R)-1)+"),dp(1)");
     649   def Rlp = changeord(list(list("dp",1:(nvars(R)-1)),list("dp",1:1)));
    650650   setring Rlp;
    651651   poly p = imap(R,p);
    … …  
    728728//=== vector space dim(basering/J)=deg(F),
    729729//=== F a poly in Q[T] such that <F>=kernel(Q[T]--->basering) mapping T to r
    730 //=== if not found returns a generic (randomly choosen) r
     730//=== if not found returns a generic (randomly chosen) r
    731731   int d = vdim(J);
    732732   def R = basering;
    … …  
    919919   int i,j,p;
    920920   def R0 = basering;
    921    
     921
    922922   while(!i)
    923923   {
    … …  
    929929      }
    930930   }
    931    
     931
    932932   ringL[1] = p;
    933933   def @R = ring(ringL);
    … …  
    942942   poly testF = fetch(R0,testF);
    943943   int k = (testF == F);
    944    
     944
    945945   setring R0;
    946946   return(k);

  • Singular/LIB/atkins.lib

    r8f296a r1e1ec4  
    117117    {
    118118      D=-4*a+B;
    119       if((D<0)&&((D mod 4)!=2)&&((D mod 4)!=3)&&(absValue(D)<4*N)&&(newTest(L,D)==1))
     119      if((D<0)&&((D mod 4)!=-2)&&((D mod 4)!=-1)&&(absValue(D)<4*N)&&(newTest(L,D)==1))
    120120      {
    121121        L[size(L)+1]=D;
    … …  
    162162       if((p/2>=x(0))||(p<=x(0)))
    163163       {
    164            x(0)=-x(0) mod p;
     164           x(0)=-x(0) mod p+p;
    165165       }
    166166       a=p;
    … …  
    203203"
    204204{
    205   if((D>=0)||((D mod 4)==2)||((D mod 4)==3)||(absValue(D)>=4*p))
     205  if((D>=0)||((D mod 4)==-2)||((D mod 4)==-1)||(absValue(D)>=4*p))
    206206  {                                  // (0)[Test if assumptions well-defined]
    207207    return(0);
    … …  
    234234      {
    235235        x(0)=squareRoot(D,p);        // (3)[Compute square root]
    236         if((x(0) mod 2)!=(D mod 2))
     236        if((x(0) mod 2)!=(-(D mod 2))) // D is <0
    237237        {
    238238          x(0)=p-x(0);
    … …  
    248248        }
    249249        c=(4*p-b^2)/absValue(D);     // (5)[Test solution]
    250         if((((4*p-b^2) mod absValue(D))!=0)||(c!=intRoot(c)^2))
     250        number root_c=intRoot(c);
     251        if((((4*p-b^2) mod absValue(D))!=0)||(c!=root_c^2))
    251252        {
    252253          return(-1);
    … …  
    255256        else
    256257        {
    257           list L=b,intRoot(c);
     258          list L=b,root_c;
    258259          return(L);
    259260        }
    … …  
    294295  }
    295296  number a=random(2,2147483629);
    296   number d=gcdN(powerN(a,k,n)-1,n);
     297  number d=gcd(powerN(a,k,n)-1,n);
    297298  if((d>1)&&(d<n)){return(d);}
    298299  return(n);
    … …  
    797798  if(N==1)           {return(-1);} // (0)[Test if assumptions well-defined]
    798799  if((N==2)||(N==3)) {return(1);}
    799   if(gcdN(N,6)!=1)
    800   {
    801     if(printlevel>=1) {"gcd(N,6) = "+string(gcdN(N,6));pause();"";}
     800  if(gcd(N,6)!=1)
     801  {
     802    if(printlevel>=1) {"gcd(N,6) = "+string(gcd(N,6));pause();"";}
    802803    return(-1);
    803804  }

  • Singular/LIB/binresol.lib

    r8f296a r1e1ec4  
    1818 inidata(K,k);                      verifies input data, a binomial ideal K of k generators
    1919 identifyvar();                     identifies status of variables
    20  data(K,k,n);                       transforms data on lists of lenght n
     20 data(K,k,n);                       transforms data on lists of length n
    2121 Edatalist(Coef,Exp,k,n,flag);      gives the E-order of each term in Exp
    2222 EOrdlist(Coef,Exp,k,n,flag);       computes the E-order of an ideal (giving in the language of lists)

  • Singular/LIB/classify.lib

    r8f296a r1e1ec4  
    27442744    s2 = Typ + s3 +"]";
    27452745  }  // es kommt mindestens ein komma vor...
    2746   //----------------------- more than 1 paramater -----------------------------
     2746  //----------------------- more than 1 parameter -----------------------------
    27472747  else {
    27482748    b  = find(s1, ",");

  • Singular/LIB/crypto.lib

    r8f296a r1e1ec4  
    109109  number x=a mod n;
    110110  if(x==0){return(list(0,1,n))}
     111  if (x<0) { x=x+n;}
    111112  list l=exgcdN(n,x);
    112113  return(list(l[2],l[1]-(a-x)*l[2]/n,l[3]))
    … …  
    183184"
    184185{
    185    if(d==0){return(1)}
     186   if(d==0){return(1);}
    186187   int i;
    187188   if(n==0)
    … …  
    195196   for(i=12;i>=2;i--)
    196197   {
    197       if((d mod i)==0){return(powerN(m,d/i,n)^i mod n);}
     198      if((d mod i)==0)
     199      {
     200        number rr=powerN(m,d/i,n)^i mod n;
     201        if (rr<0) { rr=rr+n;}
     202        return(rr);
     203      }
    198204   }
    199205   return(m*powerN(m,d-1,n) mod n);
    … …  
    348354"
    349355{
    350    return((numerator(x)-(numerator(x) mod denominator(x)))/denominator(x));
     356  if (x>=0)
     357  {
     358    return((numerator(x)-(numerator(x) mod denominator(x)))/denominator(x));
     359  }
     360  else
     361  {
     362    return((numerator(x)-(numerator(x) mod denominator(x)+denominator(x)))/denominator(x));
     363  }
    351364}
    352365example
    … …  
    367380   while(((m>t)&&(m>s))||((m<t)&&(m<s)))
    368381   {
    369       x=intPart(x/2+m/(2*x));  //Newton step
    370       t=x^2;
     382      if (x==0) { t=0; }
     383      else
     384      {
     385        x=intPart(x/2+m/(2*x));  //Newton step
     386        t=x^2;
     387      }
    371388      if(t>m)
    372389      {
    … …  
    397414   if(p==2){return(a);}
    398415   if((a mod p)==0){return(0);}
     416   if (a<0) { a=a+p; }
    399417   if(powerN(a,p-1,p)!=1)
    400418   {
    … …  
    628646   number a=(t*u/d) mod (p-1);
    629647
     648   number pn=powerN(b,a,p);
     649   if (pn<0) { pn=pn+p;}
    630650   while(powerN(b,a,p)!=y)
    631651   {
    632652      a=(a+(p-1)/d) mod (p-1);
     653      if (a<0) { a=a+p-1; }
    633654   }
    634655   return(a);
    … …  
    713734   if((n mod 2)==0){return(0);}
    714735
    715    number a;
     736   number a,pn,jn;
    716737   int i;
    717738   while(i<k)
    … …  
    719740      i++;
    720741      a=random(2,2147483629) mod n; if(a==0){a=3;}
    721       if(gcdN(a,n)!=1){return(0);}
    722       if(powerN(a,(n-1)/2,n)!=(Jacobi(a,n) mod n)){return(0);}
     742      if(gcd(a,n)!=1){return(0);}
     743      pn=powerN(a,(n-1)/2,n);
     744      if (pn<0) { pn=pn+n;}
     745      jn=Jacobi(a,n) mod n;
     746      if (jn<0) { jn=jn+n;}
     747      if(pn!=jn){return(0);}
    723748   }
    724749   return(1);
    … …  
    804829         a=a+1;
    805830         if(powerN(a,N-1,N)!=1){return("not prime");}
    806          g=gcdN(powerN(a,(N-1)/PA[i],N),N);
     831         g=gcd(powerN(a,(N-1)/PA[i],N),N);
    807832         if(g==1)
    808833         {
    … …  
    884909          y=powerN(y,2,n); y=(y+a) mod n;
    885910          y=powerN(y,2,n); y=(y+a) mod n;
    886           d=gcdN(x-y,n);
     911          d=gcd(x-y,n);
    887912          if(d>1)
    888913          {
    … …  
    900925          }
    901926       }
    902 
    903927    }
    904928    if(allFactors)      //want to obtain all prime factors
    … …  
    957981   }
    958982   number a=random(2,2147483629);
    959    number d=gcdN(powerN(a,k,n)-1,n);
     983   number d=gcd(powerN(a,k,n)-1,n);
    960984   if((d>1)&&(d<n)){return(d);}
    961985   return(n);
    … …  
    10931117         y=y*B[l]^(v[l] div 2) mod n;
    10941118      }
    1095       d=gcdN(x-y,n);
     1119      d=gcd(x-y,n);
    10961120      if((d>1)&&(d<n)){return(d);}
    10971121   }
    … …  
    11411165   for(i=1;i<=3;i++)
    11421166   {
    1143       P[i]=P[i] mod N;
    1144       Q[i]=Q[i] mod N;
     1167      P[i]=P[i] mod N; if (P[i]<0) { P[i]=P[i]+N:}
     1168      Q[i]=Q[i] mod N; if (Q[i]<0) { Q[i]=Q[i]+N;}
    11451169   }
    11461170   list Resu;
    … …  
    11521176   //test for ellictic curve
    11531177   number D=4*a^3+27*b^2;
    1154    number g=gcdN(D,N);
     1178   number g=gcd(D,N);
    11551179   if(g==N){return(Error);}
    11561180   if(g!=1)
    … …  
    12101234   }
    12111235   Resu[1]=(L^2-P[1]-Q[1]) mod N;
     1236   if (Resu[1]<0) { Resu[1]=Resu[1]+N; }
    12121237   Resu[2]=(L*(P[1]-Resu[1])-P[2]) mod N;
     1238   if (Resu[2]<0) { Resu[2]=Resu[2]+N; }
    12131239   Resu[3]=number(1);
    12141240   return(Resu);
    … …  
    12961322     //test for ellictic curve
    12971323     number D=4*a^3+27*b^4; //the constant term is b^2
    1298      number g=gcdN(D,N);
     1324     number g=gcd(D,N);
    12991325     if(g<N){return(list(a,b,g));}
    13001326   }
    … …  
    14631489         while(j<m)
    14641490         {
    1465             j=j+1;
     1491            j++;
    14661492            if((Q[1]==B[j][1])&&(Q[2]==B[j][2])&&(Q[3]==B[j][3]))
    14671493            {
    … …  
    14741500                  if(K[size(K)][2]!=K[size(K)-1][2])
    14751501                  {
    1476                      d=gcdN(K[size(K)][2],K[size(K)-1][2]);
     1502                     d=gcd(K[size(K)][2],K[size(K)-1][2]);
    14771503                     if(ellipticMult(q,a,b,K[size(K)],d)[3]==0)
    14781504                     {
    … …  
    15301556//test for ellictic curve
    15311557   number D=4*a^3+27*b^2;
    1532    number G=gcdN(D,N);
     1558   number G=gcd(D,N);
    15331559   if(G==N){ERROR("not an elliptic curve");}
    15341560   if(G!=1){ERROR("not a prime");}
    … …  
    19391965   {
    19401966      P=ellipticRandomPoint(N,L[1],L[2]);  //a random point on C
     1967      "P=",P;
    19411968      if(ellipticMult(N,L[1],L[2],P,m)[3]!=0){"N is not prime";return(-5);}
    19421969      if(ellipticMult(N,L[1],L[2],P,u)[3]!=0)

  • Singular/LIB/curvepar.lib

    r8f296a r1e1ec4  
    1 ///////////////////////////////////////////////////////////////////////////////
     1/////////////////////////////////////////////////////////////////////////////////
    22version="$Id$";
    3 category="Singularities";
     3category="Singularity Theory";
    44info="
    5 LIBRARY:  space_curve.lib
    6 
    7 AUTHOR:   Viazovska Maryna, email: viazovsk@mathematik.uni-kl.de
     5LIBRARY:  curvepar.lib  Resolution of space curve singularities, semi-group
     6
     7AUTHOR:     Gerhard Pfister    email: pfister@mathematik.uni-kl.de
     8            Nil Sahin          email: e150916@metu.edu.tr
     9            Maryna Viazovska   email: viazovsk@mathematik.uni-kl.de
     10
     11SEE ALSO: spcurve_lib
    812
    913PROCEDURES:
    10 BlowingUp(f,I,l);          BlowingUp of V(I) at the point 0;
    11 CurveRes(I);               Resolution of V(I)
    12 CurveParam(I);             Parametrization of algebraic branches of V(I)
    13 WSemigroup(X,b);           Weierstrass semigroup of the curve
     14BlowingUp(f,I,l);             BlowingUp of V(I) at the point 0;
     15CurveRes(I);                  Resolution of V(I)
     16CurveParam(I);                Parametrization of algebraic branches of V(I)
     17WSemigroup(X,b);              Weierstrass semigroup of the curve
     18primparam(x,y,c);             HN matrix of parametrization(x(t),y(t))
     19MultiplicitySequence(I);      Multiplicity sequences of the branches of plane curve V(I)
     20CharacteristicExponents(I);   Characteristic exponents of the branches of plane curve V(I)
     21IntersectionMatrix(I);        Intersection Matrix of the branches of plane curve V(I)
     22ContactMatrix(I);             Contact Matrix of the branches of plane curve V(I)
     23plainInvariants(I);           Invariants of the branches of plane curve V(I)
    1424";
    1525
    … …  
    1727LIB "primdec.lib";
    1828LIB "linalg.lib";
     29LIB "ring.lib";
     30LIB "alexpoly.lib";
     31LIB "matrix.lib";
    1932
    2033//////////////////////////////////////////////////////////////
    … …  
    2235//////////////////////////////////////////////////////////////
    2336
    24 proc BlowingUp(poly f,ideal I,list l)
     37proc BlowingUp(poly f,ideal I,list l,list #)
    2538"USAGE:   BlowingUp(f,I,l);
    2639          f=poly
    … …  
    3548
    3649RETURN:   list C of charts.
    37           Each chart C[i] is a list of size 5
     50          Each chart C[i] is a list of size 5 (reps. 6 in case of plane curves)
    3851          C[i][1] is an integer j. It shows, which standard chart do we consider.
    39           C[i][2] is an irreducible polynomial g in k[a]. It is a minimal polynomial for the new parameter.
     52          C[i][2] is an irreducible poly g in k[a]. It is a minimal polynomial
     53                  for the new parameter.
    4054          C[i][3] is an ideal H in k[a].
    4155                  c_i=F_i(a_new) for i=1..n,
    4256                  a_old=H[n+1](a_new).
    43           C[i][4] is a map teta:k[x(1)..x(n),a]-->k[x(1)..x(n),a] from the new curve to the one one.
     57          C[i][4] is a map teta:k[x(1)..x(n),a]-->k[x(1)..x(n),a] from the new
     58                  curve to the old one.
    4459                  x(1)-->x(j)*x(1)
    4560                  . . .
    … …  
    4863                  x(n)-->x(j)*(c_n+x(n))
    4964          C[i][5] is an ideal J of a new curve. J=teta(I).
     65          C[i][6] is the list of exceptional divisors in the chart
    5066
    5167EXAMPLE: example BlowingUp; shows an example"
    … …  
    5975  ideal locI=std(I);
    6076  ideal J=tangentcone(I);
    61 
    6277  setring r;
    6378  ideal J=imap(r1,J);
    6479  ideal locI=imap(r1,locI);
    65   int i,j,k,ind;
    66   list C,C1,Z,Z1;
     80  int j;
     81  int i;
     82  list C,E;
     83  list C1;
    6784  ideal B;
    68   poly g,b;
    69   ideal F,D,D1,I1,I2;
    70   map teta,teta1;
    71 
     85  poly g;
     86  ideal F;
     87  poly b,p;
     88  list Z;
     89  list Z1;
     90  ideal D;
     91  map teta;
     92  ideal D1;
     93  map teta1;
     94  int k,e;
     95  ideal I1;
     96  ideal I2;
     97  int ind;
    7298  list w=mlist(l,n);
    73 
    7499  for(j=1;j<=n;j++)
    75100  {
    76101    B=J;
    77     for(i=1;i<j;i++)
    78     {
    79       B=B+x(w[i]);
    80     }
     102    for(i=1;i<j;i++) {B=B+x(w[i]);}
    81103    B=B+(x(w[j])-1);
    82104    B=B+f;
    … …  
    104126        while(ind==1)
    105127        {
    106           if(gcd(x(w[j]),I1[i])==x(w[j]))
    107           {
    108             I1[i]=I1[i]/x(w[j]);
    109           }
    110           else
    111           {ind=0;}
     128          if(gcd(x(w[j]),I1[i])==x(w[j])){I1[i]=I1[i]/x(w[j]);}
     129          else{ind=0;}
    112130        }
    113131      }
    114132    }
    115 
    116133    for(k=1;k<=size(Z);k++)
    117134    {
    … …  
    123140      C1[3]=F;
    124141      for(i=1;i<j;i++)
    125       {
    126         D[w[i]]=x(w[j])*x(w[i]);
    127       }
     142      {D[w[i]]=x(w[j])*x(w[i]);}
    128143      D[w[j]]=x(w[j]);
    129144      for(i=j+1;i<=n;i++)
    130       {
    131         D[w[i]]=x(w[j])*(F[w[i]]+x(w[i]));
    132       }
     145      {D[w[i]]=x(w[j])*(F[w[i]]+x(w[i]));}
    133146      D[n+1]=Z[k][2][n+1];
    134147      teta=r,D;
    135148      C1[4]=D;
    136149      for(i=1;i<=j;i++)
    137       {
    138         D1[w[i]]=x(w[i]);
    139       }
     150      {D1[w[i]]=x(w[i]);}
    140151      for(i=j+1;i<=n;i++)
    141       {
    142         D1[w[i]]=F[w[i]]+x(w[i]);
    143       }
     152      {D1[w[i]]=F[w[i]]+x(w[i]);}
    144153      D1[n+1]=a;
    145154      teta1=r,D1;
    146       I2=teta1(I1);
    147       I2=reduce(I2,std(g));
     155      if(nvars(basering)==3)
     156      {
     157        I2=quickSubst(I1[1],teta1[1],teta1[2],std(g));
     158      }
     159      else
     160      {
     161        I2=teta1(I1);
     162        I2=reduce(I2,std(g));
     163      }
    148164      C1[5]=I2;
     165      if(nvars(basering)==3)
     166      {
     167        if(size(#)>0)
     168        {
     169          E=#;
     170          E=teta(E);
     171          for(e=1;e<=size(E);e++)
     172          {
     173            p=E[e];
     174            while(subst(p,x(w[j]),0)==0)
     175            {
     176              p=p/x(w[j]);
     177            }
     178            if((deg(E[e])>0)&&(deg(p)==0))
     179            {
     180              E[e]=size(E);
     181            }
     182            else
     183            {
     184              E[e]=p;
     185            }
     186          }
     187          E[size(E)+1]=x(w[j]);
     188          C1[6]=E;
     189        }
     190        else
     191        {
     192          C1[6]=list(x(w[j]));
     193        }
     194      }
    149195      C=insert(C,C1);
    150196    }
    … …  
    153199}
    154200example
    155 { "EXAMPLE:"; echo=2;
     201{
     202  "EXAMPLE:";echo = 2;
    156203  ring r=0,(x(1..3),a),dp;
    157204  poly f=a2+1;
    … …  
    161208  B;
    162209}
    163 
    164 
     210//============= ACHTUNG ZeroIdeal ueberarbeiten / minAssGTZ rein  ========================
    165211//////////////////////////////////////////////////////////////////////////////////////////
    166212static proc ZeroIdeal(ideal J)
     213
    167214"USAGE:   ZeroIdeal(J);
    168215          J=ideal
     216
    169217ASSUME:   J is a zero-dimensional ideal in k[x(1),...,x(n)].
     218
    170219COMPUTE:  Primary decomposition of radical(J). Each prime ideal J[i] has the form:
    171220          x(1)-f[1](b),...,x(n)-f[n](b),
    … …  
    178227          Z[k][2] is an ideal H, H[n]=f[n],
    179228          Z[k][3] is a poly x(1)*a(1)+...+x(n)*a(n)
    180 "
    181 {
     229
     230EXAMPLE:"
     231{
     232  intvec opt = option(get);
    182233  def r=basering;
    183234  int n=nvars(r);
    184235  if(dim(std(J))!=0){return(0);}
    185236  ring s=0,(x(1..n)),lp;
    186   ideal A,S;
    187   int i,j,ind,q;
    188   for(i=1;i<=n;i++)
    189   {A[i]=x(i);}
     237  ideal A; ideal S; int i; int j;
     238  for(i=1;i<=n;i++) {A[i]=x(i);}
    190239  map phi=r,A;
    191240  ideal J=phi(J);
    192241  ideal I=radical(J);
    193242  list D=zerodec(I);
    194   list Z,u;
    195   ideal H,T,Di;
    196   intvec w,v;
    197   map tau;
    198   poly h;
    199 
     243  list Z; ideal H; intvec w; intvec v; int ind; ideal T; map tau; int q; list u;
     244  ideal Di; poly h;
    200245  for(i=1;i<=size(D);i++)
    201246  {
    … …  
    271316  list Z;
    272317  for(i=1;i<=n;i++)
    273   {
    274     A[i]=var(i);
    275   }
     318  {A[i]=var(i);}
    276319  map psi=s,A;
    277320  Z=psi(Z);
     321  option(set, opt);
    278322  return(Z);
    279323}
    280 
    281324/////////////////////////////////////////////////////////////////////////////////////////////
    282325//assume that the basering is k[x(1),...,x(n),a]
    283326
    284 static proc main(ideal I,ideal Psi,poly f,list m,list l,list HN)
     327static proc main(ideal I,ideal Psi,poly f,list m,list l,list HN,intvec v,list HI,list #)
    285328{
    286329  def s=basering;
    287   int i,j;
    288   list C,C1,C2,C3,l1,m1,HN1;
     330  int i,z;
     331  int j;
     332  list C,E,resTree;
     333  list C1;
     334  list C2;
     335  list C3;
     336  list l1;
     337  C2[8]=HI;
     338  list m1;
     339  list HN1;
    289340  ideal J;
    290341  map psi;
    291 
    292   if(SmoothTest(I,f)==1)
     342  intvec w;
     343  z=(SmoothTest(I,f)==1);
     344  if((nvars(basering)==3)&&z&&(size(#)>0))
     345  {
     346    z=transversalTest(I,f,#);
     347  }
     348  if(z)
    293349  {
    294350    C2[1]=I;
    … …  
    298354    C2[5]=l;
    299355    C2[6]=HN;
     356    if(nvars(basering)==3)
     357    {
     358      if(size(#)>0)
     359      {
     360        C2[9]=#;
     361      }
     362      C2[7]=v;
     363    }
     364    //C2[8][size(C2[8])+1]=list(C2[7],C2[9]);
    300365    C[1]=C2;
    301366  }
    302   if(SmoothTest(I,f)==0)
     367  if(!z)
    303368  {
    304369    int mm=mmult(I,f);
    305370    m1=insert(m,mm,size(m));
    306     C1=BlowingUp(f,I,l);
     371    if(nvars(basering)==3)
     372    {
     373      if(size(#)>0)
     374      {
     375        E=#;
     376        C1=BlowingUp(f,I,l,E);
     377      }
     378      else
     379      {
     380        C1=BlowingUp(f,I,l);
     381      }
     382    }
     383    else
     384    {
     385      C1=BlowingUp(f,I,l);
     386    }
    307387    for(j=1;j<=size(C1);j++)
    308388    {
    … …  
    318398      HN1=insert(HN1,C1[j][3],size(HN1)-1);
    319399      C2[6]=HN1;
    320       C3=main(C2[1],C2[2],C2[3],C2[4],C2[5],C2[6]);
    321       C=C+C3;
     400      if(deg(C2[3])>1)
     401      {
     402        w=v,-j;
     403      }
     404      else
     405      {
     406        w=v,j;
     407      }
     408      C2[7]=w;
     409      if(nvars(basering)==3)
     410      {
     411        C2[9]=C1[j][6];
     412        C2[8][size(C2[8])+1]=list(C2[7],C2[9]);
     413        C3=main(C2[1],C2[2],C2[3],C2[4],C2[5],C2[6],C2[7],C2[8],C2[9]);
     414        C=C+C3;
     415      }
     416      else
     417      {
     418        C3=main(C2[1],C2[2],C2[3],C2[4],C2[5],C2[6],C2[7],C2[8]);
     419        C=C+C3;
     420      }
    322421    }
    323422  }
    324423  return(C);
    325424}
    326 
    327425////////////////////////////////////////////////////////////////////////////////////////////////
    328426
     427static proc transversalTest(ideal I,poly f,list L)
     428{
     429  def r=basering;
     430  int n=nvars(r)-1;
     431  int i;
     432  ring r1=(0,a),(x(1..n)),ds;
     433  number f=leadcoef(imap(r,f));
     434  minpoly=f;
     435  ideal I=imap(r,I);
     436  list L=imap(r,L);
     437  ideal K=jet(L[size(L)],deg(lead(L[size(L)])));
     438  ideal T=1;
     439  if(size(L)>1)
     440  {
     441    for(i=1;i<=size(L)-1;i++)
     442    {
     443      if(subst(L[i],x(1),0,x(2),0)==0) break;
     444    }
     445    if(i<=size(L)-1)
     446    {
     447      T=jet(L[i],deg(lead(L[i])));
     448    }
     449  }
     450  ideal J=jet(I[1],deg(lead(I[1])));
     451  setring r;
     452  ideal J=imap(r1,J);
     453  ideal K=imap(r1,K);
     454  ideal T=imap(r1,T);
     455  int m=size(reduce(J,std(K)))+size(reduce(K,std(J)));
     456  if(m)
     457  {
     458    m=size(reduce(J+K+T,std(ideal(x(1),x(2)))));
     459  }
     460  return(m);
     461}
     462////////////////////////////////////////////////////////////////////////////////////////////////
    329463static proc SmoothTest(ideal I,poly f)
    330464//Assume I is a radical ideal of dimension 1 in a ring k[x(1..n),a]
    331465//Returns 1 if a curve V(I) is smooth at point 0 and returns 0 otherwise
    332466{
    333   int ind,l;
     467  int ind;
     468  int l;
    334469  def t=basering;
    335470  int n=nvars(t)-1;
    … …  
    343478  return(ind);
    344479}
    345 
    346 
    347480////////////////////////////////////////////////////////////////////////////////////////////////
    348 
    349481proc CurveRes(ideal I)
    350482"USAGE:   CurveRes(I);
    351483          I ideal
    352 
    353484ASSUME:   The basering is r=0,(x(1..n))
    354485          V(I) is a curve with a singular point 0.
    355 
    356486COMPUTE:  Resolution of the curve V(I).
    357 
    358487RETURN:   a ring R=basering+k[a]
    359488          Ring R contains a list Resolve
    360489          Resolve is a list of charts
    361490          Each Resolve[i] is a list of size 6
    362 
    363491          Resolve[i][1] is an ideal J of a new curve. J=teta(I).
    364492          Resolve[i][2] ideal which represents the map
    365493                        teta:k[x(1)..x(n),a]-->k[x(1)..x(n),a] from the
    366494                        new curve to the old one.
    367           Resolve[i][3] is an irreducible polynomial g in k[a]. It is a minimal polynomial for the new parameter a.
    368           Resolve[i][4] sequence of multiplicities
     495          Resolve[i][3] is an irreducible poly g in k[a]. It is a minimal polynomial for the
     496                        new parameter a. deg(g) gives the number of branches in Resolve[i]
     497          Resolve[i][4] sequence of multiplicities (sum over all branches in Resolve as long as
     498                        they intersect each other !)
    369499          Resolve[i][5] is a list of integers l. It shows, which standard charts we considered.
    370500          Resolve[i][6] HN matrix
     501          Resolve[i][7] (only for plane curves) the development of exceptional divisors
     502                        the entries correspond to the i-th blowing up. The first entry is an
     503                        intvec. The first negative entry gives the splitting of the (over Q
     504                        irreducible) branches. The second entry is a list of the exceptional
     505                        divisors. If the entry is an integer i, it says that the divisor is not
     506                        visible in this chart after the i-th blowing up.
    371507
    372508EXAMPLE: example CurveRes; shows an example"
    … …  
    376512  ring s=0,(x(1..n),a),dp;
    377513  ideal A;
    378   int i,j;
    379   for(i=1;i<=n;i++)
    380   {A[i]=x(i);}
     514  int i;
     515  int j;
     516  for(i=1;i<=n;i++){A[i]=x(i);}
    381517  map phi=r,A;
    382518  ideal I=phi(I);
    383519  poly f=a;
    384   list l,m;
     520  list l;
     521  list m;
    385522  list HN=x(1);
    386523  ideal psi;
    387   for(i=1;i<=n;i++)
    388   {psi[i]=x(i);}
     524  for(i=1;i<=n;i++){psi[i]=x(i);}
    389525  psi[n+1]=a;
    390   list Resolve=main(I,psi,f,m,l,HN);
    391   for(i=1;i<=size(Resolve);i++)
    392   {
    393     Resolve[i][6]=delete(Resolve[i][6],size(Resolve[i][6]));
     526  intvec v;
     527  list L,Resolve;
     528  if(n==2)
     529  {
     530    ideal J=factorize(I[1],1);
     531    list resolve;
     532    for(int k=1;k<=size(J);k++)
     533    {
     534      I=J[k];
     535      resolve=main(I,psi,f,m,l,HN,v,L);
     536      for(i=1;i<=size(resolve);i++)
     537      {
     538        resolve[i][6]=delete(resolve[i][6],size(resolve[i][6]));
     539        if(size(resolve[i])>=9){resolve[i]=delete(resolve[i],9);}
     540        resolve[i]=delete(resolve[i],7);
     541      }
     542      if(k==1){Resolve=resolve;}
     543      else{Resolve=Resolve+resolve;}
     544    }
     545  }
     546  else
     547  {
     548    Resolve=main(I,psi,f,m,l,HN,v,L);
     549    for(i=1;i<=size(Resolve);i++)
     550    {
     551      Resolve[i][6]=delete(Resolve[i][6],size(Resolve[i][6]));
     552      Resolve[i]=delete(Resolve[i],8);
     553    }
    394554  }
    395555  export(Resolve);
    … …  
    397557}
    398558example
    399 {"EXAMPLE:"; echo=2;
     559{
     560  "EXAMPLE:"; echo=2;
    400561  ring r=0,(x,y,z),dp;
    401562  ideal i=x2-y3,z2-y5;
    … …  
    404565  Resolve;
    405566}
    406 
    407567//////////////////////////////////////////////////////////////////
    408 
    409568static proc mlist(list l,int n)
    410569{
    411   list N,M;
    412   int i,j;
    413   for(i=1;i<=n;i++)
    414   {
    415     M[i]=i;
    416   }
     570  list N;
     571  list M;
     572  int i;
     573  int j;
     574  for(i=1;i<=n;i++) {M[i]=i;}
    417575  N=l+M;
    418576  for(i=1;i<=size(N)-1;i++)
    … …  
    422580    {
    423581      if(N[i]==N[j]){N=delete(N,j);}
    424       else{j++;}
     582      else
     583      {j++;}
    425584    }
    426585  }
    427586  return(N);
    428587}
    429 
    430588/////////////////////////////////////////////////////////////////////
    431589//Assume that the basering is k[x(1..n),a]
    … …  
    442600  return(m);
    443601}
    444 
    445602//////////////////////////////////////////////////////////////
    446603//--------Parametrization of smooth curve-------------------//
    447604//////////////////////////////////////////////////////////////
    448 
    449605
    450606////////////////////////////////////////////////////////////////////////
    … …  
    457613  int k=size(I);
    458614  matrix M[k][n];
    459   int i,j,l;
     615  int i;
     616  int j;
     617  int l;
    460618  for(i=1;i<=k;i++)
    461619  {
    … …  
    468626  return(M);
    469627}
    470 
    471 
    472628//////////////////////////////////////////////////////////
    473 
    474629static proc mmi(matrix M,int n)
    475630{
    476631  ideal l;
    477632  int k=nrows(M);
    478   int i,j;
     633  int i;
     634  int j;
    479635  for(i=1;i<=k;i++)
    480636  {
    … …  
    510666  return(lmi);
    511667}
    512 
    513 
    514668//////////////////////////////////////////////////////////
    515 
    516669static proc mC(matrix Mi,int n)
    517670{
    … …  
    530683  return(c);
    531684}
    532 
    533685//////////////////////////////////////////////////////////
    534 
    535 proc mmF(ideal C, matrix Mi,int n,int k)
     686static proc mmF(ideal C, matrix Mi,int n,int k)
    536687{
    537688  int s=size(C);
    … …  
    539690  int p=0;
    540691  int t=0;
    541   int i, j;
     692  int i;
     693  int j;
    542694  int v=0;
    543695  for(i=s;i>0;i--)
    … …  
    562714  return(mmf);
    563715}
    564 
    565 
    566716/////////////////////////////////////////////////////
    567 
    568717static proc cparam(ideal I,poly f,int n,int m,int N)
    569718{
    … …  
    589738  matrix B=imap(q,B);
    590739  matrix C=inverse(B);
    591 
    592   int i,j,l;
     740  int i;
     741  int j;
    593742  ideal P;
    594743  for(i=1;i<mi;i++){P[i]=x(i);}
    … …  
    597746  map phi=s,P;
    598747  ideal I1=phi(I);
     748  if(nvars(basering)==2)
     749  {
     750    setring r;
     751    ideal I1=imap(s,I1);
     752    matrix C=imap(s,C);
     753    list X;
     754    matrix d[n-1][1];
     755    matrix b[n-1][1];
     756    ideal Q;
     757    map psi;
     758    int l;
     759    for(i=1;i<=N;i++)
     760    {
     761      for(j=1;j<=n-1;j++)
     762      {
     763        d[j,1]=diff(I1[mf[j]],x(n));
     764        for(l=1;l<=n;l++)
     765        {
     766          d[j,1]=subst(d[j,1],x(l),0);
     767        }
     768      }
     769      b=-C*d;
     770      I1=jet(I1,N-i+2);
     771      X[i]=b;
     772      for(j=1;j<=n-1;j++){Q[j]=x(n)*(b[j,1]+x(j));}
     773      Q[n]=x(n);
     774      I1[1]=quickSubst(I1[1],Q[1],Q[2],std(f));
     775      I1=I1/x(n);
     776    }
     777    list Y=X,mi;
     778    return(Y);
     779  }
    599780  list X;
    600781  matrix d[n-1][1];
    … …  
    602783  ideal Q;
    603784  map psi;
     785  int l;
    604786  for(i=1;i<=N;i++)
    605787  {
    … …  
    623805  return(Y);
    624806}
    625 
    626807//////////////////////////////////////////////////////////////
    627808//--------Parametrization of singular curve-----------------//
    628809//////////////////////////////////////////////////////////////
    629 
    630 proc CurveParam (ideal I)
     810proc CurveParam (list #)
    631811"USAGE:   CurveParam(I);
    632812          I ideal
    633 
    634813ASSUME:   I is an ideal of a curve C with a singular point 0.
    635 
    636814COMPUTE:  Parametrization for algebraic branches of the curve C.
    637 
    638815RETURN:   list L of size 1.
    639816          L[1] is a ring ring rt=0,(t,a),ds;
    … …  
    641818          Param is a list of algebraic branches
    642819          Each Param[i] is a list of size 3
    643 
    644820          Param[i][1] is a list of polynomials
    645           Param[i][2] is an irredusible polynomial f\in k[a].It is a minimal polynomial for the parameter a.
     821          Param[i][2] is an irredusible polynomial f\in k[a].It is a minimal polynomial for
     822                      the parameter a.
    646823          Param[i][3] is an integer b--upper bound for the conductor of Weierstrass semigroup
    647824
    648 
    649825EXAMPLE: example curveparam; shows an example"
    650826{
    651   int i, j,mi,b,k;
    652   def s=CurveRes(I);
     827  int i;
     828  int j;
     829  if(typeof(#[1])=="ideal")
     830  {
     831    int d=deg(#[1][1]);
     832    def s=CurveRes(#[1]);
     833  }
     834  else
     835  {
     836    def s=#[1];
     837  }
    653838  setring s;
    654839  int n=nvars(s)-1;
    655   list Param,l;
     840  list Param;
     841  list l;
    656842  ideal D,P,Q,T;
    657843  poly f;
    658844  map tau;
    659   list Z, Y,X;
     845  list Z;
     846  list Y;
     847  list X;
     848  int mi;
     849  int b;
     850  int k;
     851  int dd;
    660852  for(j=1;j<=size(Resolve);j++)
    661853  {
    … …  
    665857      b=b+Resolve[j][4][k]*(Resolve[j][4][k]+1);
    666858    }
    667     if(b==0){b=1;}
     859    if((n==2)&&(size(Resolve[j][4])==0))
     860    {b=d;}
    668861    Y=cparam(Resolve[j][1],Resolve[j][3],n,1,b);
    669862    X=Y[1];
    … …  
    684877    tau=s,P;
    685878    T=Resolve[j][2];
    686     Q=tau(T);
     879//HERE A TEST FOR dd
     880    if(nvars(basering)==3)
     881    {
     882      dd=boundparam(P[2]);
     883      if(dd==1){dd=boundparam(P[1]);}
     884      P[1]=jet(P[1],dd);
     885      P[2]=jet(P[2],dd);
     886      Q[1]=quickSubst(T[1],P[1],P[2],std(f));
     887      Q[2]=quickSubst(T[2],P[1],P[2],std(f));
     888      Q[3]=a;
     889    }
     890    else
     891    {
     892      Q=tau(T);
     893    }
    687894    for(i=1;i<=n;i++){Z[i]=jet(reduce(Q[i],std(f)),b+1);}
    688895    l[1]=Z;
    … …  
    700907  return(rt);
    701908}
    702 
    703909example
    704 {"EXAMPLE:";echo=2;
     910{
     911  "EXAMPLE:";echo=2;
    705912  ring r=0,(x,y,z),dp;
    706913  ideal i=x2-y3,z2-y5;
    … …  
    708915  setring s;
    709916  Param;
    710   ring r=0,(x,y,z),dp;
    711   ideal i=x2-y3,z2-y5;
    712   def s=CurveParam(i);
    713   setring s;
    714   Param;
    715 }
    716 
     917}
    717918///////////////////////////////////////////////////////////////////////////////////////////
    718919//----------Computation of Weierstrass Semigroup from parametrization--------------------//
    719920///////////////////////////////////////////////////////////////////////////////////////////
    720 
    721 proc Semi(intvec G,int b)
     921static proc Semi(intvec G,int b)
    722922"USAGE:   Semi(G,b);
    723923          G=intvec
    724924          b=int
    725 
    726925ASSUME:   G[1]<=G[2]<=...<=G[k],
    727 
    728926COMPUTE:  elements of semigroup S generated by the enteries of G till the bound b.
    729927          For each element i of S computes the list of integer vectors v of dimension
    … …  
    731929          conductor  of semigroup S  c<b-n, where n is minimal element of G, then
    732930          computes also c+n.
    733 
    734931RETURN:   list M of size 2.
    735932          L=M[1] is a list  of size min(b,c+n).
    … …  
    738935          M[2] is an integer =min(b,c+n)
    739936          M[3] minimal generators of S
    740 "
    741 {
    742   list L,l;
    743   int i,j,t,q;
     937EXAMPLE:"
     938{
     939  list L;
     940  list l;
     941  int i;
    744942  for(i=1;i<=b;i++){L[i]=l;}
    745943  int k=size(G);
    746944  int n=G[1];
     945  int j;
     946  int t;
     947  int q;
    747948  int c=0;
    748   intvec w,v;
     949  intvec w;
     950  intvec v;
    749951  for(i=1;i<=k;i++)
    750952  {
    … …  
    791993  for(j=1;j<=k;j++)
    792994  {
    793     if(size(L[G[j]])==size(L1[G[j]]))
     995    if(size(L[G[j]])==size(L1[G[j]]) && G[j]<i)
    794996    {
    795997      Gmin[jmin]=G[j];
    … …  
    8011003  return(M);
    8021004}
    803 
    8041005///////////////////////////////////////////////////////////////////////////////////////////
    8051006static proc AddElem(list L,int b,int k,int g,int n)
    … …  
    8091010          g=int new generator
    8101011          n=int. minimal generator
    811 
    8121012RETURN: list M
    8131013        M[1]=new L;
    … …  
    8571057  return(M);
    8581058}
    859 
    8601059///////////////////////////////////////////////////////////////////////////////////////////
    8611060proc WSemigroup(list X,int b0)
    … …  
    8631062           X a list of polinomials in one vaiable, say t.
    8641063           b0 an integer
    865 
    866 COMPUTE:   Weierstrass semigroup of space curve C,which is given by a parametrization X[1](t),...,X[k](t), till the bound b0.
     1064COMPUTE:   Weierstrass semigroup of space curve C,which is given by a parametrization
     1065          X[1](t),...,X[k](t), till the bound b0.
    8671066
    8681067ASSUME:    b0 is greater then conductor
    869 
    8701068RETURN:    list M of size 5.
    8711069           M[1]= list of integers, which are minimal generators set of the Weierstrass semigroup.
    … …  
    8731071           M[3]=intvec, all elements of the Weierstrass semigroup till some bound b,
    8741072           which is greather than conductor.
    875 
    876 
    877 
    878 
    8791073WARNING:   works only over the ring with one variable with ordering ds
    880 
    8811074EXAMPLE: example WSemigroup; shows an example"
     1075
    8821076{
    8831077  int k=size(X);
    8841078  intvec G;
    885   int i,i2,g;
     1079  int i,i2;
    8861080  poly t=var(1);
    8871081  poly h;
     1082  int g;
    8881083  for(i=1;i<=k;i++)
    889   {
    890     G[i]=ord(X[i]);
    891   }
     1084  {G[i]=ord(X[i]);}
    8921085  for(i=1;i<k;i++)
    8931086  {
    … …  
    9061099  G=U[3];
    9071100  int k1=size(G);
    908   list N, l;
     1101  list N;
     1102  list l;
    9091103  for(i=1;i<=b;i++){N[i]=l;}
    9101104  int j;
    9111105  for(j=b0;j>b;j--){L=delete(L,j);}
    9121106  poly p;
    913   int s, e;
     1107  int s;
     1108  int e;
    9141109  for(i=1;i<=b;i++)
    9151110  {
    … …  
    9281123    }
    9291124  }
    930   int j1, j2,m,b1,q,i1;
     1125  int j1;
     1126  int j2;
    9311127  list M;
    932   poly c1, c2, f;
     1128  poly c1;
     1129  poly c2;
     1130  poly f;
     1131  int m;
     1132  int b1;
    9331133  ideal I;
    934   matrix C, C1;
     1134  matrix C;
     1135  matrix C1;
     1136  int q;
     1137  int i1;
    9351138  i=1;
    9361139  while(i<=b)
    … …  
    9961199}
    9971200example
    998 {"EXAMPLE:";echo=2;
     1201{
     1202  "EXAMPLE:";echo=2;
    9991203  ring r=0,(t),ds;
    10001204  list X=t4,t5+t11,t9+2*t7;
    … …  
    10021206  L;
    10031207}
    1004 
    10051208////////////////////////////////////////////////////////////////////////////////////////////
    1006 /*
    1007 LIB"spacecurve.lib";
    1008 
    1009 ring r=0,(x,y,z),dp;
    1010 ideal i=x2-y3,z2-y5;
    1011 def s=CurveParam(i);
    1012 setring s;
    1013 Param;
    1014 
    1015 ring r=0,(t),ds;
    1016 list X=t4,t5+t11,t9+2*t7;
    1017 list L=WSemigroup(X,30);
    1018 L;
    1019 
    1020 
    1021 ring r=0,(x,y,z,t),dp;
    1022 ideal I=x-t4,y-t5-t11,z-t9-2t7;
    1023 I=eliminate(I,t);
    1024 ring r1=0,(x,y,z),dp;
    1025 ideal I=imap(r,I);
    1026 def s=CurveParam(I);
    1027 setring s;
    1028 Param;
    1029 WSemigroup(Param[1][1],30);
    1030 
    1031 ring r=0,(x,y,z),dp;
    1032 ideal I=(x5-y11)*(x2+(y+1)^3),z;
    1033 def s=CurveParam(I);
    1034 setring s;
    1035 Param;
    1036 
    1037 
    1038 ring r=0,(x,y,z),dp;
    1039 ideal I=y2-x3-x2,z;
    1040 def s=CurveParam(I);
    1041 setring s;
    1042 Param;
    1043 setring r;
    1044 I=intersect(I,ideal(y2-x3,z));
    1045 s=CurveParam(I);
    1046 setring s;
    1047 Param;
    1048 
    1049 ring r=0,(x,y,z),dp;
    1050 ideal I=y2+x3+x2,z;
    1051 def s=CurveParam(I);
    1052 setring s;
    1053 Param;
    1054 */
     1209
     1210static proc quickSubst(poly h, poly r, poly s,ideal I)
     1211{
     1212//=== computes h(r,s) mod I for h in Q[x(1),x(2),a]
     1213  attrib(I,"isSB",1);
     1214  if((r==x(1))&&(s==x(2))){return(reduce(h,I));}
     1215  poly q1 = 1;
     1216  poly q2 = 1;
     1217  poly q3 = 1;
     1218  int i,j,e1,e2,e3;
     1219  list L,L1,L2,L3;
     1220  if(r==x(1))
     1221  {
     1222    matrix M=coeffs(h,x(2));
     1223    L[1]=1;
     1224    for(i=2;i<=nrows(M);i++)
     1225    {
     1226      q2 = reduce(q2*s,I);
     1227      L[i]=q2;
     1228    }
     1229    i=1;
     1230    h=0;
     1231    while(i <= nrows(M))
     1232    {
     1233      if(M[i,1]!=0)
     1234      {
     1235         h=h+M[i,1]*L[i];
     1236      }
     1237      i++;
     1238    }
     1239    h=reduce(h,I);
     1240    return(h);
     1241  }
     1242  if(s==x(2))
     1243  {
     1244    matrix M=coeffs(h,x(1));
     1245    L[1]=1;
     1246    for(i=2;i<=nrows(M);i++)
     1247    {
     1248      q1 = reduce(q1*r,I);
     1249      L[i]=q1;
     1250    }
     1251    i=1;
     1252    h=0;
     1253    while(i <= nrows(M))
     1254    {
     1255      if(M[i,1]!=0)
     1256      {
     1257        h=h+M[i,1]*L[i];
     1258      }
     1259      i++;
     1260    }
     1261    h=reduce(h,I);
     1262    return(h);
     1263  }
     1264  for(i=1;i<=size(h);i++)
     1265  {
     1266    if(leadexp(h[i])[1]>e1){e1=leadexp(h[i])[1];}
     1267    if(leadexp(h[i])[2]>e2){e2=leadexp(h[i])[2];}
     1268    if(leadexp(h[i])[3]>e3){e3=leadexp(h[i])[3];}
     1269  }
     1270  for(i = 1; i <= size(h); i++)
     1271  {
     1272    L[i] = list(leadcoef(h[i]),leadexp(h[i]));
     1273  }
     1274  L1[1]=1;
     1275  L2[1]=1;
     1276  L3[1]=1;
     1277  for(i=1;i<=e1;i++)
     1278  {
     1279    q1 = reduce(q1*r,I);
     1280    L1[i+1]=q1;
     1281  }
     1282  for(i=1;i<=e2;i++)
     1283  {
     1284    q2 = reduce(q2*s,I);
     1285    L2[i+1]=q2;
     1286  }
     1287  for(i=1;i<=e3;i++)
     1288  {
     1289    q3 = reduce(q3*var(3),I);
     1290    L3[i+1]=q3;
     1291  }
     1292  int m=size(L);
     1293  i = 1;
     1294  h = 0;
     1295  while(i <= m)
     1296  {
     1297    h=h+L[i][1]*L1[L[i][2][1]+1]*L2[L[i][2][2]+1]*L3[L[i][2][3]+1];
     1298    i++;
     1299  }
     1300  h=reduce(h,I);
     1301  return(h);
     1302}
     1303
     1304static proc semi2char(intvec v)
     1305{
     1306  intvec k=v[1..2];
     1307  intvec w=v[1];
     1308  int i,j,p,q;
     1309  for(i=2;i<size(v);i++)
     1310  {
     1311    w[i]=gcd(w[i-1],v[i]);
     1312  }
     1313  for(i=3;i<=size(v);i++)
     1314  {
     1315    k[i]=v[i];
     1316    for(j=2;j<i;j++)
     1317    {
     1318      k[i]=k[i]-(w[j-1] div w[j]-1)*v[j];
     1319    }
     1320  }
     1321  return(k);
     1322}
     1323/////////////////////////////////////////////////////////////////////////////////////////////////
     1324proc primparam(poly x,poly y,int c)
     1325"USAGE:  MultiplicitySequence(x,y,c);  x poly, y poly, c integer
     1326ASSUME:  x and y are polynomials in k(a)[t] such that (x,y) is a primitive parametrization of
     1327         a plane curve branch and ord(x)<ord(y).
     1328RETURN:  Hamburger-Noether Matrix of the curve branch given parametrically by (x,y).
     1329EXAMPLE: example primparam;  shows an example
     1330"
     1331{
     1332  int i,h;
     1333  poly F,z;
     1334  list L;
     1335  while(ord(x)>1)
     1336  {
     1337    list v;
     1338    while(ord(y)>=ord(x))
     1339    {
     1340      F=divide(y,x,c);
     1341      if(ord(F)==0)
     1342      {
     1343        v=insert(v,subst(F,t,0),size(v));
     1344        y=F-subst(F,t,0);
     1345      }
     1346      else
     1347      {
     1348        v=insert(v,0,size(v));
     1349        y=F;
     1350      }
     1351    }
     1352    v=insert(v,t,size(v));
     1353    L=insert(L,transform(v),size(L));
     1354    z=x;
     1355    x=y;
     1356    y=z;
     1357    kill v;
     1358  }
     1359  if(ord(x)==1)
     1360  {
     1361    list v;
     1362    while(i<c)
     1363    {
     1364      F=divide(y,x,c);
     1365      if(ord(F)==0)
     1366      {
     1367        v=insert(v,subst(F,t,0),size(v));
     1368        y=F-subst(F,t,0);
     1369      }
     1370      else
     1371      {
     1372        v=insert(v,0,size(v));
     1373        y=F;
     1374      }
     1375      if(y==0)
     1376      {
     1377        v=insert(v,t,size(v));
     1378        break;
     1379      }
     1380      i++;
     1381    }
     1382    L=insert(L,transform(v),size(L));
     1383  }
     1384  return(compose(L));
     1385}
     1386example
     1387{
     1388  "EXAMPLE:"; echo=2;
     1389  ring r=(0,a),t,ds;
     1390  poly x=t6;
     1391  poly y=t8+t11;
     1392  int c=15;
     1393  primparam(x,y,c);
     1394}
     1395
     1396//////////////////////////////////////
     1397//L is a list of polynomials
     1398//////////////////////////////////////
     1399static proc transform(list L)
     1400{
     1401  matrix m2;
     1402  matrix m1=matrix(L[1]);
     1403  for(int i=2;i<=size(L);i++)
     1404  {
     1405    m2=matrix(L[i]);
     1406    m1=concat(m1,m2);
     1407  }
     1408  return(m1);
     1409}
     1410/////////////////////////////////////////////////////////////////
     1411//L is a list of matrices
     1412///////////////////////////////////////
     1413static proc compose(list L)
     1414{
     1415  matrix M[ncols(L[1])][1]=transpose(L[1]);
     1416  for(int i=2;i<=size(L);i++)
     1417  {
     1418    M=concat(M,transpose(L[i]));
     1419  }
     1420  return(transpose(M));
     1421}
     1422//////////////////////////////////////////////////////////
     1423static proc rduce(poly p)
     1424{
     1425  int n=ord(p);
     1426  poly q=p/(t^n);
     1427  return(q);
     1428}
     1429////////////////////////////////////////////////////
     1430static proc divide(poly p,poly q,int c)i
     1431{
     1432  int n=ord(p);
     1433  int m=ord(q);
     1434  poly p'=rduce(p);
     1435  poly q'=rduce(q);
     1436  poly r=t^(n-m)*p'*jet(1,q',c);
     1437  return(jet(r,c));
     1438}
     1439/////////////////////////////////////////////////////
     1440static proc contact(list L)
     1441{
     1442  def M=L[1];
     1443  intvec v=L[2];
     1444  int s,j,i;
     1445  for(i=1;i<=size(v);i++)
     1446  {
     1447    if(v[i]<0){v[i]=-1-v[i];}
     1448    for(j=1;j<=v[i];j++)
     1449    {
     1450      s=s+1;
     1451      if(find(string(M[i,j]),"a")!=0){return(s);}
     1452    }
     1453  }
     1454}
     1455///////////////////////////////////////////////////////////
     1456static proc converter(list L)
     1457{
     1458def s=basering;
     1459list D;
     1460int i,c;
     1461for(i=1;i<=size(L);i++)
     1462   {D=insert(D,deg(L[i][2]),size(D));}
     1463ring r=(0,a),(t),ds;
     1464list L=imap(s,L);
     1465poly x,y,z;
     1466list A,B;
     1467for(i=1;i<=size(L);i++)
     1468   {A[5]=D[i];
     1469   x=L[i][1][1];
     1470   y=L[i][1][2];
     1471   if(ord(x)<=ord(y)){A[3]=0;}
     1472   else{A[3]=1;
     1473       z=x;
     1474       x=y;
     1475       y=z;
     1476       }
     1477   c=bound(x,y);
     1478   if(c==-1){ERROR("Bound is not enough");}
     1479   A[1]=primparam(x,y,c);
     1480   A[2]=lengths(A[1]);
     1481   A[4]=0;
     1482   B=insert(B,A,size(B));
     1483   A=list();
     1484   }
     1485ring r1=(0,a),(x,y),ds;
     1486list hne=fetch(r,B);
     1487export(hne);
     1488return(r1);
     1489}
     1490//////////////////////////////////////////////////////////
     1491static proc intermat(list L)
     1492{
     1493    int s=size(L);
     1494    intvec v=L[1][5];
     1495    intvec w1;
     1496    int i,j,d,b,l,k,c,o,p;
     1497    for(i=2;i<=s;i++)
     1498    {v=v,L[i][5];}
     1499    intvec u=v[1];
     1500    for(i=2;i<=s;i++)
     1501    {
     1502      l=u[size(u)]+v[i];
     1503      u=u,l;
     1504    }
     1505    int m=u[size(u)];
     1506    intmat M[m][m];
     1507    for(i=1;i<=m;i++)
     1508    {
     1509      for(j=i+1;j<=m;j++)
     1510      {
     1511        d=sorting(u,i);
     1512        b=sorting(u,j);
     1513        if(d==b){k=contact(L[d]);
     1514        w1=multsequence(L[d]);
     1515        if(size(w1)<k){for(p=size(w1)+1;p<=k;p++)
     1516        {w1=w1,1;} }
     1517        for(o=1;o<=k;o++){c=c+w1[o]*w1[o];}
     1518        M[i,j]=c;
     1519        c=0;
     1520      }
     1521      else
     1522      {M[i,j]=intersection(L[d],L[b]);}
     1523    }
     1524  }
     1525  return(M);
     1526}
     1527/////////////////////////////////////////////////////////////////
     1528static proc lengths(matrix M)
     1529{
     1530  intvec v;
     1531  int s,i,j;
     1532  for(i=1;i<=nrows(M);i++)
     1533  {
     1534    for(j=1;j<=ncols(M);j++)
     1535    {
     1536      if(M[i,j]==t)
     1537      {
     1538        v[i]=j-1;
     1539        if(i==nrows(M)){s=1;}
     1540        break;
     1541      }
     1542    }
     1543  }
     1544  if(s==0){v[nrows(M)]=-j;}
     1545  return(v);
     1546}
     1547//////////////////////////////////////////////////////////////////////
     1548static proc sorting(intvec u,int k)
     1549{
     1550  int s=size(u);
     1551  int i;
     1552  for(i=1;i<=s;i++)
     1553  {
     1554    if(u[i]>=k){break;}
     1555  }
     1556  return(i);
     1557}
     1558//////////////////////////////////////////////////////////////////////
     1559proc MultiplicitySequence(ideal i)
     1560"USAGE:  MultiplicitySequence(i); i ideal
     1561ASSUME:  i is the defining ideal of a (reducible) plane curve singularity.
     1562RETURN:  list X of charts. Each chart contains the multiplicity sequence of
     1563         the corresponding branch.
     1564EXAMPLE: example MultiplicitySequence;  shows an example
     1565"
     1566{
     1567  def s=CurveParam(i);
     1568  setring s;
     1569  int j,k;
     1570  def r1=converter(Param);
     1571  setring r1;
     1572  list Y=hne;
     1573  list X;
     1574  for(j=1;j<=size(Y);j++)
     1575  {
     1576    for(k=1;k<=Y[j][5];k++)
     1577    {
     1578      X=insert(X,multsequence(Y[j]),size(X));
     1579    }
     1580  }
     1581  return(X);
     1582}
     1583example
     1584{
     1585  "EXAMPLE:"; echo = 2;
     1586  ring r=0,(x,y),ds;
     1587  ideal i=x14-x4y7-y11;
     1588  MultiplicitySequence(i);
     1589}
     1590/////////////////////////////////////////////////////////////////////////
     1591proc IntersectionMatrix(ideal i)
     1592"USAGE:  IntersectionMatrix(i); i ideal
     1593ASSUME:  i is the defining ideal of a (reducible) plane curve singularity.
     1594RETURN:  intmat of the intersection multiplicities of the branches.
     1595EXAMPLE: example IntersectionMatrix;  shows an example
     1596"
     1597{
     1598  def s=CurveParam(i);
     1599  setring s;
     1600  int j,k;
     1601  def r1=converter(Param);
     1602  setring r1;
     1603  list Y=hne;
     1604  return(intermat(Y));
     1605}
     1606example
     1607{
     1608  "EXAMPLE:"; echo = 2;
     1609  ring r=0,(x,y),ds;
     1610  ideal i=x14-x4y7-y11;
     1611  IntersectionMatrix(i);
     1612}
     1613///////////////////////////////////////////////////////////////////////////
     1614proc CharacteristicExponents(ideal i)
     1615"USAGE:  CharacteristicExponents(i); i ideal
     1616ASSUME:  i is the defining ideal of a (reducible) plane curve singularity.
     1617RETURN:  list X of charts. Each chart contains the characteristic exponents
     1618         of the corresponding branch.
     1619EXAMPLE: example CharacteristicExponents;  shows an example
     1620"
     1621{
     1622  def s=CurveParam(i);
     1623  setring s;
     1624  int j,k;
     1625  def r1=converter(Param);
     1626  setring r1;
     1627  list X;
     1628  list Y=hne;
     1629  for(j=1;j<=size(Y);j++)
     1630  {
     1631    for(k=1;k<=Y[j][5];k++)
     1632    {
     1633      X=insert(X,invariants(Y[j])[1],size(X));
     1634    }
     1635  }
     1636  return(X);
     1637}
     1638example
     1639{
     1640  "EXAMPLE:"; echo = 2;
     1641  ring r=0,(x,y),ds;
     1642  ideal i=x14-x4y7-y11;
     1643  CharacteristicExponents(i);
     1644}
     1645/////////////////////////////////////////////////////////////////////////////
     1646static proc contactNumber(int a,intvec v1,intvec v2)
     1647{
     1648//====  a is the intersection multiplicity of the branches
     1649//====  v1,v2 are the multiplicity sequences
     1650  int i,c,d;
     1651  if(size(v1)>size(v2))
     1652  {
     1653    for(i=size(v2)+1;i<=size(v1);i++)
     1654    {
     1655      v2[i]=1;
     1656    }
     1657  }
     1658  if(size(v1)<size(v2))
     1659  {
     1660    for(i=size(v1)+1;i<=size(v2);i++)
     1661    {
     1662      v1[i]=1;
     1663    }
     1664  }
     1665  while((c<a)&&(d<size(v1)))
     1666  {
     1667    d++;
     1668    c=c+v1[d]*v2[d];
     1669  }
     1670  if(c<a)
     1671  {
     1672    d=d+a-c;
     1673  }
     1674   return(d);
     1675}
     1676////////////////////////////////////////////////////////////////////////////
     1677proc ContactMatrix(ideal I)
     1678"USAGE:  ContactMatrix(I); I ideal
     1679ASSUME:  i is the defining ideal of a (reducible) plane curve singularity.
     1680RETURN:  intmat N of the contact matrix of the braches of the curve.
     1681EXAMPLE: example ContactMatrix;  shows an example
     1682"
     1683{
     1684  def s=CurveParam(I);
     1685  setring s;
     1686  int j,k,i;
     1687  def r1=converter(Param);
     1688  setring r1;
     1689  list Y=hne;
     1690  list L;
     1691  for(j=1;j<=size(Y);j++)
     1692  {
     1693    for(k=1;k<=Y[j][5];k++)
     1694    {
     1695      L=insert(L,multsequence(Y[j]),size(L));
     1696     }
     1697  }
     1698  intmat M=intermat(Y);
     1699  intmat N[nrows(M)][ncols(M)];
     1700  for(i=1;i<=nrows(M);i++)
     1701  {
     1702  for(j=i+1;j<=ncols(M);j++)
     1703  {
     1704    N[i,j]=contactNumber(M[i,j],L[i],L[j]);
     1705    N[j,i]=N[i,j];}
     1706  }
     1707  return(N);
     1708}
     1709example
     1710{
     1711  "EXAMPLE:"; echo = 2;
     1712  ring r=0,(x,y),ds;
     1713  ideal i=x14-x4y7-y11;
     1714  ContactMatrix(i);
     1715}
     1716///////////////////////////////////////////////////////////////////////////
     1717proc plainInvariants(ideal i)
     1718"USAGE:  plainInvariants(i); i ideal
     1719ASSUME:  i is the defining ideal of a (reducible) plane curve singularity.
     1720RETURN:  list L of charts. L[j] is the invariants of the jth branch and the last entry
     1721         of L is a list containing the intersection matrix,contact matrix,resolution
     1722         graph of the curve.
     1723         L[j][1]:     intvec (characteristic exponents of the jth branch)
     1724         L[j][2]:    intvec (generators of the semigroup of the jth branch)
     1725         L[j][3]:    intvec (first components of the puiseux pairs of the jth branch)
     1726         L[j][4]:    intvec (second components of the puiseux pairs of the jth branch)
     1727         L[j][5]:    int    (degree of conductor of the jth branch)
     1728         L[j][6]:    intvec (multiplicity sequence of the jth branch.
     1729         L[last][1]: intmat (intersection matrix of the branches)
     1730         L[last][2]: intmat (contact matrix of the branches)
     1731         L[last][3]: intmat (resolution graph of the curve)
     1732SEE ALSO: MultiplicitySequence, CharacteristicExponents, IntersectionMatrix,
     1733          ContactMatrix
     1734EXAMPLE: example plainInvariants;  shows an example
     1735"
     1736{
     1737  def s=CurveParam(i);
     1738  setring s;
     1739  int j,k;
     1740  def r1=converter(Param);
     1741  setring r1;
     1742  list Y=hne;
     1743  list L,X,Q;
     1744  for(j=1;j<=size(Y);j++)
     1745  {
     1746    for(k=1;k<=Y[j][5];k++)
     1747    {
     1748      L=insert(L,invariants(Y[j]),size(L));   //computes the same thing again
     1749      X=insert(X,invariants(Y[j])[1],size(X));
     1750    }
     1751  }
     1752  L=insert(L,list(),size(L));
     1753  L[size(L)]=insert(L[size(L)],intermat(Y),size(L[size(L)]));
     1754  intmat N[nrows(intermat(Y))][ncols(intermat(Y))];
     1755  for(k=1;k<=nrows(intermat(Y));k++)
     1756  {
     1757    for(j=k+1;j<=ncols(intermat(Y));j++)
     1758    {
     1759      N[k,j]=contactNumber(intermat(Y)[k,j],L[k][6],L[j][6]);
     1760      N[j,k]=N[k,j];
     1761    }
     1762  }
     1763  L[size(L)]=insert(L[size(L)],N,size(L[size(L)]));
     1764  Q=L[size(L)][2],X;
     1765  L[size(L)]=insert(L[size(L)],resolutiongraph(Q),size(L[size(L)]));
     1766  return(L);
     1767}
     1768example
     1769{
     1770  "EXAMPLE:"; echo = 2;
     1771  ring r=0,(x,y),ds;
     1772  ideal i=x14-x4y7-y11;
     1773  plainInvariants(i);
     1774}
     1775////////////////////////////////////////////////////////////////////////////////////
     1776static proc bound(poly x,poly y)
     1777{
     1778  poly z=x+y;
     1779  int m=ord(z);
     1780  int i;
     1781  int c=-1;
     1782  for(i=2;i<=size(z);i++)
     1783  {
     1784    if(gcd(m,leadexp(z[i])[1])==1)
     1785    {
     1786      c=2*leadexp(z[i])[1];
     1787      break;
     1788    }
     1789    else
     1790    {
     1791      m=gcd(m,leadexp(z[i])[1]);
     1792    }
     1793  }
     1794  return(c);
     1795}
     1796/////////////////////////////////////////////////////////////////////////////////////
     1797static proc boundparam(poly f)
     1798{
     1799  int i;
     1800  int l=size(f);
     1801  int m=leadexp(f[l])[1];
     1802  for(i=l-1;i>=1;i--)
     1803  {
     1804    if(gcd(m,leadexp(f[i])[1])==1)
     1805    {
     1806      i=i-1;
     1807      break;
     1808    }
     1809    else
     1810    {
     1811      m=gcd(m,leadexp(f[i])[1]);
     1812    }
     1813  }
     1814  return(2*leadexp(f[i+1])[1]);
     1815}

  • Singular/LIB/decodegb.lib

    r8f296a r1e1ec4  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: decodegb.lib 15103 2012-07-11 10:00:13Z motsak $";
     2version="$Id$";
    33category="Coding theory";
    44info="

  • Singular/LIB/deform.lib

    r8f296a r1e1ec4  
    880880  int @c = ncols(A);
    881881  int i,j;
    882   string ord_str = "wp("+string(w_vec)+")";
    883882  def br = basering;
    884   def nr=changeord(ord_str);
     883  def nr=changeord(list(list("wp",w_vec)));
    885884  setring nr;
    886885  matrix A    = imap(br,A);
    … …  
    946945   A  = transpose(A);
    947946   def br = basering;
    948    string o_str = "wp("+string(d_vec)+")";
    949    def nr=changeord(o_str);
     947   def nr=changeord(list(list("wp",d_vec)));
    950948   setring nr;
    951949   module A = fetch(br,A);

  • Singular/LIB/dmod.lib

    r8f296a r1e1ec4  
    13661366  dbprint(ppl-1, @R4);
    13671367  ideal K4 = imap(@R,K2);
     1368  intvec saveopt=option(get);
    13681369  option(redSB);
    13691370  dbprint(ppl,"// -4-2- the final cosmetic std");
    … …  
    13771378  ideal LD = K4;
    13781379  export LD;
     1380  option(set,saveopt);
    13791381  return(@R4);
    13801382}
    … …  
    20312033  dbprint(ppl,"// -4-3- an operator PS found, PS*f^(s+1) = b(s)*f^s");
    20322034  dbprint(ppl-1,PS);
     2035  intvec saveopt=option(get);
    20332036  option(redSB);
    20342037  dbprint(ppl,"// -4-4- the final cosmetic std");
    … …  
    20362039  LD  = engine(LD,eng);
    20372040  export F,LD,LD0,bs,BS,PS;
     2041  option(set,saveopt);
    20382042  return(@R4);
    20392043}
    … …  
    23992403  ideal K = imap(@R,K);
    24002404  kill @R;
     2405  intvec saveopt=option(get);
    24012406  option(redSB);
    24022407  dbprint(ppl,"// -3-2- the final cosmetic std");
    … …  
    24122417  export Param;
    24132418  kill sParam;
     2419  option(set,saveopt);
    24142420  return(@R2);
    24152421}
    … …  
    24412447
    24422448proc annfsBMI(ideal F, list #)
    2443 "USAGE:  annfsBMI(F [,eng,met]);  F an ideal, eng, met optional ints
     2449"USAGE:  annfsBMI(F [,eng,met,us]);  F an ideal, eng, met, us optional ints
    24442450RETURN:  ring
    24452451PURPOSE: compute two kinds of Bernstein-Sato ideals, associated to
    … …  
    24502456@*       If eng <>0, @code{std} is used for Groebner basis computations,
    24512457@*       otherwise, and by default @code{slimgb} is used.
    2452 @*       If met <0, the B-Sigma ideal (cf. Castro and Ucha, 
    2453 @*       'On the computation of Bernstein-Sato ideals', 2005) is computed. 
    2454 @*       If 0 < met < P, then the ideal B_P (cf. the paper) is computed. 
     2458@*       If met <0, the B-Sigma ideal (cf. Castro and Ucha,
     2459@*       'On the computation of Bernstein-Sato ideals', 2005) is computed.
     2460@*       If 0 < met < P, then the ideal B_P (cf. the paper) is computed.
    24552461@*       Otherwise, and by default, the ideal B (cf. the paper) is computed.
     2462@*       If us<>0, then syzygies-driven method is used.
    24562463@*       If the output ideal happens to be principal, the list of factors
    24572464@*       with their multiplicities is returned instead of the ideal.
    … …  
    24632470  int eng = 0;
    24642471  int met = 0;
     2472  int usesyz = 0;
    24652473  if ( size(#)>0 )
    24662474  {
    … …  
    24732481      if ( typeof(#[2]) == "int" )
    24742482      {
    2475         met = int(#[2]);
     2483        met = int(#[2]);
    24762484      }
    24772485    }
     2486    if ( size(#)>2 )
     2487    {
     2488      if ( typeof(#[3]) == "int" )
     2489      {
     2490        usesyz = int(#[3]);
     2491      }
     2492    }
     2493
    24782494  }
    24792495  // returns a list with a ring and an ideal LD in it
    … …  
    24832499  int N = nvars(basering);
    24842500  //if F has some generators which are zero, int P = ncols(I);
    2485   int P = size(F); 
    2486   if (P < ncols(F)) 
     2501  int P = size(F);
     2502  if (P < ncols(F))
    24872503  {
    24882504    F = simplify(F,2);
    … …  
    25962612  }
    25972613  // -------- the ideal I is ready ----------
     2614  if (usesyz)
     2615  {
     2616    dbprint(ppl,"// -1-1-1 adding syzygy-driven elements to the ideal");
     2617    // -- form the extended Jacobian matrix; do the comps over K[x]
     2618    setring save;
     2619    module JC = module(transpose(jacob(F))); // NxP
     2620    for (j=1; j<=P; j++)
     2621    {
     2622      JC[j] = JC[j] + F[j]*gen(N+j);
     2623    }
     2624    //    matrix JAC[N+P][P];
     2625    dbprint(ppl,"// -1-1-2 computing syzygies of the extended Jacobian matrix over K[_x]");
     2626    dbprint(ppl-1, matrix(JC));
     2627    matrix SJ = transpose(syz(transpose(JC))); //or of std(syz)?
     2628    dbprint(ppl,"// -1-1-3 finished computing syzygies of the ext Jacobian matrix over K[_x]");
     2629    dbprint(ppl-1, matrix(SJ));
     2630    setring @R;
     2631    // add generators: first N comps d_i, then P comps s_j
     2632    matrix SJ = imap(save, SJ);
     2633    poly pi;
     2634    for (i=1; i<=nrows(SJ); i++)
     2635    {
     2636      pi = 0;
     2637      for (j=1; j<=N; j++)
     2638      {
     2639        pi = pi + SJ[i,j]*var(2*P+N+j);
     2640      }
     2641      for (j=1; j<=P; j++)
     2642      {
     2643        pi = pi + SJ[i,N+j]*s(j);
     2644      }
     2645      dbprint(ppl-1, "// -1-1-4 adding element:" + string(pi));
     2646      I = I, pi;
     2647    }
     2648  }
    25982649  dbprint(ppl,"// -1-2- starting the elimination of "+string(t(1..P))+" in @R");
    25992650  dbprint(ppl-1, I);
    … …  
    26812732  if ( (met>0) && (met<=ncols(F) ) )
    26822733  {
    2683         /* compute the ideal B_met */
     2734  /* compute the ideal B_met */
    26842735  //j=2;         // for example
    26852736  //K = K,F[j];  // to compute Bj (see "On the computation of Bernstein-Sato ideals"; Castro, Ucha)
    … …  
    26892740  if ( ( met == 0) || (met > ncols(F) ) )
    26902741  {
    2691     // that is met ==0 or met> ncols(F)   
     2742    // that is met ==0 or met> ncols(F)
    26922743    /* compute the ideal B */
    26932744    poly f=1;
    … …  
    27872838  dbprint(ppl-1, @R4);
    27882839  ideal K4 = imap(@R,K);
     2840  intvec saveopt=option(get);
    27892841  option(redSB);
    27902842  dbprint(ppl,"// -4-2- the final cosmetic std");
    … …  
    27982850  ideal LD = K4;
    27992851  export LD;
     2852  option(set,saveopt);
    28002853  return(@R4);
    28012854}
    … …  
    28152868  print(matrix(lead(LD))); // compact form of leading
    28162869  //  monomials from the annihilator
    2817   BS; // Bernstein-Sato B-Sigma ideal: it is principal, 
     2870  BS; // Bernstein-Sato B-Sigma ideal: it is principal,
    28182871  // so factors and multiplicities are returned
    2819   // *3* and now, let us compute B-i ideal 
    2820   setring r; 
     2872  // *3* and now, let us compute B-i ideal
     2873  setring r;
    28212874  def Bi = annfsBMI(F,0,3); // that is F[3]=x+y is taken
    28222875  setring Bi;
    … …  
    31813234  dbprint(ppl-1, @R5);
    31823235  ideal K5 = imap(@R2,K3);
     3236  intvec saveopt=option(get);
    31833237  option(redSB);
    31843238  dbprint(ppl,"// -5-2- the final cosmetic std");
    … …  
    31943248  kill @R4;
    31953249  export LD;
     3250  option(set,saveopt);
    31963251  return(@R5);
    31973252}
    … …  
    48604915  dbprint(ppl-1, @R4);
    48614916  ideal K4 = imap(@R2,K2);
     4917  intvec saveopt=option(get);
    48624918  option(redSB);
    48634919  dbprint(ppl,"// -3-2- the final cosmetic std");
    … …  
    48704926  ideal LD = K4;
    48714927  export LD;
     4928  option(set,saveopt);
    48724929  return(@R4);
    48734930}

  • Singular/LIB/dmodapp.lib

    r8f296a r1e1ec4  
    134134  17.03.11 by DA:
    135135  - bugfixes for inForm with polynomial input, typo in restrictionIdealEngine
     136
     137  06.06.12 by DA:
     138  - bugfix and documentation in deRhamCohomIdeal, output and
     139  documentation in deRhamCohom
     140  - changed charVariety: no homogenization is needed
     141  - changed inForm: code is much simpler using jet
     142
    136143*/
    137144
    … …  
    256263
    257264proc engine(def I, int i)
    258   "USAGE:  engine(I,i);  I  ideal/module/matrix, i an int
     265"USAGE:  engine(I,i);  I  ideal/module/matrix, i an int
    259266RETURN:  the same type as I
    260267PURPOSE: compute the Groebner basis of I with the algorithm, chosen via i
    … …  
    291298
    292299proc poly2list (poly f)
    293   "USAGE:  poly2list(f); f a poly
     300"USAGE:  poly2list(f); f a poly
    294301RETURN:  list of exponents and corresponding terms of f
    295302PURPOSE: converts a poly to a list of pairs consisting of intvecs (1st entry)
    … …  
    330337
    331338proc fl2poly(list L, string s)
    332   "USAGE:  fl2poly(L,s); L a list, s a string
     339"USAGE:  fl2poly(L,s); L a list, s a string
    333340RETURN:  poly
    334341PURPOSE: reconstruct a monic polynomial in one variable from its factorization
    … …  
    379386
    380387proc insertGenerator (list #)
    381   "USAGE:  insertGenerator(id,p[,k]);
     388"USAGE:  insertGenerator(id,p[,k]);
    382389@*       id an ideal/module, p a poly/vector, k an optional int
    383390RETURN:  of the same type as id
    … …  
    460467
    461468proc deleteGenerator (def id, int k)
    462   "USAGE:   deleteGenerator(id,k);  id an ideal/module, k an int
     469"USAGE:   deleteGenerator(id,k);  id an ideal/module, k an int
    463470RETURN:   of the same type as id
    464471PURPOSE:  deletes the k-th generator from the first argument and returns
    … …  
    534541
    535542proc bFactor (poly F)
    536   "USAGE:  bFactor(f);  f poly
     543"USAGE:  bFactor(f);  f poly
    537544RETURN:  list of ideal and intvec and possibly a string
    538545PURPOSE: tries to compute the roots of a univariate poly f
    … …  
    626633
    627634proc isInt (number n)
    628   "USAGE:  isInt(n); n a number
     635"USAGE:  isInt(n); n a number
    629636RETURN:  int, 1 if n is an integer or 0 otherwise
    630637PURPOSE: check whether given object of type number is actually an int
    … …  
    654661
    655662proc intRoots (list l)
    656   "USAGE:  isInt(L); L a list
     663"USAGE:  isInt(L); L a list
    657664RETURN:  list
    658665PURPOSE: extracts integer roots from a list given in @code{bFactor} format
    … …  
    713720
    714721proc sortIntvec (intvec v)
    715   "USAGE:  sortIntvec(v); v an intvec
     722"USAGE:  sortIntvec(v); v an intvec
    716723RETURN:  list of two intvecs
    717724PURPOSE: sorts an intvec
    … …  
    795802
    796803proc isFsat(ideal I, poly F)
    797   "USAGE:  isFsat(I, F);  I an ideal, F a poly
     804"USAGE:  isFsat(I, F);  I an ideal, F a poly
    798805RETURN:  int, 1  if I is F-saturated and 0 otherwise
    799806PURPOSE: checks whether the ideal I is F-saturated
    … …  
    834841
    835842proc annRat(poly g, poly f)
    836   "USAGE:   annRat(g,f);  f, g polynomials
     843"USAGE:   annRat(g,f);  f, g polynomials
    837844RETURN:   ring (a Weyl algebra) containing an ideal 'LD'
    838845PURPOSE:  compute the annihilator of the rational function g/f in the
    … …  
    960967
    961968proc annPoly(poly f)
    962   "USAGE:   annPoly(f);  f a poly
     969"USAGE:   annPoly(f);  f a poly
    963970RETURN:   ring (a Weyl algebra) containing an ideal 'LD'
    964971PURPOSE:  compute the complete annihilator ideal of f in the corresponding
    … …  
    10271034
    10281035proc DLoc(ideal I, poly F)
    1029   "USAGE:  DLoc(I, f);  I an ideal, f a poly
     1036"USAGE:  DLoc(I, f);  I an ideal, f a poly
    10301037RETURN:  list of ideal and list
    10311038ASSUME:  the basering is a Weyl algebra
    … …  
    10811088
    10821089proc DLoc0(ideal I, poly F)
    1083   "USAGE:  DLoc0(I, f);  I an ideal, f a poly
     1090"USAGE:  DLoc0(I, f);  I an ideal, f a poly
    10841091RETURN:  ring (a Weyl algebra) containing an ideal 'LD0' and a list 'BS'
    10851092PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s,
    … …  
    13231330
    13241331proc SDLoc(ideal I, poly F)
    1325   "USAGE:  SDLoc(I, f);  I an ideal, f a poly
     1332"USAGE:  SDLoc(I, f);  I an ideal, f a poly
    13261333RETURN:  ring (basering extended by a new variable) containing an ideal 'LD'
    13271334PURPOSE: compute a generic presentation of the localization of D/I w.r.t. f^s
    … …  
    15411548
    15421549proc GBWeight (ideal I, intvec u, intvec v, list #)
    1543   "USAGE:  GBWeight(I,u,v [,s,t,w]);
     1550"USAGE:  GBWeight(I,u,v [,s,t,w]);
    15441551@*       I ideal, u,v intvecs, s,t optional ints, w an optional intvec
    15451552RETURN:  ideal, Groebner basis of I w.r.t. the weights u and v
    … …  
    16961703  }
    16971704  ideal II = I;
    1698   if (simplify(II,2) == 0)
    1699   {
    1700     return(I);
    1701   }
    1702   int j,i;
    1703   bigint s,m;
    1704   list l;
     1705  int j;
    17051706  poly g;
    17061707  ideal J;
    17071708  for (j=1; j<=ncols(II); j++)
    17081709  {
    1709     l = poly2list(II[j]);
    1710     m = scalarProd(w,l[1][1]);
    1711     g = l[1][2];
    1712     for (i=2; i<=size(l); i++)
    1713     {
    1714       s = scalarProd(w,l[i][1]);
    1715       if (s == m)
    1716       {
    1717         g = g + l[i][2];
    1718       }
    1719       else
    1720       {
    1721         if (s > m)
    1722         {
    1723           m = s;
    1724           g = l[i][2];
    1725         }
    1726       }
    1727     }
    1728     J[j] = g;
     1710    g = II[j];
     1711    J[j] = g - jet(g,deg(g,w)-1,w);
    17291712  }
    17301713  if (inp1 == "ideal")
    … …  
    17551738
    17561739proc initialIdealW(ideal I, intvec u, intvec v, list #)
    1757   "USAGE:  initialIdealW(I,u,v [,s,t,w]);
     1740"USAGE:  initialIdealW(I,u,v [,s,t,w]);
    17581741@*       I ideal, u,v intvecs, s,t optional ints, w an optional intvec
    17591742RETURN:  ideal, GB of initial ideal of the input ideal wrt the weights u and v
    … …  
    18141797
    18151798proc initialMalgrange (poly f,list #)
    1816   "USAGE:  initialMalgrange(f,[,a,b,v]); f poly, a,b optional ints, v opt. intvec
     1799"USAGE:  initialMalgrange(f,[,a,b,v]); f poly, a,b optional ints, v opt. intvec
    18171800RETURN:  ring, Weyl algebra induced by basering, extended by two new vars t,Dt
    18181801PURPOSE: computes the initial Malgrange ideal of a given polynomial w.r.t. the
    … …  
    20262009    vector v = pIntersect(s,inG);
    20272010    list L = bFactor(vec2poly(v));
    2028     dbprint(ppl,"// found b-function of given ideal wrt given weight");
     2011    dbprint(ppl,"// found b-function of given ideal wrt weight " + string(w));
    20292012    dbprint(ppl-1,"// roots: "+string(L[1]));
    20302013    dbprint(ppl-1,"// multiplicities: "+string(L[2]));
    … …  
    21802163
    21812164proc restrictionModule (ideal I, intvec w, list #)
    2182   "USAGE:  restrictionModule(I,w,[,eng,m,G]);
     2165"USAGE:  restrictionModule(I,w,[,eng,m,G]);
    21832166@*       I ideal, w intvec, eng and m optional ints, G optional ideal
    21842167RETURN:  ring (a Weyl algebra) containing a module 'resMod'
    … …  
    23022285
    23032286proc restrictionIdeal (ideal I, intvec w, list #)
    2304   "USAGE:  restrictionIdeal(I,w,[,eng,m,G]);
     2287"USAGE:  restrictionIdeal(I,w,[,eng,m,G]);
    23052288@*       I ideal, w intvec, eng and m optional ints, G optional ideal
    23062289RETURN:  ring (a Weyl algebra) containing an ideal 'resIdeal'
    … …  
    23542337
    23552338proc fourier (ideal I, list #)
    2356   "USAGE:  fourier(I[,v]); I an ideal, v an optional intvec
     2339"USAGE:  fourier(I[,v]); I an ideal, v an optional intvec
    23572340RETURN:  ideal
    23582341PURPOSE: computes the Fourier transform of an ideal in a Weyl algebra
    … …  
    24232406
    24242407proc inverseFourier (ideal I, list #)
    2425   "USAGE:  inverseFourier(I[,v]); I an ideal, v an optional intvec
     2408"USAGE:  inverseFourier(I[,v]); I an ideal, v an optional intvec
    24262409RETURN:  ideal
    24272410PURPOSE: computes the inverse Fourier transform of an ideal in a Weyl algebra
    … …  
    24922475
    24932476proc integralModule (ideal I, intvec w, list #)
    2494   "USAGE:  integralModule(I,w,[,eng,m,G]);
     2477"USAGE:  integralModule(I,w,[,eng,m,G]);
    24952478@*       I ideal, w intvec, eng and m optional ints, G optional ideal
    24962479RETURN:  ring (a Weyl algebra) containing a module 'intMod'
    … …  
    26062589
    26072590proc integralIdeal (ideal I, intvec w, list #)
    2608   "USAGE:  integralIdeal(I,w,[,eng,m,G]);
     2591"USAGE:  integralIdeal(I,w,[,eng,m,G]);
    26092592@*       I ideal, w intvec, eng and m optional ints, G optional ideal
    26102593RETURN:  ring (a Weyl algebra) containing an ideal 'intIdeal'
    … …  
    26532636
    26542637proc deRhamCohomIdeal (ideal I, list #)
    2655   "USAGE:  deRhamCohomIdeal (I[,w,eng,k,G]);
     2638"USAGE:  deRhamCohomIdeal (I[,w,eng,k,G]);
    26562639@*       I ideal, w optional intvec, eng and k optional ints, G optional ideal
    26572640RETURN:  ideal
    … …  
    26682651PURPOSE: computes a basis of the n-th de Rham cohomology group of the complement
    26692652@*       of the hypersurface defined by f
    2670 NOTE:    If I does not satisfy the assumptions described above, the result might
     2653NOTE:    The elements of the basis are of the form f^m*p, where p runs over the
     2654@*       entries of the returned ideal.
     2655@*       If I does not satisfy the assumptions described above, the result might
    26712656@*       have no meaning. Note that I can be computed with @code{annfs}.
    26722657@*       If w is an intvec with exactly n strictly positive entries, w is used
    … …  
    27932778      DR[size(DR)+1] = B[i]*Dt;
    27942779      j=1;
    2795       while (p<N[j])
     2780      while ((j<size(N)) && (p<N[j]))
    27962781      {
    27972782        j++;
    … …  
    28162801
    28172802proc deRhamCohom (poly f, list #)
    2818   "USAGE:  deRhamCohom(f[,eng,m]);  f poly, eng and m optional ints
    2819 RETURN:  ring (a Weyl Algebra) containing an ideal 'DR'
    2820 ASSUME:  Basering is a commutative and over a field of characteristic 0.
     2803"USAGE:  deRhamCohom(f[,w,eng,m]);  f poly, w optional intvec,
     2804                                    eng and m optional ints
     2805RETURN:  ring (a Weyl Algebra) containing a list 'DR' of ideal and int
     2806ASSUME:  Basering is commutative and over a field of characteristic 0.
    28212807PURPOSE: computes a basis of the n-th de Rham cohomology group of the complement
    28222808@*       of the hypersurface defined by f, where n denotes the number of
    28232809@*       variables of the basering
    2824 NOTE:    The output ring is the n-th Weyl algebra. It contains an ideal 'DR'
    2825 @*       being a basis of the n-th de Rham cohomology group of the complement of
    2826 @*       the hypersurface defined by f.
     2810NOTE:    The output ring is the n-th Weyl algebra. It contains a list 'DR' with
     2811@*       two entries (ideal J and int m) such that {f^m*J[i] : i=1..size(I)} is
     2812@*       a basis of the n-th de Rham cohomology group of the complement of the
     2813@*       hypersurface defined by f.
     2814@*       If w is an intvec with exactly n strictly positive entries, w is used
     2815@*       in the computation. Otherwise, and by default, w is set to (1,...,1).
    28272816@*       If eng<>0, @code{std} is used for Groebner basis computations,
    28282817@*       otherwise, and by default, @code{slimgb} is used.
    … …  
    28392828{
    28402829  int ppl = printlevel - voice + 2;
     2830  def save = basering;
     2831  int n = nvars(save);
     2832  intvec w = 1:n;
    28412833  int eng,l0,l0given;
    28422834  if (size(#)>0)
    28432835  {
    2844     if(intLike(#[1]))
    2845     {
    2846       eng = int(#[1]);
     2836    if (typeof(#[1])=="intvec")
     2837    {
     2838      w = #[1];
    28472839    }
    28482840    if (size(#)>1)
    … …  
    28502842      if(intLike(#[2]))
    28512843      {
    2852         l0 = int(#[2]);
    2853         l0given = 1;
     2844        eng = int(#[2]);
    28542845      }
     2846      if (size(#)>2)
     2847      {
     2848        if(intLike(#[3]))
     2849        {
     2850          l0 = int(#[3]);
     2851          l0given = 1;
     2852        }
     2853      }
    28552854    }
    28562855  }
    … …  
    28602859  }
    28612860  int i;
    2862   def save = basering;
    2863   int n = nvars(save);
    28642861  dbprint(ppl,"// Computing s-parametric annihilator Ann(f^s)...");
    28652862  def A = Sannfs(f);
    … …  
    29082905  dbprint(ppl-1,"//    Got: " + string(LD));
    29092906  kill A;
    2910   intvec w = 1:n;
    2911   ideal DR = deRhamCohomIdeal(LD,w,eng);
     2907  ideal DRJ = deRhamCohomIdeal(LD,w,eng);
     2908  list DR = DRJ,l0;
    29122909  export(DR);
    29132910  setring save;
    … …  
    29272924
    29282925proc appelF1()
    2929   "USAGE:  appelF1();
     2926"USAGE:  appelF1();
    29302927RETURN:  ring (a parametric Weyl algebra) containing an ideal 'IAppel1'
    29312928PURPOSE: defines the ideal in a parametric Weyl algebra,
    … …  
    29582955
    29592956proc appelF2()
    2960   "USAGE:  appelF2();
     2957"USAGE:  appelF2();
    29612958RETURN:  ring (a parametric Weyl algebra) containing an ideal 'IAppel2'
    29622959PURPOSE: defines the ideal in a parametric Weyl algebra,
    … …  
    29882985
    29892986proc appelF4()
    2990   "USAGE:  appelF4();
     2987"USAGE:  appelF4();
    29912988RETURN:  ring (a parametric Weyl algebra) containing an ideal 'IAppel4'
    29922989PURPOSE: defines the ideal in a parametric Weyl algebra,
    … …  
    30213018
    30223019proc charVariety(ideal I, list #)
    3023   "USAGE:  charVariety(I [,eng]); I an ideal, eng an optional int
     3020"USAGE:  charVariety(I [,eng]); I an ideal, eng an optional int
    30243021RETURN:  ring (commutative) containing an ideal 'charVar'
    30253022PURPOSE: computes an ideal whose zero set is the characteristic variety of I in
    … …  
    30513048  def save = basering;
    30523049  int n = nvars(save) div 2;
    3053   intvec u = 0:n;
    3054   intvec v = 1:n;
    3055   dbprint(ppl,"// Computing Groebner basis wrt weights...");
    3056   I = GBWeight(I,u,v,eng);
    3057   dbprint(ppl,"// ...finished");
     3050  intvec uv = (0:n),(1:n);
     3051  list RL = ringlist(save);
     3052  list L = RL[3];
     3053  L = insert(L,list("a",uv));
     3054  RL[3] = L;
     3055  // TODO printlevel
     3056  def Ra = ring(RL);
     3057  setring Ra;
     3058  dbprint(ppl,"// Starting Groebner basis computation...");
     3059  ideal I = imap(save,I);
     3060  I = engine(I,eng);
     3061  dbprint(ppl,"// ... done.");
    30583062  dbprint(ppl-1,"//    Got: " + string(I));
    3059   list RL = ringlist(save);
     3063  setring save;
     3064  RL = ringlist(save);
    30603065  RL = RL[1..4];
    30613066  def newR = ring(RL);
    30623067  setring newR;
    3063   ideal charVar = imap(save,I);
    3064   intvec uv = u,v;
     3068  ideal charVar = imap(Ra,I);
    30653069  charVar = inForm(charVar,uv);
    3066   charVar = groebner(charVar);
     3070  // charVar = groebner(charVar);
    30673071  export(charVar);
    30683072  setring save;
    … …  
    30813085  setring CA; CA; // commutative ring
    30823086  charVar;
    3083   dim(charVar);   // hence I is holonomic
     3087  dim(std(charVar));   // hence I is holonomic
    30843088}
    30853089
    30863090proc charInfo(ideal I)
    3087   "USAGE:  charInfo(I);  I an ideal
     3091"USAGE:  charInfo(I);  I an ideal
    30883092RETURN:  ring (commut.) containing ideals 'charVar','singLoc' and list 'primDec'
    30893093PURPOSE: computes characteristic variety of I (in the sense of D-module theory),

  • Singular/LIB/ehv.lib

    r8f296a r1e1ec4  
    137137{
    138138  if(printlevel > 2){"Entering equiMaxEHV.";}
    139   if(ord_test(basering)!=1)
     139  if(attrib(basering,"global")==0)
    140140    {
    141141      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    188188EXAMPLE: example removeComponent; shows an example"
    189189{
    190   if(ord_test(basering)!=1)
     190  if(attrib(basering,"global")==0)
    191191    {
    192192      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    240240EXAMPLE: example AssOfDim; shows an example"
    241241{
    242   if(ord_test(basering)!=1)
     242  if(attrib(basering,"global")==0)
    243243    {
    244244      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    335335{
    336336  if(printlevel > 2){"Entering equiRadEHV.";}
    337   if(ord_test(basering)!=1)
     337  if(attrib(basering,"global")==0)
    338338    {
    339339      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    546546{
    547547  if(printlevel > 2){"Entering radEHV.";}
    548   if(ord_test(basering)!=1)
     548  if(attrib(basering,"global")==0)
    549549    {
    550550      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    589589EXAMPLE: example IntAssOfDim1; shows an example"
    590590{
    591   if(ord_test(basering)!=1)
     591  if(attrib(basering,"global")==0)
    592592    {
    593593      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    640640EXAMPLE: example IntAssOfDim2; shows an example"
    641641{
    642   if(ord_test(basering)!=1)
     642  if(attrib(basering,"global")==0)
    643643    {
    644644      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    674674{
    675675  if(printlevel > 2){"Entering decompEHV.";}
    676   if(ord_test(basering)!=1)
     676  if(attrib(basering,"global")==0)
    677677    {
    678678      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    799799{
    800800  if(printlevel > 2){"Entering idempotent.";}
    801   if(ord_test(basering)!=1)
     801  if(attrib(basering,"global")==0)
    802802    {
    803803      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    882882{
    883883  if(printlevel > 2){"Entering equiAssEHV.";}
    884   if(ord_test(basering)!=1)
     884  if(attrib(basering,"global")==0)
    885885    {
    886886      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    956956EXAMPLE: example AssEHV; shows an example"
    957957{
    958   if(ord_test(basering)!=1)
     958  if(attrib(basering,"global")==0)
    959959    {
    960960      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    10481048              //...homogenize K w.r.t. t,...
    10491049              if(printlevel > 2){"Input not homogeneous; must homogenize.";}
    1050               changeord("homoR","dp");
     1050              def homoR=changeord(list(list("dp",1:nvars(basering))));
     1051              setring homoR;
    10511052              ideal homoJ = fetch(base,K);
    10521053              homoJ = groebner(homoJ);
    … …  
    11161117EXAMPLE: example minAssEHV; shows an example"
    11171118{
    1118   if(ord_test(basering)!=1)
     1119  if(attrib(basering,"global")==0)
    11191120    {ERROR("// Not implemented for this ordering, please change to global ordering.");}
    11201121
    … …  
    11791180EXAMPLE: example localize; shows an example"
    11801181{
    1181   if(ord_test(basering)!=1)
     1182  if(attrib(basering,"global")==0)
    11821183    {
    11831184      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    12331234{
    12341235  if(printlevel > 2){"Entering componentEHV.";}
    1235   if(ord_test(basering)!=1)
     1236  if(attrib(basering,"global")==0)
    12361237    {
    12371238      ERROR("// Not implemented for this ordering, please change to global ordering.");
    … …  
    13181319{
    13191320  if(printlevel > 2){"Entering primdecEHV.";}
    1320   if(ord_test(basering)!=1)
     1321  if(attrib(basering,"global")==0)
    13211322    {
    13221323      ERROR("// Not implemented for this ordering, please change to global ordering.");

  • Singular/LIB/elim.lib

    r8f296a r1e1ec4  
    406406METHOD: elim uses elimRing to create a ring with an elimination ordering for
    407407        the variables to be eliminated and then applies std if \"std\"
    408         is given, or slimgb if \"slimgb\" is given, or a heuristically choosen
     408        is given, or slimgb if \"slimgb\" is given, or a heuristically chosen
    409409        method.
    410410@*      If the variables in the basering have weights these weights are used

  • Singular/LIB/ffsolve.lib

    r8f296a r1e1ec4  
    11////////////////////////////////////////////////////////////////////
    2 version="$Id: ffsolve.lib f0fdb24 2012-04-23 12:50:26 UTC$";
     2version="$Id$";
    33category="Symbolic-numerical solving";
    44info="
    … …  
    2121LIB "standard.lib";
    2222LIB "matrix.lib";
    23 LIB "arcpoint.lib"; //stattdessen das, wo auch varstr vorkommt
    2423
    2524////////////////////////////////////////////////////////////////////
    … …  
    321320  ideal I, linear_solution, unsolved_part, univar_part, multivar_part, unsolved_vars;
    322321  intvec unsolved_var_nums;
    323   list new_vars;
     322  string new_vars;
    324323  // check assumptions
    325324  if(npars(basering)>1){
    … …  
    389388      number_new_vars = ncols(unsolved_vars);
    390389
    391       list new_vars;
    392       for(int i = number_new_vars; i>=1; i--)
    393       {
    394         new_vars[i]=  "@y("+string(i)+")";
    395       }
     390      new_vars = "@y(1.."+string(number_new_vars)+")";
    396391      def R_new = changevar(new_vars, original_ring);
    397392      setring R_new;
    … …  
    588583
    589584  setring original_ring;
    590 //  string var_str = varstr(original_ring)+",@X,@y";
    591   list l1 = varlist(original_ring);
    592   list l2 = "@X", "@Y";
    593   list var_list = l1+l2;
     585  string var_str = varstr(original_ring)+",@X,@y";
    594586  string minpoly_str = "minpoly="+string(minpoly)+";";
    595   def ring_for_substitution = Ring::changevar(var_list, original_ring);
     587  def ring_for_substitution = Ring::changevar(var_str, original_ring);
    596588
    597589  setring ring_for_substitution;
    … …  
    715707    +get_minpoly_str(size(original_ring),parstr(original_ring,1))+";" ;
    716708  }
    717   list old_vars = varlist(original_ring);
    718   list new_vars;
    719   for(int i = 1; i<=number_of_monomials; i++){
    720     new_vars = insert(new_vars, "@y("+string(i)+")", i-1);
    721   }
    722 
    723   def ring_for_var_change = changevar( old_vars+new_vars, original_ring);
     709  string old_vars = varstr(original_ring);
     710  string new_vars = "@y(1.."+string( number_of_monomials )+")";
     711
     712  def ring_for_var_change = changevar( old_vars+","+new_vars, original_ring);
    724713  setring ring_for_var_change;
    725714  if( prime_field == 0){
    … …  
    730719  ideal I = imap(original_ring, I);
    731720  ideal C;
    732   string weights = "wp(";
    733   for(i=1; i<=nvars(original_ring); i++){
    734     weights = weights+string(1)+",";
    735   }
     721  intvec weights=1:nvars(original_ring);
     722
    736723  for(i=1; i<= number_of_monomials; i++){
    737724    C[i] = monomials[i] - @y(i);
    738     weights = weights+string(deg(monomials[i])+1)+",";
    739   }
    740   weights[size(weights)]=")";
     725    weights = weights,deg(monomials[i])+1;
     726  }
    741727  ideal linear_eqs = I;
    742728  for(i=1; i<=ncols(C); i++){
    … …  
    752738  ideal I = imap(ring_for_var_change, linear_eqs);
    753739  ideal lin_sol = solvelinearpart(I);
    754   string new_ordstr = weights+",C";
    755   def ring_for_back_change = changeord( new_ordstr, ring_for_var_change);
     740  def ring_for_back_change = changeord( list(list("wp",weights),list("C",0:1)), ring_for_var_change);
    756741
    757742  setring ring_for_back_change;
    … …  
    872857    int newvars = nvars(original_ring);
    873858  }
    874   list newvarlist;
    875   for(int i = 1; i<=newvars; i++){
    876     newvarlist = insert(newvarlist, "v("+string(i)+")", i-1);
    877   }
    878   def newring = changevar(newvarlist, original_ring);
     859  string newvarstr = "v(1.."+string(newvars)+")";
     860  def newring = changevar(newvarstr, original_ring);
    879861  setring newring;
    880862  if( prime_field ){

  • Singular/LIB/finvar.lib

    r8f296a r1e1ec4  
    45924592  int dgb=degBound;
    45934593  degBound = 0;
     4594  intvec saveopt=option(get);
    45944595  option(redSB);
    45954596  ideal sP = groebner(ideal(P));
    … …  
    48794880    }
    48804881  }
     4882  option(set,saveopt);
    48814883  return(S,matrix(compress(IS)));
    48824884}
    … …  
    57755777  int dgb=degBound;
    57765778  degBound = 0;
     5779  intvec saveopt=option(get);
    57775780  option(redSB);
    57785781  ideal sP = groebner(ideal(P));
    … …  
    60766079  { kill `ring_name`;
    60776080  }
     6081  option(set,saveopt);
    60786082  return(matrix(S),matrix(compress(IS)));
    60796083}

  • Singular/LIB/gkdim.lib

    r8f296a r1e1ec4  
    77@*         Rabelo, C.,        crabelo@ugr.es
    88
    9 SUPPORT: 'Metodos algebraicos y efectivos en grupos cuanticos', BFM2001-3141, MCYT, Jose Gomez-Torrecillas (Main researcher).
     9Support: 'Metodos algebraicos y efectivos en grupos cuanticos', BFM2001-3141, MCYT, Jose Gomez-Torrecillas (Main researcher).
     10
     11NOTE: The built-in command @code{dim}, executed for a module in @plural, computes the Gelfand-Kirillov dimension.
    1012
    1113PROCEDURES:

  • Singular/LIB/gmspoly.lib

    r8f296a r1e1ec4  
    9494  def @X=basering;
    9595  int n=nvars(@X);
    96   def @WX=changechar("(0,w(1.."+string(n)+"))");
     96  ring ext_ring=0,(w(1..n)),dp;
     97  setring @X;
     98  def @WX=changechar(ringlist(ext_ring));
    9799  setring @WX;
     100  kill ext_ring;
    98101  ideal J=jacob(imap(@X,f));
    99102  int i;

  • Singular/LIB/grobcov.lib

    r8f296a r1e1ec4  
    44info="
    55LIBRARY:  grobcov.lib   Groebner Cover for parametric ideals.
    6 PURPOSE:  Comprehensive Groebner Systems, Groebner Cover, Canonical Forms.
    7           The library contains Montes's algorithms to compute the
     6PURPOSE:  Comprehensive Groebner Systems, Groebner Cover, Canonical Forms,
     7          Parametric Polynomial Systems.
     8          The library contains Montes-Wibmer's algorithms to compute the
    89          canonical Groebner cover of a parametric ideal as described in
    910          the paper:
    1011
    1112          Montes A., Wibmer M.,
    12           Groebner Bases for Polynomial Systems with parameters.
     13          \"Groebner Bases for Polynomial Systems with parameters\".
    1314          Journal of Symbolic Computation 45 (2010) 1391-1425.
    1415
    1516          The central routine is grobcov. Given a parametric
    16           ideal, grobcov outputs its canonical Groebner cover, consisting
     17          ideal, grobcov outputs its Canonical Groebner Cover, consisting
    1718          of a set of pairs of (basis, segment). The basis (after
    1819          normalization) is the reduced Groebner basis for each point
    … …  
    2324          whole parameter space. The output is canonical, it only
    2425          depends on the given parametric ideal and the monomial order.
    25           This is much more than a simple comprehensive Groebner system.
     26          This is much more than a simple Comprehensive Groebner System.
    2627          The algorithm grobcov allows options to solve partially the
    2728          problem when the whole automatic algorithm does not finish
    … …  
    2930
    3031          grobcov uses a first algorithm cgsdr that outputs a disjoint
    31           reduced comprehensive Groebner system with constant lpp.
     32          reduced Comprehensive Groebner System with constant lpp.
     33          For this purpose, in this library, the implemented algorithm is
     34          Kapur-Sun-Wang algorithm, because it is the most efficient
     35          algorithm known for this purpose.
     36
     37          D. Kapur, Y. Sun, and D.K. Wang.
     38          \"A New Algorithm for Computing Comprehensive Groebner Systems\".
     39          Proceedings of ISSAC'2010, ACM Press, (2010), 29-36.
     40
    3241          cgsdr can be called directly if only a disjoint reduced
    33           comprehensive Groebner system is required.
    34 
    35           Two other routines: gencase1 and multigrobcov can be used
    36           in problems with basis of the generic case equal to 1
    37           (for example in automatic geometric theorem discovering)
    38           that allow to obtain partial results even when grobcov does
    39           not finish in reasonable time.
    40 
    41           For completeness, the library also contains the algorithms
    42           with similar purposes contained in the old library redcgs.lib.
    43           These algorithms are, in general, less efficient and do not
    44           ensure a canonical results, even if they are similar to the
    45           results obtained with grobcov.
    46           The old routines are no more recommended and remain in
    47           this library for didactic purposes. These are
    48           cgsdrold, grobcovold, buildtreetoMaple, cantreetoMaple.
     42          Comprehensive Groebner System (CGS) is required.
    4943
    5044AUTHORS:  Antonio Montes , Hans Schoenemann.
    … …  
    5751@*         basering Q[a][x]; (a=parameters, x=variables)
    5852@*         After defining the ring, the main routines
    59 @*         grobcov, cgsdr, gencase1, multigrobcov
     53@*         grobcov, cgsdr,
    6054@*         generate the global rings
    6155@*         @R   (Q[a][x]),
    … …  
    6660@*         create before the above rings by calling setglobalrings();
    6761@*         because most of the internal routines use these rings.
    68 @*         The call to the basic routines grobcov, cgsdr, gencase1, multigrobcov
    69 @*         or even the older grobcovold, cgsdrold will kill these rings.
     62@*         The call to the basic routines grobcov, cgsdr will
     63@*         kill these rings.
    7064
    7165PROCEDURES:
    7266
    73 grobcov(F);          Is the basic routine giving the canonical
    74                      Groebner cover of the parametric ideal F.
    75                      This routine accepts many options, that
    76                      allow to obtain results even when the canonical
    77                      computation does not finish in reasonable time.
    78 
    79 cgsdr(F);            Is the procedure for obtaining a first disjoint,
    80                      reduced comprehensive Groebner system that
    81                      is used in grobcov, but that can be used
    82                      independently if only the CGS is required.
    83                      It is a more efficient version of buildtree
    84                      that does not output the complete discussion tree
    85                      but only the terminal vertices giving the
    86                      disjoint reduced comprehensive Groebner system.
    87 
    88 gencase1(F);         Returns the segment of the generic case when his
    89                      basis is 1. This is useful for automatic discovering
    90                      of geometrical theorems, as it gives the components
    91                      where a solution exists and is much more efficient
    92                      than the complete computation of grobcov.
    93 
    94 multigrobcov(F);     In problems like automatic discovery of theorems,
    95                      when grobcov does not give the answer in reasonable
    96                      time, and the generic case is expected to
    97                      have basis 1, one can try with multigrobcov procedure
    98                      to obtain an answer over the different irreducible
    99                      components: the generic case with basis 1, and the
    100                      components not corresponding to the generic case. To
    101                      deduce from its result the true Groebner cover one
    102                      must discuss theoretically in which segment
    103                      must be located the intersecting parts in the
    104                      different irreducible components.
    105 
    106 setglobalrings();    Generates the global rings @R, @P and @PR that are
    107                      respectively the rings Q[a][x], Q[a], Q[x,a].
    108                      It is called inside each of the fundamental routines of the
    109                      library: grobcov, cgsdr, gencase1, multigrobcov, as well as
    110                      by the old routines cgsdrold, grobcovold and killed
    111                      before the output.
    112                      If the user want to use some other internal routine,
    113                      then setglobalrings() is to be called before, as
    114                      the rings @R, @P and @RP are needed in most of them.
    115                      globally, and more internal routines can be used, but
    116                      These rings are destroyed by the call to any of the basic
    117                      routines.
    118 
    119 pdivi(f,F);          Performs a pseudodivision of a parametric polynomial
    120                      by a parametric ideal.
    121 
    122 pnormalform(f,N,W);  Reduces a parametric polynomial f by a reduced-representation
    123                      (N,W) of null and non-null conditions over the parameters.
    124                      Before using it setglobalrings() must be called.
    125 
    126 Also included from the old library redcgs.lib the following routines
    127 
    128 cgsdrold(F);         Similar to cgsdr using the algorithm buildtree
    129                      of the old library.
    130 grobcovold(F);       Similar to grobcov with the algorithms of the old
    131                      library.
    132 buildtreetoMaple(T); Writes into a file the output of cgsdrold called
    133                      with option ('old',0) into a text file that is Maple
    134                      readable and can be plotted in Maple using
    135                      the tplot routine of the library dpgb.
    136 cantreetoMaple(M);   Writes into a text file the output of grobcovold called
    137                      with  option ('out',1), that is readable
    138                      in Maple and can be plotted using the routine
    139                      plotcantree of the Maple library dpgb.
     67grobcov(F);        Is the basic routine giving the canonical
     68                   Groebner cover of the parametric ideal F.
     69                   This routine accepts many options, that
     70                   allow to obtain results even when the canonical
     71                   computation does not finish in reasonable time.
     72
     73cgsdr(F);          Is the procedure for obtaining a first disjoint,
     74                   reduced Comprehensive Groebner System that
     75                   is used in grobcov, but that can be used
     76                   independently if only the CGS is required.
     77                   It is a more efficient routine than buildtree
     78                   (the own routine that is no more used). It uses
     79                   now KSW algorithm.
     80
     81setglobalrings();  Generates the global rings @R, @P and @PR that are
     82                   respectively the rings Q[a][x], Q[a], Q[x,a].
     83                   It is called inside each of the fundamental routines
     84                   of the library: grobcov, cgsdr and killed before
     85                   the output.
     86                   If the user want to use some other internal routine,
     87                   then setglobalrings() is to be called before, as
     88                   the rings @R, @P and @RP are needed in most of them.
     89                   globally, and more internal routines can be used, but
     90                   these rings are killed by the call to any of the
     91                   basic routines.
     92
     93pdivi(f,F);     Performs a pseudodivision of a parametric polynomial
     94                   by a parametric ideal.
     95
     96pnormalf(f,E,N);   Reduces a parametric polynomial f over V(E) \ V(N)
     97                   E is the null ideal and N the non-null ideal
     98                   over the parameters.
     99
     100extend(GC); When the grobcov of an ideal has been computed
     101                   with the default option ('ext',0) and the explicit
     102                   option ('rep',2) (which is not the default), then
     103                   one can call extend (GC) (and options) to obtain the
     104                   full representation of the bases. With the default
     105                   option ('ext',0) only the generic representation of
     106                   the bases are computed, and one can obtain the full
     107                   representation using extend.
     108
     109locus2d:      Special routine for determining the locus of points
     110                   of a two dimensional object. Given an ideal J with
     111                   two parameters (a,b) and so many variables as
     112                   needed, representing the system determining
     113                   the locus of points (a,b) who verify certain
     114                   geometrical properties, computing the grobcov of
     115                   J and applying to it locus2d, determines the locus.
     116
     117locus2dto:   Transforms the output of locus2d to a string that
     118                   can be reed from different computational systems.
    140119
    141120SEE ALSO: compregb_lib
    … …  
    149128// Library grobcov.lib
    150129// (Groebner cover):
     130// Release 1: (public)
     131// Initial data: 21-1-2008
     132// Final data: 3-7-2008
     133// Release 2: (private)
    151134// Initial data: 6-9-2009
    152 // Release 1:
    153 // Final data: 30-12-2010
    154 // Contains also the old redcgs.lib library that was created
    155 // Initial data: 21-1-2008
    156 // Release 1:
    157 // Final data: 3-7-2008
    158 // Given and determined polynomials and ideals are in the
     135// Final data: 25-10-2011
     136// Release 3: (this release, public)
     137// Initial data: 1-7-2012
     138// Final data: 4-9-2012
    159139// basering Q[a][x];
    160140
    … …  
    167147          defined as global variables.
    168148NOTE:     It is called internally by the fundamental routines of the
    169           library grobcov, cgsdr, gencase1, muligrobcov as well as by the
    170           old ones grobcovold,cgsdrold, and killed before the output.
     149          library grobcov, cgsdr, extend, pdivi, pnormalf, locus2d, locus2dto,
     150          and killed before the output.
    171151          The user does not need to call it, except when it is interested
    172152          in using some internal routine of the library that
    … …  
    177157EXAMPLE:  setglobalrings; shows an example"
    178158{
    179   if (defined(@P)==1)
     159  if (defined(@P))
    180160  {
    181161    kill @P; kill @R; kill @RP;
    … …  
    197177  exportto(Top,@RP);     // global ring K[x,a] with product order
    198178  setring(RR);
    199 }
     179};
    200180example
    201181{ "EXAMPLE:"; echo = 2;
    … …  
    205185  @P;
    206186  @RP;
     187ringlist(R);
     188ringlist(@P);
     189ringlist(@RP);
    207190}
    208191
    … …  
    216199//    ideal Jc (the new form of ideal J without denominators and
    217200//       normalized to content 1)
    218 static proc cld(ideal J)
     201proc cld(ideal J)
    219202{
    220203  if (size(J)==0){return(ideal(0));}
    … …  
    223206  def Ja=imap(RR,J);
    224207  ideal Jb;
    225   if (size(Ja)==0){return(ideal(0));}
     208  if (size(Ja)==0){setring(RR); return(ideal(0));}
    226209  int i;
    227210  def j=0;
    … …  
    230213  def Jc=imap(@RP,Jb);
    231214  return(Jc);
    232 }
    233 
    234 static proc memberpos(f,J)
     215};
     216
     217proc memberpos(f,J)
    235218//"USAGE:  memberpos(f,J);
    236219//         (f,J) expected (polynomial,ideal)
    … …  
    354337//  list L=(7,4,5,1,1,4,9);
    355338//  memberpos(1,L);
    356 //  >
    357339//}
    358340
    359 
    360 static proc subset(J,K)
     341proc subset(J,K)
    361342//"USAGE:   subset(J,K);
    362343//          (J,K)  expected (ideal,ideal)
    … …  
    385366
    386367// elimintfromideal: elimine the constant numbers from an ideal
    387 //     (designed for W, nonnull conditions)
     368//        (designed for W, nonnull conditions)
    388369// input: ideal J
    389 // output:ideal K with the elements of J that are non constants, in the ring @P
    390 static proc elimintfromideal(ideal J)
     370// output:ideal K with the elements of J that are non constants, in the
     371//        ring @P
     372proc elimintfromideal(ideal J)
    391373{
    392374  int i;
    … …  
    401383// input: two coeficients (or terms), that are considered as a quotient
    402384// output: the two coeficients reduced without common factors
    403 static proc simpqcoeffs(poly n,poly m)
     385proc simpqcoeffs(poly n,poly m)
    404386{
    405387  def nc=content(n);
    … …  
    410392}
    411393
    412 // pdivi : pseudodivision of a poly f by an ideal F in a parametric ideal
    413 //         Q[a][x]
     394// pdivi : pseudodivision of a poly f by a parametric ideal F in Q[a][x].
    414395// input:
    415 //   poly f0 (in the parametric ring @R)
    416 //   ideal F0 (in the parametric ring @R)
     396//   poly  f (in the parametric ring @R)
     397//   ideal F (in the parametric ring @R)
    417398// output:
    418399//   list (poly r, ideal q, poly mu)
    … …  
    423404RETURN:   A list (poly r, ideal q, poly m). r is the remainder of the
    424405          pseudodivision, q is the set of quotients, and m is the
    425           factor by which f is to be multiplied.
     406          coefficient factor by which f is to be multiplied.
    426407NOTE:     pseudodivision of a poly f by an ideal F in @R. Returns a
    427408          list (r,q,m) such that m*f=r+sum(q.G), and no lpp of a divisor
    … …  
    430411EXAMPLE:  pdivi; shows an example"
    431412{
     413  int te=0;
     414  if (defined(@P)==1){te=1;}
     415  else{setglobalrings();}
     416  def R=basering;
    432417  int i;
    433418  int j;
    … …  
    436421  def p=f;
    437422  ideal q;
    438   for (i=1; i<=size(F); i++){q[i]=0;}
     423  for (i=1; i<=size(F); i++){q[i]=0;};
    439424  ideal lpf;
    440425  ideal lcf;
    … …  
    478463  }
    479464  list res=r,q,mu;
     465  if(te==0){kill @P; kill @R; kill @RP;}
    480466  return(res);
    481467}
    … …  
    497483// @R
    498484// input:
    499 //   poly f  (given in the ring @R)
     485//   poly f (given in the ring @R)
    500486//   poly g (given in the ring @R)
    501487// output:
    502488//   list (S, red):  S is the S-poly(f,g) and red is a Boolean variable
    503 //                if red==1 then S reduces by Buchberger 1st criterion (not used)
    504 static proc pspol(poly f,poly g)
     489//                if red then S reduces by Buchberger 1st criterion
     490//                (not used)
     491proc pspol(poly f,poly g)
    505492{
    506493  def lcf=leadcoef(f);
    … …  
    523510//         Operates in the ring @P, but can be called from ring @R,
    524511//         and the ideal @P must be defined calling first setglobalrings();
    525 // input:   ideal J
    526 // output:  ideal Jc: Returns all the free-square factors of the elements
     512// input:  ideal J
     513// output: ideal Jc: Returns all the free-square factors of the elements
    527514//         of ideal J (non repeated). Integer factors are ignored,
    528 //         even 0 is ignored. It can be called from ideal @R, but
    529 //         the given ideal J must only contain poynomials in the
    530 //         parameters.
    531 static proc facvar(ideal J)
     515//         even 0 is ignored. It can be called from ideal @R.
     516proc facvar(ideal J)
    532517//"USAGE:   facvar(J);
    533518//          J: an ideal in the parameters
    … …  
    546531  setring(@P);
    547532  def Ja=imap(RR,J);
    548   if(size(Ja)==0){return(ideal(0));}
     533  if(size(Ja)==0){setring(RR); return(ideal(0));}
    549534  Ja=elimintfromideal(Ja); // also in ideal @P
    550535  ideal Jb;
    … …  
    569554//}
    570555
    571 // Wred: eliminate the factors in the polynom f that are in W
    572 //       in ring @RP
     556// Ered: eliminates the factors in the polynom f that are non-null.
     557//       In ring @R
    573558// input:
    574559//   poly f:
    575 //   ideal W  of non-null conditions (already supposed that it is facvar)
     560//   ideal E  of null-conditions
     561//   ideal N  of non-null conditions
     562//        (E,N) represents V(E)\V(N),
     563//        Ered eliminates the non-null factors of f in V(E)\V(N)
    576564// output:
    577 //   poly f2  where the non-null conditions in W have been dropped from f
    578 static proc Wred(poly f, ideal W)
    579 {
    580   if (f==0){return(f);}
     565//   poly f2  where the non-null conditions have been dropped from f
     566proc Ered(poly f,ideal E, ideal N)
     567{
    581568  def RR=basering;
    582   setring(@RP);
    583   def ff=imap(RR,f);
    584   def RPW=imap(RR,W);
    585   def l=factorize(ff,2);
     569  setring(@R);
     570  poly ff=imap(RR,f);
     571  ideal EE=imap(RR,E);
     572  ideal NN=imap(RR,N);
     573  if((ff==0) or (equalideals(NN,ideal(1)))){setring(RR); return(f);}
     574  def v=variables(ff);
    586575  int i;
    587   poly f1=1;
    588   for(i=1;i<=size(l[1]);i++)
    589   {
    590     if ((memberpos(l[1][i],RPW)[1]) or (memberpos(-l[1][i],RPW)[1])){;}
    591     else{f1=f1*((l[1][i])^(l[2][i]));}
    592   }
    593   setring(RR);
    594   def f2=imap(@RP,f1);
    595   return(f2);
    596 }
    597 
    598 // pnormalform: reduces a polynomial wrt a red-spec dividing by N and eliminating factors in W.
    599 //              called in the ring @R
    600 //              operates in the ring @RP
    601 //              both ideals must be defined calling first setglobalrings();
     576  poly X=1;
     577  for(i=1;i<=size(v);i++){X=X*v[i];}
     578  matrix M=coef(ff,X);
     579  setring(@P);
     580  def RPE=imap(@R,EE);
     581  def RPN=imap(@R,NN);
     582  matrix Mp=imap(@R,M);
     583  poly g=Mp[2,1];
     584  if (size(Mp)!=2)
     585  {
     586    for(i=2;i<=size(Mp) div 2;i++)
     587    {
     588      g=gcd(g,Mp[2,i]);
     589    }
     590  }
     591  if (g==1){setring(RR); return(f);}
     592  else
     593  {
     594    def wg=factorize(g,2);
     595    if (wg[1][1]==1){setring(RR); return(f);}
     596    else
     597    {
     598      poly simp=1;
     599      int te;
     600      for(i=1;i<=size(wg[1]);i++)
     601      {
     602        te=inconsistent(RPE+wg[1][i],RPN);
     603        if(te)
     604        {
     605          simp=simp*(wg[1][i])^(wg[2][i]);
     606        }
     607      }
     608    }
     609    if (simp==1){setring(RR); return(f);}
     610    else
     611    {
     612      setring(RR);
     613      def simp0=imap(@P,simp);
     614      def f2=f/simp0;
     615      return(f2);
     616    }
     617  }
     618}
     619
     620// pnormalf: reduces a polynomial f wrt a V(E)\V(N)
     621//           dividing by E and eliminating factors in N.
     622//           called in the ring @R,
     623//           operates in the ring @RP.
    602624// input:
    603625//         poly  f
     626//         ideal E  (depends only on the parameters)
    604627//         ideal N  (depends only on the parameters)
    605 //         ideal W  (depends only on the parameters)
    606 //                   (N,W) must be a red-spec (depends only on the parameters)
    607 // output: poly f2 reduced wrt to the red-spec (N,W)
    608 // note:   for security a lot of work is done. If (N,W) is already a red-spec
    609 //         it should be simplified
    610 proc pnormalform(poly f, ideal N, ideal W)
    611 "USAGE:   pnormalform(f,N,W);
    612           f: the polynomial to be reduced modulo (N,W) a reduced representation
     628//                  (E,N) represents V(E)\V(N)
     629//         optional:
     630// output: poly f2 reduced wrt to V(E)\V(N)
     631proc pnormalf(poly f, ideal E, ideal N)
     632"USAGE:   pnormalf(f,E,N);
     633          f: the polynomial to be reduced modulo V(E)\V(N)
    613634          of a segment in the parameters.
    614           N: the null conditions ideal
    615           W: the non-null conditions (set of irreducible polynomials)
     635          E: the null conditions ideal
     636          N: the non-null conditions
    616637RETURN:   a reduced polynomial g of f, whose coefficients are reduced
    617           modulo N and having no factor in W.
    618 NOTE:     Should be called from ring Q[a][x], and the global rings @R, @P
    619           and @RP must be defined. These rings can be created by calling
    620           previously setglobalrings();
    621           Ideals N and W must be given by polynomials
    622           in the parameters forming a reduced-representation (see
    623           definition in the paper).
     638          modulo E and having no factor in N.
     639NOTE: Should be called from ring Q[a][x].
     640          Ideals E and N must be given by polynomials
     641          in the parameters.
    624642KEYWORDS: division, pdivi, reduce
    625 EXAMPLE:  pnormalform; shows an example"
     643EXAMPLE:  pnormalf; shows an example"
    626644{
    627645    def RR=basering;
    628     setglobalrings();
     646    int te=0;
     647    if (defined(@P)){te=1;}
     648    else{setglobalrings();}
    629649    setring(@RP);
    630650    def fa=imap(RR,f);
     651    def Ea=imap(RR,E);
    631652    def Na=imap(RR,N);
    632     def Wa=imap(RR,W);
    633653    option(redSB);
    634     Na=std(Na);
    635     def r=cld(reduce(fa,Na));
    636     def f1=Wred(r[1],Wa);
     654    Ea=std(Ea);
     655    def r=cld(reduce(fa,Ea));
     656    poly f1=r[1];
     657    f1=Ered(r[1],Ea,Na);
    637658    setring(RR);
    638659    def f2=imap(@RP,f1);
     660    if(te==0){kill @R; kill @RP; kill @P;}
    639661    return(f2)
    640 }
     662};
    641663example
    642664{ "EXAMPLE:"; echo = 2;
    643665  ring R=(0,a,b,c),(x,y),dp;
    644   setglobalrings();
    645666  poly f=(b^2-1)*x^3*y+(c^2-1)*x*y^2+(c^2*b-b)*x+(a-bc)*y;
    646   ideal N=(ab-c)*(a-b),(a-bc)*(a-b);
    647   ideal W=a^2-b^2,bc;
    648   def r=redspec(N,W);
    649   pnormalform(f,r[1],r[2]);
     667  ideal E=(c-1);
     668  ideal N=a-b;
     669  pnormalf(f,E,N);
    650670}
    651671
    … …  
    655675// input: two ideals in the ring @P
    656676// output the intersection of both (is not a GB)
    657 static proc idint(ideal I, ideal J)
     677proc idint(ideal I, ideal J)
    658678{
    659679  def RR=basering;
    … …  
    672692}
    673693
    674 // redspec: generates a red-representation
    675 //          called in any ring
    676 //          it changes to the ring @P
    677 //          So the globalrings @P, @RP, @R, must be created before
    678 //          using it by calling setglobalrings();
    679 // input:
    680 //   ideal N : the ideal of null-conditions
    681 //   ideal W : set of non-null polynomials:
    682 //             if W corresponds to no non null conditions then W=ideal(0)
    683 //             otherwise it should be given as an ideal.
    684 // returns: list (Na,Wa,DGN)
    685 // the completely reduced representation:
    686 //   Na = ideal reduced and radical of the red-spec
    687 //   facvar(Wa) = ideal the reduced non-null set of polynomials of the red-spec.
    688 //             if it corresponds to no non null conditions then it is ideal(0)
    689 //             otherwise the ideal is returned.
    690 //   DGN = the list of prime ideals associated to Na (uses primASSGTZ in "primdec.lib")
    691 //   none of the polynomials in facvar(Wa) are contained in none of the ideals in DGN
    692 //   If the given conditions are not compatible, then N=ideal(1) and DGN=list(ideal(1))
    693 proc redspec(ideal Ni, ideal Wi)
    694 //"USAGE:   redspec(N,W);
    695 //          N: null conditions ideal
    696 //          W: set of non-null polynomials (ideal)
    697 //RETURN:   a list (N1,W1,L1) containing a red-representation of the segment (N,W).
    698 //          N1 is the radical reduced ideal characterizing the segment.
    699 //          V(N1) is the Zariski closure of the segment (N,W).
    700 //          The segment S=V(N1) \ V(h), where h=prod(w in W1)
    701 //          N1 is uniquely determined and no prime component of N1 contains none of
    702 //          the polynomials in W1. The polynomials in W1 are prime and reduced
    703 //          wrt N1, and are considered non-null on the segment.
    704 //          L1 contains the list of prime components of N1.
    705 //NOTE:     Called from ring @R it works in ring @P, that must be defined
    706 //          by the call to setglobalrings();
    707 //          Used in the old library redcgs.lib.
    708 //KEYWORDS: representation
    709 //EXAMPLE:  redspec; shows an example"
    710 {
    711   ideal Nc;
    712   ideal Wc;
    713   def RR=basering;
    714   setring(@P);
    715   def N=imap(RR,Ni);
    716   def W=imap(RR,Wi);
    717   ideal Wa;
    718   ideal Wb;
    719   if(size(W)==0){Wa=ideal(0);}
    720      //when there are no non-null conditions then W=ideal(0)
    721   else
    722   {
    723     Wa=facvar(W);
    724   }
    725   if (size(N)==0)
    726   {
    727     setring(RR);
    728     Wc=imap(@P,Wa);
    729     return(list(ideal(0), Wc, list(ideal(0))));
    730   }
    731   int i;
    732   list LNb;
    733   list LNa;
    734   def LN=minGTZ(N);
    735   for (i=1;i<=size(LN);i++)
    736   {
    737     option(redSB);
    738     LNa[i]=std(LN[i]);
    739   }
    740   poly h=1;
    741   if (size(Wa)!=0)
    742   {
    743     for(i=1;i<=size(Wa);i++){h=h*Wa[i];}
    744   }
    745   ideal Na;
    746   intvec save_opt=option(get);
    747   if (size(N)!=0 and (size(LNa)>0))
    748   {
    749     option(returnSB);
    750     Na=intersect(LNa[1..size(LNa)]);
    751     option(redSB);
    752     Na=std(Na);
    753     option(set,save_opt);
    754   }
    755   attrib(Na,"isSB",1);
    756   if (reduce(h,Na,1)==0)
    757   {
    758     setring(RR);
    759     Wc=imap(@P,Wa);
    760     return(list (ideal(1),Wc,list(ideal(1))));
    761   }
    762   i=1;
    763   while(i<=size(LNa))
    764   {
    765     if (reduce(h,LNa[i],1)==0){LNa=delete(LNa,i);}
    766     else{ i++;}
    767   }
    768   if (size(LNa)==0)
    769   {
    770     setring(RR);
    771     return(list(ideal(1),ideal(0),list(ideal(1))));
    772   }
    773   option(returnSB);
    774   ideal Nb=intersect(LNa[1..size(LNa)]);
    775   option(redSB);
    776   Nb=std(Nb);
    777   option(set,save_opt);
    778   if (size(Wa)==0)
    779   {
    780     setring(RR);
    781     Nc=imap(@P,Nb);
    782     Wc=imap(@P,Wa);
    783     LNb=imap(@P,LNa);
    784     return(list(Nc,Wc,LNb));
    785   }
    786   Wb=ideal(0);
    787   attrib(Nb,"isSB",1);
    788   for (i=1;i<=size(Wa);i++){Wb[i]=reduce(Wa[i],Nb);}
    789   Wb=facvar(Wb);
    790   if (size(LNa)!=0)
    791   {
    792     setring(RR);
    793     Nc=imap(@P,Nb);
    794     Wc=imap(@P,Wb);
    795     LNb=imap(@P,LNa);
    796     return(list(Nc,Wc,LNb))
    797   }
    798   else
    799   {
    800     setring(RR);
    801     Nd=imap(@P,Nb);
    802     Wc=imap(@P,Wb);
    803     kill LNb;
    804     list LNb;
    805     return(list(Nd,Wc,LNb))
    806   }
    807 }
    808 //example
    809 //{ "EXAMPLE:"; echo = 2;
    810 //  ring r=(0,a,b,c),(x,y),dp;
    811 //  setglobalrings();
    812 //  ideal N=(ab-c)*(a-b),(a-bc)*(a-b);
    813 //  ideal W=a^2-b^2,bc;
    814 //  redspec(N,W);
    815 //}
    816 
    817694// lesspol: compare two polynomials by its leading power products
    818695// input:  two polynomials f,g in the ring @R
    819696// output: 0 if f<g,  1 if f>=g
    820 static proc lesspol(poly f, poly g)
     697proc lesspol(poly f, poly g)
    821698{
    822699  if (leadmonom(f)<leadmonom(g)){return(1);}
    … …  
    830707    }
    831708  }
    832 }
     709};
    833710
    834711// delfromideal: deletes the i-th polynomial from the ideal F
    835 static proc delfromideal(ideal F, int i)
     712proc delfromideal(ideal F, int i)
    836713{
    837714  int j;
    … …  
    839716  if (size(F)<i){ERROR("delfromideal was called incorrect arguments");}
    840717  if (size(F)<=1){return(ideal(0));}
    841   if (i==0){return(F);}
     718  if (i==0){return(F)};
    842719  for (j=1;j<=size(F);j++)
    843720  {
    … …  
    850727// input: ideals I,J
    851728// output: the ideal J without the polynomials in I
    852 static proc delidfromid(ideal I, ideal J)
     729proc delidfromid(ideal I, ideal J)
    853730{
    854731  int i; list r;
    … …  
    866743
    867744// sortideal: sorts the polynomials in an ideal by lm in ascending order
    868 static proc sortideal(ideal Fi)
     745proc sortideal(ideal Fi)
    869746{
    870747  def RR=basering;
    … …  
    894771// mingb: given a basis (gb reducing) it
    895772// order the polynomials is ascending order and
    896 // eliminate the polynomials whose lpp is divisible by some
     773// eliminates the polynomials whose lpp are divisible by some
    897774// smaller one
    898 static proc mingb(ideal F)
     775proc mingb(ideal F)
    899776{
    900777  int t; int i; int j;
    … …  
    917794}
    918795
    919 // redgb: given a minimal basis (gb reducing) it
     796// redgbn: given a minimal basis (gb reducing) it
     797// reduces each polynomial wrt to V(E) \ V(N)
     798proc redgbn(ideal F, ideal E, ideal N)
     799{
     800  int te=0;
     801  if (defined(@P)==1){te=1;}
     802  ideal G=F;
     803  ideal H;
     804  int i;
     805  if (size(G)==0){return(ideal(0));}
     806  for (i=1;i<=size(G);i++)
     807  {
     808    H=delfromideal(G,i);
     809    G[i]=pnormalf(pdivi(G[i],H)[1],E,N);
     810    G[i]=primepartZ(G[i]);
     811  }
     812  if(te==1){setglobalrings();}
     813  return(G);
     814};
     815
     816// eliminates repeated elements form an ideal
     817proc elimrepeated(ideal F)
     818{
     819  int i;
     820  ideal FF;
     821  FF[1]=F[1];
     822  for (i=2;i<=ncols(F);i++)
     823  {
     824    if (not(memberpos(F[i],FF)[1]))
     825    {
     826      FF[size(FF)+1]=F[i];
     827    }
     828  }
     829  return(FF);
     830}
     831
     832// equalideals
     833// input: 2 ideals F and G;
     834// output: 1 if they are identical (the same polynomials in the same order)
     835//         0 else
     836proc equalideals(ideal F, ideal G)
     837{
     838  int i=1; int t=1;
     839  if (size(F)!=size(G)){return(0);}
     840  while ((i<=size(F)) and (t))
     841  {
     842    if (F[i]!=G[i]){t=0;}
     843    i++;
     844  }
     845  return(t);
     846}
     847
     848// delintvec
     849// input: intvec V
     850//        int i
     851// output:
     852//        intvec W (equal to V but the coordinate i is deleted
     853proc delintvec(intvec V, int i)
     854{
     855  int j;
     856  intvec W;
     857  for (j=1;j<i;j++){W[j]=V[j];}
     858  for (j=i+1;j<=size(V);j++){W[j-1]=V[j];}
     859  return(W);
     860}
     861
     862//**************Begin homogenizing************************
     863
     864// ishomog:
     865// Purpose: test if a polynomial is homogeneous in the variables or not
     866// input:  poly f
     867// output  1 if f is homogeneous, 0 if not
     868proc ishomog(f)
     869{
     870  int i; poly r; int d; int dr;
     871  if (f==0){return(1);}
     872  d=deg(f); dr=d; r=f;
     873  while ((d==dr) and (r!=0))
     874  {
     875    r=r-lead(r);
     876    dr=deg(r);
     877  }
     878  if (r==0){return(1);}
     879  else{return(0);}
     880}
     881
     882// postredgb: given a minimal basis (gb reducing) it
    920883// reduces each polynomial wrt to the others
    921 static proc redgb(ideal F, ideal N, ideal W)
    922 {
     884proc postredgb(ideal F)
     885{
     886  int te=0;
     887  if(defined(@P)==1){te=1;}
    923888  ideal G;
    924889  ideal H;
    … …  
    928893  {
    929894    H=delfromideal(F,i);
    930     G[i]=pnormalform(pdivi(F[i],H)[1],N,W);
    931   }
     895    G[i]=pdivi(F[i],H)[1];
     896  }
     897  if(te==1){setglobalrings();}
    932898  return(G);
    933899}
    934900
    935 //********************Main routines for buildtree******************
    936 
    937 // splitspec: a new leading coefficient f is given to a red-spec
    938 //            then splitspec computes the two new red-spec by
    939 //            considering it null, and non null.
    940 // in ring @P
    941 // given f, and the red-spec (N,W)
    942 //     it outputs the null and the non-null red-spec adding f.
    943 //     if some of the output representations has N=1 then
    944 //     there must be no split and buildtree must continue on
    945 //     the compatible red-spec
    946 // input:  poly f coefficient to split if needed
    947 //         list r=(N,W,LN) redspec
    948 // output: list L = list(ideal N0, ideal W0), list(ideal N1, ideal W1), cond
    949 static proc splitspec(poly fi, list ri)
    950 {
    951   def RR=basering;
    952   def Ni=ri[1];
    953   def Wi=ri[2];
    954   setring(@P);
    955   def f=imap(RR,fi);
    956   def N=imap(RR,Ni);
    957   def W=imap(RR,Wi);
    958   f=Wred(f,W);
    959   def N0=N;
    960   def W1=W;
    961   N0[size(N0)+1]=f;
    962   def r0=redspec(N0,W);
    963   W1[size(W1)+1]=f;
    964   def r1=redspec(N,W1);
    965   setring(RR);
    966   def ra0=imap(@P,r0);
    967   def ra1=imap(@P,r1);
    968   def cond=imap(@P,f);
    969   return (list(ra0,ra1,cond));
    970 }
    971 
    972 // redcoefs
    973 // 15/09/2010
    974 static proc redcoefs(poly f, ideal N)
    975 {
    976   def f1=f; int test0=1; poly lc; poly lm;
    977   poly lc1;
    978   def RR=basering;
    979   setring(@P);
    980   poly lcp;
    981   def Np=imap(RR,N);
    982   attrib(Np,"isSB",1);
    983   setring(RR);
    984   while((test0==1) and (f1<>0))
    985   {
    986     lc=leadcoef(f1);
    987     lm=leadmonom(f1);
    988     setring(@P);
    989     lcp=imap(RR,lc);
    990     lcp=reduce(lcp,Np);
    991     setring(RR);
    992     lc1=imap(@P,lcp);
    993     if(lc1<>0){test0=0;}
    994     f1=f1+(lc1-lc)*lm;
    995   }
    996   return(f1);
    997 }
    998 
    999 // discusspolys: given a basis B and a red-spec (N,W), it analyzes the
    1000 //               leadcoef of the polynomials in B until it finds
    1001 //               that one of them can be either null or non null.
    1002 //               If at the end only the non null option is compatible
    1003 //               then the reduced B has all the leadcoef non null.
    1004 //               Else recbuildtree must split.
    1005 // ring @R
    1006 // input:  ideal B
    1007 //         ideal N
    1008 //         ideal W (a reduced-representation)
    1009 // output: list of ((N0,W0,LN0),(N1,W1,LN1),Br,cond)
    1010 //         cond is the condition to branch
    1011 static proc discusspolys(ideal B, list r)
    1012 {
    1013   poly f;     poly f1;    poly f2;
    1014   poly cond;
    1015   def N=r[1]; def W=r[2]; def LN=r[3];
    1016   def Ba=B;   def F=B;
    1017   ideal N0=1; def W0=W;   list LN0=ideal(1);
    1018   def N1=N;   def W1=W;   def LN1=LN;
    1019   list L;
    1020   list M;     list M0;    list M1;
    1021   list rr;
    1022   if (size(B)==0)
    1023   {
    1024     M0=N0,W0,LN0; // incompatible
    1025     M1=N1,W1,LN1;
    1026     M=M0,M1,B,poly(1);
    1027     return(M);
    1028   }
    1029   while ((size(F)!=0) and ((N0[1]==1) or (N1[1]==1)))
    1030   {
    1031     f=F[1];
    1032     F=delfromideal(F,1);
    1033     f1=pnormalform(f,N,W);
    1034     rr=memberpos(f,Ba);
    1035     if (f1!=0)
    1036     {
    1037       Ba[rr[2]]=f1;
    1038       if (pardeg(leadcoef(f1))!=0)
    1039       {
    1040         f2=Wred(leadcoef(f1),W);
    1041         L=splitspec(f2,list(N,W,LN));
    1042         N0=L[1][1]; W0=L[1][2]; LN0=L[1][3]; N1=L[2][1]; W1=L[2][2]; LN1=L[2][3];
    1043         cond=L[3];
    1044       }
    1045     }
    1046     else
    1047     {
    1048       Ba=delfromideal(Ba,rr[2]);
    1049       N0=ideal(1); //F=ideal(0);
    1050     }
    1051   }
    1052   M0=N0,W0,LN0;
    1053   M1=N1,W1,LN1;
    1054   M=M0,M1,Ba,cond;
    1055   return(M);
    1056 }
    1057 
    1058 // discussSpolys: given a basis B and a red-spec (N,W), it analyzes the
    1059 //                leadcoef of the polynomials in B until it finds
    1060 //                that one of them can be either null or non null.
    1061 //                If at the end only the non null option is compatible
    1062 //                then the reduced B has all the leadcoef non null.
    1063 //                Else recbuildtree must split.
    1064 // ring @R
    1065 // input:  ideal B
    1066 //         ideal N
    1067 //         ideal W (a reduced-representation)
    1068 //         list  P current set of pairs of polynomials from B to be tested.
    1069 // output: list of (N0,W0,LN0),(N1,W1,LN1),Br,Pr,cond]
    1070 //         list Pr the not checked list of pairs.
    1071 static proc discussSpolys(ideal B, list r, list P)
    1072 {
    1073   int i; int j; int k;
    1074   int npols; int nSpols; int tt;
    1075   poly cond=1;
    1076   poly lm; poly lpf; poly lpg;
    1077   def F=B; def Pa=P; list Pa0;
    1078   def N=r[1]; def W=r[2]; def LN=r[3];
    1079   ideal N0=1; def W0=W; list LN0=ideal(1);
    1080   def N1=N; def W1=W; def LN1=LN;
    1081   ideal Bw;
    1082   poly S;
    1083   list L; list L0; list L1;
    1084   list M; list M0; list M1;
    1085   list pair;
    1086   list KK; int loc;
    1087   int crit;
    1088   poly h;
    1089   if (size(B)==0)
    1090   {
    1091     M0=N0,W0,LN0;
    1092     M1=N1,W1,LN1;
    1093     M=M0,M1,ideal(0),Pa,cond;
    1094     return(M);
    1095   }
    1096   tt=1;
    1097   i=1;
    1098   while ((tt) and (i<=size(B)))
    1099   {
    1100     h=B[i];
    1101     for (j=1;j<=npars(@R);j++)
    1102     {
    1103       h=subst(h,par(j),0);
    1104     }
    1105     if (h!=B[i]){tt=0;}
    1106     i++;
    1107   }
    1108   if (tt)
    1109   {
    1110     //"T_ a non parametric system occurred";
    1111     def RR=basering;
    1112     def RL=ringlist(RR);
    1113     RL[1]=0;
    1114     def LRR=ring(RL);
    1115     setring(LRR);
    1116     def BP=imap(RR,B);
    1117     option(redSB);
    1118     BP=std(BP);
    1119     setring(RR);
    1120     B=imap(LRR,BP);
    1121     M0=ideal(1),W0,LN0;
    1122     M1=N1,W1,LN1;
    1123     M=M0,M1,B,list(),cond;
    1124     return(M);
    1125   }
    1126   if (size(Pa)==0){npols=size(B); Pa=orderingpairs(F); nSpols=size(Pa);}
    1127   while ((size(Pa)!=0) and (N0[1]==1) or (N1[1]==1))
    1128   {
    1129     pair=Pa[1];
    1130     i=pair[1];
    1131     j=pair[2];
    1132     Pa=delete(Pa,1);
    1133     // Buchberger 1st criterion (not needed here, it is already eliminated
    1134     // when creating the list of pairs
    1135     for (k=1;k<=size(Pa);k++){Pa0[k]=delete(Pa[k],3);}
    1136     crit=0;
    1137     if (not(crit))
    1138     {
    1139       S=pspol(F[i],F[j]);
    1140       KK=pdivi(S,F);
    1141       S=KK[1];
    1142       if (S!=0)
    1143       {
    1144         S=pnormalform(S,N,W);
    1145         if (S!=0)
    1146         {
    1147           L=discusspolys(ideal(S),list(N,W,LN));
    1148           N0=L[1][1];
    1149           W0=L[1][2];
    1150           LN0=L[1][3];
    1151           N1=L[2][1];
    1152           W1=L[2][2];
    1153           LN1=L[2][3];
    1154           S=L[3][1];
    1155           cond=L[4];
    1156           if (S==1)
    1157           {
    1158             M0=ideal(1),W0,list(ideal(1));
    1159             M1=N1,W1,LN1;
    1160             M=M0,M1,ideal(1),list(),cond;
    1161             return(M);
    1162           }
    1163           if (S!=0)
    1164           {
    1165             F[size(F)+1]=S;
    1166             npols=size(F);
    1167             for (k=1;k<size(F);k++)
    1168             {
    1169               lm=lcmlmonoms(F[k],S);
    1170               // Buchberger 1st criterion
    1171               lpf=leadmonom(F[k]);
    1172               lpg=leadmonom(S);
    1173               if (lpf*lpg!=lm)
    1174               {
    1175                 pair=k,size(F),lm;
    1176                 Pa=placepairinlist(pair,Pa);
    1177                 nSpols=size(Pa);
    1178               }
    1179             }
    1180             if (N0[1]==1){N=N1; W=W1; LN=LN1;}
    1181           }
    1182         }
    1183       }
    1184     }
    1185   }
    1186   M0=N0,W0,LN0;
    1187   M1=N1,W1,LN1;
    1188   M=M0,M1,F,Pa,cond;
    1189   return(M);
    1190 }
    1191 
    1192 // lcmlmonoms: computes the lcm of the leading monomials
    1193 //             of the polynomils f and g
    1194 // ring @R
    1195 static proc lcmlmonoms(poly f,poly g)
    1196 {
    1197   def lf=leadmonom(f);
    1198   def lg=leadmonom(g);
    1199   def gls=gcd(lf,lg);
    1200   return((lf*lg)/gls);
    1201 }
    1202 
    1203 // placepairinlist
    1204 // 15/09/2010
    1205 // input:  given a new pair of the form (i,j,lmij)
    1206 //         and a list of pairs of the same form
    1207 // ring @R
    1208 // output: it inserts the new pair in ascending order of lmij
    1209 static proc placepairinlist(list pair,list P)
    1210 {
    1211   list Pr;
    1212   if (size(P)==0){Pr=insert(P,pair); return(Pr);}
    1213   if (pair[3]<P[1][3]){Pr=insert(P,pair); return(Pr);}
    1214   if (pair[3]>=P[size(P)][3]){Pr=insert(P,pair,size(P)); return(Pr);}
    1215   kill Pr;
    1216   list Pr;
    1217   int j;
    1218   int i=1;
    1219   int loc=0;
    1220   while((i<=size(P)) and (loc==0))
    1221   {
    1222     if (pair[3]>=P[i][3]){j=i; i++;}
    1223     else{loc=1; j=i-1;}
    1224   }
    1225   Pr=insert(P,pair,j);
    1226   return(Pr);
    1227 }
    1228 
    1229 // orderingpairs:
    1230 // input:  ideal F
    1231 // output: list of ordered pairs (i,j,lcmij) of F in ascending order of lcmij
    1232 //         if a pair verifies Buchberger 1st criterion it is not stored
    1233 // ring @R
    1234 static proc orderingpairs(ideal F)
    1235 {
    1236   int i;
    1237   int j;
    1238   poly lm;
    1239   poly lpf;
    1240   poly lpg;
    1241   list P;
    1242   list pair;
    1243   if (size(F)<=1){return(P);}
    1244   for (i=1;i<=size(F)-1;i++)
    1245   {
    1246     for (j=i+1;j<=size(F);j++)
    1247     {
    1248       lm=lcmlmonoms(F[i],F[j]);
    1249       // Buchberger 1st criterion
    1250       lpf=leadmonom(F[i]);
    1251       lpg=leadmonom(F[j]);
    1252       if (lpf*lpg!=lm)
    1253       {
    1254         pair=(i,j,lm);
    1255         P=placepairinlist(pair,P);
    1256       }
    1257     }
    1258   }
    1259   return(P);
    1260 }
    1261 
    1262 // Buchberger 2nd criterion
    1263 // input:  integers i,j
    1264 //         list P of pairs of the form (i,j) not yet verified
    1265 // ring @R
    1266 //         not used (it increases time)
    1267 static proc criterion(int i, int j, list P, ideal B)
    1268 {
    1269   def lcmij=lcmlmonoms(B[i],B[j]);
    1270   int crit=0;
    1271   int k=1;
    1272   list ik; list jk;
    1273   while ((k<=size(B)) and (crit==0))
    1274   {
    1275     if ((k!=i) and (k!=j))
    1276     {
    1277       if (i<k){ik=i,k;} else{ik=k,i;}
    1278       if (j<k){jk=i,k;} else{jk=k,j;}
    1279       if (not((memberpos(ik,P)[1]) or (memberpos(jk,P)[1])))
    1280       {
    1281         if ((lcmij)/leadmonom(B[k])!=0){crit=1;}
    1282       }
    1283     }
    1284     k++;
    1285   }
    1286   return(crit);
    1287 }
    1288 
    1289 // buildtree: Basic routine of the old redcgs.lib generating a
    1290 //     first reduced CGS
    1291 //     it will define the rings @R, @P and @RP as global rings
    1292 //     and the list @T a global list that will be killed at the output
    1293 // input:  ideal F in ring K[a][x];
    1294 // output: list T of lists whose list elements are of the form
    1295 //         T[i]=list(list lab, boolean terminal, ideal B, ideal N, ideal W, list of ideals decomp of N,
    1296 //              ideal of monomials lpp);
    1297 // all the ideals are in the ring K[a][x];
    1298 static proc buildtree(ideal F, list #)
    1299 //"USAGE:   buildtree(F);
    1300 //          F: ideal in Q[a][x] (parameters and variables) to be discussed.
    1301 //          It outputs the whole discussion tree to construct the
    1302 //          first disjoint reduced CGS. It is the old version of the new
    1303 //          cgsdr routine. It remains on the library for didactic purposes
    1304 //          and is, in general, less efficient.
    1305 //          Also, for some problems where cgsdr does stack, sometimes
    1306 //          buildtree is able to obtain the result.
    1307 //          The output of buildtree contains the whole information about the discussion
    1308 //          process (the whole tree discussion) and can be reduced to
    1309 //          somewhat similar to the output of cgsdr after calling
    1310 //          setglobalrings(); then applying finalcases and then groupsegments to the
    1311 //          output of buidtree. This is automatically done by the routine
    1312 //          cgsdrold also contained in the library, that outputs only the
    1313 //          CGS like the new cgsdr.
    1314 //
    1315 //RETURN:   Returns a list T describing the complete discussion tree
    1316 //          for obtaining a reduced and disjoint comprehensive
    1317 //          Groebner system (CGS) of the ideal F of Q[a][x] with
    1318 //          constant leading power products (lpp) of the reduced Groebner
    1319 //          basis.
    1320 //          The first element of the list is the root, and contains
    1321 //            [1] label: intvec(-1)
    1322 //            [2] number of children : int
    1323 //            [3] the ideal F
    1324 //            [4], [5], [6] the red-representation of the segment
    1325 //                (null, non-null conditions, prime components of the null
    1326 //                conditions) given (as option).
    1327 //                ideal (0), ideal (1), list(ideal(0)) is assumed if
    1328 //                no optional conditions are given.
    1329 //            [7] the set of lpp of ideal F
    1330 //            [8] condition that was taken to reach the vertex
    1331 //                (poly 1, for the root).
    1332 //          The remaining elements of the list represent vertices of the tree:
    1333 //          with the same structure:
    1334 //            [1] label: intvec (1,0,0,1,...) gives its position in the tree:
    1335 //                first branch condition is taken non-null, second null,...
    1336 //            [2] number of children (0 if it is a terminal vertex)
    1337 //            [3] the specialized ideal with the previous assumed conditions
    1338 //                to reach the vertex
    1339 //            [4],[5],[6] the red-representation of the segment corresponding
    1340 //                to the previous assumed conditions to reach the vertex
    1341 //            [7] the set of lpp of the specialized ideal at this stage
    1342 //            [8] condition that was taken to reach the vertex from the
    1343 //                father's vertex (that was taken non-null if the last
    1344 //                integer in the label is 1, and null if it is 0)
    1345 //          The terminal vertices form a disjoint partition of the parameter space
    1346 //          whose bases specialize to the reduced Groebner basis of the
    1347 //          specialized ideal on each point of the segment and preserve
    1348 //          the lpp. So they form a disjoint reduced CGS.
    1349 //NOTE:     The basering R, must be of the form Q[a][x], a=parameters,
    1350 //          x=variables, and should be defined previously. The ideal must
    1351 //          be defined on R.
    1352 //          The call of finalcases applied to the output of buildtree
    1353 //          selects the terminal vertices forming the disjoint and reduced
    1354 //          CGS. To obtain the output similar
    1355 //          to that of the new cgsdr procedure one can call instead
    1356 //          cgsdrold.
    1357 //
    1358 //          The content of buildtree can be written in a file that is readable
    1359 //          by Maple in order to plot its content using buildtreetoMaple;
    1360 //          The file written by buildtreetoMaple when is read in a Maple
    1361 //          worksheet can be plotted using the dbgb routine tplot;
    1362 //
    1363 //KEYWORDS: CGS, disjoint, reduced, comprehensive Groebner system
    1364 //EXAMPLE:  buildtree; shows an example"
    1365 {
    1366   list @T;
    1367   exportto(Top,@T);
    1368   setglobalrings();
    1369   int i;
    1370   int j;
    1371   poly f;
    1372   poly cond=1;
    1373   list LN;
    1374   LN[1]=ideal(0);
    1375   def N=ideal(0);
    1376   def W=ideal(1);
    1377   int comment=0;
    1378   list L=#;
    1379   for(i=1;i<=size(L) div 2;i++)
    1380   {
    1381     if(L[2*i-1]=="null"){N=L[2*i];}
    1382     else
    1383     {
    1384       if(L[2*i-1]=="nonnull"){W=L[2*i];}
    1385       else
    1386       {
    1387         if(L[2*i-1]=="comment"){comment=L[2*i];}
    1388       }
    1389     }
    1390   }
    1391   ideal B;
    1392   if(equalideals(N,ideal(0))==0)
    1393   {
    1394     def LL=redspec(N,W);
    1395     N=LL[1];
    1396     W=LL[2];
    1397     LN=LL[3];
    1398     for (i=1;i<=size(F);i++)
    1399     {
    1400       f=pnormalform(F[i],N,W);
    1401       if (f!=0){B[size(B)+1]=f;}
    1402     }
    1403   }
    1404   else {B=F;}
    1405   def lpp=ideal(0);
    1406   if (size(B)==0){lpp=ideal(0);}
    1407   else
    1408   {
    1409      for (i=1;i<=size(B);i++){lpp[i]=leadmonom(B[i]);}
    1410     // lpp=ideal of lead power product of the polys in B
    1411   }
    1412   intvec lab=-1;
    1413   int term=0;
    1414   list root;
    1415   root[1]=lab;
    1416   root[2]=term;
    1417   root[3]=B;
    1418   root[4]=N;
    1419   root[5]=W;
    1420   root[6]=LN;
    1421   root[7]=lpp;
    1422   root[8]=cond;
    1423   @T[1]=root;
    1424   list P;
    1425   recbuildtree(root,P);
    1426   def T=@T;
    1427   kill @T;
    1428   kill @P; kill @RP; kill @R;
    1429   return(T)
    1430 }
    1431 //example
    1432 //{ "EXAMPLE:"; echo = 2;
    1433 //  ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
    1434 //  "Casas conjecture for degree 4";
    1435 //  ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
    1436 //          x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
    1437 //          x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
    1438 //          x2^2+(2*a3)*x2+(a2),
    1439 //          x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
    1440 //          x3+(a3);
    1441 //  def T=buildtree(F); "buildtree(F)="; T;
    1442 //  setglobalrings();
    1443 //  def FC=finalcases(T);
    1444 //  "finalcases(buildtree(F))="; FC;
    1445 //  "groupsegments(finalcases(buildtree(F)))=";
    1446 //  groupsegments(FC);
    1447 //  buildtreetoMaple(T,"Tb","Tb.txt"); " ";
    1448 //  "Compare with cgsdrold"; " ";
    1449 //  def CDR=cgsdrold(F);
    1450 //  "cgsdrold(F)="; CDR;
    1451 //}
    1452 
    1453 // recbuildtree: auxilliary recursive routine called by buildtree
    1454 static proc recbuildtree(list v, list P)
    1455 {
    1456   def vertex=v;
    1457   int i;
    1458   int j;
    1459   int pos;
    1460   list P0;
    1461   list P1;
    1462   poly f;
    1463   def lab=vertex[1];
    1464   if ((size(lab)>1) and (lab[1]==-1))
    1465   {lab=lab[2..size(lab)];}
    1466   def term=vertex[2];
    1467   def B=vertex[3];
    1468   def N=vertex[4];
    1469   def W=vertex[5];
    1470   def LN=vertex[6];
    1471   def lpp=vertex[7];
    1472   def cond=vertex[8];
    1473   def lab0=lab;
    1474   def lab1=lab;
    1475   if ((size(lab)==1) and (lab[1]==-1))
    1476   {
    1477     lab0=0;
    1478     lab1=1;
    1479   }
    1480   else
    1481   {
    1482     lab0[size(lab)+1]=0;
    1483     lab1[size(lab)+1]=1;
    1484   }
    1485   list vertex0;
    1486   list vertex1;
    1487   ideal B0;
    1488   ideal lpp0;
    1489   ideal lpp1;
    1490   ideal N0=1;
    1491   def W0=ideal(0);
    1492   list LN0=ideal(1);
    1493   def B1=B;
    1494   def N1=N;
    1495   def W1=W;
    1496   list LN1=LN;
    1497   list L;
    1498   if (size(P)==0)
    1499   {
    1500     L=discusspolys(B,list(N,W,LN));
    1501     N0=L[1][1];
    1502     W0=L[1][2];
    1503     LN0=L[1][3];
    1504     N1=L[2][1];
    1505     W1=L[2][2];
    1506     LN1=L[2][3];
    1507     B1=L[3];
    1508     cond=L[4];
    1509   }
    1510   if ((size(B1)!=0) and (N0[1]==1))
    1511   {
    1512     L=discussSpolys(B1,list(N1,W1,LN1),P);
    1513     N0=L[1][1];
    1514     W0=L[1][2];
    1515     LN0=L[1][3];
    1516     N1=L[2][1];
    1517     W1=L[2][2];
    1518     LN1=L[2][3];
    1519     B1=L[3];
    1520     P1=L[4];
    1521     cond=L[5];
    1522     lpp=ideal(0);
    1523     for (i=1;i<=size(B1);i++){lpp[i]=leadmonom(B1[i]);}
    1524   }
    1525   vertex[3]=B1;
    1526   vertex[4]=N1; // unnecessary
    1527   vertex[5]=W1; // unnecessary
    1528   vertex[6]=LN1;// unnecessary
    1529   vertex[7]=lpp;
    1530   vertex[8]=cond;
    1531   if (size(@T)>0)
    1532   {
    1533     pos=size(@T)+1;
    1534     @T[pos]=vertex;
    1535   }
    1536   if ((N0[1]!=1) and (N1[1]!=1))
    1537   {
    1538     vertex1[1]=lab1;
    1539     vertex1[2]=0;
    1540     vertex1[3]=B1;
    1541     vertex1[4]=N1;
    1542     vertex1[5]=W1;
    1543     vertex1[6]=LN1;
    1544     vertex1[7]=lpp1;
    1545     vertex1[8]=cond;
    1546     if (size(B1)==0){B0=ideal(0); lpp0=ideal(0);}
    1547     else
    1548     {
    1549       j=1;
    1550       lpp0=ideal(0);
    1551       for (i=1;i<=size(B1);i++)
    1552       {
    1553         f=pnormalform(B1[i],N0,W0);
    1554         if (f!=0){B0[j]=f; lpp0[j]=leadmonom(f);j++;}
    1555       }
    1556     }
    1557     vertex0[1]=lab0;
    1558     vertex0[2]=0;
    1559     vertex0[3]=B0;
    1560     vertex0[4]=N0;
    1561     vertex0[5]=W0;
    1562     vertex0[6]=LN0;
    1563     vertex0[7]=lpp0;
    1564     vertex0[8]=cond;
    1565     recbuildtree(vertex0,P0);
    1566     recbuildtree(vertex1,P1);
    1567   }
    1568   else
    1569   {
    1570     if (equalideals(N1,ideal(1))==0)
    1571     {
    1572       vertex[2]=1;
    1573       B1=mingb(B1);
    1574       vertex[3]=redgb(B1,N1,W1);
    1575       vertex[4]=N1;
    1576       vertex[5]=W1;
    1577       vertex[6]=LN1;
    1578       lpp=ideal(0);
    1579       for (i=1;i<=size(vertex[3]);i++){lpp[i]=leadmonom(vertex[3][i]);}
    1580       vertex[7]=lpp;
    1581       vertex[8]=cond;
    1582       @T[pos]=vertex;
    1583       //print(vertex);
    1584     }
    1585   }
    1586 }
    1587 
    1588 // RtoPrep
    1589 // Computes the P-representaion of a R-representaion (N,W,L) of a set
    1590 // input:
    1591 //    ideal N (null conditions, must be radical)
    1592 //    ideal W (non-null conditions ideal)
    1593 //    list L  must contain the radical decomposition of N.
    1594 // output:
    1595 //    the ((p_1,(p_11,..,p_1k_1)),..,(p_r,(p_r1,..,p_rk_r)));
    1596 //    the Prep of V(N) \ V(h), where h=prod(w in W).
    1597 static proc RtoPrep(ideal N, ideal W, list L)
    1598 {
    1599   int i; int j; list L0;
    1600   if (N[1]==1)
    1601   {
    1602     L0[1]=list(ideal(1),list(ideal(1)));
    1603     return(L0);
    1604   }
    1605   def RR=basering;
    1606   setring(@P);
    1607   ideal Np=imap(RR,N);
    1608   ideal Wp=imap(RR,W);
    1609   list Lp=imap(RR,L);
    1610   poly h=1;
    1611   for (i=1;i<=size(Wp);i++){h=h*Wp[i];}
    1612   list r; list Ti; list LL;
    1613   for (i=1;i<=size(Lp);i++)
    1614   {
    1615     Ti=minGTZ(Lp[i]+h);
    1616     for(j=1;j<=size(Ti);j++)
    1617     {
    1618       option(redSB);
    1619       Ti[j]=std(Ti[j]);
    1620     }
    1621     //list LL[i];
    1622     LL[i]=list(Lp[i],Ti);
    1623   }
    1624   setring(RR);
    1625   return(imap(@P,LL));
    1626 }
    1627 
    1628 // groupRtoPrep
    1629 // input:  L (list) is the output of groupsegments
    1630 // output: LL (list) the same list but the segments are expressed
    1631 //                   in canonical representations:
    1632 //  ( (lpp, (lab BuildTree, basis,
    1633 //             ((P_1),(P_{11},...,P_{1t1}))
    1634 //             ...
    1635 //             ((P_j),(P_{j1},...,P_{jtj}))
    1636 //          )
    1637 //          ...
    1638 //          (lab BuildTree, basis,
    1639 //             ((P_1),(P_{11},...,P_{1t1}))
    1640 //             ...
    1641 //             ((P_j),(P_{j1},...,P_{jtj}))
    1642 //          )
    1643 //    )
    1644 //    ...
    1645 //    (lpp, (lab BuildTree, basis,
    1646 //             ((P_1),(P_{11},...,P_{1t1}))
    1647 //             ...
    1648 //             ((P_j),(P_{j1},...,P_{jtj}))
    1649 //          )
    1650 //          ...
    1651 //          (lab BuildTree, basis,
    1652 //             ((P_1),(P_{11},...,P_{1t1}))
    1653 //             ...
    1654 //             ((P_j),(P_{j1},...,P_{jtj}))
    1655 //          )
    1656 //    )
    1657 //  )
    1658 static proc groupRtoPrep(list L)
    1659 {
    1660   int i; int j;
    1661   list LL; list ct;
    1662   // size(L)=number of lpp-segments
    1663   for (i=1;i<=size(L);i++)
    1664   {
    1665     LL[i]=list();
    1666     LL[i][1]=L[i][1];
    1667     // L[i][1]=lpp
    1668     LL[i][2]=list();
    1669     for (j=1;j<=size(L[i][2]);j++)
    1670     {
    1671       ct=RtoPrep(L[i][2][j][3],L[i][2][j][4],L[i][2][j][5]);
    1672       LL[i][2][j]=list();
    1673       LL[i][2][j][1]=L[i][2][j][1];
    1674       // L[i][2][j][1]=label
    1675       LL[i][2][j][2]=L[i][2][j][2];
    1676       // L[i][2][j][2]=basis
    1677       LL[i][2][j][3]=ct;
    1678     }
    1679   }
    1680   return(LL);
    1681 }
    1682 
    1683 // NEW
    1684 // input:  L (list) is the output of groupsegments
    1685 // output: LL (list) the same list but the segments are expressed
    1686 //                   in canonical representations:
    1687 //  ( (lpp, (lab BuildTree, basis,
    1688 //             ((1,u1),(lab,child,P_1)),
    1689 //             ((1,1,1),(lab,child,P_{11})),
    1690 //             ...
    1691 //             ((1,1,t1),(lab,child,P_{1t1})),
    1692 //             ...
    1693 //             ((1,u1),(lab,child,P_u1)),
    1694 //             ((1,u1,1),(lab,child,P_{u1,1})),
    1695 //             ...
    1696 //             ((1,u1,tu),(lab,child,P_{u1,tu})),
    1697 //          (lab BuildTree, basis,
    1698 //             ((1,u2),(lab,child,P_2)),
    1699 //             ((1,u1+1,1),(lab,child,P_{21})),
    1700 //             ...
    1701 //             ((1,u1+1,t2),(lab,child,P_{2,t2})),
    1702 //             ...
    1703 //             ((1,u1+..+ut),(lab,child,P_ut)),
    1704 //             ((1,u1+..+ut,1),(lab,child,P_{ut,1})),
    1705 //             ...
    1706 //             ((1,u1+..+ut,tu),(lab,child,P_{ut,tu})),
    1707 // ...
    1708 static proc groupredtocan(list L)
    1709 {
    1710   int i; int j;
    1711   list LL; list ct;
    1712   for (i=1;i<=size(L);i++)
    1713   {
    1714     LL[i]=list();
    1715     LL[i][1]=L[i][1];
    1716     LL[i][2]=list();
    1717     for (j=1;j<=size(L[i][2]);j++)
    1718     {
    1719       ct=redtocanspec(intvec(i),j-1,list(L[i][2][j][3],L[i][2][j][4],L[i][2][j][5]));
    1720       LL[i][2][j]=list();
    1721       LL[i][2][j][1]=L[i][2][j][1];
    1722       LL[i][2][j][2]=L[i][2][j][2];
    1723       LL[i][2][j][3]=ct;
    1724     }
    1725   }
    1726   return(LL);
    1727 }
    1728 
    1729 //****************End of BuildTree*************************************
    1730 
    1731 //****************Begin BuildTree To Maple*****************************
    1732 
    1733 // buildtreetoMaple: writes the list provided by buildtree to a file
    1734 //    containing the table representing it in Maple
    1735 
    1736 // writes the list L=buildtree(F) to a file "writefile" that
    1737 // is readable by Maple whith name T
    1738 // input:
    1739 //   L: the list output by buildtree
    1740 //   T: the name (string) of the output table in Maple
    1741 //   writefile: the name of the datafile where the output is to be stored
    1742 // output:
    1743 //   the result is written on the datafile "writefile" containig
    1744 //   the assignement to the table with name "T"
    1745 proc buildtreetoMaple(list L, string T, string writefile)
    1746 "USAGE:   buildtreetoMaple(T, TM, writefile);
    1747           T: is the list provided by grobcovold called with option "old",0;
    1748           TM: is the name (string) of the table variable in Maple that will represent
    1749           the output of cgsdrold;
    1750           writefile: is the name (string) of the file whereas to write the
    1751           content.
    1752 RETURN:   writes the list provided by grobcovold called with option "old",0,
    1753           (old buildtree) to a file containing the table representing it in
    1754           Maple.
    1755 KEYWORDS: cgsdrold, buildtree, Maple
    1756 EXAMPLE:  buildtreetoMaple; shows an example"
    1757 {
    1758   def R=basering;
    1759   if(size(T[1])!=8)
    1760   {
    1761     "  'Warning!' cgsdrold must be called with option 'old' set to 0 to be operative";
    1762     return();
    1763   }
    1764   short=0;
    1765   poly cond;
    1766   int i;
    1767   link LLw=":w "+writefile;
    1768   string La=string("table(",T,");");
    1769   write(LLw, La);
    1770   close(LLw);
    1771   link LLa=":a "+writefile;
    1772   def RL=ringlist(R);
    1773   list p=RL[1][2];
    1774   string param=string(p[1]);
    1775   if (size(p)>1)
    1776   {
    1777     for(i=2;i<=size(p);i++){param=string(param,",",p[i]);}
    1778   }
    1779   list v=RL[2];
    1780   string vars=string(v[1]);
    1781   if (size(v)>1)
    1782   {
    1783     for(i=2;i<=size(v);i++){vars=string(vars,",",v[i]);}
    1784   }
    1785   list xord;
    1786   list pord;
    1787   if (RL[1][3][1][1]=="dp"){pord=string("tdeg(",param);}
    1788   if (RL[1][3][1][1]=="lp"){pord=string("plex(",param);}
    1789   if (RL[3][1][1]=="dp"){xord=string("tdeg(",vars);}
    1790   if (RL[3][1][1]=="lp"){xord=string("plex(",vars);}
    1791   write(LLa,string(T,"[[9]]:=",xord,");"));
    1792   write(LLa,string(T,"[[10]]:=",pord,");"));
    1793   write(LLa,string(T,"[[11]]:=true; "));
    1794   list S;
    1795   for (i=1;i<=size(L);i++)
    1796   {
    1797     if (L[i][2]==0)
    1798     {
    1799       cond=L[i][8];
    1800       S=btcond(T,L[i],cond);
    1801       write(LLa,S[1]);
    1802       write(LLa,S[2]);
    1803     }
    1804     S=btbasis(T,L[i]);
    1805     write(LLa,S);
    1806     S=btN(T,L[i]);
    1807     write(LLa,S);
    1808     S=btW(T,L[i]);
    1809     write(LLa,S);
    1810     if (L[i][2]==1) {S=btterminal(T,L[i]); write(LLa,S);}
    1811     S=btlpp(T,L[i]);
    1812     write(LLa,S);
    1813   }
    1814   close(LLa);
    1815 }
    1816 example
    1817 { "EXAMPLE:"; echo = 2;
    1818   ring R=(0,a1,a2,a3,a4),(x1,x2,x3,x4),dp;
    1819   ideal F=x4-a4+a2,
    1820    x1+x2+x3+x4-a1-a3-a4,
    1821    x1*x3*x4-a1*a3*a4,
    1822    x1*x3+x1*x4+x2*x3+x3*x4-a1*a4-a1*a3-a3*a4;
    1823   def T=cgsdrold(F,"old",0); "T="; T;
    1824   buildtreetoMaple(T,"Tb","Tb.txt");
    1825 }
    1826 
    1827 // auxiliary routine called by buildtreetoMaple
    1828 // input:
    1829 //   list L: element i of the list of buildtree(F)
    1830 // output:
    1831 //   the string of T[[lab,1]]:=label; in Maple
    1832 static proc btterminal(string T, list L)
    1833 {
    1834   int i;
    1835   string Li;
    1836   string term;
    1837   string coma=",";
    1838   if (L[2]==0){term="false";} else {term="true";}
    1839   def lab=L[1];
    1840   string slab;
    1841   if ((size(lab)==1) and lab[1]==-1)
    1842   {slab="";coma="";} //if (size(lab)==0)
    1843   else
    1844   {
    1845     slab=string(lab[1]);
    1846     if (size(lab)>=1)
    1847     {
    1848       for (i=2;i<=size(lab);i++){slab=string(slab,",",lab[i]);}
    1849     }
    1850   }
    1851   Li=string(T,"[[",slab,coma,"6]]:=",term,"; ");
    1852   return(Li);
    1853 }
    1854 
    1855 // auxiliary routine called by buildtreetoMaple
    1856 // input:
    1857 //   list L: element i of the list of buildtree(F)
    1858 // output:
    1859 //   the string of T[[lab,3]] (basis); in Maple
    1860 static proc btbasis(string T, list L)
    1861 {
    1862   int i;
    1863   string Li;
    1864   string coma=",";
    1865   def lab=L[1];
    1866   string slab;
    1867   if ((size(lab)==1) and lab[1]==-1)
    1868   {slab="";coma="";} //if (size(lab)==0)
    1869   else
    1870   {
    1871     slab=string(lab[1]);
    1872     if (size(lab)>=1)
    1873     {
    1874       for (i=2;i<=size(lab);i++){slab=string(slab,",",lab[i]);}
    1875     }
    1876   }
    1877   Li=string(T,"[[",slab,coma,"3]]:=[",L[3],"]; ");
    1878   return(Li);
    1879 }
    1880 
    1881 // auxiliary routine called by buildtreetoMaple
    1882 // input:
    1883 //   list L: element i of the list of buildtree(F)
    1884 // output:
    1885 //   the string of T[[lab,4]] (null conditions ideal); in Maple
    1886 static proc btN(string T, list L)
    1887 {
    1888   int i;
    1889   string Li;
    1890   string coma=",";
    1891   def lab=L[1];
    1892   string slab;
    1893   if ((size(lab)==1) and lab[1]==-1)
    1894   {slab=""; coma="";}
    1895   else
    1896   {
    1897     slab=string(lab[1]);
    1898     if (size(lab)>=1)
    1899     {
    1900       for (i=2;i<=size(lab);i++){slab=string(slab,",",lab[i]);}
    1901     }
    1902   }
    1903   if ((size(lab)==1) and lab[1]==-1)
    1904     {Li=string(T,"[[",slab,coma,"4]]:=[ ]; ");}
    1905   else
    1906     {Li=string(T,"[[",slab,coma,"4]]:=[",L[4],"]; ");}
    1907   return(Li);
    1908 }
    1909 
    1910 // auxiliary routine called by buildtreetoMaple
    1911 // input:
    1912 //   list L: element i of the list of buildtree(F)
    1913 // output:
    1914 //   the string of T[[lab,5]] (null conditions ideal); in Maple
    1915 static proc btW(string T, list L)
    1916 {
    1917   int i;
    1918   string Li;
    1919   string coma=",";
    1920   def lab=L[1];
    1921   string slab;
    1922   if ((size(lab)==1) and lab[1]==-1)
    1923   {slab=""; coma="";}
    1924   else
    1925   {
    1926     slab=string(lab[1]);
    1927     if (size(lab)>=1)
    1928     {
    1929       for (i=2;i<=size(lab);i++){slab=string(slab,",",lab[i]);}
    1930     }
    1931   }
    1932   if (size(L[5])==0)
    1933     {Li=string(T,"[[",slab,coma,"5]]:={ }; ");}
    1934   else
    1935     {Li=string(T,"[[",slab,coma,"5]]:={",L[5],"}; ");}
    1936   return(Li);
    1937 }
    1938 
    1939 // auxiliary routine called by buildtreetoMaple
    1940 // input:
    1941 //   list L: element i of the list of buildtree(F)
    1942 // output:
    1943 //   the string of T[[lab,12]] (lpp); in Maple
    1944 static proc btlpp(string T, list L)
    1945 {
    1946   int i;
    1947   string Li;
    1948   string coma=",";;
    1949   def lab=L[1];
    1950   string slab;
    1951   if ((size(lab)==1) and lab[1]==-1)
    1952   {slab=""; coma="";}
    1953   else
    1954   {
    1955     slab=string(lab[1]);
    1956     if (size(lab)>=1)
    1957     {
    1958       for (i=2;i<=size(lab);i++){slab=string(slab,",",lab[i]);}
    1959     }
    1960   }
    1961   if (size(L[7])==0)
    1962   {
    1963     Li=string(T,"[[",slab,coma,"12]]:=[ ]; ");
    1964   }
    1965   else
    1966   {
    1967     Li=string(T,"[[",slab,coma,"12]]:=[",L[7],"]; ");
    1968   }
    1969   return(Li);
    1970 }
    1971 
    1972 // auxiliary routine called by buildtreetoMaple
    1973 // input:
    1974 //   list L: element i of the list of buildtree(F)
    1975 // output:
    1976 //   the list of strings of (T[[lab,0]]=0,T[[lab,1]]<>0); in Maple
    1977 static proc btcond(string T, list L, poly cond)
    1978 {
    1979   int i;
    1980   string Li1;
    1981   string Li2;
    1982   def lab=L[1];
    1983   string slab;
    1984   string coma=",";;
    1985     if ((size(lab)==1) and lab[1]==-1)
    1986     {slab=""; coma="";}
    1987   else
    1988   {
    1989     slab=string(lab[1]);
    1990     if (size(lab)>=1)
    1991     {
    1992       for (i=2;i<=size(lab);i++){slab=string(slab,",",lab[i]);}
    1993     }
    1994   }
    1995   Li1=string(T,"[[",slab+coma,"0]]:=",L[8],"=0; ");
    1996   Li2=string(T,"[[",slab+coma,"1]]:=",L[8],"<>0; ");
    1997   return(list(Li1,Li2));
    1998 }
    1999 
    2000 //*****************End of BuildtreetoMaple*********************
    2001 
    2002 //*****************Begin of Selectcases************************
    2003 
    2004 // given an intvec with sum=n
    2005 // it returns the list of intvect with the sum=n+1
    2006 static proc comp1(intvec l)
    2007 {
    2008   list L;
    2009   int p=size(l);
    2010   int i;
    2011   if (p==0){return(l);}
    2012   if (p==1){return(list(intvec(l[1]+1)));}
    2013   L[1]=intvec((l[1]+1),l[2..p]);
    2014   L[p]=intvec(l[1..p-1],(l[p]+1));
    2015   for (i=2;i<p;i++)
    2016   {
    2017     L[i]=intvec(l[1..(i-1)],(l[i]+1),l[(i+1)..p]);
    2018   }
    2019   return(L);
    2020 }
    2021 
    2022 // comp: p-compositions of n
    2023 // input
    2024 //   int n;
    2025 //   int p;
    2026 // return
    2027 //   the list of all intvec (p-composition of n)
    2028 static proc comp(int n,int p)
    2029 {
    2030   if (n<0){ERROR("comp was called with negative argument");}
    2031   if (n==0){return(list(0:p));}
    2032   int i;
    2033   int k;
    2034   list L1=comp(n-1,p);
    2035   list L=comp1(L1[1]);
    2036   list l;
    2037   list la;
    2038   for (i=2; i<=size(L1);i++)
    2039   {
    2040     l=comp1(L1[i]);
    2041     for (k=1;k<=size(l);k++)
    2042     {
    2043       if(not(memberpos(l[k],L)[1]))
    2044       {L[size(L)+1]=l[k];}
    2045     }
    2046   }
    2047   return(L);
    2048 }
    2049 
    2050 // given the matrices of coefficients and monomials m amd m1 of
    2051 // two polynomials (the first one contains all the terms of f
    2052 // and the second only those of f
    2053 // it returns the list with the comon monomials and the list of coefficients
    2054 // of the polynomial f with zeroes if necessary.
    2055 static proc adaptcoef(matrix m, matrix m1)
    2056 {
    2057   int i;
    2058   int j;
    2059   int ncm=ncols(m);
    2060   int ncm1=ncols(m1);
    2061   ideal T;
    2062   for (i=1;i<=ncm;i++){T[i]=m[1,i];}
    2063   ideal C;
    2064   for (i=1;i<=ncm;i++){C[i]=0;}
    2065   for (i=1;i<=ncm;i++)
    2066   {
    2067     j=1;
    2068     while((j<ncm1) and (m1[1,j]>m[1,i])){j++;}
    2069     if (m1[1,j]==m[1,i]){C[i]=m1[2,j];}
    2070   }
    2071   return(list(T,C));
    2072 }
    2073 
    2074 // given teh ideal of non-null conditions and an intvec lambda
    2075 // with the exponents of each w in W
    2076 // it returns the polynomial prod (w_i)^(lambda_i).
    2077 static proc WW(ideal W, intvec lambda)
    2078 {
    2079   if (size(W)==0){return(poly(1));}
    2080   poly w=1;
    2081   int i;
    2082   for (i=1;i<=ncols(W);i++)
    2083   {
    2084     w=w*(W[i])^(lambda[i]);
    2085   }
    2086   return(w);
    2087 }
    2088 
    2089 // given a polynomial f and the non-null conditions W
    2090 // WPred eliminates the factors in f that are in W
    2091 // ring @PAB
    2092 // input:
    2093 //   poly f:
    2094 //   ideal W  of non-null conditions (already supposed that it is facvar)
    2095 // output:
    2096 //   poly f2  where the non-null conditions in W have been dropped from f
    2097 static proc WPred(poly f, ideal W)
    2098 {
    2099   if (f==0){return(f);}
    2100   def l=factorize(f,2);
    2101   int i;
    2102   poly f1=1;
    2103   for(i=1;i<=size(l[1]);i++)
    2104   {
    2105     if (memberpos(l[1][i],W)[1]){;}
    2106     else{f1=f1*((l[1][i])^(l[2][i]));}
    2107   }
    2108   return(f1);
    2109 }
    2110 
    2111 //genimage
    2112 // ring @R
    2113 //input:
    2114 //   poly f1, idel N1,ideal W1,poly f2, ideal N2, ideal W2
    2115 //   corresponding to two polynomials having the same lpp
    2116 //   f1 in the redspec given by N1,W1,  f2 in the redspec given by N2,W2
    2117 //output:
    2118 //   the list of (ideal GG, list(list r1, list r2))
    2119 //   where g an ideal whose elements have the same lpp as f1 and f2
    2120 //   that specialize well to f1 in N1,W1 and to f2 in N2,W2.
    2121 //   If it doesn't exist a genimage, then g=ideal(0).
    2122 static proc genimage(poly f1, ideal N1, ideal W1, poly f2, ideal N2, ideal W2)
    2123 {
    2124   int i; ideal W12;  poly ff1; poly g1=0; ideal GG;
    2125   int tt=1;
    2126   // detect weather f1 reduces to 0 on segment 2
    2127   ff1=pnormalform(f1,N2,W2);
    2128   if (ff1==0)
    2129   {
    2130     // detect weather N1 is included in N2
    2131     def RR=basering;
    2132     setring @P;
    2133     def NP1=imap(RR,N1);
    2134     def NP2=imap(RR,N2);
    2135     attrib(NP2,"isSB",1);
    2136     poly nr;
    2137     i=1;
    2138     while ((tt) and (i<=size(NP1)))
    2139     {
    2140       nr=reduce(NP1[i],NP2);
    2141       if (nr!=0){tt=0;}
    2142       i++;
    2143     }
    2144     setring(RR);
    2145   }
    2146   else{tt=0;}
    2147   if (tt==1)
    2148   {
    2149     // detect weather W1 intersect W2 is non-empty
    2150     for (i=1;i<=size(W1);i++)
    2151     {
    2152       if (memberpos(W1[i],W2)[1])
    2153       {
    2154         W12[size(W12)+1]=W1[i];
    2155       }
    2156       else
    2157       {
    2158         if (nonnull(W1[i],N2,W2))
    2159         {
    2160           W12[size(W12)+1]=W1[i];
    2161         }
    2162       }
    2163     }
    2164     for (i=1;i<=size(W2);i++)
    2165     {
    2166       if (not(memberpos(W2[i],W12)[1]))
    2167       {
    2168         W12[size(W12)+1]=W2[i];
    2169       }
    2170     }
    2171   }
    2172   if (tt==1){g1=extendpoly(f1,N1,W12);}
    2173   if (g1!=0)
    2174   {
    2175     if (pnormalform(g1,N1,W1)==0)
    2176     {
    2177       GG=f1,g1;
    2178     }
    2179     else
    2180     {
    2181       GG=g1;
    2182     }
    2183     return(GG);
    2184   }
    2185 
    2186   // begins the second step;
    2187   int bound=6;
    2188   // in ring @R
    2189   int j; int g=0; int alpha; int r1; int s1=1; int s2=1;
    2190   poly G;
    2191   matrix qT;
    2192   matrix T;
    2193   ideal N10;
    2194   poly GT;
    2195   ideal N12=N1,N2;
    2196   def varx=maxideal(1);
    2197   int nx=size(varx);
    2198   poly pvarx=1;
    2199   for (i=1;i<=nx;i++){pvarx=pvarx*varx[i];}
    2200   def m=coef(43*f1+157*f2,pvarx);
    2201   def m1=coef(f1,pvarx);
    2202   def m2=coef(f2,pvarx);
    2203   list L1=adaptcoef(m,m1);
    2204   list L2=adaptcoef(m,m2);
    2205   ideal Tm=L1[1];
    2206   ideal c1=L1[2];
    2207   ideal c2=L2[2];
    2208   poly ww1;
    2209   poly ww2;
    2210   poly cA1;
    2211   poly cB1;
    2212   matrix TT;
    2213   poly H;
    2214   list r;
    2215   ideal q;
    2216   poly mu;
    2217   ideal N;
    2218 
    2219   // in ring @PAB
    2220   list Px=ringlist(@P);
    2221   list v="@A","@B";
    2222   Px[2]=Px[2]+v;
    2223   def npx=size(Px[3][1][2]);
    2224   Px[3][1][2]=1:(npx+size(v));
    2225   def @PAB=ring(Px);
    2226   setring(@PAB);
    2227 
    2228   poly PH;
    2229   ideal NP;
    2230   list rP;
    2231   def PN1=imap(@R,N1);
    2232   def PW1=imap(@R,W1);
    2233   def PN2=imap(@R,N2);
    2234   def PW2=imap(@R,W2);
    2235   def a1=imap(@R,c1);
    2236   def a2=imap(@R,c2);
    2237   matrix PT;
    2238   ideal PN;
    2239   ideal PN12=PN1,PN2;
    2240   PN=liftstd(PN12,PT);
    2241   list compos1;
    2242   list compos2;
    2243   list compos0;
    2244   intvec comp0;
    2245   poly w1=0;
    2246   poly w2=0;
    2247   poly h;
    2248   poly cA=0;
    2249   poly cB=0;
    2250   int t=0;
    2251   list l;
    2252   poly h1;
    2253   g=0;
    2254   while ((g<=bound) and not(t))
    2255   {
    2256     compos0=comp(g,2);
    2257     r1=1;
    2258     while ((r1<=size(compos0)) and not(t))
    2259     {
    2260       comp0=compos0[r1];
    2261       if (comp0[1]<=bound div 2)
    2262       {
    2263         compos1=comp(comp0[1],ncols(PW1));
    2264         s1=1;
    2265         while ((s1<=size(compos1)) and not(t))
    2266         {
    2267           if (comp0[2]<=bound div 2)
    2268           {
    2269             compos2=comp(comp0[2],ncols(PW2));
    2270             s2=1;
    2271             while ((s2<=size(compos2)) and not(t))
    2272             {
    2273               w1=WW(PW1,compos1[s1]);
    2274               w2=WW(PW2,compos2[s2]);
    2275               h=@A*w1*a1[1]-@B*w2*a2[1];
    2276               h=reduce(h,PN);
    2277               if (h==0){cA=1;cB=-1;}
    2278               else
    2279               {
    2280                 l=factorize(h,2);
    2281                 h1=1;
    2282                 for(i=1;i<=size(l[1]);i++)
    2283                 {
    2284                   if ((memberpos(@A,variables(l[1][i]))[1]) or  (memberpos(@B,variables(l[1][i]))[1]))
    2285                   {h1=h1*l[1][i];}
    2286                 }
    2287                 cA=diff(h1,@B);
    2288                 cB=diff(h1,@A);
    2289               }
    2290               if ((cA!=0) and (cB!=0) and (jet(cA,0)==cA) and (jet(cB,0)==cB))
    2291               {
    2292                 t=1;
    2293                 alpha=1;
    2294                 while((t) and (alpha<=ncols(a1)))
    2295                 {
    2296                   h=cA*w1*a1[alpha]+cB*w2*a2[alpha];
    2297                   if (not(reduce(h,PN,1)==0)){t=0;}
    2298                   alpha++;
    2299                 }
    2300               }
    2301               else{t=0;}
    2302               s2++;
    2303             }
    2304           }
    2305           s1++;
    2306         }
    2307       }
    2308       r1++;
    2309     }
    2310     g++;
    2311   }
    2312   setring(@R);
    2313   ww1=imap(@PAB,w1);
    2314   ww2=imap(@PAB,w2);
    2315   T=imap(@PAB,PT);
    2316   N=imap(@PAB,PN);
    2317   cA1=imap(@PAB,cA);
    2318   cB1=imap(@PAB,cB);
    2319   if (t)
    2320   {
    2321     G=0;
    2322     for (alpha=1;alpha<=ncols(Tm);alpha++)
    2323     {
    2324       H=cA1*ww1*c1[alpha]+cB1*ww2*c2[alpha];
    2325       setring(@PAB);
    2326       PH=imap(@R,H);
    2327       PN=imap(@R,N);
    2328       rP=division(PH,PN);
    2329       setring(@R);
    2330       r=imap(@PAB,rP);
    2331       if (r[2][1]!=0){ERROR("the division is not null and it should be");}
    2332       q=r[1];
    2333       qT=transpose(matrix(q));
    2334       N10=N12;
    2335       for (i=size(N1)+1;i<=size(N1)+size(N2);i++){N10[i]=0;}
    2336       G=G+(cA1*ww1*c1[alpha]-(matrix(N10)*T*qT)[1,1])*Tm[alpha];
    2337     }
    2338     GG=ideal(G);
    2339   }
    2340   else{GG=ideal(0);}
    2341   return(GG);
    2342 }
    2343 
    2344 // purpose: given a polynomial f (in the reduced basis)
    2345 //          the null-conditions ideal N in the segment
    2346 //          end the set of non-null polynomials common to the segment and
    2347 //          a new segment,
    2348 //          to obtain an equivalent polynomial with a leading coefficient
    2349 //          that is non-null in the second segment.
    2350 // input:
    2351 // poly f:    a polynomials of the reduced basis in the segment (N,W)
    2352 // ideal N:   the null-conditions ideal in the segment
    2353 // ideal W12: the set of non-null polynomials common to the segment and
    2354 //            a second segment
    2355 static proc extendpoly(poly f, ideal N, ideal W12)
    2356 {
    2357   int bound=4;
    2358   ideal cfs;
    2359   ideal cfsn;
    2360   ideal ppfs;
    2361   poly p=f;
    2362   poly fn;
    2363   poly lm; poly lc;
    2364   int tt=0;
    2365   int i;
    2366   while (p!=0)
    2367   {
    2368     lm=leadmonom(p);
    2369     lc=leadcoef(p);
    2370     cfs[size(cfs)+1]=lc;
    2371     ppfs[size(ppfs)+1]=lm;
    2372     p=p-lc*lm;
    2373   }
    2374   def lcf=cfs[1];
    2375   int r1=0; int s1;
    2376   def RR=basering;
    2377   setring @P;
    2378   list compos1;
    2379   poly w1;
    2380   ideal q;
    2381   def lcfp=imap(RR,lcf);
    2382   def W=imap(RR,W12);
    2383   def Np=imap(RR,N);
    2384   def cfsp=imap(RR,cfs);
    2385   ideal cfspn;
    2386   matrix T;
    2387   ideal H=lcfp,Np;
    2388   def G=liftstd(H,T);
    2389   list r;
    2390   while ((r1<=bound) and not(tt))
    2391   {
    2392     compos1=comp(r1,ncols(W));
    2393     s1=1;
    2394     while ((s1<=size(compos1)) and not(tt))
    2395     {
    2396       w1=WW(W,compos1[s1]);
    2397       cfspn=ideal(0);
    2398       cfspn[1]=w1;
    2399       tt=1;
    2400       i=2;
    2401       while ((i<=size(cfsp)) and (tt))
    2402       {
    2403         r=division(w1*cfsp[i],G);
    2404         if (r[2][1]!=0){tt=0;}
    2405         else
    2406         {
    2407           q=r[1];
    2408           cfspn[i]=(T*transpose(matrix(q)))[1,1];
    2409         }
    2410         i++;
    2411       }
    2412       s1++;
    2413     }
    2414     r1++;
    2415   }
    2416   setring RR;
    2417   if (tt)
    2418   {
    2419     cfsn=imap(@P,cfspn);
    2420     fn=0;
    2421     for (i=1;i<=size(ppfs);i++)
    2422     {
    2423       fn=fn+cfsn[i]*ppfs[i];
    2424     }
    2425   }
    2426   else{fn=0;}
    2427   return(fn);
    2428 }
    2429 
    2430 // nonnull
    2431 // ring @P (or @R)
    2432 // input:
    2433 //   poly f
    2434 //   ideal N
    2435 //   ideal W
    2436 // output:
    2437 //   1 if f is nonnull in the segment (N,W)
    2438 //   0 if it can be zero
    2439 static proc nonnull(poly f, ideal N, ideal W)
    2440 {
    2441   int tt;
    2442   ideal N0=N;
    2443   N0[size(N0)+1]=f;
    2444   poly h=1;
    2445   int i;
    2446   for (i=1;i<=size(W);i++){h=h*W[i];}
    2447   def RR=basering;
    2448   setring(@P);
    2449   list Px=ringlist(@P);
    2450   list v="@C";
    2451   Px[2]=Px[2]+v;
    2452   def npx=size(Px[3][1][2]);
    2453   Px[3][1][1]="dp";
    2454   Px[3][1][2]=1:(npx+size(v));
    2455   def @PC=ring(Px);
    2456   setring(@PC);
    2457   def N1=imap(RR,N0);
    2458   def h1=imap(RR,h);
    2459   ideal G=1-@C*h1;
    2460   G=G+N1;
    2461   option(redSB);
    2462   ideal G1=std(G);
    2463   if (G1[1]==1){tt=1;} else{tt=0;}
    2464   setring(RR);
    2465   return(tt);
    2466 }
    2467 
    2468 // decide
    2469 // input:
    2470 //   given two corresponding polynomials g1 and g2 with the same lpp
    2471 //   g1 belonging to the basis in the segment N1,W1
    2472 //   g2 belonging to the basis in the segment N2,W2
    2473 // output:
    2474 //   an ideal (with a single polynomial or more if a sheaf is needed)
    2475 //   that specializes well on both segments to g1 and g2 respectivelly.
    2476 //   If ideal(0) is output, then no such polynomial nor sheaf exists.
    2477 static proc decide(poly g1, ideal N1, ideal W1, poly g2, ideal N2, ideal W2)
    2478 {
    2479   poly S;
    2480   poly S1;
    2481   poly S2;
    2482   S=leadcoef(g2)*g1-leadcoef(g1)*g2;
    2483   def RR=basering;
    2484   setring(@RP);
    2485   def SR=imap(RR,S);
    2486   def N1R=imap(RR,N1);
    2487   def N2R=imap(RR,N2);
    2488   attrib(N1R,"isSB",1);
    2489   attrib(N2R,"isSB",1);
    2490   poly S1R=reduce(SR,N1R);
    2491   poly S2R=reduce(SR,N2R);
    2492   setring(RR);
    2493   S1=imap(@RP,S1R);
    2494   S2=imap(@RP,S2R);
    2495   if ((S2==0) and (nonnull(leadcoef(g1),N2,W2))){return(ideal(g1));}
    2496   if ((S1==0) and (nonnull(leadcoef(g2),N1,W1))){return(ideal(g2));}
    2497   if ((S1==0) and (S2==0))
    2498   {
    2499     return(ideal(g1,g2));
    2500   }
    2501   return(ideal(genimage(g1,N1,W1,g2,N2,W2)));
    2502 }
    2503 
    2504 // input:  the tree (list) from buildtree output
    2505 // output: the list of terminal vertices.
    2506 static proc finalcases(list T)
    2507 //"USAGE:   finalcases(T);
    2508 //          T is the list provided by buildtree
    2509 //RETURN:   A list with the CGS determined by buildtree.
    2510 //          Each element of the list represents one segment
    2511 //          of the terminal vertices of buildtree givieng the CGS.
    2512 //          The list elements have the following structure:
    2513 //           [1]: label (an intvec(1,0,..)) that indicates the position
    2514 //                in the buildtree but that is irrelevant for the CGS
    2515 //           [2]: 1 (integer) it is also irrelevant and indicates
    2516 //                that this was a terminal vertex in buildtree.
    2517 //           [3]: the reduced basis of the segment.
    2518 //           [4], [5], [6]: the red-representation of the segment
    2519 //                [4] are the null-conditions radical ideal N,
    2520 //                [5] are the non-null polynomials set (ideal) W,
    2521 //                [6] is the set of prime components (ideals) of N.
    2522 //           [7]: is the set of lpp
    2523 //           [8]: poly 1 (irrelevant) is the condition to branch (but no
    2524 //                more branch is necessary in the discussion, so 1 is the result.
    2525 //NOTE:     It can be called having as argument the list output by buildtree
    2526 //KEYWORDS: buildtree, buildtreetoMaple, CGS
    2527 //EXAMPLE:  finalcases; shows an example"
    2528 {
    2529   int i;
    2530   list L;
    2531   for (i=1;i<=size(T);i++)
    2532   {
    2533     if (T[i][2])
    2534     {L[size(L)+1]=T[i];}
    2535   }
    2536   return(L);
    2537 }
    2538 //example
    2539 //{ "EXAMPLE:"; echo = 2;
    2540 //  ring R=(0,a1,a2,a3,a4),(x1,x2,x3,x4),dp;
    2541 //  ideal F=x4-a4+a2, x1+x2+x3+x4-a1-a3-a4, x1*x3*x4-a1*a3*a4, x1*x3+x1*x4+x2*x3+x3*x4-a1*a4-a1*a3-a3*a4;
    2542 //  def T=buildtree(F);
    2543 //  setglobalrings();
    2544 //  finalcases(T);
    2545 //}
    2546 
    2547 // input:  the list of terminal vertices of buildtree (output of finalcases)
    2548 // output: the same terminal vertices grouped by lpp
    2549 static proc groupsegments(list T)
    2550 {
    2551   int i;
    2552   list L;
    2553   list lpp;
    2554   list lp;
    2555   list ls;
    2556   int n=size(T);
    2557   lpp[1]=T[n][7];
    2558   L[1]=list(lpp[1],list(list(T[n][1],T[n][3],T[n][4],T[n][5],T[n][6])));
    2559   if (n>1)
    2560   {
    2561     for (i=1;i<=size(T)-1;i++)
    2562     {
    2563       lp=memberpos(T[n-i][7],lpp);
    2564       if(lp[1]==1)
    2565       {
    2566         ls=L[lp[2]][2];
    2567         ls[size(ls)+1]=list(T[n-i][1],T[n-i][3],T[n-i][4],T[n-i][5],T[n-i][6]);
    2568         L[lp[2]][2]=ls;
    2569       }
    2570       else
    2571       {
    2572         lpp[size(lpp)+1]=T[n-i][7];
    2573         L[size(L)+1]=list(T[n-i][7],list(list(T[n-i][1],T[n-i][3],T[n-i][4],T[n-i][5],T[n-i][6])));
    2574       }
    2575     }
    2576   }
    2577   //"L in groupsegments="; L;
    2578   return(L);
    2579 }
    2580 
    2581 // eliminates repeated elements form an ideal
    2582 static proc elimrepeated(ideal F)
    2583 {
    2584   int i;
    2585   int j;
    2586   ideal FF;
    2587   FF[1]=F[1];
    2588   for (i=2;i<=ncols(F);i++;)
    2589   {
    2590     if (not(memberpos(F[i],FF)[1]))
    2591     {
    2592       FF[size(FF)+1]=F[i];
    2593     }
    2594   }
    2595   return(FF);
    2596 }
    2597 
    2598 // decide F is the same as decide but allows as first element a sheaf F
    2599 static proc decideF(ideal F,ideal N,ideal W, poly f2, ideal N2, ideal W2)
    2600 {
    2601   int i;
    2602   ideal G=F;
    2603   ideal g;
    2604   if (ncols(F)==1) {return(decide(F[1],N,W,f2,N2,W2));}
    2605   for (i=1;i<=ncols(F);i++)
    2606   {
    2607     G=G+decide(F[i],N,W,f2,N2,W2);
    2608   }
    2609   return(elimrepeated(G));
    2610 }
    2611 
    2612 // newredspec
    2613 // input:  two redspec in the form of N,W and Nj,Wj
    2614 // output: a redspec representing the minimal redspec segment that contains
    2615 //         both input segments.
    2616 static proc newredspec(ideal N,ideal W, ideal Nj, ideal Wj)
    2617 {
    2618   ideal nN;
    2619   ideal nW;
    2620   int u;
    2621   def RR=basering;
    2622   setring(@P);
    2623   list r;
    2624   def Np=imap(RR,N);
    2625   def Wp=imap(RR,W);
    2626   def Njp=imap(RR,Nj);
    2627   def Wjp=imap(RR,Wj);
    2628   Np=intersect(Np,Njp);
    2629   ideal WR;
    2630   for(u=1;u<=size(Wjp);u++)
    2631   {
    2632     if(nonnull(Wjp[u],Np,Wp)){WR[size(WR)+1]=Wjp[u];}
    2633   }
    2634   for(u=1;u<=size(Wp);u++)
    2635   {
    2636     if((not(memberpos(Wp[u],WR)[1])) and (nonnull(Wp[u],Njp,Wjp)))
    2637     {
    2638       WR[size(WR)+1]=Wp[u];
    2639     }
    2640   }
    2641   r=redspec(Np,WR);
    2642   option(redSB);
    2643   Np=std(r[1]);
    2644   Wp=r[2];
    2645   setring(RR);
    2646   nN=imap(@P,Np);
    2647   nW=imap(@P,Wp);
    2648   return(list(nN,nW));
    2649 }
    2650 
    2651 // selectcases
    2652 // input:
    2653 //   list bT: the list output by buildtree.
    2654 // output:
    2655 //   list L   it contins the list of segments allowing a common
    2656 //            reduced basis. The elements of L are of the form
    2657 //            list (lpp,B,list(list(N,W,L),..list(N,W,L)) )
    2658 static proc selectcases(list bT)
    2659 {
    2660   list T=groupsegments(finalcases(bT));
    2661   //NEW
    2662   //groupredtocan(T);
    2663   list T0=bT[1];
    2664              // first element of the list of buildtree
    2665   list TT0;
    2666   TT0[1]=list(T0[7],T0[3],list(list(T0[4],T0[5],T0[6])));
    2667              // first element of the output of selectcases
    2668   list T1=T; // the initial list; it is only actualized (split)
    2669              // when a segment is completly revised (all split are
    2670              // already be considered);
    2671              // ( (lpp, ((lab,B,N,W,L),.. ()) ), .. (..) )
    2672   list TT;   // the output list ( (lpp,B,((N,W,L),..()) ),.. (..) )
    2673   // case i
    2674   list S1;   // the segments in case i T1[i][2]; ( (lab,B,N,W,L),..() )
    2675   list S2;   // the segments in case i that are being summarized in
    2676              // actual segment ( (N,W,L),..() )
    2677   list S3;   // the segments in case i that cannot be summarized in
    2678              // the actual case. When the case is finished a new case
    2679              // is created with them ( (lab,B,N,W,L),..() )
    2680   list s3;   // list of integers s whose segment cannot be summarized
    2681              // in the actual case
    2682   ideal lpp; // the summarized lpp (can contain repetitions)
    2683   ideal lppi;// in process of sumarizing lpp (can contain repetitions)
    2684   ideal B;   // the summarized B (can contain polynomials with
    2685              // the same lpp (sheaves))
    2686   ideal Bi;  // in process of summarizing B (can contain polynomials with
    2687              // the same lpp (sheaves))
    2688   ideal N;   // the summarized N
    2689   ideal W;   // the summarized W
    2690   ideal F;   // the summarized poly j (can contain a sheaf instead of
    2691              // a single poly)
    2692   ideal FF;  // the same as F but it can be ideal(0)
    2693   poly lpj;
    2694   poly fj;
    2695   ideal Nj;
    2696   ideal Wj;
    2697   ideal G;
    2698   int i;     // the index of the case i in T1;
    2699   int j;     // the index of the polynomial j of the basis
    2700   int s;     // the index of the segment s in S1;
    2701   int u;
    2702   int tests; // true if al the polynomial in segment s have been generalized;
    2703   list r;
    2704   // initializing the new list
    2705   i=1;
    2706   while(i<=size(T1))
    2707   {
    2708     S1=T1[i][2]; // ((lab,B,N,W,L)..) of the segments in case i
    2709     if (size(S1)==1)
    2710     {
    2711       TT[i]=list(T1[i][1],S1[1][2],list(list(S1[1][3],S1[1][4],S1[1][5])));
    2712     }
    2713     else
    2714     {
    2715       S2=list();
    2716       S3=list(); // ((lab,B,N,W,L)..) of the segments in case i to
    2717                  // create another segment i+1
    2718       s3=list();
    2719       B=S1[1][2];
    2720       Bi=ideal(0);
    2721       lpp=T1[i][1];
    2722       j=1;
    2723       tests=1;
    2724       while (j<=size(S1[1][2]))
    2725       { // j desings the new j-th polynomial
    2726         N=S1[1][3];
    2727         W=S1[1][4];
    2728         F=ideal(S1[1][2][j]);
    2729         s=2;
    2730         while (s<=size(S1) and not(memberpos(s,s3)[1]))
    2731         { // s desings the new segment s
    2732           fj=S1[s][2][j];
    2733           Nj=S1[s][3];
    2734           Wj=S1[s][4];
    2735           FF=decideF(F,N,W,fj,Nj,Wj);
    2736           if (FF[1]==0)
    2737           {
    2738             if (@ish)
    2739             {
    2740               "Warning: Dealing with an homogeneous ideal";
    2741               "mrcgs was not able to summarize all lpp cases into a single segment";
    2742               "Please send a mail with your Problem to antonio.montes@upc.edu";
    2743               "You found a counterexample of the complete success of the actual mrcgs algorithm";
    2744               //NEW
    2745               "f1:"; F; "N1:"; N; "W1:"; W; "f2:"; fj; "N2:"; Nj; "W2:"; Wj;
    2746             }
    2747             S3[size(S3)+1]=S1[s];
    2748             s3[size(s3)+1]=s;
    2749             tests=0;
    2750           }
    2751           else
    2752           {
    2753             F=FF;
    2754             lpj=leadmonom(fj);
    2755             r=newredspec(N,W,Nj,Wj);
    2756             N=r[1];
    2757             W=r[2];
    2758           }
    2759           s++;
    2760         }
    2761         if (Bi[1]==0){Bi=FF;}
    2762         else
    2763         {
    2764           Bi=Bi+FF;
    2765         }
    2766         j++;
    2767       }
    2768       if (tests)
    2769       {
    2770         B=Bi;
    2771         lpp=ideal(0);
    2772         for (u=1;u<=size(B);u++){lpp[u]=leadmonom(B[u]);}
    2773       }
    2774       for (s=1;s<=size(T1[i][2]);s++)
    2775       {
    2776         if (not(memberpos(s,s3)[1]))
    2777         {
    2778           S2[size(S2)+1]=list(S1[s][3],S1[s][4],S1[s][5]);
    2779         }
    2780       }
    2781       TT[i]=list(lpp,B,S2);
    2782       // for (s=1;s<=size(s3);s++){S1=delete(S1,s);}
    2783       T1[i][2]=S2;
    2784       if (size(S3)>0){T1=insert(T1,list(T1[i][1],S3),i);}
    2785     }
    2786     i++;
    2787   }
    2788   for (i=1;i<=size(TT);i++){TT0[i+1]=TT[i];}
    2789   return(TT0);
    2790 }
    2791 
    2792 //*****************End of Selectcases**************************
    2793 
    2794 //*****************Begin of CanTree****************************
    2795 
    2796 // equalideals
    2797 // input: 2 ideals F and G;
    2798 // output: 1 if they are identical (the same polynomials in the same order)
    2799 //         0 else
    2800 static proc equalideals(ideal F, ideal G)
    2801 {
    2802   int i=1; int t=1;
    2803   if (size(F)!=size(G)){return(0);}
    2804   while ((i<=size(F)) and (t))
    2805   {
    2806     if (F[i]!=G[i]){t=0;}
    2807     i++;
    2808   }
    2809   return(t);
    2810 }
    2811 
    2812 // delintvec
    2813 // input: intvec V
    2814 //        int i
    2815 // output:
    2816 //        intvec W (equal to V but the coordinate i is deleted
    2817 static proc delintvec(intvec V, int i)
    2818 {
    2819   int j;
    2820   intvec W;
    2821   for (j=1;j<i;j++){W[j]=V[j];}
    2822   for (j=i+1;j<=size(V);j++){W[j-1]=V[j];}
    2823   return(W);
    2824 }
    2825 
    2826 // redtocanspec
    2827 // Computes the canonical representation of a redspec (N,W,L).
    2828 // input:
    2829 //    ideal N (null conditions, must be radical)
    2830 //    ideal W (non-null conditions ideal)
    2831 //    list L  must contain the radical decomposition of N.
    2832 // output:
    2833 //    the list of elements of the (ideal N1,list(ideal M11,..,ideal M1k))
    2834 //    determining the canonical representation of the difference of
    2835 //    V(N) \ V(h), where h=prod(w in W).
    2836 static proc redtocanspec(intvec lab, int child, list rs)
    2837 {
    2838   ideal N=rs[1]; ideal W=rs[2]; list L=rs[3];
    2839   intvec labi; intvec labij;
    2840   int childi;
    2841   int i; int j; list L0;
    2842   L0[1]=list(lab,size(L));
    2843   if (W[1]==0)
    2844   {
    2845     for (i=1;i<=size(L);i++)
    2846     {
    2847       labi=lab,child+i;
    2848       L0[size(L0)+1]=list(labi,1,L[i]);
    2849       labij=labi,1;
    2850       L0[size(L0)+1]=list(labij,0,ideal(1));
    2851     }
    2852     return(L0);
    2853   }
    2854   if (N[1]==1)
    2855   {
    2856     L0[1]=list(lab,1);
    2857     labi=lab,child+1;
    2858     L0[size(L0)+1]=list(labi,1,ideal(1));
    2859     labij=labi,1;
    2860     L0[size(L0)+1]=list(labij,0,ideal(1));
    2861   }
    2862   def RR=basering;
    2863   setring(@P);
    2864   ideal Np=imap(RR,N);
    2865   ideal Wp=imap(RR,W);
    2866   poly h=1;
    2867   for (i=1;i<=size(Wp);i++){h=h*Wp[i];}
    2868   list Lp=imap(RR,L);
    2869   list r; list Ti; list LL;
    2870   LL[1]=list(lab,size(Lp));
    2871   for (i=1;i<=size(Lp);i++)
    2872   {
    2873     Ti=minGTZ(Lp[i]+h);
    2874     for(j=1;j<=size(Ti);j++)
    2875     {
    2876       option(redSB);
    2877       Ti[j]=std(Ti[j]);
    2878     }
    2879     labi=lab,child+i;
    2880     childi=size(Ti);
    2881     LL[size(LL)+1]=list(labi,childi,Lp[i]);
    2882     for (j=1;j<=childi;j++)
    2883     {
    2884       labij=labi,j;
    2885       LL[size(LL)+1]=list(labij,0,Ti[j]);
    2886     }
    2887   }
    2888   LL[1]=list(lab,size(Lp));
    2889   setring(RR);
    2890   return(imap(@P,LL));
    2891 }
    2892 
    2893 // difftocanspec
    2894 // Computes the canonical representation of a diffspec V(N) \ V(M)
    2895 // input:
    2896 //    intvec lab: label where to hang the canspec
    2897 //    list  N ideal of null conditions.
    2898 //    ideal M ideal of the variety to be substacted
    2899 // output:
    2900 //    the list of elements determining the canonical representation of
    2901 //    the difference  V(N) \ V(M):
    2902 //      ( (intvec(i),children), ...(lab, children, prime ideal),...)
    2903 static proc difftocanspec(intvec lab, int child, ideal N, ideal M)
    2904 {
    2905   int i; int j; list LLL;
    2906   def RR=basering;
    2907   setring(@P);
    2908   ideal Np=imap(RR,N);
    2909   ideal Mp=imap(RR,M);
    2910   def L=minGTZ(Np);
    2911   for(j=1;j<=size(L);j++)
    2912   {
    2913     option(redSB);
    2914     L[j]=std(L[j]);
    2915   }
    2916   intvec labi; intvec labij;
    2917   int childi;
    2918   list LL;
    2919   if ((Mp[1]==0) or ((size(L)==1) and (L[1][1]==1)))
    2920   {
    2921     //LL[1]=list(lab,1);
    2922     //labi=lab,1;
    2923     //LL[2]=list(labi,1,ideal(1));
    2924     //labij=labi,1;
    2925     //LL[3]=list(labij,0,ideal(1));
    2926     setring(RR);
    2927     return(LLL);
    2928   }
    2929   list r; list Ti;
    2930   def k=0;
    2931   LL[1]=list(lab,0);
    2932   for (i=1;i<=size(L);i++)
    2933   {
    2934     Ti=minGTZ(L[i]+Mp);
    2935     for(j=1;j<=size(Ti);j++)
    2936     {
    2937       option(redSB);
    2938       Ti[j]=std(Ti[j]);
    2939     }
    2940     if (not((size(Ti)==1) and (equalideals(L[i],Ti[1]))))
    2941     {
    2942       k++;
    2943       labi=lab,child+k;
    2944       childi=size(Ti);
    2945       LL[size(LL)+1]=list(labi,childi,L[i]);
    2946       for (j=1;j<=childi;j++)
    2947       {
    2948         labij=labi,j;
    2949         LL[size(LL)+1]=list(labij,0,Ti[j]);
    2950       }
    2951     }
    2952     else{setring(RR); return(LLL);}
    2953   }
    2954   if (size(LL)>0)
    2955   {
    2956     LL[1]=list(lab,k);
    2957     setring(RR);
    2958     return(imap(@P,LL));
    2959   }