source: git/Singular/LIB/primdec.lib @ d2b2a7

// $Id: primdec.lib,v 1.13 1998-05-05 11:55:34 krueger Exp $
///////////////////////////////////////////////////////
// primdec.lib
// algorithms for primary decomposition based on
// the ideas of Gianni,Trager,Zacharias
// written by Gerhard Pfister
//
// algorithms for primary decomposition based on
// the ideas of Shimoyama/Yokoyama
// written by Wolfram Decker and Hans Schoenemann
//////////////////////////////////////////////////////

version="$Id: primdec.lib,v 1.13 1998-05-05 11:55:34 krueger Exp $";
info="
LIBRARY: primdec.lib: PROCEDURE FOR PRIMARY DECOMPOSITION (I)

  minAssPrimes (ideal I)
  //computes the minimal associated primes of the ideal I

  primdecGTZ (ideal I)
  // Computes a complete primary decomposition via Gianni,Trager,Zacharias

  radical(ideal I)
  //computes the radical of the ideal I

  equiRadical(ideal I)
  //computes the radical of the equidimensional part of the ideal I

  prepareAss(ideal I)
  //computes the radicals of the equidimensional parts of I

  min_ass_prim_charsets (ideal I, int choose)
  // minimal associated primes via the  characteristic set
  // package written by Michael Messollen.
  // The integer choose must be either 0 or 1.
  // If choose=0, the given ordering of the variables is used.
  // If choose=1, the system tries to find an \"optimal ordering\",
  // which in some cases may considerably speed up the algorithm

  // You may also may want to try one of the algorithms for
  // minimal associated primes in the the library
  // primdec.lib, written by Gerhard Pfister.
  // These algorithms are variants of the algorithm
  // of Gianni-Trager-Zacharias

  primdecSY (ideal I, int choose)

  // Computes a complete primary decomposition via
  // a variant of the pseudoprimary approach of
  // Shimoyama-Yokoyama.
  // The integer choose must be either 0, 1, 2 or 3.
  // If choose=0, min_ass_prim_charsets with the given
  // ordering of the variables is used.
  // If choose=1, min_ass_prim_charsets with the \"optimized\"
  // ordering of the variables is used.
  // If choose=2, minAssPrimes from primdec.lib is used
  // If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
";

LIB "general.lib";
LIB "elim.lib";
LIB "poly.lib";
LIB "random.lib";
///////////////////////////////////////////////////////////////////////////////

proc sat1 (ideal id, poly p)
"USAGE:   sat1(id,j);  id ideal, j polynomial
RETURN:  saturation of id with respect to j (= union_(k=1...) of id:j^k)
NOTE:    result is a std basis in the basering
EXAMPLE: example sat; shows an example
"
{
   int @k;
   ideal inew=std(id);
   ideal iold;
   option(returnSB);
   while(specialIdealsEqual(iold,inew)==0 )
   {
      iold=inew;
      inew=quotient(iold,p);
      @k++;
   }
   @k--;
   option(noreturnSB);
   list L =inew,p^@k;
   return (L);
}

///////////////////////////////////////////////////////////////////////////////

proc sat2 (ideal id, ideal h)
"USAGE:   sat2(id,j);  id ideal, j polynomial
RETURN:  saturation of id with respect to j (= union_(k=1...) of id:j^k)
NOTE:    result is a std basis in the basering
EXAMPLE: example sat2; shows an example
"
{
   int @k,@i;
   def @P= basering;
   if(ordstr(basering)[1,2]!="dp")
   {
      execute "ring @Phelp=("+charstr(@P)+"),("+varstr(@P)+"),(C,dp);";
      ideal inew=std(imap(@P,id));
      ideal  @h=imap(@P,h);
   }
   else
   {
      ideal @h=h;
      ideal inew=std(id);
   }
   ideal fac;

   for(@i=1;@i<=ncols(@h);@i++)
   {
     if(deg(@h[@i])>0)
     {
        fac=fac+factorize(@h[@i],1);
     }
   }
   fac=simplify(fac,4);
   poly @f=1;
   if(deg(fac[1])>0)
   {
      ideal iold;

      for(@i=1;@i<=size(fac);@i++)
      {
        @f=@f*fac[@i];
      }
      option(returnSB);
      while(specialIdealsEqual(iold,inew)==0 )
      {
         iold=inew;
         if(deg(iold[size(iold)])!=1)
         {
            inew=quotient(iold,@f);
         }
         else
         {
            inew=iold;
         }
         @k++;
      }
      option(noreturnSB);
      @k--;
   }

   if(ordstr(@P)[1,2]!="dp")
   {
      setring @P;
      ideal inew=std(imap(@Phelp,inew));
      poly @f=imap(@Phelp,@f);
   }
   list L =inew,@f^@k;
   return (L);
}

///////////////////////////////////////////////////////////////////////////////

proc minSat(ideal inew, ideal h)
{
   int i,k;
   poly f=1;
   ideal iold,fac;
   list quotM,l;

   for(i=1;i<=ncols(h);i++)
   {
      if(deg(h[i])>0)
      {
         fac=fac+factorize(h[i],1);
      }
   }
   fac=simplify(fac,4);
   if(size(fac)==0)
   {
      l=inew,1;
      return(l);
   }
   fac=sort(fac)[1];
   for(i=1;i<=size(fac);i++)
   {
      f=f*fac[i];
   }
   quotM[1]=inew;
   quotM[2]=fac;
   quotM[3]=f;
   f=1;
   option(returnSB);
   while(specialIdealsEqual(iold,quotM[1])==0)
   {
      if(k>0)
      {
         f=f*quotM[3];
      }
      iold=quotM[1];
      quotM=quotMin(quotM);
      k++;
   }
   option(noreturnSB);
   l=quotM[1],f;
   return(l);
}

proc quotMin(list tsil)
{
   int i,j,k,action;
   ideal verg;
   list l;
   poly g;

   ideal laedi=tsil[1];
   ideal fac=tsil[2];
   poly f=tsil[3];
   
   ideal star=quotient(laedi,f);
   action=1;

   while(action==1)
   {
      if(size(fac)==1)
      {
         action=0;
         break;
      }
      for(i=1;i<=size(fac);i++)
      {
        g=1;
        verg=laedi;
        
         for(j=1;j<=size(fac);j++)
         {
            if(i!=j)
            {
               g=g*fac[j];
            }
         }
         verg=quotient(laedi,g);
         
         if(specialIdealsEqual(verg,star)==1)
         {
            f=g;
            fac[i]=0;
            fac=simplify(fac,2);
            break;
         }
         if(i==size(fac))
         {
            action=0;
         }
      }
   }
   l=star,fac,f;   
   return(l);
}

////////////////////////////////////////////////////////////////////////////////
proc testFactor(list act,poly p)
{
  poly keep=p;
  
   int i;
   poly q=act[1][1]^act[2][1];
   for(i=2;i<=size(act[1]);i++)
   {
      q=q*act[1][i]^act[2][i];
   }
   q=1/leadcoef(q)*q;
   p=1/leadcoef(p)*p;
   if(p-q!=0)
   {
      "ERROR IN FACTOR";
      basering;
      
      act;
      keep;
      pause;
      
      p;
      q;
      pause;
   }
}
////////////////////////////////////////////////////////////////////////////////

proc factor(poly p)
"USAGE:   factor(p) p poly
RETURN:  list=;
NOTE:
EXAMPLE: example factor; shows an example
"
{

  ideal @i;
  list @l;
  intvec @v,@w;
  int @j,@k,@n;

  if(deg(p)<=1)
  {
     @i=ideal(p);
     @v=1;
     @l[1]=@i;
     @l[2]=@v;
     return(@l);
  }
  if (size(p)==1)
  {
     @w=leadexp(p);
     for(@j=1;@j<=nvars(basering);@j++)
     {
        if(@w[@j]!=0)
        {
           @k++;
           @v[@k]=@w[@j];
           @i=@i+ideal(var(@j));
        }
     }
     @l[1]=@i;
     @l[2]=@v;
     return(@l);
  }
//  @l=factorize(p,2);
  @l=factorize(p);
  // if(npars(basering)>0)
  // {
     for(@j=1;@j<=size(@l[1]);@j++)
     {
        if(deg(@l[1][@j])==0)
        {
           @n++;
        }
     }
     if(@n>0)
     {
        if(@n==size(@l[1]))
        {
           @l[1]=ideal(1);
           @v=1;
           @l[2]=@v;
        }
        else
        {
           @k=0;
           int pleh;
           for(@j=1;@j<=size(@l[1]);@j++)
           {
              if(deg(@l[1][@j])!=0)
              {
                 @k++;
                 @i=@i+ideal(@l[1][@j]);
                 if(size(@i)==pleh)
                 {
"factorization error";
@l;
                    @k--;
                    @v[@k]=@v[@k]+@l[2][@j];
                 }
                 else
                 {
                    pleh++;
                    @v[@k]=@l[2][@j];
                 }
              }
           }
           @l[1]=@i;
           @l[2]=@v;
        }
     }
     // }
  return(@l);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 0,(x,y,z),lp;
   poly  p = (x+y)^2*(y-z)^3;
   list  l = factor(p);
   l;
   ring r1 =(0,b,d,f,g),(a,c,e),lp;
   poly p  =(1*d)*e^2+(1*d*f^2*g);
   list  l = factor(p);
   l;
   ring r2 =(0,b,f,g),(a,c,e,d),lp;
   poly p  =(1*d)*e^2+(1*d*f^2*g);
   list  l = factor(p);
   l;

}



////////////////////////////////////////////////////////////////////////////////

proc idealsEqual( ideal k, ideal j)
{
   return(stdIdealsEqual(std(k),std(j)));
}

proc specialIdealsEqual( ideal k1, ideal k2)
{
   int j;

   if(size(k1)==size(k2))
   {
      for(j=1;j<=size(k1);j++)
      {
         if(leadexp(k1[j])!=leadexp(k2[j]))
         {
            return(0);
         }
      }
      return(1);
   }
   return(0);
}

proc stdIdealsEqual( ideal k1, ideal k2)
{
   int j;

   if(size(k1)==size(k2))
   {
      for(j=1;j<=size(k1);j++)
      {
         if(leadexp(k1[j])!=leadexp(k2[j]))
         {
            return(0);
         }
      }
      attrib(k2,"isSB",1);
      if(size(reduce(k1,k2,1))==0)
      {
         return(1);
      }
   }
   return(0);
}


////////////////////////////////////////////////////////////////////////////////

proc testPrimary(list pr, ideal k)
"USAGE:   testPrimary(pr,k) pr list, k ideal;
RETURN:  int = 1, if the intersection of the ideals in pr is k, 0 if not
NOTE:
EEXAMPLE: example testPrimary ; shows an example
"
{
   int i;
   ideal j=pr[1];
   for (i=2;i<=size(pr)/2;i++)
   {
       j=intersect(j,pr[2*i-1]);
   }
   return(idealsEqual(j,k));
}
example
{ "EXAMPLE:"; echo = 2;
   ring  s = 0,(x,y,z),lp;
   ideal i=x3-x2-x+1,xy-y;
   ideal i1=x-1;
   ideal i2=x-1;
   ideal i3=y,x2-2x+1;
   ideal i4=y,x-1;
   ideal i5=y,x+1;
   ideal i6=y,x+1;
   list pr=i1,i2,i3,i4,i5,i6;
   testPrimary(pr,i);
   pr[5]=y+1,x+1;
   testPrimary(pr,i);
}

////////////////////////////////////////////////////////////////////////////////
proc printPrimary( list l, list #)
"USAGE:   printPrimary(l) l list;
RETURN:  nothing
NOTE:
EXAMPLE: example printPrimary; shows an example
"
{
   if(size(#)>0)
   {
      "                                            ";
      "  The primary decomposition of the ideal    ";
      #[1];
      "                                            ";
      "  is:                                       ";
      "                                            ";
   }
   int k;
   for (k=1;k<=size(l)/2;k=k+1)
   {
      "                                            ";
      "primary ideal:                              ";
      l[2*k-1];
      "                                            ";
      "associated prime ideal                      ";
      l[2*k];
      "                                            ";
   }
}
example
{ "EXAMPLE:"; echo = 2;
}

////////////////////////////////////////////////////////////////////////////////


proc randomLast(int b)
"USAGE:   randomLast
RETURN:  ideal = maxideal(1) but the last variable exchanged by
         a sum of it with a linear random combination of the other
         variables
NOTE:
EXAMPLE: example randomLast; shows an example
"
{

  ideal i=maxideal(1);
  int k=size(i);
  i[k]=0;
  i=randomid(i,size(i),b);
  ideal ires=maxideal(1);
  ires[k]=i[1]+var(k);
  return(ires);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 0,(x,y,z),lp;
   ideal i = randomLast(10);
   i;
}


////////////////////////////////////////////////////////////////////////////////


proc primaryTest (ideal i, poly p)
{
   int m=1;
   int n=nvars(basering);
   int e;
   poly t;
   ideal h;

   ideal prm=p;
   attrib(prm,"isSB",1);

   while (n>1)
   {
      n=n-1;
      m=m+1;

      //search for i[m] which has a power of var(n) as leading term
      if (n==1)
      {
         m=size(i);
      }
      else
      {
         while (lead(i[m])/var(n-1)==0)
        {
            m=m+1;
         }
         m=m-1;
      }
      //check whether i[m] =(c*var(n)+h)^e modulo prm for some 
      //h in K[var(n+1),...,var(nvars(basering))], c in K
      //if not (0) is returned, else var(n)+h is added to prm

         e=deg(lead(i[m]));
        // t=hilfe1(i,prm,m,n);
         t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))/var(n)^(e-1);

         i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];

         if (reduce(i[m]-t^e,prm,1) !=0)
         {
           return(ideal(0));
         }
         h=interred(t);
         t=h[1];

      prm = prm,t;
      attrib(prm,"isSB",1);
   }
   return(prm);
}

