source: git/Singular/LIB/primdec.lib @ 7475b3f

///////////////////////////////////////////////////////////////////////////////
version="$Id: primdec.lib,v 1.82 2000-12-31 02:02:05 greuel Exp $";
category="Commutative Algebra";
info="
LIBRARY: primdec.lib   Primary Decomposition and Radical of Ideals
AUTHORS:  Gerhard Pfister, email: pfister@mathematik.uni-kl.de (GTZ)
          Wolfram Decker, email: decker@math.uni-sb.de         (SY)
          Hans Schoenemann, email: hannes@mathematik.uni-kl.de (SY)

OVERVIEW: 
 Algorithms for primary decomposition based on the ideas of 
 Gianni,Trager,Zacharias was written by Gerhard Pfister.                       
 Algorithms for primary decomposition based on the ideas of
 Shimoyama/Yokoyama was written by Wolfram Decker and Hans Schoenemann.   
 These procedures are implemented to be used in characteristic 0.
 @*They also work in positive characteristic >> 0.
 @*In small characteristic and for algebraic extensions, primdecGTZ,
 minAssGTZ, radical and equiRadical may not terminate and primdecSY and
 minAssChar may not give a complete decomposition.

PROCEDURES:
 primdecGTZ(I);    complete primary decomposition via Gianni,Trager,Zacharias
 primdecSY(I...);  complete primary decomposition via Shimoyama-Yokoyama
 minAssGTZ(I);     the minimal associated primes via Gianni,Trager,Zacharias
 minAssChar(I...); the minimal associated primes using characteristic sets
 testPrimary(L,k); tests the result of the primary decomposition
 radical(I);       computes the radical of the ideal I
 equiRadical(I);   the radical of the equidimensional part of the ideal I
 prepareAss(I);    list of radicals of the equidimensional components of I
 equidim(I);       weak equidimensional decomposition of I
 equidimMax(I);    equidimensional locus of I
 equidimMaxEHV(I); equidimensional locus of I via Eisenbud,Huneke,Vasconcelos
 zerodec(I);       zerodimensional decomposition via Monico
";

LIB "general.lib";
LIB "elim.lib";
LIB "poly.lib";
LIB "random.lib";
LIB "inout.lib";
LIB "matrix.lib";
///////////////////////////////////////////////////////////////////////////////
//
//                      Gianni/Trager/Zacharias
//
///////////////////////////////////////////////////////////////////////////////

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);

   if(specialIdealsEqual(star,laedi))
   {
      l=star,fac,f;
      return(l);      
   }
  
   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 primaryTest (ideal i, poly p)
{
   int m=1;
   int n=nvars(basering);
   int e,f;
   poly t;
   ideal h;
   list act;

   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]));
      if(char(basering)!=0)
      {
         f=1;
         if(e mod char(basering)==0)
         {
           if ( voice >=15 )
           { 
              "// WARNING: The characteristic is perhaps too small to use";
              "// the algorithm of Gianni/Trager/Zacharias.";
              "// This may result in an infinte loop";
              "// loop in primaryTest, voice:",voice;"";
           }
           while(e mod char(basering)==0)
           {
              f=f*char(basering);
              e=e/char(basering);  
           }
       
         }
         t=leadcoef(i[m])*e*var(n)^f+(i[m]-lead(i[m]))/var(n)^((e-1)*f);
         i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
         if (reduce(i[m]-t^e,prm,1) !=0)
         {
           return(ideal(0));
         }
         if(f>1)
         {       
            act=factorize(t);
            if(size(act[1])>2)
            {
              return(ideal(0));   
            }
            if(deg(act[1][2])>1)
            {
              return(ideal(0));   
            }
            t=act[1][2];
         }
      }
      else
      {
         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 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;
   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 splitCharp(list l)
{
  if((char(basering)==0)||(npars(basering)>0))
  {
     return(l);
  }
  def P=basering;
  int i,j,k,m,q,d;
  int n=nvars(basering);
  ideal s,t,u,sact;
  poly ni;
  string minp,gnir,va;
  list sa,keep;
  for(i=1;i<=size(l)/2;i++)
  {
    if(size(l[2*i])==0)
    {
       if(deg(l[2*i-1][1])==vdim(l[2*i-1]))
       {
          l[2*i]=l[2*i-1];
       }
    }
  }
  for(i=1;i<=size(l)/2;i++)
  {
    if(size(l[2*i])==0)
    {
      s=factorize(l[2*i-1][1],1);
      t=l[2*i-1];
      m=size(t);
      ni=s[1];
      if(deg(ni)>1)
      {
        va=varstr(P);
        j=size(va);
        while(va[j]!=","){j--;}
        va=va[1..j-1];
        gnir="ring RL=("+string(char(P))+","+string(var(n))+"),("+va+"),dp;";
        execute(gnir);
        minpoly=leadcoef(imap(P,ni));
        ideal act;
        ideal t=imap(P,t);
        for(k=2;k<=m;k++)
        { 
           act=factorize(t[k],1);
           if(size(act)>1)
           {

             break;
           }
        }
        setring P;
        sact=imap(RL,act);
        kill RL;
        if(size(sact)>1)
        {
           sa=sat1(l[2*i-1],sact[1]);
           keep[size(keep)+1]=std(l[2*i-1],sa[2]);
           l[2*i-1]=std(sa[1]);
           l[2*i]=primaryTest(sa[1],sa[1][1]);
        }
      }
    }
  }
  if(size(keep)>0)
  {
    for(i=1;i<=size(keep);i++)
    {
       l[size(l)+1]=keep[i];
       l[size(l)+1]=primaryTest(keep[i],keep[i][1]);
    }
  }
  return(l);
}

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

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 uytrewq;
  int nva = nvars(basering);
  int @k,@s,@n,@k1,zz;
  list primary,lres0,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);
  }

