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

///////////////////////////////////////////////////////////////////////////////
version="$Id: primdec.lib,v 1.102 2004-02-12 09:54:06 Singular Exp $";
category="Commutative Algebra";
info="
LIBRARY: primdec.lib   Primary Decomposition and Radical of Ideals
AUTHORS:  Gerhard Pfister, pfister@mathematik.uni-kl.de (GTZ)
@*        Wolfram Decker, decker@math.uni-sb.de         (SY)
@*        Hans Schoenemann, hannes@mathematik.uni-kl.de (SY)

OVERVIEW:
    Algorithms for primary decomposition based on the ideas of
    Gianni, Trager and Zacharias (implementation by Gerhard Pfister),
    respectively based on the ideas of Shimoyama and Yokoyama (implementation
    by Wolfram Decker and Hans Schoenemann).
 @* The procedures are implemented to be used in characteristic 0.
 @* They also work in positive characteristic >> 0.
 @* In small characteristic and for algebraic extensions, primdecGTZ
    may not terminate.
    Algorithms for the computation of the radical based on the ideas of
    Krick, Logar and Kemper (implementation by Gerhard Pfister).

PROCEDURES:
 Ann(M);           annihilator of R^n/M, R=basering, M in R^n
 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 I via  Krick/Logar and Kemper
 radicalEHV(I);    computes the radical of I via  Eisenbud,Huneke,Vasconcelos
 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";
LIB "triang.lib";
///////////////////////////////////////////////////////////////////////////////
//
//                      Gianni/Trager/Zacharias
//
///////////////////////////////////////////////////////////////////////////////

static 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 sat1; shows an example
"
{
   int @k;
   ideal inew=std(id);
   ideal iold;
   intvec op=option(get);
   option(returnSB);
   while(specialIdealsEqual(iold,inew)==0 )
   {
      iold=inew;
      inew=quotient(iold,p);
      @k++;
   }
   @k--;
   option(set,op);
   list L =inew,p^@k;
   return (L);
}

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

static 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];
     }
     intvec op = option(get);
     option(returnSB);
     while(specialIdealsEqual(iold,inew)==0 )
     {
         iold=inew;
         if(deg(iold[size(iold)])!=1)
         {
            inew=quotient(iold,@f);
         }
         else
         {
            inew=iold;
         }
         @k++;
      }
      option(set,op);
      @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;
   intvec op = option(get);
   option(returnSB);
   while(specialIdealsEqual(iold,quotM[1])==0)
   {
      if(k>0)
      {
         f=f*quotM[3];
      }
      iold=quotM[1];
      quotM=quotMin(quotM);
      k++;
   }
   option(set,op);
   l=quotM[1],f;
   return(l);
}

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

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

static 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, please inform the authors";
   }
}
///////////////////////////////////////////////////////////////////////////////

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

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

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

