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

////////////////////////////////////////////////////////////////////////////
version="version primdec.lib 4.4.0.0 Dec_2023 "; // $Id$
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)@*
          Santiago Laplagne, slaplagn@dm.uba.ar         (GTZ)

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, Laplagne and Kemper (implementation by Gerhard Pfister and Santiago Laplagne).
    They work in any characteristic.@*
    Baserings must have a global ordering and no quotient ideal.
    Exceptions: primdecGTZ, absPrimdecGTZ, minAssGTZ, primdecSY, minAssChar, radical accept non-global ordering.


PROCEDURES:
 Ann(M);            annihilator of R^n/M, R=basering, M in R^n
 primdecGTZ(I);     complete primary decomposition via Gianni,Trager,Zacharias
 primdecGTZE(I);    complete primary decomposition via Gianni,Trager,Zacharias. Returns empty list for the unit ideal
 primdecSY(I...);   complete primary decomposition via Shimoyama-Yokoyama
 primdecSYE(I,..);  complete primary decomposition via Shimoyama-Yokoyama. Returns empty list for the unit ideal
 minAssGTZ(I);      the minimal associated primes via Gianni,Trager,Zacharias (with modifications by Laplagne)
 minAssGTZE(I);     the minimal associated primes via Gianni,Trager,Zacharias. Returns empty list for unit ideal
 minAssChar(I...);  the minimal associated primes using characteristic sets
 minAssCharE(I..);  the minimal associated primes using characteristic sets. Returns empty list for unit ideal
 testPrimary(L,k);  tests the result of the primary decomposition
 testPrimaryE(L,k); tests the result of the primary decomposition. Handles also empty list L.
 radical(I);        computes the radical of I via Krick/Logar (with modifications by Laplagne) 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
 absPrimdecGTZ(I);  the absolute prime components of I
 absPrimdecGTZE(I); the absolute prime components of I. Assumes I is not unit ideal.
 sep(f,k);          the separabel part of f as polynomial in Fp(t1,...,tm)

SEE ALSO: primdecint_lib
KEYWORDS: primary decomposition
";

LIB "general.lib";
LIB "elim.lib";
LIB "polylib.lib";
LIB "random.lib";
LIB "inout.lib";
LIB "matrix.lib";
LIB "triang.lib";
LIB "absfact.lib";
LIB "ring.lib";
LIB "modfinduni.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
"
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  int @k;
  ideal inew=id;
  if (!attrib(id,"isSB")) { 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
"
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  int @k,@i;
  def @P= basering;
  if(ordstr(basering)[1,2]!="dp")
  {
    def @Phelp=changeordTo(basering,"dp");
    setring @Phelp;
    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,6);
  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)
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  ASSUME(0, hasGlobalOrdering(basering) ) ;
  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,6);
  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)
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  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 factor(poly p)
"USAGE:   factor(p) p poly
RETURN:  list=;
NOTE:
EXAMPLE: example factor; shows an example
"
{
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  ideal @i;
  list @l;
  intvec @v,@w;
  int @j,@k,@n;

  @l=factorize(p);
    for(@j=1;@j<=size(@l[1]);@j++)
    {
      if(leadcoef(@l[1][@j])==@l[1][@j])
      {
        @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(leadcoef(@l[1][@j])!=@l[1][@j])
          {
            @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)
{
  if (!attrib(k,"isSB")) { k=std(k);}
  if (!attrib(j,"isSB")) { j=std(j);}
  return(stdIdealsEqual(k,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,5))==0)
    {
      return(1);
    }
  }
  return(0);
}
///////////////////////////////////////////////////////////////////////////////

proc primaryTest (ideal i, poly p)
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  ASSUME(0, hasGlobalOrdering(basering) ) ;
  if(i[1]==1){return(ideal(1));}
  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--;
    m++;

    //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--;
    }
    //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 >=16 )
        {
          "// WARNING: The characteristic is perhaps too small to use";
          "// the algorithm of Gianni/Trager/Zacharias.";
          "// This may result in an infinite loop";
          "// loop in primaryTest, voice:",voice;"";
        }
        while(e mod char(basering)==0)
        {
          f=f*char(basering);
          e=e div 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,5) !=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,5) !=0)
      {
        return(ideal(0));
      }
    }

    h=interred(t);
    t=h[1];

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

///////////////////////////////////////////////////////////////////////////////
proc gcdTest(ideal act)
{
  ASSUME(0, not isQuotientRing(basering) ) ;
  ASSUME(0, hasGlobalOrdering(basering) ) ;
  int i,j;
  if(size(act)<=1)
  {
    return(0);
  }
  for (i=1;i<size(act);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)
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  int i,j,k,s,r,w;
  list keepresult,act,keepprime;
  poly @f;
  int sl=size(l);
  for(i=sl div 2;i>=1;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) div 2)
  {
    i++;
    if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],5))==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 div 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) div 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 div 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 div 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+2]=ideal(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 div 2+1]=interred(keepprime[i]+ideal(@f));
            if(homog(keepprime[s div 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) div 2;i++)
  {
    attrib(l[2*i-1],"isSB",1);

    if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],5))==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)
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  if((char(basering)==0)||(npars(basering)>0))
  {
    return(l);
  }
  def op = option(get);
  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) div 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) div 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];
        ring RL = create_ring("("+string(char(P))+","+string(var(n))+")", "("+va+")", "lp");
        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]);
          if(sa[1][1]==l[2*i-1][1])
          {
             l[2*i-1]=std(sa[1]);
             l[2*i]=primaryTest(sa[1],s[1]);
          }
          else
          {
             l[2*i-1]=std(sa[1]);
             l[2*i]=primaryTest(sa[1],factorize(sa[1][1],1)[1]);
          }
        }
        if((size(sact)==1)&&(m==2))
        {
          l[2*i]=std(l[2*i-1],sact[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(set,op);
        }
        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],factorize(keep[i][1],1)[1]);
      }
    }
  }
  l=l+keep1;
  option(set,op);
  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 computation 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
"
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  ASSUME(0, hasGlobalOrdering(basering) ) ;

  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);
  }
  intvec save=option(get);
  option(redSB);
  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],5))==0))
      {
        primary[2*@k-1]=ideal(1);
        primary[2*@k]=ideal(1);
      }
    }
    option(set,save);
    return(primary);
  }

  option(set,save);
  if(homog(j)==1)
  {
    primary[1]=j;
    if((size(ser)>0)&&(size(reduce(ser,j,5))==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]);
  }
  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(@wr!=0)
    {
      primary[1]=std(j,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>=7)&&(char(basering)<=181)&&(@wr==1))
  {
    primary=splitCharp(primary);
  }
/*
  if((@wr==2)&&(npars(basering)>0)&&(voice>=7)&&(char(basering)>0))
  {
  //the prime decomposition of Yokoyama in characteristic p
    list ke,ek;
    @k=0;
    while(@k<size(primary) div 2)
    {
      @k++;
      if(size(primary[2*@k])==0)
      {
        ek=insepDecomp_i( 1, 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>=9){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>=9))
  {
    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) div 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) div 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>=11))
      {
        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 ============================

      if (npars(basering)>1) { @qh=groebner(@qh,"par2var"); }
      else                   { @qh=groebner(@qh); }

      ser1=phi1(ser);
      @lh=zero_decomp (@qh,phi(ser1),@wr);

      kill lres0;
      list lres0;
      if((size(@lh)==2)&&(@lh[1]!=1))
      {
        helpprim=@lh[2];
        lres0[1]=primary[2*@k-1];
        attrib(lres0[1],"isSB",1);
        ser1=psi(helpprim);
        lres0[2]=psi1(ser1);
        if(size(reduce(lres0[2],lres0[1],5))==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) div 2;@n++)
        {
          lres0[2*@n-1]=groebner(lres0[2*@n-1]);
          lres0[2*@n]=groebner(lres0[2*@n]);
        }
      }
      primary=primary+lres0;

   }
  }
  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 #)
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  ASSUME(0, hasGlobalOrdering(basering) ) ;
//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);}

  ring RH = create_ring("("+string(p^ex)+",@a)", "("+varstr(R)+")", "lp");
  string mp="minpoly="+string(minpoly)+";";
  ring RL = create_ring("("+string(p)+",@a)", "("+varstr(R)+")", "lp");
  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)+";";

  ring RS = create_ring(p, "("+varstr(R)+",@a)", "lp");
  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) div 2-1)
  {
    i++;
    j=i;
    while(j<size(L) div 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

  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  ASSUME(0, hasGlobalOrdering(basering) ) ;
  def R=basering;

  //i has to be a reduced groebner basis
  ASSUME(1, dim(i)==0);
  if (hasQQCoefficient(basering) && (nvars(basering)*vdim(I)>300))
  {
    ideal F=modFinduni(i);
  }
  else
  {
    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;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;
  }
  def P=changeordTo(R,"dp");
  setring P;
  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])

  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  ASSUME(0, hasGlobalOrdering(basering) ) ;

  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_i(int patchPrimaryDecomposition, ideal i)
{
//decomposes i into special ideals
//computes the prime decomposition of the special ideals
//and transforms it back to a decomposition of i
// if patchPrimaryDecomposition=1,  drop unit ideal in the decomposition,
// since the unit ideal it is not prime!

  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  ASSUME(0, hasGlobalOrdering(basering) ) ;
  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_i(patchPrimaryDecomposition,wo[1],2);

    option(redSB);
    for(k=1;k<=size(qr) div 2;k++)
    {
      K=qr[2*k];
      K=phi(K);
      K=groebner(K);
      re[size(re)+1]=zeroRad(K);
    }
    option(set,op);
  }
  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
"
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  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 clearSBNeu (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
"
{
 ASSUME(1, hasFieldCoefficient(basering) );
 ASSUME(1, not isQuotientRing(basering) ) ;
 ASSUME(1, hasGlobalOrdering(basering) ) ;
 int k,j;
 intvec m,n,v,w;
 int c=size(i);
 w=leadexp(0);
 v[size(i)]=0;

 j=0;
 while(j<c-1)
 {
   j++;
   if(deg(i[j])>=0)
   {
      m=leadexp(i[j]);
      for(k=j+1;k<=c;k++)
      {
        n=leadexp(i[k]);
        if(n!=w)
        {
           if(((m==n)&&(#[j]>#[k]))||((teilt(n,m))&&(n!=m)))
           {
             i[j]=0;
             v[j]=1;
             break;
           }
           if(((m==n)&&(#[j]<=#[k]))||((teilt(m,n))&&(n!=m)))
           {
             i[k]=0;
             v[k]=1;
           }
        }
      }
    }
  }
  return(v);
}

static proc teilt(intvec a, intvec b)
{
  int i;
  for(i=1;i<=size(a);i++)
  {
    if(a[i]>b[i]){return(0);}
  }
  return(1);
}
///////////////////////////////////////////////////////////////////////////////

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
"
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  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,string order)
"USAGE:   prepareQuotientring(nnp, order); nnp int, order string
RETURN:  Kvar(nnp+1),...,var(nvars)[..rest ]
EXAMPLE: example prepareQuotientring; shows an example
"
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  list rl=ringlist(basering);
  if (typeof(rl[1])=="int")
  {
    int p=rl[1];
    list rl2=rl[2];
    rl[1]=list(p,
            list(rl2[nnp+1..nvars(basering)]),
            list(list("lp",1:(nvars(basering)-nnp))),
            ideal(0));
    rl[2]=list(rl2[1..nnp]);
    rl[3]=list(list(order,1:nnp),list("C",0));
  }
  else
  {
    if (typeof(rl[1])=="list")
    {
      if(size(rl[1])==4)
      {
        list rl1=rl[1];
        list rl2=rl[2];
        rl1=list(rl1[1][1],
                rl[1][2]+list(rl2[nnp+1..nvars(basering)]),
                list(list("lp",1:(size(rl[1][2])+nvars(basering)-nnp))),
                ideal(0));
        rl[1]=rl1;
        rl[2]=list(rl2[1..nnp]);
        rl[3]=list(list(order,1:nnp),list("C",0));
      }
      else /*Float(..), Complex(..) */
      {
        list p=rl[1];
        list rl2=rl[2];
        rl[1]=list(p,
            list(rl2[nnp+1..nvars(basering)]),
            list(list("lp",1:(nvars(basering)-nnp))),
            ideal(0));
        rl[2]=list(rl2[1..nnp]);
        rl[3]=list(list(order,1:nnp),list("C",0));
      }
    }
    else
    {
        ERROR("Unexpected case in prepareQuotientring. Please inform the authors");
    }
  }

  def quotring=ring(rl);
  return(quotring);
}
example
{ "EXAMPLE:"; echo = 2;
   ring s1=(0,x),(a,b,c,d,e,f,g),lp;
   def Q= prepareQuotientring(3,"lp");
   Q;
}