//with 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);
  }
  primary=splitCharp(primary);
//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);
         //hier geaendert +untersuchen!!!!!!!!!!!!!!
       }
       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 lres0;
       list lres0;
       if(size(@lh)==2)
       {
          helpprim=@lh[2];
          lres0[1]=primary[2*@k-1];
          ser1=psi(helpprim);
          lres0[2]=psi1(ser1);
          if(size(reduce(lres0[2],lres0[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);
         //       lres0[2*@n-1]=interred(primary[2*@k-1]+psi1(ser1));
         //       helpprim=@lh[2*@n];
         //      ser1=psi(helpprim);
         //       lres0[2*@n]=psi1(ser1);
         //    }
         // }
         // else
         // {
             lres1=psi(@lh);
             lres0=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,lres0);
       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;
       lres0=imap(@Phelp1,@lr);
       kill @Phelp1;
       for(@n=1;@n<=size(lres0);@n++)
       {
          lres0[@n]=clearSB(lres0[@n]);
          attrib(lres0[@n],"isSB",1);
       }

       primary[2*@k-1]=lres0[1];
       primary[2*@k]=lres0[2];
       @s=size(primary)/2;
       for(@n=1;@n<=size(lres0)/2-1;@n++)
       {
         primary[2*@s+2*@n-1]=lres0[2*@n+1];
         primary[2*@s+2*@n]=lres0[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(lcm(h));
   }
   def bsr= basering;
   string @id=string(h);
   execute("ring @r=0,("+@pa+","+varstr(bsr)+"),(C,dp);");
   execute ("ideal @i="+@id+";");
   poly @p=lcm(@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 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;
   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=indepSet(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=indepSet(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
EXAMPLE: example minAssPrimes; shows an example
"
{
   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;
   list pr= minAssPrimes(i);  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 equidim(ideal i,list #)
"USAGE:  equidim(i) or equidim(i,1) ; i ideal 
         equidim(i,1) uses the algorithm of Eisenbud, Hunecke and Vasconcelos          
 RETURN: list = list of equidimensional ideals a1,...,as such that
         i is the intersection of a1,...,as. as is the equidimensional
         locus of i, ai are the lower dimensional equidimensional loci
         with (perhaps) modified embedded components,i.e. an embedded
         component q (primary ideal) of i can be replaced in the decompo-
         sition by a primary ideal q1 with the same radical as q
 EXAMPLE:example equidim; shows an example
"
{
  def  P = basering;
  list eq;
  intvec w;
  int n,m;
  int g=size(i);
  int a=attrib(i,"isSB");
  int homo=homog(i);
  if(size(#)!=0)
  {
     m=1;
  }

  if(((homo==1)||(a==1))&&(find(ordstr(basering),"l")==0)
                                &&(find(ordstr(basering),"s")==0))
  {
     execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
                              +ordstr(basering)+");");
     ideal i=imap(P,i);
     ideal j=i;
     if(a==1)
     {
       attrib(j,"isSB",1); 
     }
     else
     {
       j=groebner(i);
     }
  }
  else
  {
     execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),dp;");
     ideal i=imap(P,i);
     ideal j=groebner(i);
  }  
  if(homo==1)
  {
     for(n=1;n<=nvars(basering);n++)
     {
        w[n]=ord(var(n));
     }
     intvec hil=hilb(j,1,w);
  }

  if ((dim(j)==-1)||(size(j)==0)||(nvars(basering)==1)
                  ||(dim(j)==0)||(dim(j)+g==nvars(basering)))
  {
    setring P;
    eq[1]=i;
    return(eq);
  }

  if(m==0)
  {
     ideal k=equidimMax(j); 
  }
  else
  {
     ideal k=equidimMaxEHV(j); 
  }
  if(size(reduce(k,j,1))==0)
  {
    setring P;
    eq[1]=i;
    kill gnir;
    return(eq);
  }
  option(returnSB);  
  j=quotient(j,k);
  option(noreturnSB);

  list equi=equidim(j);  
  if(deg(equi[size(equi)][1])<=0)
  {
      equi[size(equi)]=k;
  }
  else
  {
    equi[size(equi)+1]=k;
  }
  setring P;
  eq=imap(gnir,equi);
  kill gnir;
  return(eq);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 32003,(x,y,z),dp;
   ideal i=intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
   equidim(i);
}

///////////////////////////////////////////////////////////////////////////////
proc equidimMax(ideal i)
"USAGE:  equidimMax(i); i ideal           
 RETURN:  ideal = ideal of equidimensional locus
 EXAMPLE: example equidimMax; shows an example
"
{
  def  P = basering;
  ideal eq;
  intvec w;
  int n;
  int g=size(i);
  int a=attrib(i,"isSB");
  int homo=homog(i);
   
  if(((homo==1)||(a==1))&&(find(ordstr(basering),"l")==0)
                                &&(find(ordstr(basering),"s")==0))
  {
     execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
                              +ordstr(basering)+");");
     ideal i=imap(P,i);
     ideal j=i;
     if(a==1)
     {
       attrib(j,"isSB",1); 
     }
     else
     {
       j=groebner(i);
     }
  }
  else
  {
     execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),dp;");
     ideal i=imap(P,i);
     ideal j=groebner(i);
  }
  list indep;
  ideal equ,equi;
  if(homo==1)
  {
     for(n=1;n<=nvars(basering);n++)
     {
        w[n]=ord(var(n));
     }
     intvec hil=hilb(j,1,w);
  }
  if ((dim(j)==-1)||(size(j)==0)||(nvars(basering)==1)
                  ||(dim(j)==0)||(dim(j)+g==nvars(basering)))
  {
    setring P;
    return(i);
  }

  indep=maxIndependSet(j);

  execute("ring gnir1 = ("+charstr(basering)+"),("+indep[1][1]+"),("
                              +indep[1][2]+");");
  if(homo==1)
  {
     ideal j=std(imap(gnir,j),hil,w);
  }
  else
  {
     ideal j=groebner(imap(gnir,j));
  }
  string quotring=prepareQuotientring(nvars(basering)-indep[1][3]);
  execute(quotring);
  ideal j=imap(gnir1,j);
  kill gnir1;
  j=clearSB(j);
  ideal h;
  for(n=1;n<=size(j);n++)
  {
     h[n]=leadcoef(j[n]);
  }
  setring gnir;
  ideal h=imap(quring,h);
  kill quring;

  list l=minSat(j,h);
  
  if(deg(l[2])>0)
  {
    equ=l[1];
    attrib(equ,"isSB",1);
    j=std(j,l[2]);

    if(dim(equ)==dim(j))
    {
      equi=equidimMax(j);
      equ=interred(intersect(equ,equi));
    }
  }
  else
  {
    equ=i;
  }

  setring P;
  eq=imap(gnir,equ);
  kill gnir;
  return(eq);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 32003,(x,y,z),dp;
   ideal i=intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
   equidimMax(i);
}

///////////////////////////////////////////////////////////////////////////////
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,@w;
  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);
     }
     for(@n=1;@n<=nvars(basering);@n++)
     {
        @w[@n]=ord(var(@n));
     }
     intvec @hilb=hilb(tras,1,@w);
     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,@w);
     }
     else
     {
        ideal @j=fetch(@P,tras);
        attrib(@j,"isSB",1);
     }
  }
  else
  {
     ideal @j=groebner(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)==nvars(basering))
    {
       setring @P;
       primary[1]=i;
       primary[2]=i;
       return(primary);
    }
    if(size(fried)>0)
    {
       string newva;
       string newma; 
       for(@k=1;@k<=nvars(basering);@k++)
       {
          @n1=0;
          for(@n=1;@n<=size(fried);@n++)
          {
             if(leadmonom(fried[@n])==var(@k))
             {
                @n1=1;
                break;
             }
          } 
          if(@n1==0)
          {
            newva=newva+string(var(@k))+",";
            newma=newma+string(var(@k))+",";
          }
          else
          {
            newma=newma+string(0)+",";   
          }      
       }
       newva[size(newva)]=")";
       newma[size(newma)]=";";
       execute("ring @deirf=("+charstr(gnir)+"),("+newva+",lp;");
       execute("map @kappa=gnir,"+newma);
       ideal @j= @kappa(@j);
       @j=simplify(@j,2);
       attrib(@j,"isSB",1);
       list pr=decomp(@j);
       setring gnir;
       list pr=imap(@deirf,pr);
       for(@k=1;@k<=size(pr);@k++)
       {
         @j=pr[@k]+fried;
         pr[@k]=@j;
       }
       setring @P;
       return(imap(gnir,pr));
    }
  }

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

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

  //----------------------------------------------------------------
  //  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 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,@w);
    }

  }
  else
  {
    ideal jwork=groebner(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,@w);
        }
        else
        {
          ideal @j=groebner(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)

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

     //arrangement for quotientring K(var(nnp+1),..,var(nva))[..the rest..] and
     //map phi:K[var(1),...,var(nva)] --->K(var(nnpr+1),..,var(nva))[..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
     //in the quotientring: this is coded in saturn

     list saturn;
     ideal hpl;

     for(@n=1;@n<=size(uprimary);@n++)
     {
        uprimary[@n]=interred(uprimary[@n]); // temporary fix
        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,@w);
        }

        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))
           //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,@w);
              }
              else
              {
                ideal @j=groebner(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))[..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))[..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 toh 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))[..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;
   list pr= decomp(i);
   pr;
   testPrimary( pr, i);
}