///////////////////////////////////////////////////////////////////////////////
static 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;
}
///////////////////////////////////////////////////////////////////////////////
static proc splitCharp(list l)
{
  if((char(basering)==0)||(npars(basering)>0))
  {
     return(l);
  }
  def P=basering;
  int i,j,k,m,q,d,o;
  int n=nvars(basering);
  ideal s,t,u,sact;
  poly ni;
  string minp,gnir,va;
  list sa,keep,rp,keep1;
  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);   //vermeiden!!!
      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+"),lp;";
        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);

        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(sact)==1)&&(m==2))
        {
           l[2*i]=l[2*i-1];
           attrib(l[2*i],"isSB",1);
        }
        if((size(sact)==1)&&(m>2))
        {
           setring RL;
           option(redSB);
           t=std(t);

           list sp=zero_decomp(t,0,0);

           setring P;
           rp=imap(RL,sp);
           for(o=1;o<=size(rp);o++)
           {
              rp[o]=interred(simplify(rp[o],1)+ideal(ni));
           }
           l[2*i-1]=rp[1];
           l[2*i]=rp[2];
           rp=delete(rp,1);
           rp=delete(rp,1);
           keep1=keep1+rp;
           option(noredSB);
        }
        kill RL;
      }
    }
  }
  if(size(keep)>0)
  {
    for(i=1;i<=size(keep);i++)
    {
       if(deg(keep[i][1])>0)
       {
          l[size(l)+1]=keep[i];
          l[size(l)+1]=primaryTest(keep[i],keep[i][1]);
       }
    }
  }
  l=l+keep1;
  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/Trager/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);
  }

  if((voice>=6)&&(char(basering)<=181))
  {
     primary=splitCharp(primary);
  }

  if((@wr==2)&&(npars(basering)>0)&&(voice>=6)&&(char(basering)>0))
  {
  //the prime decomposition of Yokoyama in characteristic p
     list ke,ek;
     @k=0;
     while(@k<size(primary)/2)
     {
        @k++;
        if(size(primary[2*@k])==0)
        {
           ek=insepDecomp(primary[2*@k-1]);
           primary=delete(primary,2*@k);
           primary=delete(primary,2*@k-1);
           @k--;
        }
        ke=ke+ek;
     }
     for(@k=1;@k<=size(ke);@k++)
     {
        primary[size(primary)+1]=ke[@k];
        primary[size(primary)+1]=ke[@k];
     }
  }

  if(voice>=8){primary=extF(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.
  if((npars(basering)>0)&&(voice>=8))
  {
     poly randp;
     for(zz=1;zz<nvars(basering);zz++)
     {
        randp=randp
              +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(zz);
     }
     randp=randp+var(nvars(basering));
  }
  @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++)
       {
          attrib(primary[2*@k-1],"isSB",1);
          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];
       if((npars(basering)>0)&&(voice>=10))
       {
          jmap[size(jmap)]=randp;
       }

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

//=================== the new part ============================

       @qh=groebner(@qh);

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


       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
       {

          lres1=psi(@lh);
          lres0=psi1(lres1);
       }

//=================== the new part ============================

       primary=delete(primary,2*@k-1);
       primary=delete(primary,2*@k-1);
       @k--;
       if(size(lres0)==2)
       {
          lres0[2]=groebner(lres0[2]);
       }
       else
       {
          for(@n=1;@n<=size(lres0)/2;@n++)
          {
             if(specialIdealsEqual(lres0[2*@n-1],lres0[2*@n])==1)
             {
                lres0[2*@n-1]=groebner(lres0[2*@n-1]);
                lres0[2*@n]=lres0[2*@n-1];
                attrib(lres0[2*@n],"isSB",1);
             }
             else
             {
                lres0[2*@n-1]=groebner(lres0[2*@n-1]);
                lres0[2*@n]=groebner(lres0[2*@n]);
             }
          }
       }
       primary=primary+lres0;

//=============================================================
//       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 extF(list l,list #)
{
//zero_dimensional primary decomposition after finite field extension
   def R=basering;
   int p=char(R);

   if((p==0)||(p>13)||(npars(R)>0)){return(l);}

   int ex=3;
   if(size(#)>0){ex=#[1];}

   list peek,peek1;
   while(size(l)>0)
   {
      if(size(l[2])==0)
      {
         peek[size(peek)+1]=l[1];
      }
      else
      {
         peek1[size(peek1)+1]=l[1];
         peek1[size(peek1)+1]=l[2];
      }
      l=delete(l,1);
      l=delete(l,1);
   }
   if(size(peek)==0){return(peek1);}

   string gnir="ring RH=("+string(p)+"^"+string(ex)+",a),("+varstr(R)+"),lp;";
   execute(gnir);
   string mp="minpoly="+string(minpoly)+";";
   gnir="ring RL=("+string(p)+",a),("+varstr(R)+"),lp;";
   execute(gnir);
   execute(mp);
   list L=imap(R,peek);
   list pr, keep;
   int i;
   for(i=1;i<=size(L);i++)
   {
      attrib(L[i],"isSB",1);
      pr=zero_decomp(L[i],0,0);
      keep=keep+pr;
   }
   for(i=1;i<=size(keep);i++)
   {
      keep[i]=simplify(keep[i],1);
   }
   mp="poly pp="+string(minpoly)+";";

   string gnir1="ring RS="+string(p)+",("+varstr(R)+",a),lp;";
   execute(gnir1);
   execute(mp);
   list L=imap(RL,keep);

   for(i=1;i<=size(L);i++)
   {
      L[i]=eliminate(L[i]+ideal(pp),a);
   }
   i=0;
   int j;
   while(i<size(L)/2-1)
   {
      i++;
      j=i;
      while(j<size(L)/2)
      {
         j++;
         if(idealsEqual(L[2*i-1],L[2*j-1]))
         {
            L=delete(L,2*j-1);
            L=delete(L,2*j-1);
            j--;
         }
      }
   }
   setring R;
   list re=imap(RS,L);
   re=re+peek1;

   return(extF(re,ex+1));
}

///////////////////////////////////////////////////////////////////////////////
proc zeroSp(ideal i)
{
//preparation for the separable closure
//decomposition into ideals of special type
//i.e. the minimal polynomials of every variable mod i are irreducible
//returns a list of 2 lists: rr=pe,qe
//the ideals in pe[l] are special, their special elements are in qe[l]
//pe[l] is a dp-Groebnerbasis
//the radical of the intersection of the pe[l] is equal to the radical of i

   def R=basering;

   //i has to be a reduced groebner basis
   ideal F=finduni(i);

   int j,k,l,ready;
   list fa;
   fa[1]=factorize(F[1],1);
   poly te,ti;
   ideal tj;
   //avoid factorization of the same polynomial
   for(j=2;j<=size(F);j++)
   {
      for(k=1;k<=j-1;k++)
      {
         ti=F[k];
         te=subst(ti,var(k),var(j));
         if(te==F[j])
         {
            tj=fa[k];
            fa[j]=subst(tj,var(k),var(j));
            ready=1;
            break;
         }
      }
      if(!ready)
      {
         fa[j]=factorize(F[j],1);
      }
      ready=0;
   }
   execute( "ring P=("+charstr(R)+"),("+varstr(R)+"),(C,dp);");
   ideal i=imap(R,i);
   if(npars(basering)==0)
   {
      ideal J=fglm(R,i);
   }
   else
   {
      ideal J=groebner(i);
   }
   list fa=imap(R,fa);
   list qe=J;          //collects a dp-Groebnerbasis of the special ideals
   list keep=ideal(0); //collects the special elements

   list re,em,ke;
   ideal K,L;

   for(j=1;j<=nvars(basering);j++)
   {
      for(l=1;l<=size(qe);l++)
      {
         for(k=1;k<=size(fa[j]);k++)
         {
            L=std(qe[l],fa[j][k]);
            K=keep[l],fa[j][k];
            if(deg(L[1])>0)
            {
               re[size(re)+1]=L;
               ke[size(ke)+1]=K;
            }
         }
      }
      qe=re;
      re=em;
      keep=ke;
      ke=em;
   }

   setring R;
   list qe=imap(P,keep);
   list pe=imap(P,qe);
   for(l=1;l<=size(qe);l++)
   {
      qe[l]=simplify(qe[l],2);
    }
    list rr=pe,qe;
    return(rr);
}
///////////////////////////////////////////////////////////////////////////////

proc zeroSepClos(ideal I,ideal F)
{
//computes the separable closure of the special ideal I
//F is the set of special elements of I
//returns the separable closure sc(I) of I and an intvec v
//such that sc(I)=preimage(frobenius definde by v)
//i.e. var(i)----->var(i)^(p^v[i])

   if(homog(I)==1){return(maxideal(1));}

   //assume F[i] irreducible in I and depending only on var(i)

   def R=basering;
   int n=nvars(R);
   int p=char(R);
   intvec v;
   v[n]=0;
   int i,k;
   list l;

   for(i=1;i<=n;i++)
   {
      l[i]=sep(F[i],i);
      F[i]=l[i][1];
      if(l[i][2]>k){k=l[i][2];}
   }

   if(k==0){return(list(I,v));}        //the separable case
   ideal m;

   for(i=1;i<=n;i++)
   {
      m[i]=var(i)^(p^l[i][2]);
      v[i]=l[i][2];
   }
   map phi=R,m;
   ideal J=preimage(R,phi,I);
   return(list(J,v));
}
///////////////////////////////////////////////////////////////////////////////

proc insepDecomp(ideal i)
{
//decomposes i into special ideals
//computes the prime decomposition of the special ideals
//and transforms it back to a decomposition of i

   def R=basering;
   list pr=zeroSp(i);
   int l,k;
   list re,wo,qr;
   ideal m=maxideal(1);
   ideal K;
   map phi=R,m;
   int p=char(R);
   intvec op=option(get);

   for(l=1;l<=size(pr[1]);l++)
   {
      wo=zeroSepClos(pr[1][l],pr[2][l]);
      for(k=1;k<=nvars(basering);k++)
      {
         m[k]=var(k)^(p^wo[2][k]);
      }
      phi=R,m;
      qr=decomp(wo[1],2);

      option(redSB);
      for(k=1;k<=size(qr)/2;k++)
      {
         K=qr[2*k];
         K=phi(K);
         K=groebner(K);
         re[size(re)+1]=zeroRad(K);
      }
      option(noredSB);
   }
   option(set,op);
   return(re);
}


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

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

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

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


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

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


}

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

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

}

///////////////////////////////////////////////////////////////////////////////
static 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 minAssPrimesold(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);
}