proc hilfe(ideal i,ideal prm,int m)
{
   poly t;
   int e;

   if(size(i[m])==1)
   {
      t=var(n);
   }
   else
   {
      e=deg(lead(i[m]));

      if(deg(poly(e))>=0)
      {
           t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))/var(n)^(e-1);
           i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
      }
      else
      {
         i[m]=i[m]/leadcoef(i[m]);
         t=reduce(coef(i[m],var(n))[2,e+1],prm);
         t=var(n)+factorize(t,1)[1];
      }
   }
   return(t);
}
proc hilfe1(ideal i,ideal prm,int m,int n)
{
   poly t;
   int e;
   if(size(i[m])==1)
   {
      t=var(n);
   }
   else
   {
      e=deg(lead(i[m]));
      i[m]=i[m]/leadcoef(i[m]);
      t=reduce(coeffs(i[m],var(n))[1,1],prm);
      if(size(t)==0){return(var(n));}
      t=var(n)+factorize(t,1)[1];
   }
   return(t);
}

///////////////////////////////////////////////////////////////////////////////
proc splitPrimary(list l,ideal ser,int @wr,list sact)
{
   int i,j,k,s,r,w;
   list keepresult,act,keepprime;
   poly @f;
   int sl=size(l);

   for(i=1;i<=sl/2;i++)
   {
      if(sact[2][i]>1)
      {
         keepprime[i]=l[2*i-1]+ideal(sact[1][i]);
      }
      else
      {
         keepprime[i]=l[2*i-1];
      }
  }
   i=0;
   while(i<size(l)/2)
   {
      i++;
      if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==0))
      {
         l[2*i-1]=ideal(1);
         l[2*i]=ideal(1);
         continue;
      }


      if(size(l[2*i])==0)
      {
         if(homog(l[2*i-1])==1)
         {
            l[2*i]=maxideal(1);
            continue;
         }
         j=0;
         if(i<=sl/2)
         {
            j=1;
         }
         while(j<size(l[2*i-1]))
         {
            j++;
            act=factor(l[2*i-1][j]);
            r=size(act[1]);
            attrib(l[2*i-1],"isSB",1);
            if((r==1)&&(vdim(l[2*i-1])==deg(l[2*i-1][j])))
            {
              l[2*i]=std(l[2*i-1],act[1][1]);
              break;
            }
            if((r==1)&&(act[2][1]>1))
            {
               keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
               if(homog(keepprime[i])==1)
               {
                  l[2*i]=maxideal(1);
                  break;
               }
            }
            if(gcdTest(act[1])==1)
            {
               for(k=2;k<=r;k++)
               {
                  keepprime[size(l)/2+k-1]=interred(keepprime[i]+ideal(act[1][k]));
               }
               keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
               for(k=1;k<=r;k++)
               {
                  if(@wr==0)
                  {
                     keepresult[k]=std(l[2*i-1],act[1][k]^act[2][k]);
                  }
                  else
                  {
                     keepresult[k]=std(l[2*i-1],act[1][k]);
                  }
               }
               l[2*i-1]=keepresult[1];
               if(vdim(keepresult[1])==deg(act[1][1]))
               {
                  l[2*i]=keepresult[1];
               }
               if((homog(keepresult[1])==1)||(homog(keepprime[i])==1))
               {
                  l[2*i]=maxideal(1);
               }
               s=size(l)-2;
               for(k=2;k<=r;k++)
               {
                  l[s+2*k-1]=keepresult[k];
                  keepprime[s/2+k]=interred(keepresult[k]+ideal(act[1][k]));
                  if(vdim(keepresult[k])==deg(act[1][k]))
                  {
                     l[s+2*k]=keepresult[k];
                  }
                  else
                  {
                     l[s+2*k]=ideal(0);
                  }
                  if((homog(keepresult[k])==1)||(homog(keepprime[s/2+k])==1))
                  {
                    l[s+2*k]=maxideal(1);
                  }
               }
               i--;
               break;
            }
            if(r>=2)
            {
               s=size(l);
               @f=act[1][1];
               act=sat1(l[2*i-1],act[1][1]);
               if(deg(act[1][1])>0)
               {
                  l[s+1]=std(l[2*i-1],act[2]);
                  if(homog(l[s+1])==1)
                  {
                     l[s+2]=maxideal(1);
                  }
                  else
                  {
                     l[s+2]=ideal(0);
                  }
                  keepprime[s/2+1]=interred(keepprime[i]+ideal(@f));
                  if(homog(keepprime[s/2+1])==1)
                  {
                     l[s+2]=maxideal(1);
                  }
                  keepprime[i]=act[1];
                  l[2*i-1]=act[1];
                  attrib(l[2*i-1],"isSB",1);
                  if(homog(l[2*i-1])==1)
                  {
                     l[2*i]=maxideal(1);
                  }

                  i--;
                  break;
               }
            }
         }
      }
   }
   if(sl==size(l))
   {
      return(l);
   }
   for(i=1;i<=size(l)/2;i++)
   {
     attrib(l[2*i-1],"isSB",1);
     
     if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==0)&&(deg(l[2*i-1][1])>0))
     {
       "Achtung in split";
       
         l[2*i-1]=ideal(1);
         l[2*i]=ideal(1);
     }
     if((size(l[2*i])==0)&&(specialIdealsEqual(keepprime[i],l[2*i-1])!=1))
      {
         keepprime[i]=std(keepprime[i]);
         if(homog(keepprime[i])==1)
         {
             l[2*i]=maxideal(1);
         }
         else
         {
            act=zero_decomp(keepprime[i],ideal(0),@wr,1);
            if(size(act)==2)
            {
               l[2*i]=act[2];
            }
         }
      }
   }
   return(l);
}
example
{ "EXAMPLE:"; echo=2;
   LIB "primdec.lib";
   ring  r = 32003,(x,y,z),lp;
   ideal i1=x*(x+1),yz,(z+1)*(z-1);
   ideal i2=xy,yz,(x-2)*(x+3);
   list l=i1,ideal(0),i2,ideal(0),i2,ideal(1);
   list l1=splitPrimary(l,ideal(0),0);
   l1;
}

////////////////////////////////////////////////////////////////////////////////

proc zero_decomp (ideal j,ideal ser,int @wr,list #)
"USAGE:   zero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1
         (@wr=0 for primary decomposition, @wr=1 for computaion of associated
         primes)
RETURN:  list = list of primary ideals and their radicals (at even positions
         in the list) if the input is zero-dimensional and a standardbases
         with respect to lex-ordering
         If ser!=(0) and ser is contained in j or if j is not zero-dimen-
         sional then ideal(1),ideal(1) is returned
NOTE:    Algorithm of Gianni, Traeger, Zacharias
EXAMPLE: example zero_decomp; shows an example
"
{
  def   @P = basering;
  int nva = nvars(basering);
  int @k,@s,@n,@k1,zz;
  list primary,lres,lres1,act,@lh,@wh;
  map phi,psi,phi1,psi1;
  ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1;
  intvec @vh,@hilb;
  string @ri;
  poly @f;

  if (dim(j)>0)
  {
     primary[1]=ideal(1);
     primary[2]=ideal(1);
     return(primary);
  }
  
  j=interred(j);  
  attrib(j,"isSB",1);
  if(vdim(j)==deg(j[1]))
  {    
     act=factor(j[1]);
     for(@k=1;@k<=size(act[1]);@k++)
     {
       @qh=j;
       if(@wr==0)
       {
          @qh[1]=act[1][@k]^act[2][@k];
       }
       else
       {
          @qh[1]=act[1][@k];
       }
       primary[2*@k-1]=interred(@qh);
       @qh=j;
       @qh[1]=act[1][@k];
       primary[2*@k]=interred(@qh);
       attrib( primary[2*@k-1],"isSB",1);
       
       if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],1))==0))
       {
          primary[2*@k-1]=ideal(1);
          primary[2*@k]=ideal(1);          
       }
     }
     return(primary);
  }

  if(homog(j)==1)
  {
     primary[1]=j;
     if((size(ser)>0)&&(size(reduce(ser,j,1))==0))
     {
          primary[1]=ideal(1);
          primary[2]=ideal(1);
          return(primary);
     }
     if(dim(j)==-1)
     {
        primary[1]=ideal(1);
        primary[2]=ideal(1);
     }
     else
     {
        primary[2]=maxideal(1);
     }
     return(primary);
  }

//the first element in the standardbase is factorized
  if(deg(j[1])>0)
  {
    act=factor(j[1]);
    testFactor(act,j[1]);
  }
  else
  {
     primary[1]=ideal(1);
     primary[2]=ideal(1);
     return(primary);
  }

//withe the factors new ideals (hopefully the primary decomposition)
//are created

  if(size(act[1])>1)
  {
     if(size(#)>1)
     {
        primary[1]=ideal(1);
        primary[2]=ideal(1);
        primary[3]=ideal(1);
        primary[4]=ideal(1);
        return(primary);
     }
     for(@k=1;@k<=size(act[1]);@k++)
     {
        if(@wr==0)
        {
           primary[2*@k-1]=std(j,act[1][@k]^act[2][@k]);
        }
        else
        {
           primary[2*@k-1]=std(j,act[1][@k]);
        }
        if((act[2][@k]==1)&&(vdim(primary[2*@k-1])==deg(act[1][@k])))
        {
           primary[2*@k]   = primary[2*@k-1];
        }
        else
        {
           primary[2*@k]   = primaryTest(primary[2*@k-1],act[1][@k]);
        }
     }
  }
  else
  {  
     primary[1]=j;
     if((size(#)>0)&&(act[2][1]>1))
     {
        act[2]=1;
        primary[1]=std(primary[1],act[1][1]);
     }

     if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1])))
     {
        primary[2]=primary[1];
     }
     else
     {
        primary[2]=primaryTest(primary[1],act[1][1]);
     }
  }
  if(size(#)==0)
  {
    
     primary=splitPrimary(primary,ser,@wr,act);
     
  }

//test whether all ideals in the decomposition are primary and
//in general position
//if not after a random coordinate transformation of the last
//variable the corresponding ideal is decomposed again.

  @k=0;
  while(@k<(size(primary)/2))
  {
    @k++;
    if (size(primary[2*@k])==0)
    {
       for(zz=1;zz<size(primary[2*@k-1])-1;zz++)
       {
          if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][zz]))
          {
             primary[2*@k]=primary[2*@k-1];
          }
       }
    }
  }
    
  @k=0;
  ideal keep;
  while(@k<(size(primary)/2))
  {
    @k++;
    if (size(primary[2*@k])==0)
    {

       jmap=randomLast(100);
       jmap1=maxideal(1);
       jmap2=maxideal(1);
       @qht=primary[2*@k-1];

       for(@n=2;@n<=size(primary[2*@k-1]);@n++)
       {
          if(deg(lead(primary[2*@k-1][@n]))==1)
          {
             for(zz=1;zz<=nva;zz++)
             {
                if(lead(primary[2*@k-1][@n])/var(zz)!=0)
                {
                   jmap1[zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]
                   +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]);
                    jmap2[zz]=primary[2*@k-1][@n];
                    @qht[@n]=var(zz);
                     
                }
             }
             jmap[nva]=subst(jmap[nva],lead(primary[2*@k-1][@n]),0);
          }
       }
       if(size(subst(jmap[nva],var(1),0)-var(nva))!=0)
       {
          jmap[nva]=subst(jmap[nva],var(1),0);
       }
       phi1=@P,jmap1;
       phi=@P,jmap;
 
       for(@n=1;@n<=nva;@n++)
       {
          jmap[@n]=-(jmap[@n]-2*var(@n));
       }
       psi=@P,jmap;
       psi1=@P,jmap2;

       @qh=phi(@qht);
       
       if(npars(@P)>0)
       {
          @ri= "ring @Phelp ="
                  +string(char(@P))+",
                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
       }
       else
       {
          @ri= "ring @Phelp ="
                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
       }
       execute(@ri);
       ideal @qh=homog(imap(@P,@qht),@t);

       ideal @qh1=std(@qh);
       @hilb=hilb(@qh1,1);
       @ri= "ring @Phelp1 ="
                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
       execute(@ri);
       ideal @qh=homog(imap(@P,@qh),@t);
       kill @Phelp;
       @qh=std(@qh,@hilb);
       @qh=subst(@qh,@t,1);
       setring @P;
       @qh=imap(@Phelp1,@qh);
       kill @Phelp1;
       @qh=clearSB(@qh);
       attrib(@qh,"isSB",1);       
       ser1=phi1(ser);
       @lh=zero_decomp (@qh,phi(ser1),@wr);