///////////////////////////////////////////////////////////////////////////////
static proc cleanPrimary(list l)
{
  int i,j;
  list lh;
  for(i=1;i<=size(l) div 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 minAssPrimesoldE(ideal I, list #)
"USAGE:   minAssPrimesoldE(I); I ideal
         minAssPrimesold(I,1); I ideal  (to use also the factorizing Groebner)
RETURN:  list = the minimal associated prime ideals of I
EXAMPLE: example minAssPrimesoldE; shows an example
"
{
    return(minAssPrimesold_i(1,I,#));
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 32003,(x,y,z),lp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   list pr= minAssPrimesoldE(i);  pr;

   minAssPrimesoldE(i,1);
}

proc minAssPrimesold(ideal I, list #)
"USAGE:   minAssPrimesold(I); I ideal
         minAssPrimesold(i,1); I ideal  (to use also the factorizing Groebner)
RETURN:  list = the minimal associated prime ideals of I. In case I is unit ideal, returns list(ideal(1));
EXAMPLE: example minAssPrimesold; shows an example
"
{
    return(minAssPrimesold_i(0,I,#));
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 32003,(x,y,z),lp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   list pr= minAssPrimesold(i);  pr;

   minAssPrimesold(i,1);
}

static proc minAssPrimesold_i(int patchPrimaryDecomposition, ideal i, list #)
{
//
// parameter patchPrimaryDecomposition : if = 1, patch the decomposition( drop unit ideal in the decomposition),
//                                 : if = 0, taken no special action in case the unit ideal is in the decomposition
// for other parameters see minAssPrimesold, minAssPrimesoldE

   ASSUME(1, hasFieldCoefficient(basering) );
   ASSUME(0, not isQuotientRing(basering) ) ;
   ASSUME(0, hasGlobalOrdering(basering) ) ;
   def @P=basering;
   if(size(i)==0) { return(list(ideal(0))); }
   list qr=simplifyIdeal(i);
   map phi=@P,qr[2];
   i=qr[1];

   def gnir=ring(ring_list(@P));
   setring gnir;

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

      @res[1]=primary;

      tluser=union(@res);

      setring @P;
      if (size(tluser)>0)
      {
          list @res=imap(gnir,tluser);
          return(phi(@res));
      }
      else
      {
          return(tluser);
      }
   }
   list @res,empty;
   ideal ser;
   def op = option( get );
   option( redSB );
   list @pr=facstd(i);

   // option( noredSB );
   option( set, op );
   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;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_i(patchPrimaryDecomposition,@pr[j],2);
        }
        else
        {
           @res[j]=empty;
        }
     }
   }
   else
   {
     ser=ideal(1);
     for(j=1;j<=size(@pr);j++)
     {
//@pr[j];
//pause();
        @res[j]=decomp_i(patchPrimaryDecomposition,@pr[j],2);
//       @res[j]=decomp_i(patchPrimaryDecomposition,@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));
}


static proc primT(ideal i)
{
   ASSUME(1, hasFieldCoefficient(basering) );
   ASSUME(1, not isQuotientRing(basering) ) ;
   ASSUME(1, hasGlobalOrdering(basering) ) ;

   //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 minAssPrimesE(ideal I, list #)
"USAGE:   minAssPrimesE(I); I ideal
      Optional parameters in list #: (can be entered in any order)
      0, "facstd"   ->   uses facstd to first decompose the ideal
      1, "noFacstd" ->  does not use facstd (default)
      "SL" ->     the new algorithm is used (default)
      "GTZ" ->     the old algorithm is used
RETURN:  list = the minimal associated prime ideals of I
EXAMPLE: example minAssPrimesE; shows an example
"
{
  return(minAssPrimes_i(1,I,#));
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 32003,(x,y,z),lp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal I = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   list pr= minAssPrimesE(I);  pr;

   minAssPrimesE(I,1);
}

static proc minAssPrimes(ideal I, list #)
"USAGE:   minAssPrimes(I); I ideal
      Optional parameters in list #: (can be entered in any order)
      0, "facstd"   ->   uses facstd to first decompose the ideal
      1, "noFacstd" ->  does not use facstd (default)
      "SL" ->     the new algorithm is used (default)
      "GTZ" ->     the old algorithm is used
RETURN:  list = the minimal associated prime ideals of I. If I is the unit ideal returns list(ideal(1)) ;
EXAMPLE: example minAssPrimes; shows an example
"
{
  return(minAssPrimes_i(0,I,#));
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 32003,(x,y,z),lp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   list pr= minAssPrimes(i);  pr;

   minAssPrimes(i,1);
}


static proc minAssPrimes_i(int patchPrimaryDecomposition, ideal i, list #)
{
// parameter  patchPrimaryDecomposition:  1 to patch( remove unit ideal from the decomposition) ,
//                                        0 for no special action on unit ideal.
// for other parameters see 'minAssPrimes', 'minAssPrimesE'
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  if (size(i) == 0)  {  return(list(ideal(0)));  }

  intvec origOp = option(get);

  string algorithm;    // Algorithm to be used
  string facstdOption;    // To uses proc facstd
  int j;          // Counter
  def P0 = basering;
  list Pl=ring_list(P0);
  Pl[3]=list(list("dp",1:nvars(P0)),list("C",0));
  def P=ring(Pl);
  setring P;
  ideal i=imap(P0,i);

  // Set input parameters
  algorithm = "SL";         // Default: SL algorithm
  facstdOption = "Facstd";    // Default: facstd is not used
  int subsystem=0;
  if(size(#) > 0)
  {
    int valid;
    for(j = 1; j <= size(#); j++)
    {
      valid = 0;
      if((typeof(#[j]) == "int") or (typeof(#[j]) == "number"))
      {
        if (#[j] == 0) {facstdOption = "noFacstd"; valid = 1;}    // If #[j] == 0, facstd is not used.
        if (#[j] == 1) {facstdOption = "facstd";   valid = 1;}    // If #[j] == 1, facstd is used.
      }
      if(typeof(#[j]) == "string")
      {
        if(#[j] == "GTZ" || #[j] == "SL")
        {
          algorithm = #[j];
          valid = 1;
        }
        if(#[j] == "noFacstd" || #[j] == "facstd")
        {
          facstdOption = #[j];
          valid = 1;
        }
        if(#[j] == "noSubsystem" || #[j] == "subsystem")
        {
          subsystem = (#[j]=="subsystem");
          valid = 1;
        }
      }
      if(valid == 0)
      {
        dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
      }
    }
  }

  i=simplify(i,2);
  intvec op = option(get);
  list q = simplifyIdeal(i);
  list re = maxideal(1);
  int a;
  map phi = P,q[2];
  i=q[1];
  if (subsystem==0)
  {
    int k;

    list result;

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

    if (attrib(i,"isSB")!=1)
    {
      i=groebner(q[1]);
    }
    else
    {
      for(j=1;j<=nvars(basering);j++)
      {
        if (q[2][j]!=var(j)) {k=1;break;}
      }
      if(k)
      {
        i=groebner(q[1]);
      }
    }
    if(size(i)==1)
    {
      if ( deg(lead(i[1]))==0 ) // we have the unit ideal.
      {
        setring P0;
        option( set,origOp );
        if (patchPrimaryDecomposition==1)
        {
          return( list() );
        }
        else
        {
          return( list(ideal(1)) );
        }
      }
    }

    if( dim(i) == -1 )
    {
      option( set,op );
      setring P0;
      option( set,origOp );
      return( ideal(1) );
    }
    if( (dim(i) == 0 ) && ( npars(P) == 0) )
    {
      int di = vdim(i);
      def gnir=changeordTo(P,"lp");
      setring gnir;
      ideal J = std(imap(P,i));
      attrib(J, "isSB", 1);
      if(vdim(J) != di)
      {
        J = fglm(P, i);
      }
  //    list pr = triangMH(J,2); HIER KOENNEN verschiedene Mengen zu gleichen
  //                             asoziierten Primidealen fuehren
  // Aenderung
      list pr = triangMH(J,2);
      list qr, re;
      for(k = 1; k <= size(pr); k++)
      {
        if(primT(pr[k])&&(0))
        {
          re[size(re) + 1] = pr[k];
        }
        else
        {
          attrib(pr[k], "isSB", 1);
          // Lines changed
          if (algorithm == "GTZ")
          {
            qr = decomp_i(patchPrimaryDecomposition,pr[k], 2);
          }
          else
          {
            qr = minAssSL(pr[k]);
          }
          for(j = 1; j <= size(qr) div 2; j++)
          {
            re[size(re) + 1] = std(qr[2 * j]);
          }
        }
      }
      setring P;
      re = imap(gnir, re);
      re=phi(re);
      option(set, op);
      setring(P0);
      list re=imap(P,re);
      option( set,origOp );
      return(re);
    }

    if ((facstdOption == "noFacstd") || (dim(i) == 0))
    {
      if (algorithm == "GTZ")
      {
        re[1] = decomp_i(patchPrimaryDecomposition,i, 2);
      }
      else
      {
        re[1] = minAssSL(i);
      }
      re = union(re);
      option(set, op);
      re=phi(re);
      setring(P0);
      option( set,origOp );
      list re=imap(P,re);
      return(re);
    }
  }
  list qq;
  if (subsystem) { qq = prep_decomp(i); }
  else { qq=list(i); }
  q=list();
  for(int ll=size(qq);ll>0;ll--)
  {
    q=q+facstd(qq[ll]);
  }
  kill qq; kill ll;

  option(set,op);
  // Debug
  dbprint(printlevel - voice, "Components returned by facstd", size(q), q);
  for(j = 1; j <= size(q); j++)
  {
    if(a == 0){attrib(q[j], "isSB", 1);}
    // Debug
    dbprint(printlevel - voice, "We compute the decomp of component", j);
    // Lines changed
    if (algorithm == "GTZ")
    {
      re[j] = decomp_i(patchPrimaryDecomposition,q[j], 2);
    }
    else
    {
      re[j] = minAssSL(q[j]);
    }
    // Debug
    dbprint(printlevel - voice, "Number of components obtained for this component:", size(re[j]) div 2);
    dbprint(printlevel - voice, "re[j]:", re[j]);
  }
  re = union(re);
  re=phi(re);
  setring(P0);
  list re=imap(P,re);
  option( set,origOp );
  return(re);
}


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

  def P=basering;

  int liSize=size(li);
  int li1Size=0;
  if (size(li)>0)
  {
     li1Size=size(li[1]);
  }
  def ir=changeordTo(basering,"lp");
  setring ir;
  list l;
  if ( liSize > 0)
  {
     if (li1Size > 0)
     {
         l = fetch(P,li);
     }
     else
     {
         ASSUME(1, size(li)==1);
         l[1] = list();
     }
  }
  list @erg;

  for(k=1;k<=size(l);k++)
  {
     for(j=1;j<=size(l[k]) div 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,5))==0)
        {
           if(size(reduce(@erg[i],@erg[k],5))==0)
           {
              @erg[k]=ideal(1);
              i2=ideal(1);
           }
        }
        if(size(reduce(i2,i1,5))==0)
        {
           if(size(reduce(@erg[k],@erg[i],5))==0)
           {
              break;
           }
        }
        k++;
        if(k>size(@erg))
        {
           @wos[size(@wos)+1]=@erg[i];
        }
     }
  }
  if (size(@erg)>0)
  {
     if(deg(@erg[size(@erg)][1])!=0)
     {
        @wos[size(@wos)+1]=@erg[size(@erg)];
     }
  }
  int @wosSize = size(@wos);
  setring P;
  list @ser;
  if (@wosSize>0)
  {
      @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, except I is unit ideal.
        - a[1],...,a[s-1] the lower dimensional equidimensional loci.
         If I is the unit ideal, a list containing the unit ideal as a[1] is returned.
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
"
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  if(attrib(basering,"global")!=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))&&(find(ordstr(basering),"l")==0)
                                &&(find(ordstr(basering),"s")==0))
  {
     def gnir=ring(ring_list(basering));
     setring gnir;
     ideal i=imap(P,i);
     ideal j=i;
     if(a)
     {
       attrib(j,"isSB",1);
     }
     else
     {
       j=groebner(i);
     }
  }
  else
  {
     def gnir=changeordTo(basering,"dp");
     setring gnir;
     ideal i=imap(P,i);
     ideal j=groebner(i);
  }
  if(homo==1)
  {
     for(n=1;n<=nvars(basering);n++)
     {
        w[n]=ord(var(n));
     }
     bigintvec 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,5))==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
"
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  if(attrib(basering,"global")!=1)
  {
      ERROR(
      "Not implemented for this ordering, please change to a 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))&&(find(ordstr(basering),"l")==0)
                                &&(find(ordstr(basering),"s")==0))
  {
     def gnir=ring(ring_list(basering));
     setring gnir;
     ideal i=imap(P,i);
     ideal j=i;
     if(a)
     {
       attrib(j,"isSB",1);
     }
     else
     {
       j=groebner(i);
     }
  }
  else
  {
     def gnir=changeordTo(basering,"dp");
     setring gnir;
     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));
     }
     bigintvec 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);
  ring gnir1 = create_ring(ring_list(basering)[1],"("+indep[1][1]+")","("+indep[1][2]+")","no_minpoly");
  if(homo==1)
  {
     ideal j=std(imap(gnir,j),hil,w);
  }
  else
  {
     ideal j=groebner(imap(gnir,j));
  }
  def quotring=prepareQuotientring(nvars(basering)-indep[1][3],"lp");
  setring 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(quotring,h);
  kill quotring;

  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);
}
///////////////////////////////////////////////////////////////////////////////
//w=0: GTZ
//w=1: SY
//w=2: minAss
proc algeDecoE(ideal I, int w)
{//reduces primery decomposition over algebraic extensions to
//the other cases
    return( algeDeco_i( 1, I, w) );
}

//w=0: GTZ
//w=1: SY
//w=2: minAss
// deprecated. use 'algeDecoE()'
proc algeDeco(ideal I, int w)
{//reduces primery decomposition over algebraic extensions to
//the other cases
    return( algeDeco_i(0, I, w));
}


//w=0: GTZ
//w=1: SY
//w=2: minAss
static proc algeDeco_i(int patchPrimaryDecomposition, ideal i, int w)
{//reduces primery decomposition over algebraic extensions to
//the other cases
// if patchPrimaryDecomposition=1,  drop unit ideal in the decomposition,
// since the unit ideal it is not prime, otherwise take no special action.

   ASSUME(0, hasFieldCoefficient(basering) );
   ASSUME(0, not isQuotientRing(basering) ) ;
   ASSUME(0, hasGlobalOrdering(basering) ) ;

// the really needed things:
   ASSUME(1, typeof(ringlist(basering)[1])=="list"); // in alg. extension

//reduces primery decomposition over algebraic extensions to
//the other cases
   def R=basering;
   int n=nvars(R);

   intvec op = option(get);

//---Anfang Provisorium
   if((size(i)==2) && (w==2))
   {
      //treats a special case separately which would otherwise take a lot longer in factorization
      option( redSB );
      ideal J = std(i);
      option( set, op );
      if(size(J)==1)
      {
         if ( deg(lead(J[1]))==0 ) // we have the unit ideal
         {
            if (patchPrimaryDecomposition==1)
            {
                return( list() );
            }
            else
            {
                return( list( ideal(1) ) );
            }
         }
      }
      if ((size(J)==2)&&(deg(J[1])==1))
      {
         // minAssPrimes correspond to factorization of J[2]
         ideal keep;
         poly f;
         int j;
         for(j=1;j<=nvars(basering);j++)
         {
           f=J[2];
           while((f/var(j))*var(j)-f==0)
           {
             f=f/var(j);
             keep=keep,var(j);
           }
           J[2]=f;
         }
         ideal K=factorize(J[2],1);
         if(deg(K[1])==0){K=0;}
         K=K+std(keep);
         ideal L;
         list resu;
         for(j=1;j<=size(K);j++)
         {
            L=J[1],K[j];
            resu[j]=L;
         }
         option( set, op );
         return(resu);
      }
   }
//---Ende Provisorium
   list R_l=ringlist(R);
   ideal @p=R_l[1][4]; // minpoly
   R_l[2]=R_l[2]+R_l[1][2]; // vars
   R_l[1]=R_l[1][1];  // char
   R_l[3]=list(list("dp",1:size(R_l[2])),list("C",0)); // ord
   def RH=ring(R_l); kill R_l;setring RH;
   ideal @pp=imap(R,@p); poly @p=@pp[1];
   ideal i=imap(R,i);
   ideal I=subst(i,var(nvars(basering)),0);
   int j;
   for(j=1;j<=ncols(i);j++)
   {
     if(i[j]!=I[j]){break;}
   }
   if((j>ncols(i))&&(deg(@p)==1))
   {
     setring R;
     kill RH;
     // remove extension, set order to dp:
     list R_l=ringlist(R);
     R_l[1]=R_l[1][1]; // char
     R_l[3]=list(list("dp",1:nvars(R)),list("C",0)); // ord
     def RH=ring(R_l); kill R_l; setring RH;
     ideal i=imap(R,i);
     ideal J;
   }
   else
   {
      i=i,@p;
   }
   list pr;

   if(w==0)
   {
      pr=decomp_i(patchPrimaryDecomposition,i);
   }
   if(w==1)
   {
      pr=prim_dec_i(patchPrimaryDecomposition,i,1);
      pr=reconvList(pr);
   }
   if(w==2)
   {
      pr=minAssPrimes_i(patchPrimaryDecomposition,i);
   }

   int sizepr = size(pr);

   if(n<nvars(basering))
   {
      // remove extension, set order to dp(n),lp:
      list R_l=ringlist(basering);
      if (typeof(R_l[1])=="list") { R_l[1]=R_l[1][1]; }
      R_l[3]=list(list("dp",1:n),list("lp",1:(nvars(basering)-n)),list("C",0));
      def RS=ring(R_l); kill R_l; setring RS;
      if (sizepr>0) { list pr=imap(RH,pr); }
      ideal K;
      for(j=1;j<=sizepr;j++)
      {
         K=groebner(pr[j]);
         if (size(K)>1)
         {
             K = K[2..size(K)];
         }
         pr[j]=K;
      }
      setring R;
      if (sizepr>0) {  list pr=imap(RS,pr); }
   }
   else
   {
      setring R;
      if (sizepr>0)  { list pr=imap(RH,pr); }
   }

   list re;
   if(w==2)
   {
      re=pr;
   }
   else
   {
      re=convList(pr);
   }
   option( set, op );
   return( re );
}
///////////////////////////////////////////////////////////////////////////////
static proc prepare_absprimdec(list primary)
{
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  list resu,tempo;
  string absotto;
  resu[size(primary) div 2]=list();
  for(int ab=1;ab<=size(primary) div 2;ab++)
  {
    absotto= absFactorize(primary[2*ab][1],77);
    tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering)));
    resu[ab]=tempo;
  }
  return(resu);
}
///////////////////////////////////////////////////////////////////////////////

static proc decompE(ideal I,list #)
"USAGE:  decompE(I); I ideal  (for primary decomposition)   (resp.
         decompE(I,1);        (for the associated primes of dimension of I) )
         decompE(I,2);        (for the minimal associated primes) )
         decompE(I,3);        (for the absolute primary decomposition) )
RETURN:  list = list of primary ideals and their associated primes
         (at even positions in the list)
         (resp. a list of the minimal associated primes)
         if I is unit ideal, returns empty list
NOTE:    Algorithm of Gianni/Trager/Zacharias
EXAMPLE: example decompE; shows an example
"
{
    return(decomp_i(1,I,#));
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 32003,(x,y,z),lp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal I = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   list pr= decompE(I);
   pr;
   testPrimary( pr, I);
}

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) )
         decomp(I,3);        (for the absolute primary decomposition) )
RETURN:  list = list of primary ideals and their associated primes
         (at even positions in the list)
         (resp. a list of the minimal associated primes)
         if I is unit ideal, returns list(ideal(1),ideal(1)) ( resp. list(ideal(1)))
EXAMPLE: example decomp; shows an example
"
{
    return(decomp_i(0,I,#));
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 32003,(x,y,z),lp;
   poly  p = z2+1;
   poly  q = z4+2;
   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
   list pr= decomp(i);
   pr;
   testPrimary( pr, i);
}

static proc decomp_i(int patchPrimaryDecomposition, ideal i,list #)
{
// if patchPrimaryDecomposition=1,  drop unit ideal in the decomposition,
// since the unit ideal it is not prime, otherwise take no special action.
// for other parameters see 'decomp' or 'decompE'
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  intvec op,@vv;
  list primary,indep,ltras;
  intvec @vh,isat,@w;
  int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi,abspri,ab,nn;
  ideal peek=i;
  ideal ser,tras;

  if(size(#)>0)
  {
    if (typeof(#[1])=="int")
    {
      if((#[1]==0)||(#[1]==1)||(#[1]==2)||(#[1]==3))
      {
        if (#[1]!=0) { @wr=#[1];}
        if(@wr==3){abspri=1;@wr=0;}
        if(size(#)>1)
        {
          seri=1;
          peek=#[2];
          ser=#[3];
        }
      }
    }
  }
  if(abspri)
  {
    list absprimary,abskeep,absprimarytmp,abskeeptmp;
  }
  //----------------------------------------------------------------
  //i is the zero-ideal ?
  //----------------------------------------------------------------
  if(size(i)==0)
  {
    primary=ideal(0),ideal(0);
    if (abspri) { return(prepare_absprimdec(primary));}
    return(primary);
  }
  intvec initialOp = option(get);
  def  @P = basering;
  int isS=attrib(i,"isSB");

  //----------------------------------------------------------------
  //i is homgeneous ?
  //----------------------------------------------------------------
  homo=homog(i);
  if(homo)
  {
    if(isS/*attrib(i,"isSB")*/!=1)
    {
      //ltras=mstd(i);
      tras=groebner(i);
    }
    else
    {
      tras=i;
    }
    ltras = tras,tras;
    attrib( ltras[1], "isSB", 1);

    if (size(ltras[1])>0)
    {
      if ( deg(lead(ltras[1]))==0 ) // we have the unit ideal.
      {
        option(set,initialOp);
        if (patchPrimaryDecomposition==1)
        {
          if (abspri) { return(prepare_absprimdec(list()));  }
          return( list() );
        }
        else
        {
          primary[1]=ideal(1);
          primary[2]=ideal(1);
          if (abspri) { return(prepare_absprimdec(primary));}
          return( primary );
        }
      }
    }

    tras=ltras[1];
    attrib(tras,"isSB",1);
    if((dim(tras)==0) && (!abspri))
    {
      primary[1]=ltras[2];
      primary[2]=maxideal(1);
      option(set,initialOp);
      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));
    }
    bigintvec @hilb=hilb(tras,1,@w);
    bigintvec keephilb=@hilb;
  }

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

  int lp=islp();

  def gnir=changeordTo(basering,"lp");
  setring gnir;
  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));
    }
    if(size(@j)==1)
    {
      if ( deg( lead(@j[1]) )==0 ) // we have the unit ideal.
      {
        setring @P;
        option(set,initialOp);
        if (patchPrimaryDecomposition==1)
        {
          return( list() );
        }
        else
        {
          return( list(ideal(1),ideal(1)) );
        }
      }
    }
  }
  option(set,op);
  if(seri==1)
  {
    ideal peek=fetch(@P,peek);
    attrib(peek,"isSB",1);
  }
  else
  {
    ideal peek=@j;
  }
  if((size(ser)==0)&&(!abspri))
  {
    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;
      option(set,initialOp);
      primary[1]=i;
      primary[2]=i;
      if (abspri) { return(prepare_absprimdec(primary));}
      return(primary);
    }
    if(size(fried)>0)
    {
      string newva;
      string newma;
      poly f;
      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)+",";
          fried[@n]=fried[@n]/leadcoef(fried[@n]);
          f=fried[@n]-lead(fried[@n]);
          @j=subst(@j,var(@k),-f);
        }
      }
      newva[size(newva)]=")";
      newma[size(newma)]=";";
      ring @deirf = create_ring(ring_list(gnir)[1],"("+newva,"lp","no_minpoly");
      execute("map @kappa=gnir,"+newma);
      ideal @j= @kappa(@j);
      @j=std(@j);

      list pr=decomp_i(patchPrimaryDecomposition,@j,#);
      if (size(pr)==0)
      {
        setring @P;
        option(set,initialOp);
        if (abspri) { return(prepare_absprimdec(list()));}
        return(list());
      }

      setring gnir;
      list pr=imap(@deirf,pr);
      for(@k=1;@k<=size(pr);@k++)
      {
        @j=pr[@k]+fried;
        pr[@k]=@j;
      }
      setring @P;
      option(set,initialOp);
      primary=imap(gnir,pr);
      if (abspri) { return(prepare_absprimdec(primary));}
      return(primary);
    }
  }
  //----------------------------------------------------------------
  //j is the ring
  //----------------------------------------------------------------

  if (dim(@j)==-1)
  {
    setring @P;
    option(set,initialOp);
    primary=ideal(1),ideal(1);
    if (abspri) { return(prepare_absprimdec(primary));}
    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;
    option(set,initialOp);
    primary=fetch(gnir,gprimary);

    if (abspri) { return(prepare_absprimdec(primary));}
    return(primary);
  }

 //------------------------------------------------------------------
 //the zero-dimensional case
 //------------------------------------------------------------------
  if (dim(@j)==0)
  {
    op=option(get);
    option(redSB);
    list gprimary= zero_decomp(@j,ser,@wr);

    setring @P;
    primary=fetch(gnir,gprimary);

    if(size(ser)>0)
    {
      primary=cleanPrimary(primary);
    }
    if(abspri)
    {
      setring gnir;
      list primary=imap(@P,primary);
      list resu,tempo;
      string absotto;
      map sigma,invsigma;
      ideal II,jmap;
      nn=nvars(basering);
      for(ab=1;ab<=size(primary) div 2;ab++)
      {
        II=primary[2*ab];
        attrib(II,"isSB",1);
        if(deg(II[1])==vdim(II))
        {
          absotto= absFactorize(primary[2*ab][1],77);
          tempo=
            primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering)));
        }
        else
        {
          invsigma=basering,maxideal(1);
          jmap=randomLast(50);
          sigma=basering,jmap;
          jmap[nn]=2*var(nn)-jmap[nn];
          invsigma=basering,jmap;
          II=groebner(sigma(II));
          absotto = absFactorize(II[1],77);
          II=var(nn);
          tempo= primary[2*ab-1],primary[2*ab],absotto,string(invsigma(II));
        }
        resu[ab]=tempo;
      }
      primary=resu;
      setring @P;
      primary=imap(gnir,primary);
    }
    option(set,initialOp);
    return(primary);
  }

  poly @gs,@gh,@p;
  string @va;
  def 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")&&(size(ring_list(@P)[3])==2)) // weighted ordering
  {
    list gnir_l=ring_list(gnir);
    list @P_l=ring_list(@P);
    gnir_l[3]=@P_l[3]; // ord
    def @Phelp=ring(gnir_l);
    kill gnir_l,@P_l;
    setring @Phelp;
  }
  else
  {
    def @Phelp=changeordTo(gnir,"dp");
    setring @Phelp;
  }

  if(homo==1)
  {
    if(((ordstr(@P)[3]=="d")||(ordstr(@P)[1]=="d")||(ordstr(@P)[1]=="w")
       ||(ordstr(@P)[3]=="w"))&&(size(ring_list(@P)[3])==2))
    {
      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;
  ideal @Ptest=1;
  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
    {
      def gnir1=ring(ring_list(basering));
      setring gnir1;
      ideal @j=fetch(gnir,@j);
      attrib(@j,"isSB",1);
      ideal ser=fetch(gnir,ser);
    }
    else
    {
      @va=string(maxideal(1));
      if(@m==1)
      {
        @j=fetch(@P,i);
      }
      ring gnir1 = create_ring(ring_list(basering)[1],"("+indep[@m][1]+")","("+indep[@m][2]+")","no_minpoly");
      execute("map phi=gnir,"+@va+";");
      op=option(get);
      option(redSB);
      ideal @j=groebner(phi(@j));
      ideal ser=phi(ser);

      option(set,op);
    }
    if((deg(@j[1])==0)||(dim(@j)<jdim))
    {
      setring gnir;
      kill gnir1;
      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],"lp");

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

    ideal @jj=lead(@j);               //!! vorn vereinbaren
    setring quotring;

    ideal @jj=imap(gnir1,@jj);
    @vv=clearSBNeu(@jj,fett);  //!! vorn vereinbaren
    setring gnir1;
    @k=size(@j);
    for (lauf=1;lauf<=@k;lauf++)
    {
      if(@vv[lauf]==1)
      {
        @j[lauf]=0;
      }
    }
    @j=simplify(@j,2);
    setring 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

    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);
      if(abspri)
      {
        ideal II;
        ideal jmap;
        map sigma;
        nn=nvars(basering);
        map invsigma=basering,maxideal(1);
        for(ab=1;ab<=size(uprimary) div 2;ab++)
        {
          II=uprimary[2*ab];
          attrib(II,"isSB",1);
          if(deg(II[1])!=vdim(II))
          {
            jmap=randomLast(50);
            sigma=basering,jmap;
            jmap[nn]=2*var(nn)-jmap[nn];
            invsigma=basering,jmap;
            II=groebner(sigma(II));
          }
          absprimarytmp[ab]= absFactorize(II[1],77);
          II=var(nn);
          abskeeptmp[ab]=string(invsigma(II));
          invsigma=basering,maxideal(1);
        }
      }
      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(quotring,uprimary);
    lsau=imap(quotring,saturn);
    @h=imap(quotring,@h);

    kill quotring;
    def quotring;

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

    for(@n1=1;@n1<=size(collectprimary) div 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];
        if(abspri)
        {
          absprimary[@n2 div 2]=absprimarytmp[@n1];
          abskeep[@n2 div 2]=abskeeptmp[@n1];
        }
      }
    }

    //here the intersection with the polynomialring
    //mentioned above is really computed
    for(@n=@n3 div 2+1;@n<=@n2 div 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)
      {
        if(defined(@LL)){kill @LL;}
        list @LL=minSat(jwork,@h);
        @Ptest=intersect(@Ptest,@LL[1]);
        @q=@LL[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,6);
        @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;
    quprimary[1]=ideal(1);
    quprimary[2]=ideal(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)
      {
        //HIER STATT DURCHSCHNITT SATURIEREN!
        ideal htest=@Ptest;
      }
      else
      {
        ideal htest=hquprimary[2];

        for (@n1=2;@n1<=size(hquprimary) div 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,5))!=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
        {
          def gnir1=ring(ring_list(basering));
          setring gnir1;
          ideal @j=fetch(gnir,jkeep);
          attrib(@j,"isSB",1);
        }
        else
        {
          @va=string(maxideal(1));
          ring gnir1 = create_ring(ring_list(basering)[1],"("+restindep[@m][1]+")","("+restindep[@m][2]+")","no_minpoly");
          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..]
        //------------------------------------------------------------------
        if (defined(quotring)==voice) {kill quotring;}
        def quotring=prepareQuotientring(nvars(basering)-restindep[@m][3],"lp");

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

        setring 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);
        if(abspri)
        {
          ideal II;
          ideal jmap;
          map sigma;
          nn=nvars(basering);
          map invsigma=basering,maxideal(1);
          for(ab=1;ab<=size(uprimary) div 2;ab++)
          {
            II=uprimary[2*ab];
            attrib(II,"isSB",1);
            if(deg(II[1])!=vdim(II))
            {
              jmap=randomLast(50);
              sigma=basering,jmap;
              jmap[nn]=2*var(nn)-jmap[nn];
              invsigma=basering,jmap;
              II=groebner(sigma(II));
            }
            absprimarytmp[ab]= absFactorize(II[1],77);
            II=var(nn);
            abskeeptmp[ab]=string(invsigma(II));
            invsigma=basering,maxideal(1);
          }
        }
        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(quotring,uprimary);
        lsau=imap(quotring,saturn);
        @h=imap(quotring,@h);

        kill quotring;

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

        for(@n1=1;@n1<=size(collectprimary) div 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];
            if(abspri)
            {
              absprimary[@n2 div 2]=absprimarytmp[@n1];
              abskeep[@n2 div 2]=abskeeptmp[@n1];
            }
          }
        }


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

        for(@n=@n3 div 2+1;@n<=@n2 div 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 div 2+1;@n<=@n2 div 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(abspri)
      {
        list resu,tempo;
        for(ab=1;ab<=size(quprimary) div 2;ab++)
        {
          if (deg(quprimary[2*ab][1])!=0)
          {
            tempo=quprimary[2*ab-1],quprimary[2*ab],
                         absprimary[ab],abskeep[ab];
            resu[ab]=tempo;
          }
        }
        quprimary=resu;
        @wr=3;
      }
      if(size(reduce(ser,peek,5))!=0)
      {
        htprimary=decomp_i(patchPrimaryDecomposition, @j,@wr,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];
        }
      }
    }
  }
  else
  {
    if(abspri)
    {
      list resu,tempo;
      for(ab=1;ab<=size(quprimary) div 2;ab++)
      {
        tempo=quprimary[2*ab-1],quprimary[2*ab],
                   absprimary[ab],abskeep[ab];
        resu[ab]=tempo;
      }
      quprimary=resu;
    }
  }
  //---------------------------------------------------------------------------
  //back to the ring we started with
  //the final result: primary
  //---------------------------------------------------------------------------
  setring @P;
  option(set,initialOp);
  primary=imap(gnir,quprimary);
  if(!abspri)
  {
    primary=cleanPrimary(primary);
  }
  if (size(primary)>0)
  {
    if (abspri && (typeof(primary[1][1])=="poly"))
    { return(prepare_absprimdec(primary));}
  }
  return(primary);
}




///////////////////////////////////////////////////////////////////////////////
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
{
   ASSUME(1, hasFieldCoefficient(basering) );
   ASSUME(1, not isQuotientRing(basering) ) ;
   ASSUME(1, hasGlobalOrdering(basering) ) ;

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

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 considered
        as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k))
 EXAMPLE: example sep; shows an example
{
   ASSUME(0, hasFieldCoefficient(basering) );
   ASSUME(0, not isQuotientRing(basering) ) ;
   ASSUME(0, hasGlobalOrdering(basering) ) ;

   def R=basering;
   int k;
   if(size(#)>0){k=#[1];}


   poly h=gcd(f,diff(f,var(i)));
   if((reduce(f,std(h),5)!=0)||(reduce(diff(f,var(i)),std(h),5)!=0))
   {
      ERROR("FEHLER IN GCD");
   }
   poly g1=f;
   if (h!=1)
   {
     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"
{
   ASSUME(0, hasFieldCoefficient(basering) );
   ASSUME(0, not isQuotientRing(basering) ) ;
   ASSUME(0, hasGlobalOrdering(basering) ) ;

   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 d=vdim(I);
   int i,k;
   list l;
   if(((p==0)||(p>d))&&(d==deg(I[1])))
   {
     intvec e=leadexp(I[1]);
     for(i=1;i<=nvars(basering);i++)
     {
       if(e[i]!=0) break;
     }
     I[1]=sep(I[1],i)[1];
     return(interred(I));
   }
   intvec op=option(get);

   option(redSB);
   ASSUME(1, dim(I)==0);
   if (hasQQCoefficient(basering) && (nvars(basering)*vdim(I)>300))
   {
     ideal F=modFinduni(I); //F[i] generates I intersected with K[var(i)]
   }
   else
   {
     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)) //the separable case
   {
     intvec save=option(get);
     option(redSB);
     I=std(I+F);
     option(set,save);
     return(I);
   }
   //I=simplify(I,1);

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

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

   list l2;
   for (int zz = 1; zz <= m; zz++)
   {
    l2[size(l2)+1] = "@y("+string(zz)+")";
   }
   ring @S = create_ring(p, "("+string(l2)+","+varstr(R)+","+parstr(R)+")", "dp");

   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*var(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 algeRad(ideal i)
{
   ASSUME(0, hasFieldCoefficient(basering) );
   ASSUME(0, not isQuotientRing(basering) ) ;
   ASSUME(0, hasGlobalOrdering(basering) ) ;

//reduces radical computation 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;
   ideal ra=radical(i);
   gnir="ring RS="+string(char(R))+",("+varstr(RH)
                +"),(dp("+string(n)+"),lp);";
   execute(gnir);
   ideal K=imap(RH,ra);
   K=std(K);
   if (ncols(K)>1) { K=K[2..size(K)];}
   setring R;
   ideal ra=imap(RS,K);

   return( ra );
}


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

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
"
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  if(attrib(basering,"global")!=1)
  {
    ERROR(
      "// Not implemented for this ordering, please change to global ordering."
    );
  }

  if ( size(i)==0 ) { return(ideal(0)); }
  if(hasAlgExtensionCoefficient(basering))
  {
    return(algeRad(i));
  }
  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);

  if ( dim(m[1])<0 ) { return(ideal(1)); }

  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],5))==0)
    {
      return(I);
    }

    I2=sat(I,radI1);

    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)
"USAGE:   Ann(M);  M module
RETURN:  ideal, the annihilator of coker(M)
NOTE:    The output is the ideal of all elements a of the basering R such that
         a * R^m is contained in M  (m=number of rows of M).
EXAMPLE: example Ann; shows an example
"
{

  M=prune(M);  //to obtain a small embedding
  ideal ann=quotient1(M,freemodule(nrows(M)));
  return(ann);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 0,(x,y,z),lp;
   module M = x2-y2,z3;
   Ann(M);
   M = [1,x2],[y,x];
   Ann(M);
   qring Q=std(xy-1);
   module M=imap(r,M);
   Ann(M);
}

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

//computes the equidimensional part of the ideal i of codimension e
static proc int_ass_primary_e(ideal i, int e)
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  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
// !! borrowed correct code from 'ehv.lib::AnnExtEHV' by Kai Dehmann !! duplicate code!! (jk)

static proc AnnExt_R(int n,list re)
"USAGE:   AnnExt_R(n,re); n integer, re resolution
RETURN:  ideal, the annihilator of Ext^n(R/I,R) with given
         resolution re of I
"
{

  if(n < 0)
  {
    return(ideal(1));
  }
  int l = size(re);

  if(n < l)
  {
    matrix f = transpose(re[n+1]);
    if(n == 0)
    {
      matrix g = matrix(0,1,ncols(f));
    }
    else
    {
      matrix g = transpose(re[n]);
    }
    module k = syz(f);
    return(quotient1(g,k));
  }

  if(n == l)
  {
    return(Ann(transpose(re[n])));
  }

  return(ideal(1));
}
///////////////////////////////////////////////////////////////////////////////

static proc analyze(list pr)
{
   ASSUME(1, hasFieldCoefficient(basering) );
   ASSUME(1, not isQuotientRing(basering) ) ;
   ASSUME(1, hasGlobalOrdering(basering) ) ;

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

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

///////////////////////////////////////////////////////////////////////////////
//
//                  Shimoyama-Yokoyama
//
///////////////////////////////////////////////////////////////////////////////
static proc simplifyIdeal(ideal i)
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  def r=basering;

  ideal iwork=i;
  ideal imap2=maxideal(1);

  int j,k;
  map phi;
  poly p;
  ideal imap1=maxideal(1);
  // first try: very simple substitutions
  intvec tested=0:nvars(r);
  for(j=1;j<=nvars(r);j++)
  {
    for(k=1;k<=ncols(i);k++)
    {
      if(deg(iwork[k]/var(j))==0)
      {
        p=-1/leadcoef(iwork[k]/var(j))*iwork[k];
        if(size(p)<=2)
        {
          tested[j]=1;
          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;
        }
      }
    }
  }
  // second try: substitutions not so simple
  for(j=1;j<=nvars(r);j++)
  {
    if (tested[j]==0)
    {
      for(k=1;k<=ncols(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--)
  {
    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_i (int patchPrimaryDecomposition, ideal PS, int cho)
{
// if patchPrimaryDecomposition=1,  drop unit ideal in the decomposition,
// since the unit ideal it is not prime, otherwise take no special action.
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  ASSUME(1, not isQuotientRing(basering) ) ;

  if((cho<0) or (cho>1))
  {
    ERROR("<int> must be 0 or 1");
  }
  intvec saveopt=option(get);
  option(notWarnSB);
  list L;
  if(cho==0)
  {
    L=min_ass_prim_charsets0_i(patchPrimaryDecomposition,PS);
  }
  else
  {
    L=min_ass_prim_charsets1_i(patchPrimaryDecomposition,PS);
  }
  option(set,saveopt);
  return(L);
}
///////////////////////////////////////////////////////
// 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_i (int patchPrimaryDecomposition, ideal PS)
{
// if patchPrimaryDecomposition=1,  drop unit ideal in the decomposition,
// since the unit ideal it is not prime, otherwise take no special action.
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  if (size(PS)==0) { return( list(ideal(0))); }
  intvec op;
  matrix m=char_series(PS);  // We compute an irreducible
                             // characteristic series
  if ((nrows(m)==1)
  && (ncols(m)==1)
  && (m[1,1]==1)) // in case of an empty series: min_ass_prim_charsets1
  {
    return (min_ass_prim_charsets1_i(patchPrimaryDecomposition,PS));
  }
  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--)
    {
      I=I+ini_mod(PHI[i][j]);
    }
    JS=std(PHI[i]);
    sizeJS=size(JS);
    for(j=size(I);j>0;j--)
    {
      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);
    }
  }
  if(size(PSI)==1)
  {
    if (idealsEqual( PSI[1], ideal(1) ))
    {
      if (patchPrimaryDecomposition==1)
      {
        return( list() );
      }
      else
      {
        return( list(ideal(1)) );
      }
    }
  }
  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_i (int patchPrimaryDecomposition, ideal PS)
{
// if patchPrimaryDecomposition=1,  drop unit ideal in the decomposition,
// since the unit ideal it is not prime, otherwise take no special action.
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  if (size(PS)==0) { return( list(ideal(0))); }
  intvec op;
  def oldring=basering;
  string n=system("neworder",PS);
  ring r = create_ring(ring_list(oldring)[1],"("+n+")","dp","no_minpoly");
  ideal PS=imap(oldring,PS);
  matrix m=char_series(PS);  // We compute an irreducible
                             // characteristic series
                             // this series may be empty (1x1: 1)
  int i,j,k,cnt;
  while ((cnt<nvars(oldring))
  && (nrows(m)==1)
  && (ncols(m)==1)
  && (m[1,1]==1)) // in case of an empty series: permute the variables
  {
    cnt++;
    n=string(var(nvars(oldring)));
    for(i=1;i<nvars(oldring);i++) { n=n+","+string(var(i)); }
    kill r;
    ring r = create_ring(ring_list(oldring)[1],"("+n+")","dp","no_minpoly");
    ideal PS=imap(oldring,PS);
    matrix m=char_series(PS);
  }
  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]=simplify(ideal(m[i,1..ncols(m)]),2);
    I=0;
    for(j=ncols(PHI[i]);j>0;j--)
    {
      I=I,ini_mod(PHI[i][j]);
    }
    if (ncols(I)>1)
    {
      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--)
    {
      I=I+ITPHI[i][j];
    }
    JS=std(PHI[i]);
    sizeJS=size(JS);
    for(j=size(I);j>0;j--)
    {
      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);
    }
  }
  if(size(PSI)==1)
  {
    if (idealsEqual( PSI[1], ideal(1) ))
    {
      if (patchPrimaryDecomposition==1)
      {
        return( list() );
      }
      else
      {
        return( list(ideal(1)) );
      }
    }
  }

  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_i(int patchPrimaryDecomposition, ideal I, int choose)
{
// if patchPrimaryDecomposition=1,  drop unit ideal in the decomposition,
// since the unit ideal it is not prime, otherwise take no special action.
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  if((choose<0) or (choose>3))
  {
    ERROR("ERROR: <int> must be 0 or 1 or 2 or 3");
  }
  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))
  {
          ASSUME(1, ncols(SI)==1);
          if (patchPrimaryDecomposition==1)
          {
             return( list() );
          }
          else
          {
               return( list(list(ideal(1),ideal(1))) );
          }
  }
  intvec save=option(get);
  option(notWarnSB);
  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
        {
          option(set,save);
          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_i(patchPrimaryDecomposition,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(patchPrimaryDecomposition,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));
        }
      }
    }
  }
  option(set,save);
}

//////////////////////////////////////////////////////////////////////////
// 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_i(int patchPrimaryDecomposition, ideal I, ideal SI, int choo)
{
// if patchPrimaryDecomposition=1,  drop the unit ideal in the decomposition,
// since the unit ideal it is not prime, otherwise take no special action.
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  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(patchPrimaryDecomposition,I,choo);
  }
  else
  {
    if(choo==2)
    {
      L=minAssPrimes_i(patchPrimaryDecomposition,I);
    }
    else
    {
      L=minAssPrimes_i(patchPrimaryDecomposition,I,1);
    }
    for(int i=size(L);i>=1;i--)
    {
      L[i]=std(L[i]);
    }
  }
  return (pseudo_prim_dec_i_i(patchPrimaryDecomposition,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_i (int patchPrimaryDecomposition, ideal SI,list V6, int choo)
{
// if patchPrimaryDecomposition=1,  drop the unit ideal in the decomposition,
// since the unit ideal it is not prime, otherwise take no special action.
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  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_i(patchPrimaryDecomposition,SP,choo);  // a list of SB
        }
        else
        {
          if(choo==2)
          {
            m=minAssPrimes_i(patchPrimaryDecomposition,SP);
          }
          else
          {
            m=minAssPrimes_i(patchPrimaryDecomposition,SP,1);
          }
          for(j=size(m);j>=1;j--)
            {
              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_i(patchPrimaryDecomposition,SI,L));
}


////////////////////////////////////////////////////////////////
// proc pseudo_prim_dec_i_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_i (int patchPrimaryDecomposition, ideal SI, list L)
{
// if patchPrimaryDecomposition=1,  drop the unit ideal in the decomposition,
// since the unit ideal it is not prime, otherwise take no special action.
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;
  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+2);
    // 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)
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  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)
      ring RAU = create_ring(ring_list(basering)[1],"("+A+U,"(dp("+string(a)+"),dp)","no_minpoly");
      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)]=")";
      ring Rloc = create_ring("("+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)
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  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 separately 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)
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  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 separately 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)
{
   ASSUME(1, hasFieldCoefficient(basering) );
   ASSUME(1, not isQuotientRing(basering) ) ;
   ASSUME(1, hasGlobalOrdering(basering) ) ;

   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) div 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=size(l);i>0;i--)
   {
      re[2*i-1]=l[i][1];
      re[2*i]=l[i][2];
   }
   return(re);
}