static proc primT(ideal i)
{
   //assumes that all generators of i are irreducible
   //i is standard basis

   attrib(i,"isSB",1);
   int j=size(i);
   int k;
   while(j>0)
   {
     if(deg(i[j])>1){break;}
     j--;
   }
   if(j==0){return(1);}
   if(deg(i[j])==vdim(i)){return(1);}
   return(0);
}


static 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 q=simplifyIdeal(i);
   list re=maxideal(1);
   int j,a,k;
   intvec op=option(get);
   map phi=P,q[2];

   if(npars(P)==0){option(redSB);}

   i=std(q[1]);
   if(dim(i)==-1){re=ideal(1);return(re);}
   if((dim(i)==0)&&(npars(P)==0))
   {
      int di=vdim(i);
      execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;");
      ideal J=interred(imap(P,i));
      attrib(J,"isSB",1);
      if(vdim(J)!=di)
      {
         J=fglm(P,i);
      }
      list pr=triangMH(J,2);
      list qr,re;

      for(k=1;k<=size(pr);k++)
      {
         if(primT(pr[k]))
         {
            re[size(re)+1]=pr[k];
         }
         else
         {
            attrib(pr[k],"isSB",1);
            qr=decomp(pr[k],2);
            for(j=1;j<=size(qr)/2;j++)
            {
               re[size(re)+1]=qr[2*j];
            }
         }
      }
      setring P;
      re=imap(gnir,re);
      option(set,op);
      return(phi(re));
   }

   if((size(#)==0)||(dim(i)==0))
   {
      re[1]=decomp(i,2);
      re=union(re);
      option(set,op);
      return(phi(re));
   }

   q=facstd(i);

   if((size(q)==1)&&(dim(i)>1))
   {
      execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;");

      list p=facstd(fetch(P,i));
      if(size(p)>1)
      {
         a=1;
         setring P;
         q=fetch(gnir,p);
      }
      else
      {
         setring P;
      }
      kill gnir;
   }

   option(set,op);
   for(j=1;j<=size(q);j++)
   {
      if(a==0){attrib(q[j],"isSB",1);}
      re[j]=decomp(q[j],2);
   }
   re=union(re);
   return(phi(re));
}
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);
}