//       @lh=zero_decomp (@qh,psi(ser),@wr);
       
       
       kill lres;
       list lres;
       if(size(@lh)==2)
       {
          helpprim=@lh[2];
          lres[1]=primary[2*@k-1];
          ser1=psi(helpprim);
          lres[2]=psi1(ser1);
          if(size(reduce(lres[2],lres[1],1))==0)
          {
             primary[2*@k]=primary[2*@k-1];
             continue;
          }
       }
       else
       {
          act=factor(@qh[1]);
          if(2*size(act[1])==size(@lh))
          {
             for(@n=1;@n<=size(act[1]);@n++)
             {
                @f=act[1][@n]^act[2][@n];
                ser1=psi(@f);
                lres[2*@n-1]=interred(primary[2*@k-1]+psi1(ser1));
                helpprim=@lh[2*@n];
                ser1=psi(helpprim);
                lres[2*@n]=psi1(ser1);
             }
          }
          else
          {
             lres1=psi(@lh);
             lres=psi1(lres1);
          }
       }
       if(npars(@P)>0)
       {
          @ri= "ring @Phelp ="
                  +string(char(@P))+",
                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
       }
       else
       {
          @ri= "ring @Phelp ="
                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
       }
       execute(@ri);
       list @lvec;
       list @lr=imap(@P,lres);
       ideal @lr1;

       if(size(@lr)==2)
       {
          @lr[2]=homog(@lr[2],@t);
          @lr1=std(@lr[2]);
          @lvec[2]=hilb(@lr1,1);
       }
       else
       {
          for(@n=1;@n<=size(@lr)/2;@n++)
          {
             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
             {
                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
                @lr1=std(@lr[2*@n-1]);
                @lvec[2*@n-1]=hilb(@lr1,1);
                @lvec[2*@n]=@lvec[2*@n-1];
             }
             else
             {
                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
                @lr1=std(@lr[2*@n-1]);
                @lvec[2*@n-1]=hilb(@lr1,1);
                @lr[2*@n]=homog(@lr[2*@n],@t);
                @lr1=std(@lr[2*@n]);
                @lvec[2*@n]=hilb(@lr1,1);

             }
         }
       }
       @ri= "ring @Phelp1 ="
                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
       execute(@ri);
       list @lr=imap(@Phelp,@lr);

       kill @Phelp;
       if(size(@lr)==2)
       {
          @lr[2]=std(@lr[2],@lvec[2]);
          @lr[2]=subst(@lr[2],@t,1);

       }
       else
       {
          for(@n=1;@n<=size(@lr)/2;@n++)
          {
             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
             {
                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
                @lr[2*@n]=@lr[2*@n-1];
                attrib(@lr[2*@n],"isSB",1);
             }
             else
             {
                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
                @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]);
                @lr[2*@n]=subst(@lr[2*@n],@t,1);
             }
          }
       }
       kill @lvec;
       setring @P;
       lres=imap(@Phelp1,@lr);
       kill @Phelp1;
       for(@n=1;@n<=size(lres);@n++)
       {
          lres[@n]=clearSB(lres[@n]);
          attrib(lres[@n],"isSB",1);
       }

       primary[2*@k-1]=lres[1];
       primary[2*@k]=lres[2];
       @s=size(primary)/2;
       for(@n=1;@n<=size(lres)/2-1;@n++)
       {
         primary[2*@s+2*@n-1]=lres[2*@n+1];
         primary[2*@s+2*@n]=lres[2*@n+2];
       }
       @k--;
     }
  }
  return(primary);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 0,(x,y,z),lp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   i=std(i);
   list  pr= zero_decomp(i,ideal(0),0);
   pr;
}

////////////////////////////////////////////////////////////////////////////////

proc ggt (ideal i)
"USAGE:   ggt(i); i list of polynomials
RETURN:  poly = ggt(i[1],...,i[size(i)])
NOTE:
EXAMPLE: example ggt; shows an example
"
{
  int k;
  poly p=i[1];
  if(deg(p)==0)
  {
    return(1);
  }


  for (k=2;k<=size(i);k++)
  {
     if(deg(i[k])==0)
     {
        return(1)
     }
     p=GCD(p,i[k]);
     if(deg(p)==0)
     {
        return(1);
     }
  }
  return(p);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 0,(x,y,z),lp;
   poly  p = (x+y)*(y+z);
   poly  q = (z4+2)*(y+z);
   ideal l=p,q;
   poly  pr= ggt(l);
   pr;
}
///////////////////////////////////////////////////////////////////////////////
proc gcdTest(ideal act)
{
  int i,j;
  if(size(act)<=1)
  {
     return(0);
  }
  for (i=1;i<=size(act)-1;i++)
  {
     for(j=i+1;j<=size(act);j++)
     {
        if(deg(std(ideal(act[i],act[j]))[1])>0)
        {
           return(0);
        }
     }
  }
  return(1);
}

///////////////////////////////////////////////////////////////////////////////
proc coeffLcm(ideal h)
{
   string @pa=parstr(basering);
   if(size(@pa)==0)
   {
      return(lcmP(h));
   }
   def bsr= basering;
   string @id=string(h);
   execute "ring @r=0,("+@pa+","+varstr(bsr)+"),(C,dp);";
   execute "ideal @i="+@id+";";
   poly @p=lcmP(@i);
   string @ps=string(@p);
   setring bsr;
   execute "poly @p="+@ps+";";
   return(@p);
}
example
{
   "EXAMPLE:"; echo = 2;
   ring  r =( 0,a,b),(x,y,z),lp;
   poly  p = (a+b)*(y-z);
   poly  q = (a+b)*(y+z);
   ideal l=p,q;
   poly  pr= coeffLcm(l);
   pr;
}

///////////////////////////////////////////////////////////////////////////////

proc lcmP(ideal i)
"USAGE:   lcm(i); i list of polynomials
RETURN:  poly = lcm(i[1],...,i[size(i)])
NOTE:
EXAMPLE: example lcm; shows an example
"
{
  int k,j;
   poly p,q;
  i=simplify(i,10);
  for(j=1;j<=size(i);j++)
  {
    if(deg(i[j])>0)
    {
      p=i[j];
      break;
    }
  }
  if(deg(p)==-1)
  {
    return(1);
  }
  for (k=j+1;k<=size(i);k++)
  {
     if(deg(i[k])!=0)
     {
        q=GCD(p,i[k]);
        if(deg(q)==0)
        {
           p=p*i[k];
        }
        else
        {
           p=p/q;
           p=p*i[k];
        }
     }
   }
  return(p);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 0,(x,y,z),lp;
   poly  p = (x+y)*(y+z);
   poly  q = (z4+2)*(y+z);
   ideal l=p,q;
   poly  pr= lcmP(l);
   pr;
   l=1,-1,p,1,-1,q,1;
   pr=lcmP(l);
   pr;
}