///////////////////////////////////////////////////////////////////////////////
proc radicalKL (list m,ideal ser,list #)
{
   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((find(ordstr(basering),"w")!=0)||(find(ordstr(basering),"W")!=0)||(find(ordstr(basering),"a")!=0))
   {
      homo=0;
   }
   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)
   {
      if((npars(basering)>0)&&(char(basering)>0)&&(char(basering)<=181))
      {
         "        !  Warning  !  ";
         "    The field is not perfect.";
         "    The result may be wrong.";
      }
      @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, 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=groebner(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))[..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))[..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))[..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);
}

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

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)
       &&((char(basering)==0)||(char(basering)>181)))
   {
      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 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=nres(f,2)[2];            //the kernel
     matrix g=transpose(re[n]);        //the image of Hom(_,R)

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

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;
         }
      }
   }
}

///////////////////////////////////////////////////////////////////////////////
//
//                  Shimoyama-Yokoyama
//
///////////////////////////////////////////////////////////////////////////////

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);
  if(SI[1]==1)  // primdecSY(ideal(1))
  {
    return(list());
  }
  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=indepSet(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);
}


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

proc convList(list l)
{
   int i;
   list re,he;
   for(i=1;i<=size(l)/2;i++)
   {
      he=l[2*i-1],l[2*i];
      re[i]=he;
   }
   return(re);
}
///////////////////////////////////////////////////////////////////////////////