static 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 equidim(ideal i,list #)
"USAGE:  equidim(i) or equidim(i,1) ; i ideal
RETURN: list of equidimensional ideals a[1],...,a[s] with:
        - a[s] the equidimensional locus of i, i.e. the intersection
          of the primary ideals of dimension of i
        - a[1],...,a[s-1] the lower dimensional equidimensional loci.
NOTE:    An embedded component q (primary ideal) of i can be replaced in the
         decomposition by a primary ideal q1 with the same radical as q. @*
         @code{equidim(i,1)} uses the algorithm of Eisenbud/Huneke/Vasconcelos.

EXAMPLE:example equidim; shows an example
"
{
  if(ord_test(basering)!=1)
  {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
  }
  intvec op ;
  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);
  }
  op=option(get);
  option(returnSB);
  j=quotient(j,k);
  option(set,op);

  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 of equidimensional locus (of maximal dimension) of i.
EXAMPLE: example equidimMax; shows an example
"
{
  if(ord_test(basering)!=1)
  {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
  }
  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);
}
///////////////////////////////////////////////////////////////////////////////
static proc islp()
{
   string s=ordstr(basering);
   int n=find(s,"lp");
   if(!n){return(0);}
   int k=find(s,",");
   string t=s[k+1..size(s)];
   int l=find(t,",");
   t=s[1..k-1];
   int m=find(t,",");
   if(l+m){return(0);}
   return(1);
}
///////////////////////////////////////////////////////////////////////////////

proc algeDeco(ideal i, int w)
{
//reduces primery decomposition over algebraic extensions to
//the other cases
   def R=basering;
   int n=nvars(R);
   string mp="poly p="+string(minpoly)+";";
   string gnir="ring RH="+string(char(R))+",("+varstr(R)+","+string(par(1))
                +"),dp;";
   execute(gnir);
   execute(mp);
   ideal i=imap(R,i);
         i=i,p;
   list pr;
   int j;

   if(w==0)
   {
      pr=decomp(i);
   }
   if(w==1)
   {
      pr=prim_dec(i,1);
      pr=reconvList(pr);
   }
   if(w==2)
   {
      pr=minAssPrimes(i);
   }

   gnir="ring RS="+string(char(R))+",("+varstr(RH)
                +"),(dp("+string(n)+"),lp);";
   execute(gnir);
   list pr=imap(RH,pr);
   ideal K;
   for(j=1;j<=size(pr);j++)
   {
      K=groebner(pr[j]);
      K=K[2..size(K)];
      pr[j]=K;
   }
   setring R;
   list pr=imap(RS,pr);
      list re;
   if(w==2)
   {
      re=pr;
   }
   else
   {
      re=convList(pr);
   }
   return(re);
}

///////////////////////////////////////////////////////////////////////////////
static proc decomp(ideal i,list #)
"USAGE:  decomp(i); i ideal  (for primary decomposition)   (resp.
         decomp(i,1);        (for the associated primes of dimension of i) )
         decomp(i,2);        (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/Trager/Zacharias
EXAMPLE: example decomp; shows an example
"
{
  intvec op;
  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;
  int isS=(attrib(i,"isSB")==1);

  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;
      attrib(ltras[1],"isSB",1);
    }
    tras=ltras[1];
    attrib(tras,"isSB",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
  //----------------------------------------------------------------

  int lp=islp();

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

  ideal ser=fetch(@P,ser);

  if(homo==1)
  {
     if(!lp)
     {
       ideal @j=std(fetch(@P,i),@hilb,@w);
     }
     else
     {
        ideal @j=fetch(@P,tras);
        attrib(@j,"isSB",1);
     }
  }
  else
  {
      if(lp&&isS)
      {
        ideal @j=fetch(@P,i);
        attrib(@j,"isSB",1);
      }
      else
      {
        ideal @j=groebner(fetch(@P,i));
      }
  }
  option(set,op);
  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)
  {
    op=option(get);
    option(redSB);

    list gprimary= zero_decomp(@j,ser,@wr);
    option(set,op);
    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;
     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+";");
        op=option(get);
        option(redSB);
        if(homo==1)
        {
           ideal @j=std(phi(@j),@hilb,@w);
        }
        else
        {
          ideal @j=groebner(phi(@j));
        }
        ideal ser=phi(ser);

        option(set,op);
     }
     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..]
        op=option(get);
        option(redSB);

        list uprimary= zero_decomp(@j,ser,@wr);
        option(set,op);
     }
     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==1))
  {
     @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;

           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+";");
              op=option(get);
              option(redSB);
              if(homo==1)
              {
                 ideal @j=std(phi(jkeep),keephilb,@w);
              }
              else
              {
                ideal @j=groebner(phi(jkeep));
              }
              ideal ser=phi(ser);
              option(set,op);
           }

           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..]

           op=option(get);
           option(redSB);
           list uprimary= zero_decomp(@j,ser,@wr);
           option(set,op);

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