///////////////////////////////////////////////////////////////////////////////
//
//     The main procedures
//
///////////////////////////////////////////////////////////////////////////////
proc primdecGTZE(ideal I, list #)
"USAGE:   primdecGTZE(I); i ideal
RETURN:  a list pr of primary ideals and their associated primes for a proper ideal, and an empty list for the unit ideal.
@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!)
         - For local orderings, the result is considered in the localization
           of the polynomial ring, not in the power series ring
         - For local and mixed orderings, the decomposition in the
           corresponding global ring is returned if the string 'global'
           is specified as second argument
EXAMPLE: example primdecGTZE; shows an example
"
{
  return (primdecGTZ_i(1,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 = primdecGTZE(I);
   pr;
   ideal J = 1;
   list prempty = primdecGTZE(J);
   prempty;
}

proc primdecGTZ(ideal I, list #)
"USAGE:   primdecGTZ(I); I ideal
RETURN:  a list pr of primary ideals and their associated primes for a proper ideal I, otherwise pr = list( list( ideal(1), ideal(1) )
@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!)
         - For local orderings, the result is considered in the localization
           of the polynomial ring, not in the power series ring
         - For local and mixed orderings, the decomposition in the
           corresponding global ring is returned if the string 'global'
           is specified as second argument
EXAMPLE: example primdecGTZ; shows an example
"
{
  return (primdecGTZ_i(0, 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;
}

static proc primdecGTZ_i(int patchPrimaryDecomposition,ideal i, list #)
{
// if parameter patchPrimaryDecomposition=1,  drop the unit ideal in the decomposition,
//           since the unit ideal it is not prime, otherwise take no special action.
// For other parameters see 'primdecGTZ' or 'primdecGTZE'.
   ASSUME(0, hasFieldCoefficient(basering) );
   ASSUME(0, not isQuotientRing(basering) ) ;
   if(size(#)>0)
   {
      int keep_comp=1;
   }
   if(attrib(basering,"global")!=1)
   {
// algorithms only work in global case!
// pass to appropriate global ring
      def r=basering;
      def s=changeord(list(list("dp",1:nvars(basering))));
      setring s;
      ideal i=imap(r,i);
// decompose and go back
      list li=primdecGTZ_i(patchPrimaryDecomposition,i,#);
      int sizeli = size(li);
      setring r;
      if (sizeli==0)
      {
          return ( list() );
      }
      list li=imap(s,li);
// clean up
      if(!defined(keep_comp))
      {
         for(int k=size(li);k>=1;k--)
         {
            if(mindeg(std(lead(li[k][2]))[1])==0)
            {
// 1 contained in ideal, i.e. component does not meet origin in local ordering
               li=delete(li,k);
            }
         }
      }
      return(li);
   }

   if(minpoly!=0)
   {
      return(algeDeco_i(patchPrimaryDecomposition,i,0));
   }
   return(convList(decomp_i(patchPrimaryDecomposition,i,#)));
}

///////////////////////////////////////////////////////////////////////////////
proc absPrimdecGTZE(ideal I, list #)
"USAGE:   absPrimdecGTZE(I); I ideal
ASSUME:  Ground field has characteristic 0.
RETURN:  a ring containing two lists: @code{absolute_primes}, the absolute
         prime components of I, and @code{primary_decomp}, the output of
         @code{primdecGTZ(I)}. Will fail for unit ideal.
         The list absolute_primes has to be interpreted as follows:
         each entry describes a class of conjugated absolute primes,
@format
   absolute_primes[i][1]   the absolute prime component,
   absolute_primes[i][2]   the number of conjugates.
@end format
         The first entry of @code{absolute_primes[i][1]} is the minimal
         polynomial of a minimal finite field extension over which the
         absolute prime component is defined.
         For local orderings, the result is considered in the localization
         of the polynomial ring, not in the power series ring.
         For local and mixed orderings, the decomposition in the
         corresponding global ring is returned if the string 'global'
         is specified as second argument
NOTE:    Algorithm of Gianni/Trager/Zacharias combined with the
         @code{absFactorize} command.
SEE ALSO: primdecGTZ; absFactorize
EXAMPLE: example absPrimdecGTZE; shows an example
"
{
     return(absPrimdecGTZ_i(1,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;
   def S = absPrimdecGTZE(I);
   setring S;
   absolute_primes;
}

proc absPrimdecGTZ(ideal I, list #)
"USAGE:   absPrimdecGTZ(I); I ideal
ASSUME:  Ground field has characteristic 0.
RETURN:  a ring containing two lists: @code{absolute_primes}, the absolute
         prime components of I, and @code{primary_decomp}, the output of
         @code{primdecGTZ(I)}.
         The list absolute_primes has to be interpreted as follows:
         each entry describes a class of conjugated absolute primes,
@format
   absolute_primes[i][1]   the absolute prime component,
   absolute_primes[i][2]   the number of conjugates.
@end format
         The first entry of @code{absolute_primes[i][1]} is the minimal
         polynomial of a minimal finite field extension over which the
         absolute prime component is defined.
         For local orderings, the result is considered in the localization
         of the polynomial ring, not in the power series ring.
         For local and mixed orderings, the decomposition in the
         corresponding global ring is returned if the string 'global'
         is specified as second argument
NOTE:    Algorithm of Gianni/Trager/Zacharias combined with the
         @code{absFactorize} command.
SEE ALSO: primdecGTZ; absFactorize
EXAMPLE: example absPrimdecGTZ; shows an example
"
{

    return(absPrimdecGTZ_i(0,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;
   def S = absPrimdecGTZ(i);
   setring S;
   absolute_primes;
}



static proc absPrimdecGTZ_i(int @patchPrimaryDecomposition, ideal @I, list #)
{
// if parameter patchPrimaryDecomposition=1,  drop the unit ideal in the decomposition,
//           since the unit ideal it is not prime, otherwise take no special action.
// For other parameters see 'absPrimdecGTZ' or 'absPrimdecGTZE'.
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  if (char(basering) != 0)
  {
    ERROR("primdec.lib::absPrimdecGTZ is only implemented for "+
           +"characteristic 0");
  }

  if(size(#)>0)
  {
     int @keep_comp=1;
  }

  if(attrib(basering,"global")!=1)
  {
// algorithm automatically passes to the global case
// hence prepare to go back to an appropriate new ring
      def @r=basering;
      ideal @max_of_r=maxideal(1);
      def @s=changeord(list(list("dp",1:nvars(basering))));
      setring @s;
      def @I=imap(@r,@I);
      def @S=absPrimdecGTZ_i(@patchPrimaryDecomposition,@I);
      setring @S;
      ring @r1=char(basering),var(nvars(@r)+1),dp;
      def @rS=@r+@r1;
// move objects to appropriate ring and clean up
      setring @rS;
      def @max_of_r=imap(@r,@max_of_r);
      attrib(@max_of_r,"isSB",1);
      def absolute_primes=imap(@S,absolute_primes);
      def primary_decomp=imap(@S,primary_decomp);
      if(!defined(@keep_comp))
      {
         ideal @tempid;
         for(int @k=size(absolute_primes);@k>=1;@k--)
         {
            @tempid=absolute_primes[@k][1];
            @tempid[1]=0;                  // ignore minimal polynomial
            if(size(reduce(lead(@tempid),@max_of_r,5))!=0)
            {
// 1 contained in ideal, i.e. component does not meet origin in local ordering
               absolute_primes=delete(absolute_primes,@k);
            }
         }
         for(@k=size(primary_decomp);@k>=1;@k--)
         {
            if(mindeg(std(lead(primary_decomp[@k][2]))[1])==0)
            {
// 1 contained in ideal, i.e. component does not meet origin in local ordering
               primary_decomp=delete(primary_decomp,@k);
            }
         }
         kill @tempid;
      }
      export(primary_decomp);
      export(absolute_primes);
      return(@rS);
  }
  if(minpoly!=0)
  {
    //return(algeDeco_i(patchPrimaryDecomposition,I,0));
    ERROR(
      "Not implemented yet for algebraic extensions.Simulate the ring extension by adding the minpoly to the ideal"
    );
  }
  def @R=basering;
  int @n=nvars(@R);
  list @L=decomp_i(@patchPrimaryDecomposition,@I,3);
  if (@patchPrimaryDecomposition && size(@L)==0 )
  {
     "// will not handle case with unit ideal";
  }
  string @newvar=@L[1][3];
  int @k=find(@newvar,",",find(@newvar,",")+1);
  @newvar=@newvar[@k+1..size(@newvar)];
  list @lR=ring_list(@R);
  int @i,@de,@ii;
  intvec @vv=1:@n;

  list @orst;
  @orst[1]=list("dp",@vv);
  @orst[2]=list("dp",intvec(1));
  @orst[3]=list("C",0);
  @lR[3]=@orst;
  @lR[2][@n+1] = @newvar;
  def @Rz = ring(@lR);
  setring @Rz;
  list @L=imap(@R,@L);
  list absolute_primes,primary_decomp;

  ideal @I,@M,@N,@K,@J;

  @M=maxideal(1);
  @N=maxideal(1);
  poly @p,@q,@f,@g;
  map @phi,@psi;
  string @tvar;
  for(@i=1;@i<=size(@L);@i++)
  {
    @tvar=@L[@i][4];
    @ii=find(@tvar,"+");
    while(@ii)
    {
      @tvar=@tvar[@ii+1..size(@tvar)];
      @ii=find(@tvar,"+");
    }
    for(@ii=1;@ii<=nvars(basering);@ii++)
    {
      if(@tvar==string(var(@ii))) break;
    }
    @I=@L[@i][2];
    execute("@K="+@L[@i][3]+";");
    @p=@K[1];
    @q=@K[2];
    execute("@f="+@L[@i][4]+";");
    @g=2*var(@ii)-@f;
    @M[@ii]=@f;
    @N[@ii]=@g;
    @de=deg(@p);
    @psi=@Rz,@M;
    @phi=@Rz,@N;
    @I=@phi(@I),@p,@q;
    @I=std(@I);
    @J=@psi(@I);
    @J=sat(@J,@f);
    if(deg(@J[1])>0)
    {
       absolute_primes[size(absolute_primes)+1]=list(@J,@de);
    }
    else
    {
       @J=@psi(@I);
       absolute_primes[size(absolute_primes)+1]=list(@J,@de);
    }
    primary_decomp[@i]=list(@L[@i][1],@L[@i][2]);
  }
  export(primary_decomp);
  export(absolute_primes);
  setring @R;
  dbprint( printlevel-voice+4,"
// 'absPrimdecGTZ' created a ring, in which two lists absolute_primes (the
// absolute prime components) and primary_decomp (the primary and prime
// components over the current basering) are stored.
// To access the list of absolute prime components, type (if the name S was
// assigned to the return value):
        setring S; absolute_primes; ");

  return(@Rz);
}


///////////////////////////////////////////////////////////////////////////////
proc primdecSYE(ideal I, list #)
"USAGE:   primdecSYE(I, c); I ideal, c int (optional)
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
If I is the unit ideal returns an empty list.
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 (default),
   if c=2,  minAssGTZ is used,
   if c=3,  minAssGTZ and facstd are used.
@end format
         For local orderings, the result is considered in the localization
         of the polynomial ring, not in the power series ring.
         For local and mixed orderings, the decomposition in the
         corresponding global ring is returned if the string 'global'
         is specified as third argument
EXAMPLE: example primdecSY; shows an example
"
{
     return (primdecSY_i(1,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 = primdecSYE(I);
   pr;
   ideal J = x;
   list prUnit = primdecSYE(J);
   prUnit;
}

proc primdecSY( ideal I, list #)
"USAGE:   primdecSY(I, c); I ideal, c int (optional)
RETURN:  a list pr of primary ideals and their associated primes for proper ideal I, otherwise pr[1] is list( ideal(1),ideal(1) )'
@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 (default),
   if c=2,  minAssGTZ is used,
   if c=3,  minAssGTZ and facstd are used.
@end format
         For local orderings, the result is considered in the localization
         of the polynomial ring, not in the power series ring.
         For local and mixed orderings, the decomposition in the
         corresponding global ring is returned if the string 'global'
         is specified as third argument
EXAMPLE: example primdecSY; shows an example
"
{
    return (primdecSY_i(0,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 = primdecSY(i);
   pr;
}


static proc primdecSY_i(int patchPrimaryDecomposition, ideal i, list #)
{
//           if patchPrimaryDecomposition=1,  drop the unit ideal in the decomposition,
//           since the unit ideal it is not prime, otherwise take no special action.
//           For other paremetes see 'primdecSY' or 'primdecSYE'
   ASSUME(0, hasFieldCoefficient(basering) );
   ASSUME(0, not isQuotientRing(basering) ) ;
   if(size(#)>1)
   {
      int keep_comp=1;
   }
   if(attrib(basering,"global")!=1)
   {
// algorithms only work in global case!
// pass to appropriate global ring
      def r=basering;
      def s=changeord(list(list("dp",1:nvars(basering))));
      setring s;
      ideal i=imap(r,i);
// decompose and go back
      list li=primdecSY_i(patchPrimaryDecomposition,i);
      int sizeli = size(li);
      setring r;
      if (sizeli==0)  { return ( list() ); }
      def li=imap(s,li);
// clean up
      if(!defined(keep_comp))
      {
        for(int k=size(li);k>=1;k--)
        {
          if(mindeg(std(lead(li[k][2]))[1])==0)
          {
// 1 contained in ideal, i.e. component does not meet origin in local ordering
            li=delete(li,k);
          }
        }
      }
      return(li);
   }
   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(patchPrimaryDecomposition,i,1));
   }
   if (size(#)!=0)
   { return(prim_dec_i(patchPrimaryDecomposition,i,#[1])); }
   else
   { return(prim_dec_i(patchPrimaryDecomposition,i,1)); }
}



proc minAssGTZE(ideal I,list #)
"USAGE:    minAssGTZE(I[, l]); I ideal, l list (optional)
   @* Optional parameters in list l (can be entered in any order):
   @* 0, \"facstd\" -> uses facstd to first decompose the ideal (default)
   @* 1, \"noFacstd\" -> does not use facstd
   @* \"GTZ\" -> the original algorithm by Gianni, Trager and Zacharias is used
   @* \"SL\" -> GTZ algorithm with modificiations by Laplagne is used (default)

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
         - For local orderings, the result is considered in the localization
           of the polynomial ring, not in the power series ring
         - For local and mixed orderings, the decomposition in the
           corresponding global ring is returned if the string 'global'
           is specified as second argument
EXAMPLE: example minAssGTZE; shows an example
"
{
  list result = minAssGTZ_i(1,I,#);
  return(result);
}
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 = minAssGTZE(I);
   pr;
   ideal J = 1;
   list prempty = minAssGTZE(J);
   prempty;

}


proc minAssGTZ(ideal I, list #)
"USAGE:    minAssGTZ(I[, l]); I ideal, l list (optional)
   @* Optional parameters in list l (can be entered in any order):
   @* 0, \"facstd\" -> uses facstd to first decompose the ideal (default)
   @* 1, \"noFacstd\" -> does not use facstd
   @* \"GTZ\" -> the original algorithm by Gianni, Trager and Zacharias is used
   @* \"SL\" -> GTZ algorithm with modificiations by Laplagne is used (default)

RETURN:  a list, the minimal associated prime ideals of proper ideal I, otherwise ideal(1)
NOTE:    - Designed for characteristic 0, works also in char k > 0 based
           on an algorithm of Yokoyama
         - For local orderings, the result is considered in the localization
           of the polynomial ring, not in the power series ring
         - For local and mixed orderings, the decomposition in the
           corresponding global ring is returned if the string 'global'
           is specified as second argument
EXAMPLE: example minAssGTZ; shows an example
"
{
  list result = minAssGTZ_i(0,I,#);
  return(result);
}
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;
}


///////////////////////////////////////////////////////////////////////////////
static proc minAssGTZ_i(int patchPrimaryDecomposition, ideal i,list #)
 {
//           if patchPrimaryDecomposition=1,  drop the unit ideal in the decomposition,
//           since the unit ideal it is not prime, otherwise take no special action.
//           For other parameters see 'minAssGTZ' or 'minAssGTZE'
   ASSUME(0, hasFieldCoefficient(basering) );
   ASSUME(0, not isQuotientRing(basering) ) ;
   if(size(#)>0)
   {
      int keep_comp=1;
   }

  if(attrib(basering,"global")!=1)
  {
  // algorithms only work in global case!
// pass to appropriate global ring
      def r=basering;
      def s=changeord(list(list("dp",1:nvars(basering))));
      setring s;
      ideal i=imap(r,i);
// decompose and go back
      list li=minAssGTZ_i(patchPrimaryDecomposition,i);
      int sizeli = size(li);
      setring r;
      if (sizeli==0) { return(list()); }
      def li=imap(s,li);
// clean up
      if(!defined(keep_comp))
      {
         for(int k=size(li);k>=1;k--)
         {
            if(mindeg(std(lead(li[k]))[1])==0)
            {
// 1 contained in ideal, i.e. component does not meet origin in local ordering
               li=delete(li,k);
            }
         }
      }
      return(li);
  }

  int j;
  string algorithm;
  string facstdOption;
  string subsystem="noSubsystem";
  int useFac;

  // Set input parameters
  algorithm = "SL";         // Default: SL algorithm
  facstdOption = "facstd";
  if(size(#) > 0)
  {
    int valid;
    for(j = 1; j <= size(#); j++)
    {
      valid = 0;
      if((typeof(#[j]) == "int") or (typeof(#[j]) == "number"))
      {
        if (#[j] == 1) {facstdOption = "noFacstd"; valid = 1;}    // If #[j] == 1, facstd is not used.
        if (#[j] == 0) {facstdOption = "facstd";   valid = 1;}    // If #[j] == 0, facstd is used.
      }
      if(typeof(#[j]) == "string")
      {
        if((#[j] == "GTZ") || (#[j] == "SL"))
        {
          algorithm = #[j];
          valid = 1;
        }
        if((#[j] == "noFacstd") || (#[j] == "facstd"))
        {
          facstdOption = #[j];
          valid = 1;
        }
        if((#[j] == "noSubsystem") || (#[j] == "subsystem"))
        {
          subsystem = #[j];
          valid = 1;
        }
      }
      if(valid == 0)
      {
        dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
      }
    }
  }

  if(minpoly!=0)
  {
    return(algeDeco_i(patchPrimaryDecomposition,i,2));
  }

  list result = minAssPrimes_i(patchPrimaryDecomposition,i, facstdOption, algorithm, subsystem);
  return(result);
}


///////////////////////////////////////////////////////////////////////////////
proc minAssCharE(ideal I, list #)
"USAGE:   minAssCharE(I[,c]); i ideal, c int (optional).
RETURN:  list, the minimal associated prime ideals of I. If I is the unit ideal returns an empty list.
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. @*
         For local orderings, the result is considered in the localization
         of the polynomial ring, not in the power series ring
         For local and mixed orderings, the decomposition in the
         corresponding global ring is returned if the string 'global'
         is specified as third argument
EXAMPLE: example minAssCharE; shows an example
"
{
    return(minAssChar_i(1,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 = minAssCharE(I);
   pr;
   ideal J = 5;
   list prempty = minAssCharE(J);
   prempty;
}

proc minAssChar(ideal I, list #)
"USAGE:   minAssChar(I[,c]); i ideal, c int (optional).
RETURN:  list, the minimal associated prime ideals of I. If I is the unit ideal returns list( ideal(1) )
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. @*
         For local orderings, the result is considered in the localization
         of the polynomial ring, not in the power series ring
         For local and mixed orderings, the decomposition in the
         corresponding global ring is returned if the string 'global'
         is specified as third argument
EXAMPLE: example minAssChar; shows an example
"
{
    return(minAssChar_i(0,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 = minAssChar(i);
   pr;
}

proc minAssChar_i(int patchPrimaryDecomposition, ideal i, list #)
{
//           if patchPrimaryDecomposition=1,  drop the unit ideal in the decomposition,
//           since the unit ideal it is not prime, otherwise take no special action.
//           For other parameters see 'minAssChar' or 'minAssCharE'
   ASSUME(0, hasFieldCoefficient(basering) );
   ASSUME(0, not isQuotientRing(basering) ) ;
   ASSUME(0,size(#)<3);
   if(size(#)>1)
   {
      int keep_comp=1;
   }
   if(attrib(basering,"global")!=1)
   {
// algorithms only work in global case!
// pass to appropriate global ring
      def r=basering;
      def s=changeord(list(list("dp",1:nvars(basering))));
      setring s;
      ideal i=imap(r,i);
// decompose and go back
      list li=minAssChar_i(patchPrimaryDecomposition,i);
      int sizeli = size(li);
      setring r;
      if (sizeli==0) { return(list()); }
      def li=imap(s,li);
// clean up
      if(!defined(keep_comp))
      {
         for(int k=size(li);k>=1;k--)
         {
            if(mindeg(std(lead(li[k]))[1])==0)
            {
// 1 contained in ideal, i.e. component does not meet origin in local ordering
               li=delete(li,k);
            }
         }
      }
      return(li);
   }
   if (size(#)>0)
   { return(min_ass_prim_charsets_i(patchPrimaryDecomposition,i,#[1])); }
   else
   { return(min_ass_prim_charsets_i(patchPrimaryDecomposition,i,1)); }
}

///////////////////////////////////////////////////////////////////////////////
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 (with modifications by Laplagne) and Kemper is used.
         Works also in positive characteristic (Kempers algorithm).
EXAMPLE: example equiRadical; shows an example
"
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  if(attrib(basering,"global")!=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[, l]); I ideal, l list (optional)
 @*  Optional parameters in list l (can be entered in any order):
 @*  0, \"fullRad\" -> full radical is computed (default)
 @*  1, \"equiRad\" -> equiRadical is computed
 @*  \"KL\" -> Krick/Logar algorithm is used
 @*  \"SL\" -> modifications by Laplagne are used (default)
 @*  \"facstd\" -> uses facstd to first decompose the ideal (default for non homogeneous ideals)
 @*  \"noFacstd\" -> does not use facstd (default for homogeneous ideals)
RETURN:  ideal, the radical of I (or the equiradical if required in the input parameters)
NOTE:    A combination of the algorithms of Krick/Logar (with modifications by Laplagne) and Kemper is used.
         Works also in positive characteristic (Kempers algorithm).
EXAMPLE: example radical; shows an example
"
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  dbprint(printlevel - voice, "Radical, version 2006.05.08");
  if(size(i) == 0){return(ideal(0));}
  if(attrib(basering,"global")!=1)
  {
// algorithms only work in global case!
// pass to appropriate global ring
      def r=basering;
      def s=changeord(list(list("dp",1:nvars(basering))));
      setring s;
      ideal i=imap(r,i);
// compute radical and go back
      def j=radical(i);
      setring r;
      def j=imap(s,j);
      return(j);
  }
  if(hasAlgExtensionCoefficient(basering))
  {
    return(algeRad(i));
  }
  int j;
  def P0 = basering;
  list Pl=ring_list(P0);
  Pl[3]=list(list("dp",1:nvars(P0)),list("C",0));
  def @P=ring(Pl); // Ring with same variables and DP ordering in all variables
  setring @P;
  ideal i=imap(P0,i);

  int il;
  string algorithm;
  int useFac,useSubsystem;

  // Set input parameters
  algorithm = "SL";                                 // Default: SL algorithm
  il = 0;                                           // Default: Full radical (not only equiRadical)
  if (homog(i) == 1)
  {   // Default: facStd is used, except if the ideal is homogeneous.
    useFac = 0;
  }
  else
  {
    useFac = 1;
  }
  if(size(#) > 0)
  {
    int valid;
    for(j = 1; j <= size(#); j++)
    {
      valid = 0;
      if((typeof(#[j]) == "int") or (typeof(#[j]) == "number"))
      {
        il = #[j];          // If il == 1, equiRadical is computed
        valid = 1;
      }
      if(typeof(#[j]) == "string")
      {
        if(#[j] == "KL")
        {
          algorithm = "KL";
          valid = 1;
        }
        if(#[j] == "SL")
        {
          algorithm = "SL";
          valid = 1;
        }
        if(#[j] == "noFacstd")
        {
          useFac = 0;
          valid = 1;
        }
        if(#[j] == "facstd")
        {
          useFac = 1;
          valid = 1;
        }
        if(#[j] == "equiRad")
        {
          il = 1;
          valid = 1;
        }
        if(#[j] == "fullRad")
        {
          il = 0;
          valid = 1;
        }
        if(#[j] == "subsystem")
        {
          useSubsystem = 1;
          valid = 1;
        }
      }
      if(valid == 0)
      {
        dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
      }
    }
  }

  ideal rad = 1;
  intvec op = option(get);
  list qr = simplifyIdeal(i);
  map phi = @P, qr[2];  // map from the same ring, to simplify the ideal
  
  
  option(redSB);
  i = groebner(qr[1]);  // We use the transformed variables, and we 
                        // will use map phi to go back to the original
                        // variables
  option(set, op);
  int di = dim(i);

  if(di == 0)
  {
    i = zeroRad(i, qr[1]);
    option(redSB);
    i=interred(phi(i));
    option(set, op);
    setring(P0);
    i=imap(@P,i);
    return(i);
  }

  option(redSB);
  list pr;
  list qq;
  if (useSubsystem)
  {
    qq=prep_decomp(i);
  }
  else
  {
    qq=i;
  }
  if(useFac == 1)
  {
    for(int ll=size(qq);ll>0;ll--)
    {
      pr=pr+facstd(qq[ll]);
    }
    kill ll;
  }
  else
  {
    pr = i;
  }
  kill qq;
  option(set, op);
  int s = size(pr);
  if(useFac == 1)
  {
    dbprint(printlevel - voice, "Number of components returned by facstd: ", s);
  }
  for(j = 1; j <= s; j++)
  {
    attrib(pr[s + 1 - j], "isSB", 1);
    if((size(reduce(rad, pr[s + 1 - j], 5)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il))
    {
      // SL Debug messages
      dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]);
      dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j]));

      if(algorithm == "KL")
      {
        rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il));
      }
      if(algorithm == "SL")
      {
        rad = intersect(rad, radicalSL(pr[s + 1 - j], il));
      }
    }
    else
    {
      // SL Debug
      dbprint(printlevel-voice, "The radical of this component is not needed.");
      dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))",
              size(reduce(rad, pr[s + 1 - j], 5)));
      dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j]));
      dbprint(printlevel-voice, "il", il);
    }
  }
    
  rad=interred(phi(rad));  // Eliminate redundant polynomials
  setring(P0);
  i=imap(@P,rad);
  return(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 = radical(i);
   pr;
}

///////////////////////////////////////////////////////////////////////////////
//
// Computes the radical of I using KL algorithm.
// The only difference with the previous implementation of KL algorithm is
// that now it uses block dp instead of lp ordering for the reduction to the
// zerodimensional case.
// The reduction step has been moved to the new routine radicalReduction, so that it can be
// used also by radicalSL procedure.
//
static proc radicalKL(ideal I, ideal ser, list #)
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

// ideal I     The ideal for which the radical is computed
// ideal ser   Used to reduce components already obtained
// list #      If #[1] = 1, equiradical is computed.

  // I needs to be a Groebner basis.
  if (attrib(I, "isSB") != 1)
  {
    I = groebner(I);
  }

  ideal rad;                                // The radical
  int allIndep = 1;                // All max independent sets are used

  list result = radicalReduction(I, ser, allIndep, #);
  int done = result[3];
  rad = result[1];
  if (done == 0)
  {
    rad = intersect(rad, radicalKL(result[2], ideal(1), #));
  }
  return(rad);
}


///////////////////////////////////////////////////////////////////////////////
//
// Computes the radical of I via Laplagne algorithm, using zerodimensional radical in
// the zero dimensional case.
// For the reduction to the zerodimensional case, it uses the procedure
// radical, with some modifications to avoid the recursion.
//
static proc radicalSL(ideal I, list #)
// Input = I, ideal
//         #, list. If #[1] = 1, then computes only the equiradical.
// Output = (P, primaryDec) where P = rad(I) and primaryDec is the list of the radicals
// obtained in intermediate steps.
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  ideal rad = 1;
  ideal equiRad = 1;
  list primes;
  int k;                        // Counter
  int il;                 // If il = 1, only the equiradical is required.
  int iDim;                // The dimension of I
  int stop = 0;   // Checks if the radical has been obtained

  if (attrib(I, "isSB") != 1)
  {
    I = groebner(I);
  }
  iDim = dim(I);

  // Checks if only equiradical is required
  if (size(#) > 0)
  {
    il = #[1];
  }

  while(stop == 0)
  {
    dbprint (printlevel-voice, "// We call radLoopR to find new prime ideals.");
    primes = radicalSLIteration(I, rad);                         // A list of primes or intersections of primes, not included in P
    dbprint (printlevel - voice, "// Output of Iteration Step:");
    dbprint (printlevel - voice, primes);
    if (size(primes) > 0)
    {
      dbprint (printlevel - voice, "// We intersect P with the ideal just obtained.");
      for(k = 1; k <= size(primes); k++)
      {
        rad = intersect(rad, primes[k]);
        if (il == 1)
        {
          if (attrib(primes[k], "isSB") != 1)
          {
            primes[k] = groebner(primes[k]);
          }
          if (iDim == dim(primes[k]))
          {
            equiRad = intersect(equiRad, primes[k]);
          }
        }
      }
    }
    else
    {
      stop = 1;
    }
  }
  if (il == 0)
  {
    return(rad);
  }
  else
  {
    return(equiRad);
  }
}

//////////////////////////////////////////////////////////////////////////
// Based on radicalKL.
// It contains all of old version of proc radicalKL except the recursion call.
//
// Output:
// #1 -> output ideal, the part of the radical that has been computed
// #2 -> complementary ideal, the part of the ideal I whose radical remains to be computed
//       = (I, h) in KL algorithm
//       This is not used in the new algorithm. It is part of KL algorithm
// #3 -> done, 1: output = radical, there is no need to continue
//                   0: radical = output \cap \sqrt{complementary ideal}
//       This is not used in the new algorithm. It is part of KL algorithm

static proc radicalReduction(ideal I, ideal ser, int allIndep, list #)
{
// allMaximal      1 -> Indicates that the reduction to the zerodim case
//                    must be done for all indep set of the leading terms ideal
//                 0 -> Otherwise
// ideal ser       Only for radicalKL. (Same as in radicalKL)
// list #          Only for radicalKL (If #[1] = 1,
//                    only equiradical is required.
//                    It is used to set the value of done.)
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  attrib(I, "isSB", 1);   // I needs to be a reduced standard basis
  list indep, fett;
  intvec op;
  int @wr, @n, @m, lauf, di;
  ideal fac, @h, collectrad, lsau;
  poly @q;
  string @va;

  def @P = basering;
  int jdim = dim(I);               // Computes the dimension of I
  ideal rad = ideal(1);            // The unit ideal
  ideal te = ser;
  if(size(#) > 0)
  {
    @wr = #[1];
  }

  // SL 2006.04.11 1 Debug messages
  dbprint(printlevel-voice, "//Computes the radical of the ideal:", I);
  // SL 2006.04.11 2 Debug messages

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

  if (jdim==-1)
  {
    return(ideal(1), ideal(1), 1);
  }

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

  if (jdim==0)
  {
    return(zeroRad(I), ideal(1), 1);
  }

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

  // SL 2006-04-24 1   If allIndep = 0, then it only computes one maximal
  //                     independent set.
  //                     This looks better for the new algorithm but not for KL
  //                     algorithm
  list parameters = allIndep;
  indep = newMaxIndependSetDp(I, parameters);
  // SL 2006-04-24 2

  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
    {
      def gnir1=ring(ring_list(basering));
      setring gnir1;
      ideal @j = fetch(@P, I);
      attrib(@j, "isSB", 1);
    }
    else
    {
      @va = string(maxideal(1));

      ring gnir1 = create_ring(ring_list(basering)[1],"(" + indep[@m][1] + ")","(" + indep[@m][2] + ")","no_minpoly");
      execute("map phi = @P," + @va + ";");
      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 there 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 LCM, 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..]
    //------------------------------------------------------------------------
    def quotring = prepareQuotientring(nvars(basering) - indep[@m][3],"dp");
    //------------------------------------------------------------------------
    // We pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
    //------------------------------------------------------------------------

    setring 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 LCM(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 = std(@j);  //to obtain a reduced standardbasis
      option(set, op);

      // SL 1 Debug messages
      dbprint(printlevel - voice, "zero_rad", basering, @j, dim(groebner(@j)));
      ideal zero_rad = zeroRad(@j);
      dbprint(printlevel - voice, "zero_rad passed");
      // SL 2
      if(!attrib(zero_rad,"isSB")) { zero_rad = std(zero_rad);}
    }
    else
    {
      ideal zero_rad = ideal(1);
    }

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


    ideal hpl;

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

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

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

    kill quotring;

    // 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, 6);
      @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, 5), 2);
    }
    else
    {
      if(deg(collectrad[1]) > 0)
      {
        rad = collectrad;
        te = intersect(te, collectrad);
        te = simplify(reduce(te, I, 5), 2);
      }
    }

    if((dim(I) < jdim)||(size(te) == 0))
    {
      break;
    }
  }

  // SL 2006.04.11 1 Debug messages
  dbprint (printlevel-voice, "// Part of the Radical already computed:", rad);
  dbprint (printlevel-voice, "// Dimension:", dim(groebner(rad)));
  // SL 2006.04.11 2 Debug messages

  // SL 2006.04.21 1    New variable "done".
  //                      It tells if the radical is already computed or
  //                      if it still has to be computed the radical of the new ideal I
  int done;
  if(((@wr == 1) && (dim(I)<jdim)) || (deg(I[1])==0) || (size(te) == 0))
  {
    done = 1;
  }
  else
  {
    done = 0;
  }
  // SL 2006.04.21 2

  // SL 2006.04.21 1     See details of the output at the beginning of this proc.
  list result = rad, I, done;
  return(result);
  // SL 2006.04.21 2
}

///////////////////////////////////////////////////////////////////////////////
// Given an ideal I and an ideal P (intersection of some minimal prime ideals
// associated to I), it calculates the intersection of new minimal prime ideals
// associated to I which where not used to calculate P. P = 1 represents empty intersection.
// This version uses ZD Radical in the zerodimensional case.
static proc radicalSLIteration (ideal I, ideal P);
// Input: I, ideal. The ideal from which new prime components will be obtained.
//        P, ideal. Intersection of some prime ideals of I.
// Output: ideal. Intersection of some primes of I different from the ones in P.
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  int k = 1;                     // Counter
  int good  = 0;                 // Checks if an element of P is in rad(I)

  dbprint (printlevel-voice, "// We search for an element in P - sqrt(I).");
  while ((k <= size(P)) and (good == 0))
  {
    dbprint (printlevel-voice, "// We try with:", P[k]);
    good = 1 - rad_con(P[k], I);
    k++;
  }
  k--;
  if (good == 0)
  {
    dbprint (printlevel-voice, "// No element was found, P = sqrt(I).");
    list emptyList = list();
    return (emptyList);
  }
  dbprint(printlevel - voice, "// That one was good!");
  dbprint(printlevel - voice, "// We saturate I with respect to this element.");
  if (P[k] != 1)
  {
    intvec oo=option(get);
    option(redSB);
    ideal J = sat(I, P[k]);
    option(set,oo);

  }
  else
  {
    dbprint(printlevel - voice, "// The polynomial is 1, the saturation in not actually computed.");
    ideal J = I;
  }

  // We now call proc radicalNew;
  dbprint(printlevel - voice, "// We do the reduction to the zerodimensional case, via radical.");
  dbprint(printlevel - voice, "// The ideal is ", J);
  dbprint(printlevel - voice, "// The dimension is ", dim(groebner(J)));

  int allMaximal = 0;   // Compute the zerodim reduction for only one indep set.
  ideal re = 1;         // No reduction is need,
                        //    there are not redundant components.
  list emptyList = list();   // Look for primes of any dimension,
                             //   not only of max dimension.
  list result = radicalReduction(J, re, allMaximal, emptyList);

  return(result[1]);
}

///////////////////////////////////////////////////////////////////////////////////
// Based on maxIndependSet
// Added list # as parameter
// If the first element of # is 0, the output is only 1 max indep set.
// If no list is specified or #[1] = 1, the output is all the max indep set of the
// leading terms ideal. This is the original output of maxIndependSet

// The ordering given in the output has been changed to block dp instead of lp.

proc newMaxIndependSetDp(ideal j, list #)
"USAGE:   newMaxIndependentSetDp(I); I ideal (returns all maximal independent sets of the corresponding leading terms ideal)
          newMaxIndependentSetDp(I, 0); I ideal (returns only one maximal independent set)
RETURN:  list = #1. new varstring with the maximal independent set at the end,
                #2. ordstring with the corresponding dp block ordering,
                #3. the number of independent variables
NOTE:
EXAMPLE: example newMaxIndependentSetDp; shows an example
"
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  ASSUME(0, hasGlobalOrdering(basering) ) ;

  int n, k, di;
  list resu, hilf;
  string var1, var2;
  list v = indepSet(j, 0);

  // SL 2006.04.21 1 Lines modified to use only one independent Set
  int allMaximal;
  if (size(#) > 0)
  {
    allMaximal = #[1];
  }
  else
  {
    allMaximal = 1;
  }

  int nMax;
  if (allMaximal == 1)
  {
    nMax = size(v);
  }
  else
  {
    nMax = 1;
  }

  for(n = 1; n <= nMax; n++)
  // SL 2006.04.21 2
  {
    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) - 2];                         // The "- 2" removes the trailer comma
      hilf[1] = var1;
      // SL 2006.21.04 1 The order is now block dp instead of lp
      hilf[2] = "dp(" + string(nvars(basering) - di) + "), dp(" + string(di) + ")";
      // SL 2006.21.04 2
      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 = newMaxIndependSetDp(i);
   l;
   i = i, g;
   l = newMaxIndependSetDp(i);
   l;

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


///////////////////////////////////////////////////////////////////////////////
proc prepareAss(ideal i)
"USAGE:   prepareAss(I); I ideal
RETURN:  list of radicals of the equidimensional components of I
NOTE:    Uses algorithm of Eisenbud/Huneke/Vasconcelos.
EXAMPLE: example prepareAss; shows an example
"
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  if(attrib(basering,"global")!=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)
  {
     resolution re=sres(j,0);                   //the resolution
     re=minres(re);                       //minimized resolution
  }
  else
  {
    list re=mres(i,0);
  }
  for(e=nvars(basering);e>=cod;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
"
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  if(attrib(basering,"global")!=1)
  {
      ERROR(
      "// Not implemented for this ordering, please change to global ordering."
      );
  }

  ideal j=groebner(i);
  int cod=nvars(basering)-dim(j);


  if(cod > nvars(basering))
    {
      dbprint(printlevel,"//If I is the entire ring...");
      dbprint(printlevel,"//...then return the ideal generated by 1.");
      return(ideal(1));
    }

  int e;
  ideal ann;
  if(homog(i)==1)
  {
     resolution re=sres(j,0);                   //the resolution
     re=minres(re);                       //minimized resolution
  }
  else
  {
    resolution re=mres(j,0);
  }
  ann = AnnExt_R(cod,re);
  if( nvars(basering)-dim(std(ann) ) != cod)
  {
     return( ideal(1) );
  }

  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 testPrimaryE(list pr, ideal k)
"USAGE:   testPrimaryE(pr,k); pr a list, k an ideal.
ASSUME:  pr is the result of a primary decomposition and may be empty ( for the unit ideal)
RETURN:  int, 1 if the intersection of the ideals in pr is k, 0 if not
EXAMPLE: example testPrimaryE; shows an example
"
{
    return(testPrimary_i(1,pr,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);
   testPrimaryE(pr,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
"
{
    return(testPrimary_i(0,pr,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);
}


static proc testPrimary_i(int patchPrimaryDecomposition,list pr, ideal k)
{
//           if patchPrimaryDecomposition=1,  handle the case of an empty decomposition list.
//           For other parameters see 'testPrimary'
   ASSUME(0, hasFieldCoefficient(basering) );
   ASSUME(0, not isQuotientRing(basering) ) ;

   int i;
   pr=reconvList(pr);
   if (patchPrimaryDecomposition==1)
   {
       if (idealsEqual( k, ideal(1)) )
       {
          return( size(pr)==0 ); //list expected to be empty.
       }
   }
   ideal j=pr[1];


   for (i=2;i<=size(pr) div 2;i++)
   {
       j=intersect(j,pr[2*i-1]);
   }
   return(idealsEqual(j,k));
}


///////////////////////////////////////////////////////////////////////////////
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
"
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  ASSUME(2, dim(groebner(I))==0 );
  if(attrib(basering,"global")!=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 --------------
  ring P1 = create_ring(ring_list(R)[1],"T","dp","no_minpoly");
  matrix n=imap(R,n);
  time = timer;
  poly charpol=det(n-var(1)*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);
    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.
    ring P = create_ring(ring_list(R)[1],"(T,"+varstr(R)+")","(dp(1),dp)","no_minpoly");
    list rfac=imap(P1,fac);
    intvec ov=option(get);;
    option(redSB);
    list re1;
    ideal new = var(1)-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;
}


static proc newDecompStepE(ideal I, list #)
{
   return(newDecompStep_i(1,I,#));
}

static proc newDecompStep(ideal I, list #)
{
    return(newDecompStep_i(0,I,#));
}

///////////////////////////////////////////////////////////////////////////////
static proc newDecompStep_i(int patchPrimaryDecomposition, ideal i, list #)
"USAGE:  newDecompStep_i(patchPrimaryDecomposition, I); I ideal  (for primary decomposition)
         newDecompStep_i(patchPrimaryDecomposition, I,1);        (for the associated primes of dimension of i)
         newDecompStep_i(patchPrimaryDecomposition, I,2);        (for the minimal associated primes)
         newDecompStep_i(patchPrimaryDecomposition, I,3);        (for the absolute primary decomposition (not tested!))
         "oneIndep";        (for using only one max indep set)
         "intersect";        (returns also the intersection of the components founded)

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
         if patchPrimaryDecomposition=1,  drop the unit ideal in the decomposition,
         since the unit ideal it is not prime, otherwise take no special action.
EXAMPLE: example newDecompStep; shows an example
"
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  intvec op@P, op,@vv;
  def  @P = basering;
  list primary,indep,ltras;
  intvec @vh,isat,@w;
  int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi,abspri,ab,nn;
  ideal peek=i;
  ideal ser,tras;
  list data;
  list result;
  bigintvec @hilb;
  int isS=attrib(i,"isSB");

  // Debug
  dbprint(printlevel - voice, "newDecompStep, v2.0");

  string indepOption = "allIndep";
  string intersectOption = "noIntersect";

  if(size(#)>0)
  {
    int count = 1;
    if(typeof(#[count]) == "string")
    {
      if ((#[count] == "oneIndep") or (#[count] == "allIndep"))
      {
        indepOption = #[count];
        count++;
      }
    }
    if(typeof(#[count]) == "string")
    {
      if ((#[count] == "intersect") or (#[count] == "noIntersect"))
      {
        intersectOption = #[count];
        count++;
      }
    }
    if((typeof(#[count]) == "int") or (typeof(#[count]) == "number"))
    {
      if ((#[count]==1)||(#[count]==2)||(#[count]==3))
      {
        @wr=#[count];
        if(@wr==3){abspri = 1; @wr = 0;}
        count++;
      }
    }
    if(size(#)>count)
    {
      seri=1;
      peek=#[count + 1];
      ser=#[count + 2];
    }
  }
  if(abspri)
  {
    list absprimary,abskeep,absprimarytmp,abskeeptmp;
  }
  homo=homog(i);
  if(homo==1)
  {
    if(attrib(i,"isSB"))
    {
      ltras=i,i;
    }
    else
    {
      //ltras=mstd(i);
      tras=groebner(i);
      ltras=tras,tras;
    }
    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[2]=maxideal(1);
        l[1]=maxideal(1);
        if (intersectOption == "intersect")
        {
          return(list(l, maxideal(1)));
        }
        else
        {
          return(l);
        }
      }
      if (intersectOption == "intersect")
      {
        return(list(primary, primary[1]));
      }
      else
      {
        return(primary);
      }
    }
    for(@n=1;@n<=nvars(basering);@n++)
    {
      @w[@n]=ord(var(@n));
    }
    @hilb=hilb(tras,1,@w);
    bigintvec keephilb=@hilb;
  }

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

  if(size(i)==0)
  {
    primary=i,i;
    if (intersectOption == "intersect")
    {
      return(list(primary, i));
    }
    else
    {
      return(primary);
    }
  }

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

  int lp=islp();

  op@P = option(get);
  def gnir=changeordTo(basering,"lp");
  setring gnir;

  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)&&(!abspri))
  {
    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;
      option(set,op@P);
      primary[1]=i;
      primary[2]=i;
      if (intersectOption == "intersect")
      {
        return(list(primary, i));
      }
      else
      {
        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)]=";";
      ring @deirf = create_ring(ring_list(gnir)[1],"("+newva,"lp","no_minpoly");
      execute("map @kappa=gnir,"+newma);
      ideal @j= @kappa(@j);
      @j=simplify(@j, 2);
      attrib(@j,"isSB",1);
      result = newDecompStep_i(patchPrimaryDecomposition, @j, indepOption, intersectOption, @wr);
      if (intersectOption == "intersect")
      {
        list pr = result[1];
        ideal intersection = result[2];
      }
      else
      {
        list pr = result;
      }

      setring gnir;
      list pr=imap(@deirf,pr);
      for(@k=1;@k<=size(pr);@k++)
      {
        @j=pr[@k]+fried;
        pr[@k]=@j;
      }
      if (intersectOption == "intersect")
      {
        ideal intersection = imap(@deirf, intersection);
        @j = intersection + fried;
        intersection = @j;
      }
      setring @P;
      option(set,op@P);
      if (intersectOption == "intersect")
      {
        return(list(imap(gnir,pr), imap(gnir,intersection)));
      }
      else
      {
        return(imap(gnir,pr));
      }
    }
  }
  //----------------------------------------------------------------
  //j is the ring
  //----------------------------------------------------------------

  if (dim(@j)==-1)
  {
    setring @P;
    option(set,op@P);
    primary=ideal(1),ideal(1);
    if (intersectOption == "intersect")
    {
      return(list(primary, ideal(1)));
    }
    else
    {
      return(primary);
    }
  }

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

  if(nvars(basering)==1)
  {
    list fac=factor(@j[1]);
    list gprimary;
    poly generator=1;
    ideal gIntersection;
    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]);
      }
      if (intersectOption == "intersect")
      {
        generator = generator * fac[1][@k];
      }
    }
    if (intersectOption == "intersect")
    {
      gIntersection = generator;
    }
    setring @P;
    primary=fetch(gnir,gprimary);
    if (intersectOption == "intersect")
    {
      ideal intersection = fetch(gnir,gIntersection);
    }

    if(abspri)
    {
      list resu,tempo;
      string absotto;
      for(ab=1;ab<=size(primary) div 2;ab++)
      {
        absotto= absFactorize(primary[2*ab][1],77);
        tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering)));
        resu[ab]=tempo;
      }
      primary=resu;
      intersection = 1;
      for(ab=1;ab<=size(primary);ab++)
      {
        intersection = intersect(intersection, primary[ab][2]);
      }
    }
    if (intersectOption == "intersect")
    {
      return(list(primary, intersection));
    }
    else
    {
      return(primary);
    }
  }

 //------------------------------------------------------------------
 //the zero-dimensional case
 //------------------------------------------------------------------
  if (dim(@j)==0)
  {
    op=option(get);
    option(redSB);
    list gprimary= newZero_decomp(@j,ser,@wr);

    setring @P;
    primary=fetch(gnir,gprimary);

    if(size(ser)>0)
    {
      primary=cleanPrimary(primary);
    }
    if(abspri)
    {
      list resu,tempo;
      string absotto;
      for(ab=1;ab<=size(primary) div 2;ab++)
      {
        absotto= absFactorize(primary[2*ab][1],77);
        tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering)));
        resu[ab]=tempo;
      }
      primary=resu;
    }
    option(set,op@P);
    if (intersectOption == "intersect")
    {
      return(list(primary, fetch(gnir,@j)));
    }
    else
    {
      return(primary);
    }
  }

  poly @gs,@gh,@p;
  string @va;
  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 = newMaxIndependSetLp(@j, indepOption);
    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 = newMaxIndependSetLp(@j, indepOption);
  }

  ideal jkeep=@j;
  if((ordstr(@P)[1]=="w")&&(size(ring_list(@P)[3])==2)) // weighted ordering
  {
    def @Phelp=ring(ring_list(gnir));
    setring @Phelp;
  }
  else
  {
    def @Phelp=changeordTo(gnir,"dp");
    setring @Phelp;
  }

  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;
//Aenderung================
  ideal @Ptest=1;
//=========================
  di=dim(jwork);
  keepdi=di;

  ser = 1;

  setring gnir;
  for(@m=1; @m<=size(indep); @m++)
  {
    data[1] = indep[@m];
    result = newReduction(@j, ser, @hilb, @w, jdim, abspri, @wr, data);
    quprimary = quprimary + result[1];
    if(abspri)
    {
      absprimary = absprimary + result[2];
      abskeep = abskeep + result[3];
    }
    @h = result[5];
    ser = result[4];
    if(size(@h)>0)
    {
      //---------------------------------------------------------------
      //we change to @Phelp to have the ordering dp for saturation
      //---------------------------------------------------------------

      setring @Phelp;
      @h=imap(gnir,@h);
//Aenderung==================================
      if(defined(@LL)){kill @LL;}
      list @LL=minSat(jwork,@h);
      @Ptest=intersect(@Ptest,@LL[1]);
      ser = intersect(ser, @LL[1]);
//===========================================

      if(@wr!=1)
      {
//Aenderung==================================
        @q=@LL[2];
//===========================================
        //@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,6);
        @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);
        ser = imap(@Phelp, ser);
        break;
      }
      if(homo == 1)
      {
        @hilb = hilb(jwork, 1, @w);
      }

      setring gnir;
      ser = imap(@Phelp, ser);
      @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;
    quprimary[1]=ideal(1);
    quprimary[2]=ideal(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)
        {
//Aenderung====================================================
//HIER STATT DURCHSCHNITT SATURIEREN!
           ideal htest=@Ptest;
/*
           ideal htest=hquprimary[1];
           for (@n1=2;@n1<=size(hquprimary) div 2;@n1++)
           {
              htest=intersect(htest,hquprimary[2*@n1-1]);
           }
*/
//=============================================================
        }
        else
        {
           ideal htest=hquprimary[2];

           for (@n1=2;@n1<=size(hquprimary) div 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,5))!=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
          {
             def gnir1=ring(ring_list(basering));
             setring gnir1;
             ideal @j=fetch(gnir,jkeep);
             attrib(@j,"isSB",1);
          }
          else
          {
             @va=string(maxideal(1));
             ring gnir1 = create_ring(ring_list(basering)[1],"(" +restindep[@m][1]+")","(" +restindep[@m][2]+")","no_minpoly");
             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..]
          //------------------------------------------------------------------

          def quotring=prepareQuotientring(nvars(basering)-restindep[@m][3],"lp");

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

          setring 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= newZero_decomp(@j,ser,@wr);
          if(abspri)
          {
             ideal II;
             ideal jmap;
             map sigma;
             nn=nvars(basering);
             map invsigma=basering,maxideal(1);
             for(ab=1;ab<=size(uprimary) div 2;ab++)
             {
                II=uprimary[2*ab];
                attrib(II,"isSB",1);
                if(deg(II[1])!=vdim(II))
                {
                   jmap=randomLast(50);
                   sigma=basering,jmap;
                   jmap[nn]=2*var(nn)-jmap[nn];
                   invsigma=basering,jmap;
                   II=groebner(sigma(II));
                 }
                 absprimarytmp[ab]= absFactorize(II[1],77);
                 II=var(nn);
                 abskeeptmp[ab]=string(invsigma(II));
                 invsigma=basering,maxideal(1);
             }
          }
          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(quotring,uprimary);
          lsau=imap(quotring,saturn);
          @h=imap(quotring,@h);

          kill quotring;


          @n2=size(quprimary);
//================NEU=========================================
          if(deg(quprimary[1][1])<=0){ @n2=0; }
//============================================================

          @n3=@n2;

          for(@n1=1;@n1<=size(collectprimary) div 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];
                if(abspri)
                {
                  absprimary[@n2 div 2]=absprimarytmp[@n1];
                  abskeep[@n2 div 2]=abskeeptmp[@n1];
                }
             }
          }


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

          for(@n=@n3 div 2+1;@n<=@n2 div 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 div 2+1;@n<=@n2 div 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(abspri)
       {
         list resu,tempo;
         for(ab=1;ab<=size(quprimary) div 2;ab++)
         {
           if (deg(quprimary[2*ab][1])!=0)
           {
             tempo=quprimary[2*ab-1],quprimary[2*ab],
                        absprimary[ab],abskeep[ab];
             resu[ab]=tempo;
           }
         }
         quprimary=resu;
         @wr=3;
       }
       if(size(reduce(ser,peek,5))!=0)
       {
         if(@wr>0)
         {
            // The following line was dropped to avoid the recursion step:
            //htprimary=newDecompStep_i(patchPrimaryDecomposition, @j,@wr,peek,ser);
            htprimary = list();
         }
         else
         {
            // The following line was dropped to avoid the recursion step:
            //htprimary=newDecompStep_i(patchPrimaryDecomposition,@j,peek,ser);
            htprimary = list();
         }
         // here we collect now both results primary(sat(j,gh))
         // and primary(j,gh^n)
         @n=size(quprimary);
         if (deg(quprimary[1][1])<=0) { @n=0; }
         for (@k=1;@k<=size(htprimary);@k++)
         {
           quprimary[@n+@k]=htprimary[@k];
         }
       }
     }
   }
   else
   {
     if(abspri)
     {
       list resu,tempo;
       for(ab=1;ab<=size(quprimary) div 2;ab++)
       {
         tempo=quprimary[2*ab-1],quprimary[2*ab],
                   absprimary[ab],abskeep[ab];
         resu[ab]=tempo;
       }
       quprimary=resu;
     }
   }
  //---------------------------------------------------------------------------
  //back to the ring we started with
  //the final result: primary
  //---------------------------------------------------------------------------

  setring @P;
  option(set,op@P);
  primary=imap(gnir,quprimary);

  if (intersectOption == "intersect")
  {
    return(list(primary, imap(gnir, ser)));
  }
  else
  {
    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;
   int patchDecomposition = 1;
   list pr = newDecompStep_i(patchDecomposition, I);
   pr;
   testPrimary( pr, I);
}

// This was part of proc decomp.
// In proc newDecompStep, used for the computation of the minimal associated primes,
// this part was separated as a soubrutine to make the code more clear.
// Also, since the reduction is performed twice in proc newDecompStep, it should use both times this routine.
// This is not yet implemented, since the reduction is not exactly the same and some changes should be made.
static proc newReduction(ideal @j, ideal ser, bigintvec @hilb, intvec @w, int jdim, int abspri, int @wr, list data)
{
   ASSUME(1, hasFieldCoefficient(basering) );
   ASSUME(1, not isQuotientRing(basering) ) ;
   ASSUME(1, hasGlobalOrdering(basering) ) ;


   string @va;
   def quotring;
   intvec op;
   intvec @vv;
   def gnir = basering;
   ideal isat=0;
   int @n;
   int @n1 = 0;
   int @n2 = 0;
   int @n3 = 0;
   int homo = homog(@j);
   int lauf;
   int @k;
   list fett;
   int keepdi;
   list collectprimary;
   list lsau;
   list lnew;
   ideal @h;

   list indepInfo = data[1];
   list quprimary = list();

   //if(abspri)
   //{
     int ab;
     list absprimarytmp,abskeeptmp;
     list absprimary, abskeep;
   //}
   // Debug
   dbprint(printlevel - voice, "newReduction, v2.0");

   if((indepInfo[1]==varstr(basering)))  // &&(@m==1)
   //this is the good case, nothing to do, just to have the same notations
   //change the ring
   {
     def gnir1=ring(ring_list(basering));
     setring gnir1;
     ideal @j = fetch(gnir, @j);
     attrib(@j,"isSB",1);
     ideal ser = fetch(gnir, ser);
   }
   else
   {
     @va=string(maxideal(1));
//Aenderung==============
     //if(@m==1)
     //{
     //  @j=fetch(@P,i);
     //}
//=======================
     ring gnir1 = create_ring(ring_list(basering)[1],"("+indepInfo[1]+")","(" +indepInfo[2]+")","no_minpoly");
     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)-indepInfo[3],"lp");

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

   ideal @jj=lead(@j);               //!! vorn vereinbaren
   setring quotring;

   ideal @jj=imap(gnir1,@jj);
   @vv=clearSBNeu(@jj,fett);  //!! vorn vereinbaren
   setring gnir1;
   @k=size(@j);
   for (lauf=1;lauf<=@k;lauf++)
   {
     if(@vv[lauf]==1)
     {
       @j[lauf]=0;
     }
   }
   @j=simplify(@j,2);
   setring 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

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

     int zeroMinAss = @wr;
     if (@wr == 2) {zeroMinAss = 1;}
     list uprimary= newZero_decomp(@j, ser, zeroMinAss);

//HIER
     if(abspri)
     {
       ideal II;
       ideal jmap;
       map sigma;
       nn=nvars(basering);
       map invsigma=basering,maxideal(1);
       for(ab=1;ab<=size(uprimary) div 2;ab++)
       {
         II=uprimary[2*ab];
         attrib(II,"isSB",1);
         if(deg(II[1])!=vdim(II))
         {
           jmap=randomLast(50);
           sigma=basering,jmap;
           jmap[nn]=2*var(nn)-jmap[nn];
           invsigma=basering,jmap;
           II=groebner(sigma(II));
         }
         absprimarytmp[ab]= absFactorize(II[1],77);
         II=var(nn);
         abskeeptmp[ab]=string(invsigma(II));
         invsigma=basering,maxideal(1);
       }
     }
     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(quotring,uprimary);
   lsau=imap(quotring,saturn);
   @h=imap(quotring,@h);

   kill quotring;

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

   for(@n1=1;@n1<=size(collectprimary) div 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];
       if(abspri)
       {
         absprimary[@n2 div 2]=absprimarytmp[@n1];
         abskeep[@n2 div 2]=abskeeptmp[@n1];
       }
     }
   }

   //here the intersection with the polynomialring
   //mentioned above is really computed
   for(@n=@n3 div 2+1;@n<=@n2 div 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];
     }
   }

   return(quprimary, absprimary, abskeep, ser, @h);
}


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




///////////////////////////////////////////////////////////////////////////////
// Based on minAssGTZ
proc minAssE(ideal I,list #)
"USAGE:   minAssE(I[, l]); I ideal, l list (optional) of parameters, same as minAssGTZ
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 minAssE; shows an example
"
{
 return(minAss_i(1,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 = minAssE(i);
   pr;
   ideal j = 1;
   list prempty = minAssE(j);
   prempty;
}

proc minAss(ideal I,list #)
"USAGE:   minAss(I[, l]); I ideal, l list (optional) of parameters, same as minAssGTZ
RETURN:  a list, the minimal associated prime ideals of I. If I is the unit ideal, returns list(ideal(1));
NOTE:    Designed for characteristic 0, works also in char k > 0 based
         on an algorithm of Yokoyama
EXAMPLE: example minAss; shows an example
"
{
    return(minAss_i(0,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 = minAss(i);
   pr;
}

static proc minAss_i(int patchPrimaryDecomposition,ideal I,list #)
{
//         if patchPrimaryDecomposition=1,  drop the unit ideal in the decomposition,
//         since the unit ideal it is not prime, otherwise take no special action.
//         For other parameters see 'minAss' or 'minAssE'
  return(minAssGTZ_i(patchPrimaryDecomposition,I,#));
}



///////////////////////////////////////////////////////////////////////////////
//
// Computes the minimal associated primes of I via Laplagne algorithm,
// using primary decomposition in the zero dimensional case.
// For reduction to the zerodimensional case, it uses the procedure
// decomp, with some modifications to avoid the recursion.
//

static proc minAssSL(ideal I)
// Input = I, ideal
// Output = primaryDec where primaryDec is the list of the minimal
// associated primes and the primary components corresponding to these primes.
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  ideal P = 1;
  list pd = list();
  int k;
  int stop = 0;
  list primaryDec = list();

  while (stop == 0)
  {
    // Debug
    dbprint(printlevel - voice, "// We call minAssSLIteration to find new prime ideals!");
    pd = minAssSLIteration(I, P);
    // Debug
    dbprint(printlevel - voice, "// Output of minAssSLIteration:");
    dbprint(printlevel - voice, pd);
    if (size(pd[1]) > 0)
    {
      primaryDec = primaryDec + pd[1];
      // Debug
      dbprint(printlevel - voice, "// We intersect the prime ideals obtained.");
      P = intersect(P, pd[2]);
      // Debug
      dbprint(printlevel - voice, "// Intersection finished.");
    }
    else
    {
      stop = 1;
    }
  }

  // Returns only the primary components, not the radical.
  return(primaryDec);
}

///////////////////////////////////////////////////////////////////////////////
// Given an ideal I and an ideal P (intersection of some minimal prime ideals
// associated to I), it calculates new minimal prime ideals associated to I
// which were not used to calculate P. P = 1 represents empty intersetion.
// This version uses Primary Decomposition in the zerodimensional case.
static proc minAssSLIteration(ideal I, ideal P);
{
  ASSUME(1, hasFieldCoefficient(basering) );
  ASSUME(1, not isQuotientRing(basering) ) ;
  ASSUME(1, hasGlobalOrdering(basering) ) ;

  int k = 1;
  int good  = 0;
  list primaryDec = list();
  // Debug
  dbprint (printlevel-voice, "// We search for an element in P - sqrt(I).");
  while ((k <= size(P)) and (good == 0))
  {
    good = 1 - rad_con(P[k], I);
    k++;
  }
  k--;
  if (good == 0)
  {
    // Debug
    dbprint (printlevel - voice, "// No element was found, P = sqrt(I).");
    return (list(primaryDec, ideal(0)));
  }
  // Debug
  dbprint (printlevel - voice, "// We found h = ", P[k]);
  dbprint (printlevel - voice, "// We calculate the saturation of I with respect to the element just founded.");
  ideal J = sat(I, P[k]);

  // Uses decomp from primdec, modified to avoid the recursion.
  // Debug
  dbprint(printlevel - voice, "// We do the reduction to the zerodimensional case, via decomp.");

  primaryDec = newDecompStep_i( 1, J, "oneIndep", "intersect", 2);
  // Debug
  dbprint(printlevel - voice, "// Proc decomp has found", size(primaryDec) div 2, "new primary components.");

  return(primaryDec);
}



///////////////////////////////////////////////////////////////////////////////////
// Based on maxIndependSet
// Added list # as parameter
// If the first element of # is 0, the output is only 1 max indep set.
// If no list is specified or #[1] = 1, the output is all the max indep set of the
// leading terms ideal. This is the original output of maxIndependSet

proc newMaxIndependSetLp(ideal j, list #)
"USAGE:   newMaxIndependentSetLp(i); i ideal (returns all maximal independent sets of the corresponding leading terms ideal)
          newMaxIndependentSetLp(i, 0); i ideal (returns only one maximal independent set)
RETURN:  list = #1. new varstring with the maximal independent set at the end,
                #2. ordstring with the lp ordering,
                #3. the number of independent variables
NOTE:
EXAMPLE: example newMaxIndependentSetLp; shows an example
"
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  ASSUME(0, hasGlobalOrdering(basering) ) ;

  int n, k, di;
  list resu, hilf;
  string var1, var2;
  list v = indepSet(j, 0);

  // SL 2006.04.21 1 Lines modified to use only one independent Set
  string indepOption;
  if (size(#) > 0)
  {
    indepOption = #[1];
  }
  else
  {
    indepOption = "allIndep";
  }

  int nMax;
  if (indepOption == "allIndep")
  {
    nMax = size(v);
  }
  else
  {
    nMax = 1;
  }

  for(n = 1; n <= nMax; n++)
  // SL 2006.04.21 2
  {
    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) - 2];       // The "- 2" removes the trailer comma
      hilf[1] = var1;
      // SL 2006.21.04 1 The order is now block dp instead of lp
      //hilf[2] = "dp(" + string(nvars(basering) - di) + "), dp(" + string(di) + ")";
      // SL 2006.21.04 2
      // For decomp, lp ordering is needed. Nothing is changed.
      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 = newMaxIndependSetLp(i);
   l;
   i = i, g;
   l = newMaxIndependSetLp(i);
   l;

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


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

proc newZero_decomp (ideal j, ideal ser, int @wr, list #)
"USAGE:   newZero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1
         (@wr=0 for primary decomposition, @wr=1 for computation of associated
         primes)
         if #[1] = "nest", then #[2] indicates the nest level (number of recursive calls)
         When the nest level is high it indicates that the computation is difficult,
         and different methods are applied.
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 newZero_decomp; shows an example
"
{
  ASSUME(0, hasFieldCoefficient(basering) );
  ASSUME(0, not isQuotientRing(basering) ) ;
  ASSUME(0, hasGlobalOrdering(basering) ) ;

  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;
  bigintvec @vh,@hilb;
  string @ri;
  poly @f;

  // Debug
  dbprint(printlevel - voice, "proc newZero_decomp");

  if (dim(j)>0)
  {
    ERROR("dim(j)>0 . Please send the example to the authors");
  }
  j=interred(j);

  attrib(j,"isSB",1);

  int nestLevel = 0;
  if (size(#) > 0)
  {
    if (typeof(#[1]) == "string")
    {
      if (#[1] == "nest")
      {
        nestLevel = #[2];
      }
      # = list();
    }
  }

  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],5))==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,5))==0))
    {
      ERROR("dim(ser/j)==-1 unexpected. Please send the example to the authors");
    }
    if(dim(j)==-1)
    {
      ERROR("dim(j)==-1 unexpected. Please send the example to the authors");
    }
    else
    {
      primary[2]=maxideal(1);
    }
    return(primary);
  }

//the first element in the standardbase is factorized
  if(deg(j[1])>0)
  {
    act=factor(j[1]);
  }
  else
  {
    ERROR("failure in newZero_decomp. Please send the example to the authors");
  }

//with the factors new ideals (hopefully the primary decomposition)
//are created
  if(size(act[1])>1)
  {
    if(size(#)>1)
    {
      ERROR("failure in newZero_decomp. Please send the example to the authors");
    }
    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(@wr!=0)
    {
      primary[1]=std(j,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>=7)&&(char(basering)<=181))
  {
    primary=splitCharp(primary);
  }

  if((@wr==2)&&(npars(basering)>0)&&(voice>=7)&&(char(basering)>0))
  {
  //the prime decomposition of Yokoyama in characteristic p
    list ke,ek;
    @k=0;
    while(@k<size(primary) div 2)
    {
      @k++;
      if(size(primary[2*@k])==0)
      {
        ek=insepDecomp_i(1, 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(nestLevel > 1){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)&&(nestLevel > 1))
  {
    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) div 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) div 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)&&(nestLevel > 1))
      {
        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 ============================

      if (npars(basering)>1) { @qh=groebner(@qh,"par2var"); }
      else                   { @qh=groebner(@qh); }

      ser1=phi1(ser);
      @lh=newZero_decomp (@qh,phi(ser1),@wr, list("nest", nestLevel + 1));

      kill lres0;
      list lres0;
      if(size(@lh)==2)
      {
        helpprim=@lh[2];
        lres0[1]=primary[2*@k-1];
        attrib(lres0[1],"isSB",1);
        ser1=psi(helpprim);
        lres0[2]=psi1(ser1);
        if(size(reduce(lres0[2],lres0[1],5))==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)
      {
        if (npars(basering)>1) { lres0[2]=groebner(lres0[2],"par2var"); }
        else                   { lres0[2]=groebner(lres0[2]); }
      }
      else
      {
        for(@n=1;@n<=size(lres0) div 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;

    }
  }
  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= newZero_decomp(i,ideal(0),0);
   pr;
}
///////////////////////////////////////////////////////////////////////////////
static proc check_variables(alias ideal I, int NN)
{
  int i=1;
  ideal J=0;
  ideal JJ;
  poly p;
  J=0;
  i=1;
  while(i<ncols(I))
  {
    p=I[i];
    if (size(variables(p))<NN)
    {
      JJ=J+ideal(p);
      if (size(variables(JJ))<NN)
      { J=JJ; }
    }
    i++;
  }
  if (size(J)>1)
  {"found subsystem:",size(J)," of ",size(I)," vars:",size(variables(J)),size(variables(I));
    return(J);
  }
  return(J);
}

static proc prep_decomp(alias ideal I)
{
  int N=size(variables(I));
  int NN=2;
  ideal J;
  list L;
  while(NN<=N)
  {
    J=check_variables(I,NN);
    if (size(J)>0)
    {
      L=minAssPrimes_i(0,J);
      if(size(L)>1)
      {
        "prep_decomp:",size(L)," NN=",NN," N=",N;
        for(int i=size(L);i>0;i--)
        {
          L[i]=simplify(L[i]+I,4);
        }
        return(L);
      }
    }
    NN++;
  }
  return(list(I));
}

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

*/