proc reconvList(list l)
{
   int i;
   list re;
   for(i=1;i<=size(l);i++)
   {
      re[2*i-1]=l[i][1];
      re[2*i]=l[i][2];
   }
   return(re);
}

///////////////////////////////////////////////////////////////////////////////
//
//     The main procedures
//
///////////////////////////////////////////////////////////////////////////////

proc primdecGTZ(ideal i)
"USAGE:   primdecGTZ(i); i ideal
RETURN:  a list, say pr, of primary ideals and their associated primes
         pr[i][1], resp. pr[i][2] is the i-th primary resp. prime component 
NOTE:    Algorithm of Gianni, Traeger, Zacharias 
         designed for characteristic 0, works also in char k > 0, if it 
         terminates (may result in an infinite loop in small characteristic!)
EXAMPLE: example primdecGTZ; shows an example
"
{
   return(convList(decomp(i)));
}
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;
   list pr = primdecGTZ(i);
   pr;
}
///////////////////////////////////////////////////////////////////////////////

proc primdecSY(ideal i, list #))
"USAGE:   primdecSY(i); i ideal, c int
         if c=0, the given ordering of the variables is used.
         if c=1, minAssChar tries to use an optimal ordering,
         if c=2, minAssGTZ is used
         if c=3, minAssGTZ and facstd is used
RETURN:  a list, say pr, of primary ideals and their associated primes
         pr[i][1], resp. pr[i][2] is the i-th primary resp. prime component 
NOTE:    Algorithm of Shimoyama-Yokoyama
         implemented for characteristic 0, works also in char k > 0,
         the result may be not completely decomposed in small characteristic
EXAMPLE: example primdecSY; shows an example
"
{
   if (size(#)==1)
   { return(prim_dec(i,#[1])); }
   else
   { return(prim_dec(i,1)); }
}
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;
   list pr = primdecSY(i);
   pr;
}
///////////////////////////////////////////////////////////////////////////////
proc minAssGTZ(ideal i)
"USAGE:   minAssGTZ(i); i ideal
RETURN:  list = the minimal associated prime ideals of i
NOTE:    designed for characteristic 0, works also in char k > 0 if it
         terminates, may result in an infinite loop in small characteristic
EXAMPLE: example minAssGTZ; shows an example
"
{
   return(minAssPrimes(i,1));
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 32003,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   list pr= minAssGTZ(i);
   pr;
}

///////////////////////////////////////////////////////////////////////////////
proc minAssChar(ideal i, list #)
"USAGE:   minAssChar(i[,c]); i ideal,
         if c=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
RETURN:  list = the minimal associated prime ideals of i
NOTE:    implemented for characteristic 0, works also in char k >> 0,
         the result may be not compltely decomposed in small characteristic
EXAMPLE: example minAssChar; shows an example
"
{
  if (size(#)==1)
  { return(min_ass_prim_charsets(i,#[1])); }
  else
  { return(min_ass_prim_charsets(i,1)); }
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 32003,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   list pr= minAssChar(i);
   pr;
}
///////////////////////////////////////////////////////////////////////////////
proc equiRadical(ideal i)
"USAGE:   equiRadical(i); i ideal
RETURN:  ideal, intersection of associated primes of i of maximal  dimension 
NOTE:    designed for characteristic 0, works also in char k > 0 if it 
         terminates, may result in an infinite loop in small characteristic
EXAMPLE: example equiRadical; shows an example
"
{
   return(radical(i,1));
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 32003,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   ideal pr= equiRadical(i);
   pr;
}
///////////////////////////////////////////////////////////////////////////////
proc radical(ideal i,list #)
"USAGE:   radical(i); i ideal
RETURN:  ideal = the radical of i
NOTE:    a combination of the algorithms of Krick/Logar and 
         Eisenbud/Huneke/Vasconcelos
         designed for characteristic 0, works also in char k > 0 if it 
         terminates, may result in an infinite loop in small characteristic
EXAMPLE: example radical; shows an example
"
{
   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(pr[1],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[s+1-j]==0)
      {
         re=intersect(re,radicalEHV(pr[s+1-j],re,il));
      }
   }
   re=interred(re);
   return(phi(re));
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 32003,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   ideal pr= radical(i);
   pr;
}
///////////////////////////////////////////////////////////////////////////////
proc prepareAss(ideal i)
"USAGE:   prepareAss(i); i ideal
RETURN:  list = the radicals of the maximal dimensional components of i
NOTE:    uses algorithm of Eisenbud, Huneke and Vasconcelos
EXAMPLE: example prepareAss; shows an example
"
{
  ideal j=std(i);
  int cod=nvars(basering)-dim(j);
  int e;
  list er;
  ideal ann;
  if(homog(i)==1)
  {
     list re=sres(j,0);                   //the resolution
     re=minres(re);                       //minimized resolution
  }
  else
  {
    list re=mres(i,0);
  }
  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);
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 32003,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   list pr= prepareAss(i);
   pr;
}
///////////////////////////////////////////////////////////////////////////////
proc equidimMaxEHV(ideal i)
"USAGE:  equidimMaxEHV(i); i ideal
RETURN:  ideal = equidimensional componente of i
NOTE:    uses algorithm of Eisenbud, Huneke and Vasconcelos
EXAMPLE: example equidimMaxEHV; shows an example
"
{
  ideal j=groebner(i);
  int cod=nvars(basering)-dim(j);
  int e;
  ideal ann;
  if(homog(i)==1)
  {
     list re=sres(j,0);                   //the resolution
     re=minres(re);                       //minimized resolution
  }
  else
  {
    list re=mres(i,0);
  }
  ann=AnnExt_R(cod,re);
  return(ann);
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 32003,(x,y,z),dp;
   ideal i=intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
   equidimMaxEHV(i);
}

proc testPrimary(list pr, ideal k)
"USAGE:   testPrimary(pr,k); pr a list, result of primdecGTZ(k) or primdecSY(k)
RETURN:  int, 1 if intersection of the primary ideals in pr is k, 0 if not
EXAMPLE: example testPrimary; shows an example
"
{
   int i;
   pr=reconvList(pr);
   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  r = 32003,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   list pr = primdecGTZ(i);
   testPrimary(pr,i);
}

///////////////////////////////////////////////////////////////////////////////
proc zerodec(ideal I)
"USAGE:   zerodec(I); I ideal
RETURN:  a list of primary ideals, the zero-dimensional decomposition of I
ASSUME:  I is zero-dimensional, the characterisitic of the ground field is 0
NOTE:    the algorithm (of C. Monico), works well only for small total
         number of solutions (i.e. vdim(std(I)) should be < 100) 
         and without parameters. In practice, it works also in big 
         characteristic p>0 but may fail for small p.
         If printlevel > 0 (default = 0) additional information is displayed
EXAMPLE: example zerodec; shows an example
"
{
   def R=basering;
   poly q;
   int j,time;
   matrix m;
   list re;
   poly va=var(1);
   ideal J=groebner(I);
   ideal ba=kbase(J);
   int d=vdim(J);
   dbprint(printlevel-voice+2,"// multiplicity of ideal : "+ string(d));
//------ compute matrix of multiplication on R/I with generic element p -----
   int e=nvars(basering);
   poly p=randomLast(100)[e]+random(-50,50);     //the generic element
   matrix n[d][d];
   time = timer;
   for(j=2;j<=e;j++)
   {
      va=va*var(j);
   }
   for(j=1;j<=d;j++)
   {
      q=reduce(p*ba[j],J);
      m=coeffs(q,ba,va);
      n[j,1..d]=m[1..d,1];
   }
   dbprint(printlevel-voice+2,
      "// time for computing multiplication matrix (with generic element) : "+ 
      string(timer-time));
//---------------- compute characteristic polynomial of matrix --------------
   execute("ring P1=("+charstr(R)+"),T,dp;");
   matrix n=imap(R,n);
   time = timer;
   poly charpol=det(n-T*freemodule(d));
   dbprint(printlevel-voice+2,"// time for computing char poly: "+ 
           string(timer-time));
//------------------- factorize characteristic polynomial -------------------
//check first if constant term of charpoly is != 0 (which is true for 
//sufficiently generic element)
   if(charpol[size(charpol)]!=0)
   {
     time = timer;
     list fac=factor(charpol);
     testFactor(fac,charpol);
     dbprint(printlevel-voice+2,"// time for factorizing char poly: "+ 
             string(timer-time));
     int f=size(fac[1]);
//--------------------------- the irreducible case --------------------------
     if(f==1)
     {
       setring R;
       re=I;
       return(re);
     }
//---------------------------- the reducible case ---------------------------
//if f_i are the irreducible factors of charpoly, mult=ri, then <I,g_i^ri> 
//are the primary components where g_i = f_i(p). However, substituting p in 
//f_i may result in a huge object although the final result may be small. 
//Hence it is better to simultaneously reduce with I. For this we need a new 
//ring.
     execute("ring P=("+charstr(R)+"),(T,"+varstr(R)+"),(dp(1),dp);");
     list rfac=imap(P1,fac);
     intvec ov=option(get);;
     option(redSB);
     list re1;
     ideal new = T-imap(R,p),imap(R,J);
     attrib(new, "isSB",1);    //we know that new is a standard basis
     for(j=1;j<=f;j++)
     {
        re1[j]=reduce(rfac[1][j]^rfac[2][j],new); 
     }
     setring R;
     re = imap(P,re1);
     for(j=1;j<=f;j++)
     {
        J=I,re[j];
        re[j]=interred(J); 
     }
     option(set,ov);
     return(re);
  }
  else
//------------------- choice of generic element failed -------------------
  {
     dbprint(printlevel-voice+2,"// try new generic element!");
     setring R;
     return(zerodec(I));
  }
}
example
{ "EXAMPLE:";  echo = 2;
   ring r=0,(x,y),dp;
   ideal i=x2-2,y2-2;
   list pr=zerodec(i);
   pr;
}
////////////////////////////////////////////////////////////////////////////
/*
//Beispiele Wenk-Dipl (in ~/Texfiles/Diplom/Wenk/Examples/)
//Zeiten: Singular/Singular/Singular -r123456789 -v :wilde13 (PentiumPro200)
//Singular for HPUX-9 version 1-3-8  (2000060214)  Jun  2 2000 15:31:26
//(wilde13)

//1. vdim=20, 3  Komponenten
//zerodec-time:2(1)  (matrix:1 charpoly:0 factor:1)
//primdecGTZ-time: 1(0)
ring rs= 0,(a,b,c),dp;
poly f1= a^2*b*c + a*b^2*c + a*b*c^2 + a*b*c + a*b + a*c + b*c;
poly f2= a^2*b^2*c + a*b^2*c^2 + a^2*b*c + a*b*c + b*c + a + c;
poly f3= a^2*b^2*c^2 + a^2*b^2*c + a*b^2*c + a*b*c + a*c + c + 1;
ideal gls=f1,f2,f3;
int time=timer;
printlevel =1;
time=timer; list pr1=zerodec(gls); timer-time;size(pr1);
time=timer; list pr =primdecGTZ(gls); timer-time;size(pr);

//2.cyclic5  vdim=70, 20 Komponenten
//zerodec-time:36(28)  (matrix:1(0) charpoly:18(19) factor:17(9)
//primdecGTZ-time: 28(5)
ring rs= 0,(a,b,c,d,e),dp;
poly f0= a + b + c + d + e + 1;
poly f1= a + b + c + d + e;
poly f2= a*b + b*c + c*d + a*e + d*e;
poly f3= a*b*c + b*c*d + a*b*e + a*d*e + c*d*e;
poly f4= a*b*c*d + a*b*c*e + a*b*d*e + a*c*d*e + b*c*d*e;
poly f5= a*b*c*d*e - 1;
ideal gls= f1,f2,f3,f4,f5;

//3. random vdim=40, 1 Komponente
//zerodec-time:126(304)  (matrix:1 charpoly:115(298) factor:10(5))
//primdecGTZ-time:17 (11) 
ring rs=0,(x,y,z),dp;
poly f1=2*x^2 + 4*x + 3*y^2 + 7*x*z + 9*y*z + 5*z^2;
poly f2=7*x^3 + 8*x*y + 12*y^2 + 18*x*z + 3*y^4*z + 10*z^3 + 12;
poly f3=3*x^4 + 1*x*y*z + 6*y^3 + 3*x*z^2 + 2*y*z^2 + 4*z^2 + 5;
ideal gls=f1,f2,f3;

//4. introduction into resultants, sturmfels, vdim=28, 1 Komponente
//zerodec-time:4  (matrix:0 charpoly:0 factor:4)
//primdecGTZ-time:1 
ring rs=0,(x,y),dp;
poly f1= x4+y4-1;
poly f2= x5y2-4x3y3+x2y5-1;
ideal gls=f1,f2;

//5. 3 quadratic equations with random coeffs, vdim=8, 1 Komponente
//zerodec-time:0(0)  (matrix:0 charpoly:0 factor:0)
//primdecGTZ-time:1(0) 
ring rs=0,(x,y,z),dp;
poly f1=2*x^2 + 4*x*y + 3*y^2 + 7*x*z + 9*y*z + 5*z^2 + 2;
poly f2=7*x^2 + 8*x*y + 12*y^2 + 18*x*z + 3*y*z + 10*z^2 + 12;
poly f3=3*x^2 + 1*x*y + 6*y^2 + 3*x*z + 2*y*z + 4*z^2 + 5;
ideal gls=f1,f2,f3;

//6. 3 polys    vdim=24, 1 Komponente
// run("ex14",2);
//zerodec-time:5(4)  (matrix:0 charpoly:3(3) factor:2(1))
//primdecGTZ-time:4 (2)
ring rs=0,(x1,x2,x3,x4),dp;
poly f1=16*x1^2 + 3*x2^2 + 5*x3^4 - 1 - 4*x4 + x4^3;
poly f2=5*x1^3 + 3*x2^2 + 4*x3^2*x4 + 2*x1*x4 - 1 + x4 + 4*x1 + x2 + x3 + x4;
poly f3=-4*x1^2 + x2^2 + x3^2 - 3 + x4^2 + 4*x1 + x2 + x3 + x4;
poly f4=-4*x1 + x2 + x3 + x4;
ideal gls=f1,f2,f3,f4;

//7. ex43, PoSSo, caprasse   vdim=56, 16 Komponenten
//zerodec-time:23(15)  (matrix:0 charpoly:16(13) factor:3(2)) 
//primdecGTZ-time:3 (2)
ring rs= 0,(y,z,x,t),dp;
ideal gls=y^2*z+2*y*x*t-z-2*x,
4*y^2*z*x-z*x^3+2*y^3*t+4*y*x^2*t-10*y^2+4*z*x+4*x^2-10*y*t+2,
2*y*z*t+x*t^2-2*z-x,
-z^3*x+4*y*z^2*t+4*z*x*t^2+2*y*t^3+4*z^2+4*z*x-10*y*t-10*t^2+2;

//8. Arnborg-System, n=6 (II),    vdim=156, 90 Komponenten
//zerodec-time (char32003):127(45)(matrix:2(0) charpoly:106(37) factor:16(7))
//primdecGTZ-time(char32003) :81 (18)
//ring rs= 0,(a,b,c,d,x,f),dp;
ring rs= 32003,(a,b,c,d,x,f),dp;
ideal gls=a+b+c+d+x+f, ab+bc+cd+dx+xf+af, abc+bcd+cdx+d*xf+axf+abf, 
abcd+bcdx+cd*xf+ad*xf+abxf+abcf, abcdx+bcd*xf+acd*xf+abd*xf+abcxf+abcdf, 
abcd*xf-1;

//9. ex42, PoSSo, Methan6_1, vdim=27, 2 Komponenten
//zerodec-time:610  (matrix:10 charpoly:557 factor:26)
//primdecGTZ-time: 118
//zerodec-time(char32003):2
//primdecGTZ-time(char32003):4
//ring rs= 0,(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10),dp;
ring rs= 32003,(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10),dp;
ideal gls=64*x2*x7-10*x1*x8+10*x7*x9+11*x7*x10-320000*x1,
-32*x2*x7-5*x2*x8-5*x2*x10+160000*x1-5000*x2,
-x3*x8+x6*x8+x9*x10+210*x6+1300000,
-x4*x8+700000,
x10^2-2*x5,
-x6*x8+x7*x9-210*x6,
-64*x2*x7-10*x7*x9-11*x7*x10+320000*x1-16*x7+7000000,
-10*x1*x8-10*x2*x8-10*x3*x8-10*x4*x8-10*x6*x8+10*x2*x10+11*x7*x10
    +20000*x2+14*x5,
x4*x8-x7*x9-x9*x10-410*x9,
10*x2*x8+10*x3*x8+10*x6*x8+10*x7*x9-10*x2*x10-11*x7*x10-10*x9*x10
    -10*x10^2+1400*x6-4200*x10;

//10. ex33, walk-s7, Diplomarbeit von Tim, vdim=114
//zerfaellt in unterschiedlich viele Komponenten in versch. Charkteristiken:
//char32003:30, char0:3(2xdeg1,1xdeg112!), char181:4(2xdeg1,1xdeg28,1xdeg84)
//char 0: zerodec-time:10075 (ca 3h) (matrix:3 charpoly:9367, factor:680
//        + 24 sec fuer Normalform (anstatt einsetzen), total [29623k])
//        primdecGTZ-time: 214
//char 32003:zerodec-time:197(68) (matrix:2(1) charpoly:173(60) factor:15(6))
//        primdecGTZ-time:14 (5)
//char 181:zerodec-time:(87) (matrix:(1) charpoly:(58) factor:(25))
//        primdecGTZ-time:(2)
//in char181 stimmen Ergebnisse von zerodec und primdecGTZ ueberein (gecheckt)

//ring rs= 0,(a,b,c,d,e,f,g),dp;
ring rs= 32003,(a,b,c,d,e,f,g),dp;
poly f1= 2gb + 2fc + 2ed + a2 + a;
poly f2= 2gc + 2fd + e2 + 2ba + b;
poly f3= 2gd + 2fe + 2ca + c + b2;
poly f4= 2ge + f2 + 2da + d + 2cb;
poly f5= 2fg + 2ea + e + 2db + c2;
poly f6= g2 + 2fa + f + 2eb + 2dc;
poly f7= 2ga + g + 2fb + 2ec + d2;
ideal gls= f1,f2,f3,f4,f5,f6,f7;

~/Singular/Singular/Singular -r123456789 -v
LIB"./primdec.lib";
timer=1;
int time=timer;
printlevel =1;
option(prot,mem);
time=timer; list pr1=zerodec(gls); timer-time;

time=timer; list pr =primdecGTZ(gls); timer-time;
time=timer; list pr =primdecSY(gls); timer-time;
*/