///////////////////////////////////////////////////////////////////////////////
static proc powerCoeffs(poly f,int e)
//computes a polynomial with the same monomials as f but coefficients
//the p^e th power of the coefficients of f
{
   int i;
   poly g;
   int ex=char(basering)^e;
   for(i=1;i<=size(f);i++)
   {
      g=g+leadcoef(f[i])^ex*leadmonom(f[i]);
   }
   return(g);
}
///////////////////////////////////////////////////////////////////////////////

static proc sep(poly f,int i, list #)
"USAGE:  input: a polynomial f depending on the i-th variable and optional
         an integer k considering the polynomial f defined over Fp(t1,...,tm)
         as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k))
 RETURN: the separabel part of f as polynomial in Fp(t1,...,tm)
        and an integer k to indicate that f should be considerd
        as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k))
 EXAMPLE: example sep; shows an example
{
   def R=basering;
   int k;
   if(size(#)>0){k=#[1];}

   poly h=gcd(f,diff(f,var(i)));
   poly g1=lift(h,f)[1][1];    //  f/h
   poly h1;

   while(h!=h1)
   {
      h1=h;
      h=gcd(h,diff(h,var(i)));
   }

   if(deg(h1)==0){return(list(g1,k));} //in characteristic 0 we return here

   k++;

   ideal ma=maxideal(1);
   ma[i]=var(i)^char(R);
   map phi=R,ma;
   ideal hh=h;    //this is technical because preimage works only for ideals

   poly u=preimage(R,phi,hh)[1]; //h=u(x(i)^p)

   list g2=sep(u,i,k);           //we consider u(t(1)^(p^-1),...,t(m)^(p^-1))
   g1=powerCoeffs(g1,g2[2]-k+1); //to have g1 over the same field as g2[1]

   list g3=sep(g1*g2[1],i,g2[2]);
   return(g3);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R=(5,t,s),(x,y,z),dp;
   poly f=(x^25-t*x^5+t)*(x^3+s);
   sep(f,1);
}

///////////////////////////////////////////////////////////////////////////////
 proc zeroRad(ideal I,list #)
"USAGE:  zeroRad(I) , I a zero-dimensional ideal
 RETURN: the radical of I
 NOTE:  Algorithm of Kemper
 EXAMPLE: example zeroRad; shows an example
{
   if(homog(I)==1){return(maxideal(1));}
   //I needs to be a reduced standard basis
   def R=basering;
   int m=npars(R);
   int n=nvars(R);
   int p=char(R);
   int i,k;
   list l;
   intvec op=option(get);

   option(redSB);
   ideal F=finduni(I);//F[i] generates I intersected with K[var(i)]

   option(set,op);
   if(size(#)>0){I=#[1];}

   for(i=1;i<=n;i++)
   {
      l[i]=sep(F[i],i);
      F[i]=l[i][1];
      if(l[i][2]>k){k=l[i][2];}  //computation of the maximal k
   }

   if((k==0)||(m==0)){return(interred(I+F));}        //the separable case

   for(i=1;i<=n;i++)             //consider all polynomials over
   {                             //Fp(t(1)^(p^-k),...,t(m)^(p^-k))
      F[i]=powerCoeffs(F[i],k-l[i][2]);
   }

   string cR="ring @R="+string(p)+",("+parstr(R)+","+varstr(R)+"),dp;";
   execute(cR);
   ideal F=imap(R,F);

   string nR="ring @S="+string(p)+",(y(1..m),"+varstr(R)+","+parstr(R)+"),dp;";
   execute(nR);

   ideal G=fetch(@R,F);    //G[i](t(1)^(p^-k),...,t(m)^(p^-k),x(i))=sep(F[i])

   ideal I=imap(R,I);
   ideal J=I+G;
   poly el=1;
   k=p^k;
   for(i=1;i<=m;i++)
   {
      J=J,var(i)^k-var(m+n+i);
      el=el*y(i);
   }

   J=eliminate(J,el);
   setring R;
   ideal J=imap(@S,J);
   return(J);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R=(5,t),(x,y),dp;
   ideal I=x^5-t,y^5-t;
   zeroRad(I);
}

///////////////////////////////////////////////////////////////////////////////
static proc radicalKL (ideal i,ideal ser,list #)
{
   attrib(i,"isSB",1);   // i needs to be a reduced standard basis
   list indep,fett;
   intvec @w,@hilb,op;
   int @wr,@n,@m,lauf,di;
   ideal fac,@h,collectrad,lsau;
   poly @q;
   string @va,quotring;

   def  @P = basering;
   int jdim=dim(i);
   int  homo=homog(i);
   ideal rad=ideal(1);
   ideal te=ser;
   if(size(#)>0)
   {
      @wr=#[1];
   }
   if(homo==1)
   {
      for(@n=1;@n<=nvars(basering);@n++)
      {
         @w[@n]=ord(var(@n));
      }
      @hilb=hilb(i,1,@w);
   }


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

   if (jdim==-1)
   {

      return(ideal(1));
   }

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

   if (jdim==0)
   {
      return(zeroRad(i));
   }
   //-------------------------------------------------------------------------
   //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(i);

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

     //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]);
         }
         op=option(get);
         option(redSB);
         @j=interred(@j);  //to obtain a reduced standardbasis
         attrib(@j,"isSB",1);
         option(set,op);

         ideal zero_rad= zeroRad(@j);
      }
      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

      zero_rad=std(zero_rad);

      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];
         }
         op=option(get);
         option(returnSB);
         option(redSB);
         i=quotient(i+ideal(@q),rad);
         attrib(i,"isSB",1);
         option(set,op);

      }
      if((deg(rad[1])>0)&&(deg(collectrad[1])>0))
      {
         rad=intersect(rad,collectrad);
         te=intersect(te,collectrad);
         te=simplify(reduce(te,i,1),2);
      }
      else
      {
         if(deg(collectrad[1])>0)
         {
            rad=collectrad;
            te=intersect(te,collectrad);
            te=simplify(reduce(te,i,1),2);
         }
      }

      if((dim(i)<jdim)||(size(te)==0))
      {
         break;
      }
      if(homo==1)
      {
         @hilb=hilb(i,1,@w);
      }
   }
   if(((@wr==1)&&(dim(i)<jdim))||(deg(i[1])==0)||(size(te)==0))
   {
      return(rad);
   }
   rad=intersect(rad,radicalKL(i,ideal(1),@wr));
   return(rad);
}

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

proc radicalEHV(ideal i)
"USAGE:   radicalEHV(i); i ideal.
RETURN:  ideal, the radical of i.
NOTE:    Uses the algorithm of Eisenbud/Huneke/Vasconcelos, which
         reduces the computation to the complete intersection case,
         by taking, in the general case, a generic linear combination
         of the input.
         Works only in characteristic 0 or p large.
EXAMPLE: example radicalEHV; shows an example
"
{
   if(ord_test(basering)!=1)
   {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
   }
   if((char(basering)<100)&&(char(basering)!=0))
   {
      "WARNING: The characteristic is too small, the result may be wrong";
   }
   ideal J,I,I0,radI0,L,radI1,I2,radI2;
   int l,n;
   intvec op=option(get);
   matrix M;

   option(redSB);
   list m=mstd(i);
        I=m[2];
   option(set,op);

   int cod=nvars(basering)-dim(m[1]);
   //-------------------complete intersection case:----------------------
   if(cod==size(m[2]))
   {
     J=minor(jacob(I),cod);
     return(quotient(I,J));
   }
   //-----first codim elements of I are a complete intersection:---------
   for(l=1;l<=cod;l++)
   {
      I0[l]=I[l];
   }
   n=dim(std(I0))+cod-nvars(basering);
   //-----last codim elements of I are a complete intersection:----------
   if(n!=0)
   {
      for(l=1;l<=cod;l++)
      {
         I0[l]=I[size(I)-l+1];
      }
      n=dim(std(I0))+cod-nvars(basering);
   }
   //-----taking a generic linear combination of the input:--------------
   if(n!=0)
   {
      M=transpose(sparsetriag(size(m[2]),cod,95,1));
      I0=ideal(M*transpose(I));
      n=dim(std(I0))+cod-nvars(basering);
   }
   //-----taking a more generic linear combination of the input:---------
   if(n!=0)
   {
      M=transpose(sparsetriag(size(m[2]),cod,0,100));
      I0=ideal(M*transpose(I));
      n=dim(std(I0))+cod-nvars(basering);
   }
   if(n==0)
   {
      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);
      }

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

      if(deg(I2[1])<=0)
      {
         return(radI1);
      }
      return(intersect(radI1,radicalEHV(I2)));
   }
   //---------------------general case-------------------------------------
   return(radical(I));
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 0,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z3+2;
   ideal i = p*q^2,y-z2;
   ideal pr= radicalEHV(i);
   pr;
}

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

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
static 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
static 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);
}
///////////////////////////////////////////////////////////////////////////////