///////////////////////////////////////////////////////////////////////////////
proc clearSB (ideal i,list #)
"USAGE:   clearSB(i); i ideal which is SB ordered by monomial ordering
RETURN:  ideal = minimal SB
NOTE:
EXAMPLE: example clearSB; shows an example
"
{
  int k,j;
  poly m;
  int c=size(i);

  if(size(#)==0)
  {
    for(j=1;j<c;j++)
    {
      if(deg(i[j])==0)
      {
        i=ideal(1);
        return(i);
      }
      if(deg(i[j])>0)
      {
        m=lead(i[j]);
        for(k=j+1;k<=c;k++)
        {
           if(size(lead(i[k])/m)>0)
           {
             i[k]=0;
           }
        }
      }
    }
  }
  else
  {
    j=0;
    while(j<c-1)
    {
      j++;
      if(deg(i[j])==0)
      {
        i=ideal(1);
        return(i);
      }
      if(deg(i[j])>0)
      {
        m=lead(i[j]);
        for(k=j+1;k<=c;k++)
        {
           if(size(lead(i[k])/m)>0)
           {
             if((leadexp(m)!=leadexp(i[k]))||(#[j]<=#[k]))
             {
                i[k]=0;
             }
             else
             {
                i[j]=0;
                break;
             }
           }
        }
      }
    }
  }
  return(simplify(i,2));
}
example
{ "EXAMPLE:"; echo = 2;
   LIB "primdec.lib";
   ring  r = (0,a,b),(x,y,z),dp;
   ideal i=ax2+y,a2x+y,bx;
   list l=1,2,1;
   ideal j=clearSB(i,l);
   j;
}

///////////////////////////////////////////////////////////////////////////////

proc independSet (ideal j)
"USAGE:   independentSet(i); i ideal
RETURN:  list = new varstring with the independent set at the end,
                ordstring with the corresponding block ordering,
                the integer where the independent set starts in the varstring
NOTE:
EXAMPLE: example independentSet; shows an example
"
{
   int n,k,di;
   list resu,hilf;
   string var1,var2;
   list v=system("indsetall",j,1);

   for(n=1;n<=size(v);n++)
   {
      di=0;
      var1="";
      var2="";
      for(k=1;k<=size(v[n]);k++)
      {
         if(v[n][k]!=0)
         {
            di++;
            var2=var2+"var("+string(k)+"),";
         }
         else
         {
            var1=var1+"var("+string(k)+"),";
         }
      }
      if(di>0)
      {
         var1=var1+var2;
         var1=var1[1..size(var1)-1];
         hilf[1]=var1;
         hilf[2]="lp";
         //"lp("+string(nvars(basering)-di)+"),dp("+string(di)+")";
         hilf[3]=di;
         resu[n]=hilf;
      }
      else
      {
         resu[n]=varstr(basering),ordstr(basering),0;
      }
   }
   return(resu);
}
example
{ "EXAMPLE:"; echo = 2;
   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
   i=std(i);
   list  l=independSet(i);
   l;
   i=i,g;
   l=independSet(i);
   l;

   ring s=0,(x,y,z),lp;
   ideal i=z,yx;
   list l=independSet(i);
   l;


}
///////////////////////////////////////////////////////////////////////////////

proc maxIndependSet (ideal j)
"USAGE:   maxIndependentSet(i); i ideal
RETURN:  list = new varstring with the maximal independent set at the end,
                ordstring with the corresponding block ordering,
                the integer where the independent set starts in the varstring
NOTE:
EXAMPLE: example maxIndependentSet; shows an example
"
{
   int n,k,di;
   list resu,hilf;
   string var1,var2;
   list v=system("indsetall",j,0);

   for(n=1;n<=size(v);n++)
   {
      di=0;
      var1="";
      var2="";
      for(k=1;k<=size(v[n]);k++)
      {
         if(v[n][k]!=0)
         {
            di++;
            var2=var2+"var("+string(k)+"),";
         }
         else
         {
            var1=var1+"var("+string(k)+"),";
         }
      }
      if(di>0)
      {
         var1=var1+var2;
         var1=var1[1..size(var1)-1];
         hilf[1]=var1;
         hilf[2]="lp";
         hilf[3]=di;
         resu[n]=hilf;
      }
      else
      {
         resu[n]=varstr(basering),ordstr(basering),0;
      }
   }
   return(resu);
}
example
{ "EXAMPLE:"; echo = 2;
   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
   i=std(i);
   list  l=maxIndependSet(i);
   l;
   i=i,g;
   l=maxIndependSet(i);
   l;

   ring s=0,(x,y,z),lp;
   ideal i=z,yx;
   list l=maxIndependSet(i);
   l;


}

///////////////////////////////////////////////////////////////////////////////

proc prepareQuotientring (int nnp)
"USAGE:   prepareQuotientring(nnp); nnp int
RETURN:  string = to define Kvar(nnp+1),...,var(nvars)[..rest ]
NOTE:
EXAMPLE: example independentSet; shows an example
"
{
  ideal @ih,@jh;
  int npar=npars(basering);
  int @n;

  string quotring= "ring quring = ("+charstr(basering);
  for(@n=nnp+1;@n<=nvars(basering);@n++)
  {
     quotring=quotring+",var("+string(@n)+")";
     @ih=@ih+var(@n);
  }

  quotring=quotring+"),(var(1)";
  @jh=@jh+var(1);
  for(@n=2;@n<=nnp;@n++)
  {
    quotring=quotring+",var("+string(@n)+")";
    @jh=@jh+var(@n);
  }
  quotring=quotring+"),(C,lp);";

  return(quotring);

}
example
{ "EXAMPLE:"; echo = 2;
   ring s1=(0,x),(a,b,c,d,e,f,g),lp;
   def @Q=basering;
   list l= prepareQuotientring(3);
   l;
   execute l[1];
   execute l[2];
   basering;
   phi;
   setring @Q;

}

///////////////////////////////////////////////////////////////////////

proc projdim(list l)
{
   int i=size(l)+1;

   while(i>2)
   {
      i--;
      if((size(l[i])>0)&&(deg(l[i][1])>0))
      {
         return(i);
      }
   }
}

///////////////////////////////////////////////////////////////////////////////
proc cleanPrimary(list l)
{
   int i,j;
   list lh;
   for(i=1;i<=size(l)/2;i++)
   {
      if(deg(l[2*i-1][1])>0)
      {
         j++;
         lh[j]=l[2*i-1];
         j++;
         lh[j]=l[2*i];
      }
   }
   return(lh);
}
///////////////////////////////////////////////////////////////////////////////

proc minAssPrimes(ideal i, list #)
"USAGE:   minAssPrimes(i); i ideal
         minAssPrimes(i,1); i ideal  (to use also the factorizing Groebner)
RETURN:  list = the minimal associated prime ideals of i
NOTE:
EXAMPLE: example minAssPrimes; shows an example
"
{
   #[1]=1;
   def @P=basering;
   list qr=simplifyIdeal(i);
   map phi=@P,qr[2];
   i=qr[1];
   
   execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
             +ordstr(basering)+");";


   ideal i=fetch(@P,i);
   if(size(#)==0)
   {
      int @wr;
      list tluser,@res;
      list primary=decomp(i,2);

      @res[1]=primary;

      tluser=union(@res);
      setring @P;
      list @res=imap(gnir,tluser);
      return(phi(@res));
   }
   list @res,empty;
   ideal ser;
   option(redSB);
   list @pr=facstd(i);
   if(size(@pr)==1)
   {
      attrib(@pr[1],"isSB",1);
      if((dim(@pr[1])==0)&&(homog(@pr[1])==1))
      {
         setring @P;
         list @res=maxideal(1);
         return(phi(@res));
      }
      if(dim(@pr[1])>1)
      {
         setring @P;
        // kill gnir;
         execute "ring gnir1 = ("+charstr(basering)+"),
                              ("+varstr(basering)+"),(C,lp);";
         ideal i=fetch(@P,i);
         list @pr=facstd(i);
        // ideal ser;
         setring gnir;
         @pr=fetch(gnir1,@pr);
         kill gnir1;
      }
   }
    option(noredSB);
   int j,k,odim,ndim,count;
   attrib(@pr[1],"isSB",1);
   if(#[1]==77)
   {
     odim=dim(@pr[1]);
     count=1;
     intvec pos;
     pos[size(@pr)]=0;
     for(j=2;j<=size(@pr);j++)
     {
        attrib(@pr[j],"isSB",1);
        ndim=dim(@pr[j]);
        if(ndim>odim)
        {
           for(k=count;k<=j-1;k++)
           {
              pos[k]=1;
           }
           count=j;
           odim=ndim;
        }
        if(ndim<odim)
        {
           pos[j]=1;
        }
     }
     for(j=1;j<=size(@pr);j++)
     {
        if(pos[j]!=1)
        {
            @res[j]=decomp(@pr[j],2);
        }
        else
        {
           @res[j]=empty;
        }
     }
   }
   else
   {
     ser=ideal(1);
     for(j=1;j<=size(@pr);j++)
     {
//@pr[j];
//pause;
        @res[j]=decomp(@pr[j],2);
//       @res[j]=decomp(@pr[j],2,@pr[j],ser);
//       for(k=1;k<=size(@res[j]);k++)
//       {
//          ser=intersect(ser,@res[j][k]);
//       }
     }
   }

   @res=union(@res);
   setring @P;
   list @res=imap(gnir,@res);
   return(phi(@res));
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 32003,(x,y,z),lp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   LIB "primaryDecomposition.lib";
   list pr= minAssPrimes(i);
   pr;
   pr= minAssPrimes(i,1);
}

///////////////////////////////////////////////////////////////////////////////

proc union(list li)
{
  int i,j,k;

  def P=basering;

  execute "ring ir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);";
  list l=fetch(P,li);
  list @erg;

  for(k=1;k<=size(l);k++)
  {
     for(j=1;j<=size(l[k])/2;j++)
     {
        if(deg(l[k][2*j][1])!=0)
        {
           i++;
           @erg[i]=l[k][2*j];
        }
     }
  }

  list @wos;
  i=0;
  ideal i1,i2;
  while(i<size(@erg)-1)
  {
     i++;
     k=i+1;
     i1=lead(@erg[i]);
      attrib(i1,"isSB",1);
      attrib(@erg[i],"isSB",1);

     while(k<=size(@erg))
     {
        if(deg(@erg[i][1])==0)
        {
           break;
        }
        i2=lead(@erg[k]);
        attrib(@erg[k],"isSB",1);
        attrib(i2,"isSB",1);

        if(size(reduce(i1,i2,1))==0)
        {
           if(size(reduce(@erg[i],@erg[k],1))==0)
           {
              @erg[k]=ideal(1);
              i2=ideal(1);
           }
        }
        if(size(reduce(i2,i1,1))==0)
        {
           if(size(reduce(@erg[k],@erg[i],1))==0)
           {
              break;
           }
        }
        k++;
        if(k>size(@erg))
        {
           @wos[size(@wos)+1]=@erg[i];
        }
     }
  }
  if(deg(@erg[size(@erg)][1])!=0)
  {
     @wos[size(@wos)+1]=@erg[size(@erg)];
  }
  setring P;
  list @ser=fetch(ir,@wos);
  return(@ser);
}
///////////////////////////////////////////////////////////////////////////////
proc radicalOld(ideal i)
{
   list pr=minAssPrimes(i,1);
   int j;
   ideal k=pr[1];
   for(j=2;j<=size(pr);j++)
   {
      k=intersect(k,pr[j]);
   }
   return(k);
}
///////////////////////////////////////////////////////////////////////////////
proc equiRadical(ideal i)
{
   return(radical(i,1));
}

///////////////////////////////////////////////////////////////////////////////
proc decomp(ideal i,list #)
"USAGE:   decomp(i); i ideal  (for primary decomposition)   (resp.
         decomp(i,1);        (for the minimal associated primes) )
RETURN:  list = list of primary ideals and their associated primes
         (at even positions in the list)
         (resp. a list of the minimal associated primes)
NOTE:    Algorithm of Gianni, Traeger, Zacharias
EXAMPLE: example decomp; shows an example
"
{
  def  @P = basering;
  list primary,indep,ltras;
  intvec @vh,isat;
  int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi;
  ideal peek=i;
  ideal ser,tras;

  if(size(#)>0)
  {
     if((#[1]==1)||(#[1]==2))
     {
        @wr=#[1];
        if(size(#)>1)
        {
          seri=1;
          peek=#[2];
          ser=#[3];
        }
      }
      else
      {
        seri=1;
        peek=#[1];
        ser=#[2];
      }
  }
  
  homo=homog(i);
  
  if(homo==1)
  {
    if(attrib(i,"isSB")!=1)
    {
      ltras=mstd(i);
      attrib(ltras[1],"isSB",1);
    }
    else
    {
      ltras=i,i;
    }
    tras=ltras[1];
    if(dim(tras)==0)
    {
        primary[1]=ltras[2];
        primary[2]=maxideal(1);
        if(@wr>0)
        {
           list l;
           l[1]=maxideal(1);
           l[2]=maxideal(1);
           return(l);
        }
        return(primary);
     }
     intvec @hilb=hilb(tras,1);
     intvec keephilb=@hilb;
  }

  //----------------------------------------------------------------
  //i is the zero-ideal
  //----------------------------------------------------------------

  if(size(i)==0)
  {
     primary=i,i;
     return(primary);
  }

  //----------------------------------------------------------------
  //pass to the lexicographical ordering and compute a standardbasis
  //----------------------------------------------------------------

  execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);";
  option(redSB);

  ideal ser=fetch(@P,ser);  

  if(homo==1)
  {
     if((ordstr(@P)[1]!="(C,lp)")&&(ordstr(@P)[3]!="(C,lp)"))
     {
       ideal @j=std(fetch(@P,i),@hilb);
     }
     else
     {
        ideal @j=fetch(@P,tras);
        attrib(@j,"isSB",1);
     }
  }
  else
  {
     ideal @j=std(fetch(@P,i));
  }
  option(noredSB);
  if(seri==1)
  {
    ideal peek=fetch(@P,peek);
    attrib(peek,"isSB",1);
  }
  else
  {
    ideal peek=@j;
  }
  if(size(ser)==0)
  {
    ideal fried;
    @n=size(@j);
    for(@k=1;@k<=@n;@k++)
    {
      if(deg(lead(@j[@k]))==1)
      {
        fried[size(fried)+1]=@j[@k];
        @j[@k]=0;
      }
    }
    if(size(fried)>0)
    {
       @j=simplify(@j,2);
       attrib(@j,"isSB",1);
       list pr=decomp(@j);
       for(@k=1;@k<=size(pr);@k++)
       {
         @j=pr[@k]+fried;
         pr[@k]=@j;
       }
       setring @P;
       return(fetch(gnir,pr));
    }
  }
 
  //----------------------------------------------------------------
  //j is the ring
  //----------------------------------------------------------------

  if (dim(@j)==-1)
  {
    setring @P;
    return(ideal(0));
  }

  //----------------------------------------------------------------
  //  the case of one variable
  //----------------------------------------------------------------

  if(nvars(basering)==1)
  {

     list fac=factor(@j[1]);
     list gprimary;
     for(@k=1;@k<=size(fac[1]);@k++)
     {
        if(@wr==0)
        {
           gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]);
           gprimary[2*@k]=ideal(fac[1][@k]);
        }
        else
        {
           gprimary[2*@k-1]=ideal(fac[1][@k]);
           gprimary[2*@k]=ideal(fac[1][@k]);
        }
     }
     setring @P;
     primary=fetch(gnir,gprimary);

     return(primary);
  }
  
 //------------------------------------------------------------------
 //the zero-dimensional case
 //------------------------------------------------------------------

  if (dim(@j)==0)
  {
    option(redSB);
    list gprimary= zero_decomp(@j,ser,@wr);
    option(noredSB);
    setring @P;
    primary=fetch(gnir,gprimary);
    if(size(ser)>0)
    {
      primary=cleanPrimary(primary);
    }
    return(primary);
  }

  poly @gs,@gh,@p;
  string @va,quotring;
  list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep;
  ideal @h;
  int jdim=dim(@j);
  list fett;
  int lauf,di,newtest;
  //------------------------------------------------------------------
  //search for a maximal independent set indep,i.e.
  //look for subring such that the intersection with the ideal is zero
  //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
  //indep[1] is the new varstring and indep[2] the string for the block-ordering
  //------------------------------------------------------------------

  if(@wr!=1)
  {
     allindep=independSet(@j);
     for(@m=1;@m<=size(allindep);@m++)
     {
        if(allindep[@m][3]==jdim)
        {
           di++;
           indep[di]=allindep[@m];
        }
        else
        {
           lauf++;
           restindep[lauf]=allindep[@m];
        }
     }
   }
   else
   {
      indep=maxIndependSet(@j);
   }
  
  ideal jkeep=@j;

  if(ordstr(@P)[1]=="w")
  {
     execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),("+ordstr(@P)+");";
  }
  else
  {
     execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),(C,dp);";
  }

  if(homo==1)
  {
    if((ordstr(@P)[3]=="d")||(ordstr(@P)[1]=="d")||(ordstr(@P)[1]=="w")
       ||(ordstr(@P)[3]=="w"))
    {
      ideal jwork=imap(@P,tras);
      attrib(jwork,"isSB",1);
    }
    else
    {
      ideal jwork=std(imap(gnir,@j),@hilb);
    }
    
  }
  else
  {
    ideal jwork=std(imap(gnir,@j));
  }
  list hquprimary;
  poly @p,@q;
  ideal @h,fac,ser;
  di=dim(jwork);
  keepdi=di;
  
  setring gnir;
  for(@m=1;@m<=size(indep);@m++)
  {
     isat=0;
     @n2=0;
     option(redSB);
     if((indep[@m][1]==varstr(basering))&&(@m==1))
     //this is the good case, nothing to do, just to have the same notations
     //change the ring
     {
        execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
                              +ordstr(basering)+");";
        ideal @j=fetch(gnir,@j);
        attrib(@j,"isSB",1);
        ideal ser=fetch(gnir,ser);
        
     }
     else
     {
        @va=string(maxideal(1));
        execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
                              +indep[@m][2]+");";
        execute "map phi=gnir,"+@va+";";
        if(homo==1)
        {
           ideal @j=std(phi(@j),@hilb);
        }
        else
        {
          ideal @j=std(phi(@j));
        }
        ideal ser=phi(ser);
        
     }
     option(noredSB);     
     if((deg(@j[1])==0)||(dim(@j)<jdim))
     {
       setring gnir;
       break;
     }
     for (lauf=1;lauf<=size(@j);lauf++)
     {
         fett[lauf]=size(@j[lauf]);
     }
     //------------------------------------------------------------------------------------
     //we have now the following situation:
     //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
     //to this quotientring, j is their still a standardbasis, the
     //leading coefficients of the polynomials  there (polynomials in
     //K[var(nnp+1),..,var(nva)]) are collected in the list h,
     //we need their ggt, gh, because of the following:
     //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
     //intersected with K[var(1),...,var(nva)] is (j:gh^n)
     //on the other hand j=(j,gh^n) intersected with (j:gh^n)

     //------------------------------------------------------------------------------------

     //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
     //and the map phi:K[var(1),...,var(nva)] ----->K(var(nnpr+1),..,var(nva))[..the rest..]
     //-------------------------------------------------------------------------------------

     quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);

     //---------------------------------------------------------------------
     //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     //---------------------------------------------------------------------

     execute quotring;

    // @j considered in the quotientring
     ideal @j=imap(gnir1,@j);
     ideal ser=imap(gnir1,ser);

     kill gnir1;

     //j is a standardbasis in the quotientring but usually not minimal
     //here it becomes minimal

     @j=clearSB(@j,fett);
     attrib(@j,"isSB",1);

     //we need later ggt(h[1],...)=gh for saturation
     ideal @h;
     if(deg(@j[1])>0)
     {
        for(@n=1;@n<=size(@j);@n++)
        {
           @h[@n]=leadcoef(@j[@n]);
        }
        //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
        option(redSB);        
        list uprimary= zero_decomp(@j,ser,@wr);
        option(noredSB);
     }
     else
     {
       list uprimary;
       uprimary[1]=ideal(1);
       uprimary[2]=ideal(1);
     }

     //we need the intersection of the ideals in the list quprimary with the
     //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
     //but fi polynomials, then the intersection of q with the polynomialring
     //is the saturation of the ideal generated by f1,...,fr with respect to
     //h which is the lcm of the leading coefficients of the fi considered in the
     //quotientring: this is coded in saturn

     list saturn;
     ideal hpl;

     for(@n=1;@n<=size(uprimary);@n++)
     {
        hpl=0;
        for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
        {
           hpl=hpl,leadcoef(uprimary[@n][@n1]);
        }
        saturn[@n]=hpl;
     }
     //--------------------------------------------------------------------
     //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     //back to the polynomialring
     //---------------------------------------------------------------------
     setring gnir;

     collectprimary=imap(quring,uprimary);
     lsau=imap(quring,saturn);
     @h=imap(quring,@h);

     kill quring;


     @n2=size(quprimary);
     @n3=@n2;

     for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
     {
        if(deg(collectprimary[2*@n1][1])>0)
        {
           @n2++;
           quprimary[@n2]=collectprimary[2*@n1-1];
           lnew[@n2]=lsau[2*@n1-1];
           @n2++;
           lnew[@n2]=lsau[2*@n1];
           quprimary[@n2]=collectprimary[2*@n1];
        }
     }

     //here the intersection with the polynomialring
     //mentioned above is really computed
     for(@n=@n3/2+1;@n<=@n2/2;@n++)
     {       
        if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
        {
           quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
           quprimary[2*@n]=quprimary[2*@n-1];
        }
        else
        {
           if(@wr==0)
           {
              quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
           }
           quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
        }
     }
    
     if(size(@h)>0)
     {
           //---------------------------------------------------------------
           //we change to @Phelp to have the ordering dp for saturation
           //---------------------------------------------------------------
        setring @Phelp;
        @h=imap(gnir,@h);
        if(@wr!=1)
        {
          @q=minSat(jwork,@h)[2];   
        }
        else
        {
            fac=ideal(0);
            for(lauf=1;lauf<=ncols(@h);lauf++)
           {
              if(deg(@h[lauf])>0)
              {
                 fac=fac+factorize(@h[lauf],1);
              }
           }
           fac=simplify(fac,4);
           @q=1;
           for(lauf=1;lauf<=size(fac);lauf++)
           {
             @q=@q*fac[lauf];
           }
        }
        jwork=std(jwork,@q);
        keepdi=dim(jwork);        
        if(keepdi<di)
        {
           setring gnir;
           @j=imap(@Phelp,jwork);
           break;
        }
        if(homo==1)
        {
              @hilb=hilb(jwork,1);
        }

        setring gnir;
        @j=imap(@Phelp,jwork);
     }
  }
  if((size(quprimary)==0)&&(@wr>0))
  {
     @j=ideal(1);
     quprimary[1]=ideal(1);
     quprimary[2]=ideal(1);
  }
  if((size(quprimary)==0))
  {
    keepdi=di-1;
  }  
  //---------------------------------------------------------------
  //notice that j=sat(j,gh) intersected with (j,gh^n)
  //we finished with sat(j,gh) and have to start with (j,gh^n)
  //---------------------------------------------------------------
  if((deg(@j[1])!=0)&&(@wr!=1))
  {
     if(size(quprimary)>0)
     {
        setring @Phelp;
        ser=imap(gnir,ser);
        hquprimary=imap(gnir,quprimary);
        if(@wr==0)
        {
           ideal htest=hquprimary[1];
           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
           {
              htest=intersect(htest,hquprimary[2*@n1-1]);
           }
        }
        else
        {
           ideal htest=hquprimary[2];

           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
           {
              htest=intersect(htest,hquprimary[2*@n1]);
           }
        }

        if(size(ser)>0)
        { 
           ser=intersect(htest,ser); 
        }
        else
        {
          ser=htest;
        } 
        setring gnir;
        ser=imap(@Phelp,ser);
     }
     if(size(reduce(ser,peek,1))!=0)
     {       
        for(@m=1;@m<=size(restindep);@m++)
        {
         // if(restindep[@m][3]>=keepdi)
         // { 
           isat=0;
           @n2=0;
           option(redSB);
           
           if(restindep[@m][1]==varstr(basering))
           //this is the good case, nothing to do, just to have the same notations
           //change the ring
           {
              execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
                              +ordstr(basering)+");";
              ideal @j=fetch(gnir,jkeep);
              attrib(@j,"isSB",1);
           }
           else
           {
              @va=string(maxideal(1));
              execute "ring gnir1 = ("+charstr(basering)+"),("+restindep[@m][1]+"),("
                              +restindep[@m][2]+");";
              execute "map phi=gnir,"+@va+";";
              if(homo==1)
              {
                 ideal @j=std(phi(jkeep),keephilb);
              }
              else
              {
                ideal @j=std(phi(jkeep));
              }
              ideal ser=phi(ser);
           }
           option(noredSB);
           
           for (lauf=1;lauf<=size(@j);lauf++)
           {
              fett[lauf]=size(@j[lauf]);
           }
           //------------------------------------------------------------------------------------
           //we have now the following situation:
           //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
           //to this quotientring, j is their still a standardbasis, the
           //leading coefficients of the polynomials  there (polynomials in
           //K[var(nnp+1),..,var(nva)]) are collected in the list h,
           //we need their ggt, gh, because of the following:
           //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
           //intersected with K[var(1),...,var(nva)] is (j:gh^n)
           //on the other hand j=(j,gh^n) intersected with (j:gh^n)

           //------------------------------------------------------------------------------------

           //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
           //and the map phi:K[var(1),...,var(nva)] ----->K(var(nnpr+1),..,var(nva))[..the rest..]
           //-------------------------------------------------------------------------------------

           quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]);

           //---------------------------------------------------------------------
           //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
           //---------------------------------------------------------------------

           execute quotring;

           // @j considered in the quotientring
           ideal @j=imap(gnir1,@j);
           ideal ser=imap(gnir1,ser);

           kill gnir1;

           //j is a standardbasis in the quotientring but usually not minimal
           //here it becomes minimal
           @j=clearSB(@j,fett);
           attrib(@j,"isSB",1);

           //we need later ggt(h[1],...)=gh for saturation
           ideal @h;

           for(@n=1;@n<=size(@j);@n++)
           {
              @h[@n]=leadcoef(@j[@n]);
           }
           //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
           
           option(redSB);
           list uprimary= zero_decomp(@j,ser,@wr);
           option(noredSB);
           
            
           //we need the intersection of the ideals in the list quprimary with the
           //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
           //but fi polynomials, then the intersection of q with the polynomialring
           //is the saturation of the ideal generated by f1,...,fr with respect to
           //h which is the lcm of the leading coefficients of the fi considered in the
           //quotientring: this is coded in saturn

           list saturn;
           ideal hpl;

           for(@n=1;@n<=size(uprimary);@n++)
           {
              hpl=0;
              for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
              {
                 hpl=hpl,leadcoef(uprimary[@n][@n1]);
              }
               saturn[@n]=hpl;
           }
           //--------------------------------------------------------------------
           //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
           //back to the polynomialring
           //---------------------------------------------------------------------
           setring gnir;

           collectprimary=imap(quring,uprimary);
           lsau=imap(quring,saturn);
           @h=imap(quring,@h);

           kill quring;


           @n2=size(quprimary);
           @n3=@n2;
           
           for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
           {
              if(deg(collectprimary[2*@n1][1])>0)
              {
                 @n2++;
                 quprimary[@n2]=collectprimary[2*@n1-1];
                 lnew[@n2]=lsau[2*@n1-1];
                 @n2++;
                 lnew[@n2]=lsau[2*@n1];
                 quprimary[@n2]=collectprimary[2*@n1];
              }
           }
           
           
           //here the intersection with the polynomialring
           //mentioned above is really computed

           for(@n=@n3/2+1;@n<=@n2/2;@n++)
           {
              if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
              {
                 quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
                 quprimary[2*@n]=quprimary[2*@n-1];
              }
              else
              {
                 if(@wr==0)
                 {
                    quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
                 }
                 quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
              }
           }
           if(@n2>=@n3+2)
           {
              setring @Phelp;
              ser=imap(gnir,ser);
              hquprimary=imap(gnir,quprimary);
              for(@n=@n3/2+1;@n<=@n2/2;@n++)
              {
                if(@wr==0)
                {
                   ser=intersect(ser,hquprimary[2*@n-1]);
                }
                else
                {
                   ser=intersect(ser,hquprimary[2*@n]);
                }
              }
              setring gnir;
              ser=imap(@Phelp,ser);
           }
           
         // }
        }             
        if(size(reduce(ser,peek,1))!=0)
        {
           if(@wr>0)
           {
              htprimary=decomp(@j,@wr,peek,ser);
           }
           else
           {
              htprimary=decomp(@j,peek,ser);
           }
           // here we collect now both results primary(sat(j,gh))
           // and primary(j,gh^n)
           @n=size(quprimary);
           for (@k=1;@k<=size(htprimary);@k++)
           {
              quprimary[@n+@k]=htprimary[@k];
           }
        }
     }
     
   }
  //------------------------------------------------------------
  //back to the ring we started with
  //the final result: primary
  //------------------------------------------------------------
  
  setring @P;
  primary=imap(gnir,quprimary);
  return(primary);
}


example
{ "EXAMPLE:"; echo = 2;
   ring  r = 32003,(x,y,z),lp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   LIB "primdec.lib";
   list pr= decomp(i);
   pr;
   testPrimary( pr, i);
}

///////////////////////////////////////////////////////////////////////////////
proc primdecGTZ(ideal i)
{
   return(decomp(i));
}
///////////////////////////////////////////////////////////////////////////////
proc primdecSY(ideal i,int j)
{
   return(prim_dec(i,j));
}
///////////////////////////////////////////////////////////////////////////////

proc radical(ideal i,list #)
{
   def @P=basering;
   int j,il;
   if(size(#)>0)
   {
     il=#[1]; 
   }
   ideal re=1;
   option(redSB);
   list qr=simplifyIdeal(i);
   map phi=@P,qr[2];
   i=qr[1];
   
   list pr=facstd(i);

   if(size(pr)==1)
   {
      attrib(pr[1],"isSB",1);
      if((dim(pr[1])==0)&&(homog(pr[1])==1))
      {
         ideal @res=maxideal(1);
         return(phi(@res));
      }
      if(dim(pr[1])>1)
      {
         execute "ring gnir = ("+charstr(basering)+"),
                              ("+varstr(basering)+"),(C,lp);";
         ideal i=fetch(@P,i);
         list @pr=facstd(i);
         setring @P;
         pr=fetch(gnir,@pr);
      }
   }
   option(noredSB);
   int s=size(pr);

   if(s==1)
   {
     i=radicalEHV(i,ideal(1),il);
     return(phi(i));
   }
   intvec pos;
   pos[s]=0;
   if(il==1)
   {
     int ndim,k;
     attrib(pr[1],"isSB",1);
     int odim=dim(pr[1]);
     int count=1;
   
     for(j=2;j<=s;j++)
     {
        attrib(pr[j],"isSB",1);
        ndim=dim(pr[j]);
        if(ndim>odim)
        {
           for(k=count;k<=j-1;k++)
           {
              pos[k]=1;
           }
           count=j;
           odim=ndim;
        }
        if(ndim<odim)
        {
           pos[j]=1;
        }
     }
   }
   for(j=1;j<=s;j++)
   {
      if(pos[j]==0)
      {
         re=intersect(re,radicalEHV(pr[s+1-j],re,il));
      }
   }
   return(phi(re));
}

proc intersect1(ideal i,ideal j)
{
   def R=basering;
   execute "ring gnir = ("+charstr(basering)+"),
                        ("+varstr(basering)+",@t),(C,dp);";
   ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))-1)*imap(R,j);
   ideal j=eliminate(i,var(nvars(basering)));
   setring R;
   map phi=gnir,maxideal(1);
   return(phi(j));
}
 

///////////////////////////////////////////////////////////////////////////////
proc radicalKL (list m,ideal ser,list #)
"USAGE:   decomp(i); i ideal  (for primary decomposition)   (resp.
         decomp(i,1);        (for the minimal associated primes) )
RETURN:  list = list of primary ideals and their associated primes
         (at even positions in the list)
         (resp. a list of the minimal associated primes)
NOTE:    Algorithm of Gianni, Traeger, Zacharias
EXAMPLE: example decomp; shows an example
"
{
   ideal i=m[2];
  //----------------------------------------------------------------
  //i is the zero-ideal
  //----------------------------------------------------------------

   if(size(i)==0)
   {
      return(ideal(0));
   }
   
   def  @P = basering;
   list indep,allindep,restindep,fett,@mu;
   intvec @vh,isat;
   int @wr,@k,@n,@m,@n1,@n2,@n3,lauf,di;
   ideal @j,@j1,fac,@h,collectrad,htrad,lsau;
   ideal rad=ideal(1);
   ideal te=ser;

   poly @p,@q;
   string @va,quotring;
   int  homo=homog(i);

   if(size(#)>0)
   {
      @wr=#[1];
   }
   @j=m[1];
   @j1=m[2];
   int jdim=dim(@j);
   if(size(reduce(ser,@j,1))==0)
   {
      return(ser);
   }
   if(homo==1)
   {
      if(jdim==0)
      {
         option(noredSB);
         return(maxideal(1));
      }
      intvec @hilb=hilb(@j,1);
   }


  //----------------------------------------------------------------
  //j is the ring
  //----------------------------------------------------------------

   if (jdim==-1)
   {
      option(noredSB);
      return(ideal(0));
   }

  //----------------------------------------------------------------
  //  the case of one variable
  //----------------------------------------------------------------

   if(nvars(basering)==1)
   {
      fac=factorize(@j[1],1);
      @p=1;
      for(@k=1;@k<=size(fac);@k++)
      {
         @p=@p*fac[@k];
      }
      option(noredSB);

      return(ideal(@p));
   }   
   //------------------------------------------------------------------
   //the case of a complete intersection
   //------------------------------------------------------------------
    if(jdim+size(@j1)==nvars(basering))
    {
      // ideal jac=minor(jacob(@j1),size(@j1));
      // return(quotient(@j1,jac));
    }

   //------------------------------------------------------------------
   //the zero-dimensional case
   //------------------------------------------------------------------

   if (jdim==0)
   {
      @j1=finduni(@j); 
      for(@k=1;@k<=size(@j1);@k++)
      {
         fac=factorize(cleardenom(@j1[@k]),1);
         @p=fac[1];
         for(@n=2;@n<=size(fac);@n++)
         {
            @p=@p*fac[@n];
         }
         @j=@j,@p;
      }
      @j=std(@j);
      option(noredSB);
      return(@j);
   }
  
   //------------------------------------------------------------------
   //search for a maximal independent set indep,i.e.
   //look for subring such that the intersection with the ideal is zero
   //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
   //indep[1] is the new varstring and indep[2] the string for the block-ordering
   //------------------------------------------------------------------

   indep=maxIndependSet(@j);
  
   di=dim(@j);

   for(@m=1;@m<=size(indep);@m++)
   {
      if((indep[@m][1]==varstr(basering))&&(@m==1))
      //this is the good case, nothing to do, just to have the same notations
      //change the ring
      {
         execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
                              +ordstr(basering)+");";
         ideal @j=fetch(@P,@j);
         attrib(@j,"isSB",1);
      }
      else
      {
         @va=string(maxideal(1));
         execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
                              +indep[@m][2]+");";
         execute "map phi=@P,"+@va+";";
         if(homo==1)
         {
            ideal @j=std(phi(@j),@hilb);
         }
         else
         {
           ideal @j=std(phi(@j));
         }
      }
      if((deg(@j[1])==0)||(dim(@j)<jdim))
      {
         setring @P;
         break;
      }
      for (lauf=1;lauf<=size(@j);lauf++)
      {
         fett[lauf]=size(@j[lauf]);
      }
     //------------------------------------------------------------------------------------
     //we have now the following situation:
     //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
     //to this quotientring, j is their still a standardbasis, the
     //leading coefficients of the polynomials  there (polynomials in
     //K[var(nnp+1),..,var(nva)]) are collected in the list h,
     //we need their ggt, gh, because of the following:
     //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
     //intersected with K[var(1),...,var(nva)] is (j:gh^n)
     //on the other hand j=(j,gh^n) intersected with (j:gh^n)

     //------------------------------------------------------------------------------------

     //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
     //and the map phi:K[var(1),...,var(nva)] ----->K(var(nnpr+1),..,var(nva))[..the rest..]
     //-------------------------------------------------------------------------------------

      quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);

     //---------------------------------------------------------------------
     //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     //---------------------------------------------------------------------

      execute quotring;

    // @j considered in the quotientring
      ideal @j=imap(gnir1,@j);

      kill gnir1;

     //j is a standardbasis in the quotientring but usually not minimal
     //here it becomes minimal

      @j=clearSB(@j,fett);
      attrib(@j,"isSB",1);

     //we need later ggt(h[1],...)=gh for saturation
      ideal @h;
      if(deg(@j[1])>0)
      {
         for(@n=1;@n<=size(@j);@n++)
         {
            @h[@n]=leadcoef(@j[@n]);
         }
        //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
         option(redSB);
         @j=interred(@j);
         attrib(@j,"isSB",1);
         list  @mo=@j,@j;
         ideal zero_rad= radicalKL(@mo,ideal(1));
      }
      else
      {
         ideal zero_rad=ideal(1);
      }

     //we need the intersection of the ideals in the list quprimary with the
     //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
     //but fi polynomials, then the intersection of q with the polynomialring
     //is the saturation of the ideal generated by f1,...,fr with respect to
     //h which is the lcm of the leading coefficients of the fi considered in the
     //quotientring: this is coded in saturn

      ideal hpl;

      for(@n=1;@n<=size(zero_rad);@n++)
      {
         hpl=hpl,leadcoef(zero_rad[@n]);
      }

     //--------------------------------------------------------------------
     //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     //back to the polynomialring
     //---------------------------------------------------------------------
      setring @P;

      collectrad=imap(quring,zero_rad);
      lsau=simplify(imap(quring,hpl),2);
      @h=imap(quring,@h);

      kill quring;


     //here the intersection with the polynomialring
     //mentioned above is really computed

      collectrad=sat2(collectrad,lsau)[1];

      if(deg(@h[1])>=0)
      {
         fac=ideal(0);
         for(lauf=1;lauf<=ncols(@h);lauf++)   
         {
            if(deg(@h[lauf])>0)
            {
                fac=fac+factorize(@h[lauf],1);
            }
         }
         fac=simplify(fac,4);
         @q=1;
         for(lauf=1;lauf<=size(fac);lauf++)
         {
            @q=@q*fac[lauf];
         }

       
         @mu=mstd(quotient(@j+ideal(@q),rad));
         @j=@mu[1];
         attrib(@j,"isSB",1);
       
      }
      if((deg(rad[1])>0)&&(deg(collectrad[1])>0))
      {
         rad=intersect(rad,collectrad);
      }
      else
      {
         if(deg(collectrad[1])>0)
         {
            rad=collectrad;
         }
      }

      te=simplify(reduce(te*rad,@j),2);
 
      if((dim(@j)<di)||(size(te)==0))
      {
         break;
      }
      if(homo==1)
      {
         @hilb=hilb(@j,1);
      }
   }

   if(((@wr==1)&&(dim(@j)<di))||(deg(@j[1])==0)||(size(te)==0))
   {
      return(rad);
   }   
  // rad=intersect(rad,radicalKL(@mu,rad,@wr));
   rad=intersect(rad,radicalKL(@mu,ideal(1),@wr));


   option(noredSB);
   return(rad);
}


example
{ "EXAMPLE:"; echo = 2;
}


proc radicalEHV(ideal i,ideal re,list #)
{
   ideal J,I,I0,radI0,L,radI1,I2,radI2;
   int l,il;
   if(size(#)>0)
   {
      il=#[1];
   }
   
   option(redSB);
   list m=mstd(i);
        I=m[2];
   option(noredSB);
   if(size(reduce(re,m[1],1))==0)
   {
      return(re);
   }
   int cod=nvars(basering)-dim(m[1]);
   if((nvars(basering)<=5)&&(size(m[2])<=5))
   {
      if(cod==size(m[2]))
      {
        J=minor(jacob(I),cod);
        return(quotient(I,J));
      }
   
      for(l=1;l<=cod;l++)
      {
         I0[l]=I[l];
      }
      if(dim(std(I0))+cod==nvars(basering))
      {
         J=minor(jacob(I0),cod);
         radI0=quotient(I0,J);
         L=quotient(radI0,I);
         radI1=quotient(radI0,L);

         if(size(reduce(radI1,m[1],1))==0)
         {
            return(I);
         }
         if(il==1)
         {
      
            return(radI1);
         }

         I2=sat(I,radI1)[1];

         if(deg(I2[1])<=0)
         {
            return(radI1);
         }
         return(intersect(radI1,radicalEHV(I2,re,il)));       
      }
   }
   return(radicalKL(m,re,il));
}

proc Ann(module M)
{
  M=prune(M);  //to obtain a small embedding
  ideal ann=quotient1(M,freemodule(nrows(M)));
  return(ann);
}

//computes the equidimensional part of the ideal i of codimension e
proc int_ass_primary_e(ideal i, int e)
{
  if(homog(i)!=1)
  {
     i=std(i);
  }
  list re=sres(i,0);                   //the resolution
  re=minres(re);                       //minimized resolution
  ideal ann=AnnExt_R(e,re);
  if(nvars(basering)-dim(std(ann))!=e)
  {
    return(ideal(1));
  }
  return(ann);
} 
 
//computes all equidimensional parts of the ideal i
proc prepareAss(ideal i)
{
  ideal j=std(i);
  int cod=nvars(basering)-dim(j);
  int e;
  list er;
  ideal ann;
  if(homog(i)==1)
  {
     list re=sres(i,0);                   //the resolution
     re=minres(re);                       //minimized resolution
  }
  else
  {
     list re=mres(i,0);             //fehler in sres
  }  
  for(e=cod;e<=nvars(basering);e++)
  {
     ann=AnnExt_R(e,re);

     if(nvars(basering)-dim(std(ann))==e)
     {
        er[size(er)+1]=equiRadical(ann);
     }
  }
  return(er);
}  

//computes the annihilator of Ext^n(R/i,R) with given resolution re
//n is not necessarily the number of variables
proc AnnExt_R(int n,list re)
{
  if(n<nvars(basering))
  {
     matrix f=transpose(re[n+1]);      //Hom(_,R)
     module k=res(f,2)[2];             //the kernel
     matrix g=transpose(re[n]);        //the image of Hom(_,R)

     ideal ann=quotient(g,k);           //the anihilator
  }
  else
  {
     ideal ann=Ann(transpose(re[n]));
  }
  return(ann);            
}

proc quotient1(module a,module b)
{
   int i;
   ideal re=quotient(a,module(b[1]));
   for(i=2;i<=size(b);i++)
   {
      re=intersect1(re,quotient(a,module(b[i])));
   }
   return(re);   
}



proc analyze(list pr)
{ 
   int ii,jj;
   for(ii=1;ii<=size(pr)/2;ii++)
   {
      dim(std(pr[2*ii]));
      idealsEqual(pr[2*ii-1],pr[2*ii]);
      "===========================";
   }

   for(ii=size(pr)/2;ii>1;ii--)
   {
      for(jj=1;jj<ii;jj++)
      {
         if(size(reduce(pr[2*jj],std(pr[2*ii],1)))==0)
         {
            "eingebette Komponente";
            jj;
            ii;
         }
      }
   }   
}


proc simplifyIdeal(ideal i)
{
  def r=basering;
  
  int j,k;
  map phi;
  poly p;
  
  ideal iwork=i;
  ideal imap1=maxideal(1);
  ideal imap2=maxideal(1);
  

  for(j=1;j<=nvars(basering);j++)
  {
    for(k=1;k<=size(i);k++)
    {
      if(deg(iwork[k]/var(j))==0)
      {
        p=-1/leadcoef(iwork[k]/var(j))*iwork[k];
        imap1[j]=p+2*var(j);
        phi=r,imap1;
        iwork=phi(iwork);
        iwork=subst(iwork,var(j),0);
        iwork[k]=var(j);
        imap1=maxideal(1);
        imap2[j]=-p;        
        break;
      }
    }
  }
  return(iwork,imap2);
}

        
///////////////////////////////////////////////////////
// ini_mod
// input: a polynomial p
// output: the initial term of p as needed
// in the context of characteristic sets
//////////////////////////////////////////////////////

proc ini_mod(poly p)
{
  if (p==0)
  {
    return(0);
  }
  int n; matrix m;
  for( n=nvars(basering); n>0; n=n-1)
  {
    m=coef(p,var(n));
    if(m[1,1]!=1)
    {
      p=m[2,1];
      break;
    }
  }
  if(deg(p)==0)
  {
    p=0;
  }
  return(p);
}
///////////////////////////////////////////////////////
// min_ass_prim_charsets
// input: generators of an ideal PS and an integer cho
// If cho=0, the given ordering of the variables is used.
// Otherwise, the system tries to find an "optimal ordering",
// which in some cases may considerably speed up the algorithm
// output: the minimal associated primes of PS
// algorithm: via characteriostic sets
//////////////////////////////////////////////////////


proc min_ass_prim_charsets (ideal PS, int cho)
{
  if((cho<0) and (cho>1))
    {
      "ERROR: <int> must be 0 or 1"
      return();
    }
  if(system("version")>933)
  {
    option(notWarnSB);
  }
  if(cho==0)
  {
    return(min_ass_prim_charsets0(PS));
  }
  else
  {
     return(min_ass_prim_charsets1(PS));
  }
}
///////////////////////////////////////////////////////
// min_ass_prim_charsets0
// input: generators of an ideal PS
// output: the minimal associated primes of PS
// algorithm: via characteristic sets
// the given ordering of the variables is used
//////////////////////////////////////////////////////


proc min_ass_prim_charsets0 (ideal PS)
{

  matrix m=char_series(PS);  // We compute an irreducible
                             // characteristic series
  int i,j,k;
  list PSI;
  list PHI;  // the ideals given by the characteristic series
  for(i=nrows(m);i>=1; i--)
  {
     PHI[i]=ideal(m[i,1..ncols(m)]);
  }
  // We compute the radical of each ideal in PHI
  ideal I,JS,II;
  int sizeJS, sizeII;
  for(i=size(PHI);i>=1; i--)
  {
     I=0;
     for(j=size(PHI[i]);j>0;j=j-1)
     {
       I=I+ini_mod(PHI[i][j]);
     }
     JS=std(PHI[i]);
     sizeJS=size(JS);
     for(j=size(I);j>0;j=j-1)
     {
       II=0;
       sizeII=0;
       k=0;
       while(k<=sizeII)                  // successive saturation
       {
         option(returnSB);
         II=quotient(JS,I[j]);
         option(noreturnSB);
//std
//         II=std(II);
         sizeII=size(II);
         if(sizeII==sizeJS)
         {
            for(k=1;k<=sizeII;k++)
            {
              if(leadexp(II[k])!=leadexp(JS[k])) break;
            }
         }
         JS=II;
         sizeJS=sizeII;
       }
    }
    PSI=insert(PSI,JS);
  }
  int sizePSI=size(PSI);
  // We eliminate redundant ideals
  for(i=1;i<sizePSI;i++)
  {
    for(j=i+1;j<=sizePSI;j++)
    {
      if(size(PSI[i])!=0)
      {
        if(size(PSI[j])!=0)
        {
          if(size(NF(PSI[i],PSI[j],1))==0)
          {
            PSI[j]=ideal(0);
          }
          else
          {
            if(size(NF(PSI[j],PSI[i],1))==0)
            {
              PSI[i]=ideal(0);
            }
          }
        }
      }
    }
  }
  for(i=sizePSI;i>=1;i--)
  {
    if(size(PSI[i])==0)
    {
      PSI=delete(PSI,i);
    }
  }
  return (PSI);
}

///////////////////////////////////////////////////////
// min_ass_prim_charsets1
// input: generators of an ideal PS
// output: the minimal associated primes of PS
// algorithm: via characteristic sets
// input: generators of an ideal PS and an integer i
// The system tries to find an "optimal ordering" of
// the variables
//////////////////////////////////////////////////////


proc min_ass_prim_charsets1 (ideal PS)
{
  def oldring=basering;
  string n=system("neworder",PS);
  execute "ring r="+charstr(oldring)+",("+n+"),dp;";
  ideal PS=imap(oldring,PS);
  matrix m=char_series(PS);  // We compute an irreducible
                             // characteristic series
  int i,j,k;
  ideal I;
  list PSI;
  list PHI;    // the ideals given by the characteristic series
  list ITPHI;  // their initial terms
  for(i=nrows(m);i>=1; i--)
  {
     PHI[i]=ideal(m[i,1..ncols(m)]);
     I=0;
     for(j=size(PHI[i]);j>0;j=j-1)
     {
       I=I,ini_mod(PHI[i][j]);
     }
      I=I[2..ncols(I)];
      ITPHI[i]=I;
  }
  setring oldring;
  matrix m=imap(r,m);
  list PHI=imap(r,PHI);
  list ITPHI=imap(r,ITPHI);
  // We compute the radical of each ideal in PHI
  ideal I,JS,II;
  int sizeJS, sizeII;
  for(i=size(PHI);i>=1; i--)
  {
     I=0;
     for(j=size(PHI[i]);j>0;j=j-1)
     {
       I=I+ITPHI[i][j];
     }
     JS=std(PHI[i]);
     sizeJS=size(JS);
     for(j=size(I);j>0;j=j-1)
     {
       II=0;
       sizeII=0;
       k=0;
       while(k<=sizeII)                  // successive iteration
       {
         option(returnSB);
         II=quotient(JS,I[j]);
         option(noreturnSB);
//std
//         II=std(II);
         sizeII=size(II);
         if(sizeII==sizeJS)
         {
            for(k=1;k<=sizeII;k++)
            {
              if(leadexp(II[k])!=leadexp(JS[k])) break;
            }
         }
         JS=II;
         sizeJS=sizeII;
       }
    }
    PSI=insert(PSI,JS);
  }
  int sizePSI=size(PSI);
  // We eliminate redundant ideals
  for(i=1;i<sizePSI;i++)
  {
    for(j=i+1;j<=sizePSI;j++)
    {
      if(size(PSI[i])!=0)
      {
        if(size(PSI[j])!=0)
        {
          if(size(NF(PSI[i],PSI[j],1))==0)
          {
            PSI[j]=ideal(0);
          }
          else
          {
            if(size(NF(PSI[j],PSI[i],1))==0)
            {
              PSI[i]=ideal(0);
            }
          }
        }
      }
    }
  }
  for(i=sizePSI;i>=1;i--)
  {
    if(size(PSI[i])==0)
    {
      PSI=delete(PSI,i);
    }
  }
  return (PSI);
}


/////////////////////////////////////////////////////
// proc prim_dec
// input:  generators of an ideal I and an integer choose
// If choose=0, min_ass_prim_charsets with the given
// ordering of the variables is used.
// If choose=1, min_ass_prim_charsets with the "optimized"
// ordering of the variables is used.
// If choose=2, minAssPrimes from primdec.lib is used
// If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
// output: a primary decomposition of I, i.e., a list
// of pairs consisting of a standard basis of a primary component
// of I and a standard basis of the corresponding associated prime.
// To compute the minimal associated primes of a given ideal
// min_ass_prim_l is called, i.e., the minimal associated primes
// are computed via characteristic sets.
// In the homogeneous case, the performance of the procedure
// will be improved if I is already given by a minimal set of
// generators. Apply minbase if necessary.
//////////////////////////////////////////////////////////


proc prim_dec(ideal I, int choose)
{
   if((choose<0) or (choose>3))
   {
     "ERROR: <int> must be 0 or 1 or 2 or 3";
      return();
   }
   if(system("version")>933)
   {
      option(notWarnSB);
    }
  ideal H=1; // The intersection of the primary components
  list U;    // the leaves of the decomposition tree, i.e.,
             // pairs consisting of a primary component of I
             // and the corresponding associated prime
  list W;    // the non-leaf vertices in the decomposition tree.
             // every entry has 6 components:
                // 1- the vertex itself , i.e., a standard bais of the
                //    given ideal I (type 1), or a standard basis of a
                //    pseudo-primary component arising from
                //    pseudo-primary decomposition (type 2), or a
                //    standard basis of a remaining component arising from
                //    pseudo-primary decomposition or extraction (type 3)
                // 2- the type of the vertex as indicated above
                // 3- the weighted_tree_depth of the vertex
                // 4- the tester of the vertex
                // 5- a standard basis of the associated prime
                //    of a vertex of type 2, or 0 otherwise
                // 6- a list of pairs consisting of a standard
                //    basis of a minimal associated prime ideal
                //    of the father of the vertex and the
                //    irreducible factors of the "minimal
                //    divisor" of the seperator or extractor
                //    corresponding to the prime ideal
                //    as computed by the procedure minsat,
                //    if the vertex is of type 3, or
                //    the empty list otherwise
  ideal SI=std(I);
  int ncolsSI=ncols(SI);
  int ncolsH=1;
  W[1]=list(I,1,0,poly(1),ideal(0),list()); // The root of the tree
  int weighted_tree_depth;
  int i,j;
  int check;
  list V;  // current vertex
  list VV; // new vertex
  list QQ;
  list WI;
  ideal Qi,SQ,SRest,fac;
  poly tester;

  while(1)
  {
    i=1;
    while(1)
    {
      while(i<=size(W)) // find vertex V of smallest weighted tree-depth
      {
        if (W[i][3]<=weighted_tree_depth) break;
        i++;
      }
      if (i<=size(W)) break;
      i=1;
      weighted_tree_depth++;
    }
    V=W[i];
    W=delete(W,i); // delete V from W

    // now proceed by type of vertex V

    if (V[2]==2)  // extraction needed
    {
      SQ,SRest,fac=extraction(V[1],V[5]);
                        // standard basis of primary component,
                        // standard basis of remaining component,
                        // irreducible factors of
                        // the "minimal divisor" of the extractor
                        // as computed by the procedure minsat,
      check=0;
      for(j=1;j<=ncolsH;j++)
      {
        if (NF(H[j],SQ,1)!=0) // Q is not redundant
        {
          check=1;
          break;
        }
      }
      if(check==1)             // Q is not redundant
      {
        QQ=list();
        QQ[1]=list(SQ,V[5]);  // primary component, associated prime,
                              // i.e., standard bases thereof
        U=U+QQ;
        H=intersect(H,SQ);
        H=std(H);
        ncolsH=ncols(H);
        check=0;
        if(ncolsH==ncolsSI)
        {
          for(j=1;j<=ncolsSI;j++)
          {
            if(leadexp(H[j])!=leadexp(SI[j]))
            {
              check=1;
              break;
            }
          }
        }
        else
        {
          check=1;
        }
        if(check==0) // H==I => U is a primary decomposition
        {
          return(U);
        }
      }
      if (SRest[1]!=1)        // the remaining component is not
                              // the whole ring
      {
        if (rad_con(V[4],SRest)==0) // the new vertex is not the
                                    // root of a redundant subtree
        {
          VV[1]=SRest;     // remaining component
          VV[2]=3;         // pseudoprimdec_special
          VV[3]=V[3]+1;    // weighted depth
          VV[4]=V[4];      // the tester did not change
          VV[5]=ideal(0);
          VV[6]=list(list(V[5],fac));
          W=insert(W,VV,size(W));
        }
      }
    }
    else
    {
      if (V[2]==3) // pseudo_prim_dec_special is needed
      {
        QQ,SRest=pseudo_prim_dec_special_charsets(V[1],V[6],choose);
                         // QQ = quadruples:
                         // standard basis of pseudo-primary component,
                         // standard basis of corresponding prime,
                         // seperator, irreducible factors of
                         // the "minimal divisor" of the seperator
                         // as computed by the procedure minsat,
                         // SRest=standard basis of remaining component
      }
      else     // V is the root, pseudo_prim_dec is needed
      {
        QQ,SRest=pseudo_prim_dec_charsets(I,SI,choose);
                         // QQ = quadruples:
                         // standard basis of pseudo-primary component,
                         // standard basis of corresponding prime,
                         // seperator, irreducible factors of
                         // the "minimal divisor" of the seperator
                         // as computed by the procedure minsat,
                         // SRest=standard basis of remaining component

      }
//check
      for(i=size(QQ);i>=1;i--)
      //for(i=1;i<=size(QQ);i++)
      {
        tester=QQ[i][3]*V[4];
        Qi=QQ[i][2];
        if(NF(tester,Qi,1)!=0)  // the new vertex is not the
                                // root of a redundant subtree
        {
          VV[1]=QQ[i][1];
          VV[2]=2;
          VV[3]=V[3]+1;
          VV[4]=tester;      // the new tester as computed above
          VV[5]=Qi;          // QQ[i][2];
          VV[6]=list();
          W=insert(W,VV,size(W));
        }
      }
      if (SRest[1]!=1)        // the remaining component is not
                              // the whole ring
      {
        if (rad_con(V[4],SRest)==0) // the vertex is not the root
                                    // of a redundant subtree
        {
          VV[1]=SRest;
          VV[2]=3;
          VV[3]=V[3]+2;
          VV[4]=V[4];      // the tester did not change
          VV[5]=ideal(0);
          WI=list();
          for(i=1;i<=size(QQ);i++)
          {
            WI=insert(WI,list(QQ[i][2],QQ[i][4]));
          }
          VV[6]=WI;
          W=insert(W,VV,size(W));
        }
      }
    }
  }
}

//////////////////////////////////////////////////////////////////////////
// proc pseudo_prim_dec_charsets
// input: Generators of an arbitrary ideal I, a standard basis SI of I,
// and an integer choo
// If choo=0, min_ass_prim_charsets with the given
// ordering of the variables is used.
// If choo=1, min_ass_prim_charsets with the "optimized"
// ordering of the variables is used.
// If choo=2, minAssPrimes from primdec.lib is used
// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
// output: a pseudo primary decomposition of I, i.e., a list
// of pseudo primary components together with a standard basis of the
// remaining component. Each pseudo primary component is
// represented by a quadrupel: A standard basis of the component,
// a standard basis of the corresponding associated prime, the
// seperator of the component, and the irreducible factors of the
// "minimal divisor" of the seperator as computed by the procedure minsat,
// calls  proc pseudo_prim_dec_i
//////////////////////////////////////////////////////////////////////////


proc pseudo_prim_dec_charsets (ideal I, ideal SI, int choo)
{
  list L;          // The list of minimal associated primes,
                   // each one given by a standard basis
  if((choo==0) or (choo==1))
    {
       L=min_ass_prim_charsets(I,choo);
    }
  else
    {
      if(choo==2)
      {
        L=minAssPrimes(I);
      }
      else
      {
        L=minAssPrimes(I,1);
      }
      for(int i=size(L);i>=1;i=i-1)
        {
          L[i]=std(L[i]);
        }
    }
  return (pseudo_prim_dec_i(SI,L));
}

////////////////////////////////////////////////////////////////
// proc pseudo_prim_dec_special_charsets
// input: a standard basis of an ideal I whose radical is the
// intersection of the radicals of ideals generated by one prime ideal
// P_i together with one polynomial f_i, the list V6 must be the list of
// pairs (standard basis of P_i, irreducible factors of f_i),
// and an integer choo
// If choo=0, min_ass_prim_charsets with the given
// ordering of the variables is used.
// If choo=1, min_ass_prim_charsets with the "optimized"
// ordering of the variables is used.
// If choo=2, minAssPrimes from primdec.lib is used
// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
// output: a pseudo primary decomposition of I, i.e., a list
// of pseudo primary components together with a standard basis of the
// remaining component. Each pseudo primary component is
// represented by a quadrupel: A standard basis of the component,
// a standard basis of the corresponding associated prime, the
// seperator of the component, and the irreducible factors of the
// "minimal divisor" of the seperator as computed by the procedure minsat,
// calls  proc pseudo_prim_dec_i
////////////////////////////////////////////////////////////////


proc pseudo_prim_dec_special_charsets (ideal SI,list V6, int choo)
{
  int i,j,l;
  list m;
  list L;
  int sizeL;
  ideal P,SP; ideal fac;
  int dimSP;
  for(l=size(V6);l>=1;l--)   // creates a list of associated primes
                             // of I, possibly redundant
  {
    P=V6[l][1];
    fac=V6[l][2];
    for(i=ncols(fac);i>=1;i--)
    {
      SP=P+fac[i];
      SP=std(SP);
      if(SP[1]!=1)
      {
        if((choo==0) or (choo==1))
        {
          m=min_ass_prim_charsets(SP,choo);  // a list of SB
        }
        else
        {
          if(choo==2)
          {
            m=minAssPrimes(SP);
          }
          else
          {
            m=minAssPrimes(SP,1);
          }
          for(j=size(m);j>=1;j=j-1)
            {
              m[j]=std(m[j]);
            }
        }
        dimSP=dim(SP);
        for(j=size(m);j>=1; j--)
        {
          if(dim(m[j])==dimSP)
          {
            L=insert(L,m[j],size(L));
          }
        }
      }
    }
  }
  sizeL=size(L);
  for(i=1;i<sizeL;i++)     // get rid of redundant primes
  {
    for(j=i+1;j<=sizeL;j++)
    {
      if(size(L[i])!=0)
      {
        if(size(L[j])!=0)
        {
          if(size(NF(L[i],L[j],1))==0)
          {
            L[j]=ideal(0);
          }
          else
          {
            if(size(NF(L[j],L[i],1))==0)
            {
              L[i]=ideal(0);
            }
          }
        }
      }
    }
  }
  for(i=sizeL;i>=1;i--)
  {
    if(size(L[i])==0)
    {
      L=delete(L,i);
    }
  }
  return (pseudo_prim_dec_i(SI,L));
}


////////////////////////////////////////////////////////////////
// proc pseudo_prim_dec_i
// input: A standard basis of an arbitrary ideal I, and standard bases
// of the minimal associated primes of I
// output: a pseudo primary decomposition of I, i.e., a list
// of pseudo primary components together with a standard basis of the
// remaining component. Each pseudo primary component is
// represented by a quadrupel: A standard basis of the component Q_i,
// a standard basis of the corresponding associated prime P_i, the
// seperator of the component, and the irreducible factors of the
// "minimal divisor" of the seperator as computed by the procedure minsat,
////////////////////////////////////////////////////////////////


proc pseudo_prim_dec_i (ideal SI, list L)
{
  list Q;
  if (size(L)==1)               // one minimal associated prime only
                                // the ideal is already pseudo primary
  {
    Q=SI,L[1],1;
    list QQ;
    QQ[1]=Q;
    return (QQ,ideal(1));
  }

  poly f0,f,g;
  ideal fac;
  int i,j,k,l;
  ideal SQi;
  ideal I'=SI;
  list QP;
  int sizeL=size(L);
  for(i=1;i<=sizeL;i++)
  {
    fac=0;
    for(j=1;j<=sizeL;j++)           // compute the seperator sep_i
                                    // of the i-th component
    {
      if (i!=j)                       // search g not in L[i], but L[j]
      {
        for(k=1;k<=ncols(L[j]);k++)
        {
          if(NF(L[j][k],L[i],1)!=0)
          {
            break;
          }
        }
        fac=fac+L[j][k];
      }
    }
    // delete superfluous polynomials
    fac=simplify(fac,8);
    // saturation
    SQi,f0,f,fac=minsat_ppd(SI,fac);
    I'=I',f;
    QP=SQi,L[i],f0,fac;
             // the quadrupel:
             // a standard basis of Q_i,
             // a standard basis of P_i,
             // sep_i,
             // irreducible factors of
             // the "minimal divisor" of the seperator
             //  as computed by the procedure minsat,
    Q[i]=QP;
  }
  I'=std(I');
  return (Q, I');
                   // I' = remaining component
}


////////////////////////////////////////////////////////////////
// proc extraction
// input: A standard basis of a pseudo primary ideal I, and a standard
// basis of the unique minimal associated prime P of I
// output: an extraction of I, i.e., a standard basis of the primary
// component Q of I with associated prime P, a standard basis of the
// remaining component, and the irreducible factors of the
// "minimal divisor" of the extractor as computed by the procedure minsat
////////////////////////////////////////////////////////////////


proc extraction (ideal SI, ideal SP)
{
  list indsets=system("indsetall",SP,0);
  poly f;
  if(size(indsets)!=0)      //check, whether dim P != 0
  {
    intvec v;               // a maximal independent set of variables
                            // modulo P
    string U;               // the independent variables
    string A;               // the dependent variables
    int j,k;
    int a;                  //  the size of A
    int degf;
    ideal g;
    list polys;
    int sizepolys;
    list newpoly;
    def R=basering;
    //intvec hv=hilb(SI,1);
    for (k=1;k<=size(indsets);k++)
    {
      v=indsets[k];
      for (j=1;j<=nvars(R);j++)
      {
        if (v[j]==1)
        {
          U=U+varstr(j)+",";
        }
        else
        {
          A=A+varstr(j)+",";
          a++;
        }
      }

      U[size(U)]=")";           // we compute the extractor of I (w.r.t. U)
      execute "ring RAU="+charstr(basering)+",("+A+U+",(dp("+string(a)+"),dp);";
      ideal I=imap(R,SI);
      //I=std(I,hv);            // the standard basis in (R[U])[A]
      I=std(I);            // the standard basis in (R[U])[A]
      A[size(A)]=")";
      execute "ring Rloc=("+charstr(basering)+","+U+",("+A+",dp;";
      ideal I=imap(RAU,I);
      //"std in lokalisierung:"+newline,I;
      ideal h;
      for(j=ncols(I);j>=1;j--)
      {
        h[j]=leadcoef(I[j]);  // consider I in (R(U))[A]
      }
      setring R;
      g=imap(Rloc,h);
      kill RAU,Rloc;
      U="";
      A="";
      a=0;
      f=lcm(g);
      newpoly[1]=f;
      polys=polys+newpoly;
      newpoly=list();
    }
    f=polys[1];
    degf=deg(f);
    sizepolys=size(polys);
    for (k=2;k<=sizepolys;k++)
    {
      if (deg(polys[k])<degf)
      {
        f=polys[k];
        degf=deg(f);
      }
    }
  }
  else
  {
    f=1;
  }
  poly f0,h0; ideal SQ; ideal fac;
  if(f!=1)
  {
    SQ,f0,h0,fac=minsat(SI,f);
    return(SQ,std(SI+h0),fac);
             // the tripel
             // a standard basis of Q,
             // a standard basis of remaining component,
             // irreducible factors of
             // the "minimal divisor" of the extractor
             // as computed by the procedure minsat
  }
  else
  {
    return(SI,ideal(1),ideal(1));
  }
}

/////////////////////////////////////////////////////
// proc minsat
// input:  a standard basis of an ideal I and a polynomial p
// output: a standard basis IS of the saturation of I w.r. to p,
// the maximal squarefree factor f0 of p,
// the "minimal divisor" f of f0 such that the saturation of
// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
// the irreducible factors of f
//////////////////////////////////////////////////////////


proc minsat(ideal SI, poly p)
{
  ideal fac=factorize(p,1);       //the irreducible factors of p
  fac=sort(fac)[1];
  int i,k;
  poly f0=1;
  for(i=ncols(fac);i>=1;i--)
  {
    f0=f0*fac[i];
  }
  poly f=1;
  ideal iold;
  list quotM;
  quotM[1]=SI;
  quotM[2]=fac;
  quotM[3]=f0;
  // we deal seperately with the first quotient;
  // factors, which do not contribute to this one,
  // are omitted
  iold=quotM[1];
  quotM=minquot(quotM);
  fac=quotM[2];
  if(quotM[3]==1)
    {
      return(quotM[1],f0,f,fac);
    }
  while(special_ideals_equal(iold,quotM[1])==0)
    {
      f=f*quotM[3];
      iold=quotM[1];
      quotM=minquot(quotM);
    }
  return(quotM[1],f0,f,fac);           // the quadrupel ((I:p),f0,f, irr. factors of f)
}

/////////////////////////////////////////////////////
// proc minsat_ppd
// input:  a standard basis of an ideal I and a polynomial p
// output: a standard basis IS of the saturation of I w.r. to p,
// the maximal squarefree factor f0 of p,
// the "minimal divisor" f of f0 such that the saturation of
// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
// the irreducible factors of f
//////////////////////////////////////////////////////////


proc minsat_ppd(ideal SI, ideal fac)
{
  fac=sort(fac)[1];
  int i,k;
  poly f0=1;
  for(i=ncols(fac);i>=1;i--)
  {
    f0=f0*fac[i];
  }
  poly f=1;
  ideal iold;
  list quotM;
  quotM[1]=SI;
  quotM[2]=fac;
  quotM[3]=f0;
  // we deal seperately with the first quotient;
  // factors, which do not contribute to this one,
  // are omitted
  iold=quotM[1];
  quotM=minquot(quotM);
  fac=quotM[2];
  if(quotM[3]==1)
    {
      return(quotM[1],f0,f,fac);
    }
  while(special_ideals_equal(iold,quotM[1])==0)
  {
    f=f*quotM[3];
    iold=quotM[1];
    quotM=minquot(quotM);
    k++;
  }
  return(quotM[1],f0,f,fac);           // the quadrupel ((I:p),f0,f, irr. factors of f)
}
/////////////////////////////////////////////////////////////////
// proc minquot
// input: a list with 3 components: a standard basis
// of an ideal I, a set of irreducible polynomials, and
// there product f0
// output: a standard basis of the ideal (I:f0), the irreducible
// factors of the "minimal divisor" f of f0 with (I:f0) = (I:f),
// the "minimal divisor" f
/////////////////////////////////////////////////////////////////

proc minquot(list tsil)
{
   int i,j,k,action;
   ideal verg;
   list l;
   poly g;
   ideal laedi=tsil[1];
   ideal fac=tsil[2];
   poly f=tsil[3];

//std
//   ideal star=quotient(laedi,f);
//   star=std(star);
   option(returnSB);
   ideal star=quotient(laedi,f);
   option(noreturnSB);
   if(special_ideals_equal(laedi,star)==1)
     {
       return(laedi,ideal(1),1);
     }
   action=1;
   while(action==1)
   {
      if(size(fac)==1)
      {
         action=0;
         break;
      }
      for(i=1;i<=size(fac);i++)
      {
        g=1;
         for(j=1;j<=size(fac);j++)
         {
            if(i!=j)
            {
               g=g*fac[j];
            }
         }
//std
//         verg=quotient(laedi,g);
//         verg=std(verg);
         option(returnSB);
         verg=quotient(laedi,g);
         option(noreturnSB);
         if(special_ideals_equal(verg,star)==1)
         {
            f=g;
            fac[i]=0;
            fac=simplify(fac,2);
            break;
         }
         if(i==size(fac))
         {
            action=0;
         }
      }
   }
   l=star,fac,f;
   return(l);
}
/////////////////////////////////////////////////
// proc special_ideals_equal
// input: standard bases of ideal k1 and k2 such that
// k1 is contained in k2, or k2 is contained ink1
// output: 1, if k1 equals k2, 0 otherwise
//////////////////////////////////////////////////

proc special_ideals_equal( ideal k1, ideal k2)
{
   int j;
   if(size(k1)==size(k2))
   {
      for(j=1;j<=size(k1);j++)
      {
         if(leadexp(k1[j])!=leadexp(k2[j]))
         {
            return(0);
         }
      }
      return(1);
   }
   return(0);
}