static 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
//
///////////////////////////////////////////////////////////////////////////////

static 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
//////////////////////////////////////////////////////

static 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
//////////////////////////////////////////////////////


static 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
//////////////////////////////////////////////////////


static proc min_ass_prim_charsets0 (ideal PS)
{
  intvec op;
  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
       {
         op=option(get);
         option(returnSB);
         II=quotient(JS,I[j]);
         option(set,op);
         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
//////////////////////////////////////////////////////


static proc min_ass_prim_charsets1 (ideal PS)
{
  intvec op;
  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
       {
         op=option(get);
         option(returnSB);
         II=quotient(JS,I[j]);
         option(set,op);
//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.
//////////////////////////////////////////////////////////


static 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
//////////////////////////////////////////////////////////////////////////


static 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
////////////////////////////////////////////////////////////////


static 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,
////////////////////////////////////////////////////////////////


static 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
////////////////////////////////////////////////////////////////


static 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
//////////////////////////////////////////////////////////


static 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
//////////////////////////////////////////////////////////


static 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
/////////////////////////////////////////////////////////////////

static proc minquot(list tsil)
{
   intvec op;
   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);
   op=option(get);
   option(returnSB);
   ideal star=quotient(laedi,f);
   option(set,op);
   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);
         op=option(get);
         option(returnSB);
         verg=quotient(laedi,g);
         option(set,op);
         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
//////////////////////////////////////////////////

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


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

static 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);
}
///////////////////////////////////////////////////////////////////////////////

static 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 pr of primary ideals and their associated primes:
@format
   pr[i][1]   the i-th primary component,
   pr[i][2]   the i-th prime component.
@end format
NOTE:    Algorithm of Gianni/Trager/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
"
{
   if(ord_test(basering)!=1)
   {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
   }
   if(minpoly!=0)
   {
      return(algeDeco(i,0));
      ERROR(
      "// Not implemented yet for algebraic extensions.Simulate the ring extension by adding the minpoly to the ideal"
      );
   }
  return(convList(decomp(i)));
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 0,(x,y,z),lp;
   poly  p = z2+1;
   poly  q = z3+2;
   ideal i = p*q^2,y-z2;
   list pr = primdecGTZ(i);
   pr;
}
///////////////////////////////////////////////////////////////////////////////

proc primdecSY(ideal i, list #)
"USAGE:   primdecSY(i); i ideal, c int
RETURN:  a list pr of primary ideals and their associated primes:
@format
   pr[i][1]   the i-th primary component,
   pr[i][2]   the i-th prime component.
@end format
NOTE:    Algorithm of Shimoyama/Yokoyama.
@format
   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 are used.
@end format
EXAMPLE: example primdecSY; shows an example
"
{
   if(ord_test(basering)!=1)
   {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
   }
   i=simplify(i,2);
   if ((i[1]==0)||(i[1]==1))
   {
     list L=list(ideal(i[1]),ideal(i[1]));
     return(list(L));
   }
   if(minpoly!=0)
   {
      return(algeDeco(i,1));
   }
   if (size(#)==1)
   { return(prim_dec(i,#[1])); }
   else
   { return(prim_dec(i,1)); }
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 0,(x,y,z),lp;
   poly  p = z2+1;
   poly  q = z3+2;
   ideal i = p*q^2,y-z2;
   list pr = primdecSY(i);
   pr;
}
///////////////////////////////////////////////////////////////////////////////
proc minAssGTZ(ideal i,list #)
"USAGE:   minAssGTZ(i); i ideal
          minAssGTZ(i,1); i ideal  does not use the factorizing Groebner
RETURN:  a list, the minimal associated prime ideals of i.
NOTE:    Designed for characteristic 0, works also in char k > 0 based
         on an algorithm of Yokoyama
EXAMPLE: example minAssGTZ; shows an example
"
{
   if(ord_test(basering)!=1)
   {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
   }
   if(minpoly!=0)
   {
      return(algeDeco(i,2));
      ERROR(
      "// Not implemented for algebraic extensions. Simulate the ring extension by adding the minpoly to the ideal"
      );
   }
  if(size(#)==0)
   {
      return(minAssPrimes(i,1));
   }
   return(minAssPrimes(i));
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 0,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z3+2;
   ideal i = p*q^2,y-z2;
   list pr = minAssGTZ(i);
   pr;
}

///////////////////////////////////////////////////////////////////////////////
proc minAssChar(ideal i, list #)
"USAGE:   minAssChar(i[,c]); i ideal, c int.
RETURN:  list, the minimal associated prime ideals of i.
NOTE:    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. @*
         Due to a bug in the factorization, the result may be not completely
         decomposed in small characteristic.
EXAMPLE: example minAssChar; shows an example
"
{
   if(ord_test(basering)!=1)
   {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
   }
   if (size(#)==1)
   { return(min_ass_prim_charsets(i,#[1])); }
   else
   { return(min_ass_prim_charsets(i,1)); }
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 0,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z3+2;
   ideal i = p*q^2,y-z2;
   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:    A combination of the algorithms of Krick/Logar and Kemper is used.
         Works also in positive characteristic (Kempers algorithm).
EXAMPLE: example equiRadical; shows an example
"
{
   if(ord_test(basering)!=1)
   {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
   }
   return(radical(i,1));
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 0,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z3+2;
   ideal i = p*q^2,y-z2;
   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 Kemper is used.
         Works also in positive characteristic (Kempers algorithm).
EXAMPLE: example radical; shows an example
"
{
   if(ord_test(basering)!=1)
   {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
   }
   def @P=basering;
   int j,il;
   if(size(#)>0){il=#[1];}
   if(size(i)==0){return(ideal(0));}
   ideal re=1;
   intvec op = option(get);
   list qr=simplifyIdeal(i);
   map phi=@P,qr[2];

   option(redSB);
   i=groebner(qr[1]);
   option(set,op);
   int di=dim(i);

   if(di==0)
   {
     i=zeroRad(i,qr[1]);
     return(interred(phi(i)));
   }

   option(redSB);
   list pr=i;
   if (!homog(i))
   {
     pr=facstd(i);
   }
   option(set,op);
   int s=size(pr);

   for(j=1;j<=s;j++)
   {
      attrib(pr[s+1-j],"isSB",1);
      if((size(reduce(re,pr[s+1-j],1))!=0)&&((dim(pr[s+1-j])==di)||!il))
      {
         re=intersect(re,radicalKL(pr[s+1-j],re,il));
      }
   }
   return(interred(phi(re)));
}
example
{ "EXAMPLE:";  echo = 2;
   ring  r = 0,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z3+2;
   ideal i = p*q^2,y-z2;
   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/Vasconcelos.
EXAMPLE: example prepareAss; shows an example
"
{
  if(ord_test(basering)!=1)
  {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
  }
  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 = 0,(x,y,z),dp;
   poly  p = z2+1;
   poly  q = z3+2;
   ideal i = p*q^2,y-z2;
   list pr = prepareAss(i);
   pr;
}
///////////////////////////////////////////////////////////////////////////////
proc equidimMaxEHV(ideal i)
"USAGE:  equidimMaxEHV(i); i ideal
RETURN:  ideal, the equidimensional component (of maximal dimension) of i.
NOTE:    Uses algorithm of Eisenbud, Huneke and Vasconcelos.
EXAMPLE: example equidimMaxEHV; shows an example
"
{
  if(ord_test(basering)!=1)
  {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
  }
  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 = 0,(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, k an ideal.
ASSUME:  pr is the result of primdecGTZ(k) or primdecSY(k).
RETURN:  int, 1 if the intersection of the 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
ASSUME:  I is zero-dimensional, the characteristic of the ground field is 0
RETURN:  list of primary ideals, the zero-dimensional decomposition of I
NOTE:    The algorithm (of Monico), works well only for a small total number
         of solutions (@code{vdim(std(I))} should be < 100) and without
         parameters. In practice, it works also in large characteristic p>0
         but may fail for small p.
@*       If printlevel > 0 (default = 0) additional information is displayed.
EXAMPLE: example zerodec; shows an example
"
{
   if(ord_test(basering)!=1)
   {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
   }
   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);
time=timer; ideal ra =radical(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)
//radical : 0
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;
time=timer; ideal ra =radical(gls); timer-time;size(pr);
LIB"all.lib";

ring R=0,(a,b,c,d,e,f),dp;
ideal I=cyclic(6);
minAssGTZ(I);


ring S=(2,a,b),(x,y),lp;
ideal I=x8-b,y4+a;
minAssGTZ(I);

ring S1=2,(x,y,a,b),lp;
ideal I=x8-b,y4+a;
minAssGTZ(I);


ring S2=(2,z),(x,y),dp;
minpoly=z2+z+1;
ideal I=y3+y+1,x4+x+1;
primdecGTZ(I);
minAssGTZ(I);

ring S3=2,(x,y,z),dp;
ideal I=y3+y+1,x4+x+1,z2+z+1;
primdecGTZ(I);
minAssGTZ(I);


ring R1=2,(x,y,z),lp;
ideal I=y6+y5+y3+y2+1,x4+x+1,z2+z+1;
primdecGTZ(I);
minAssGTZ(I);


ring R2=(2,z),(x,y),lp;
minpoly=z3+z+1;
ideal I=y2+y+(z2+z+1),x4+x+1;
primdecGTZ(I);
minAssGTZ(I);

*/