Changeset e801fe in git for Singular/LIB/primdec.lib

Timestamp:

Feb 16, 1998, 1:07:04 PM (26 years ago)

Author:

Gerhard Pfister <pfister@…>

Branches:

(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'fc741b6502fd8a97288eaa3eba6e5220f3c3df87')

Children:

1e7d8d9fbdaf1094cc869ca3a6be98b9e1fbfc9a

Parents:

437e2c21158af6efe6402cc1da62ee40760c03e0

Message:

* updated invar.lib normal.lib, merged prim_dec.lib with primdec.lib


git-svn-id: file:///usr/local/Singular/svn/trunk@1145 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• Singular/LIB/primdec.lib (modified) (108 diffs)

Singular/LIB/primdec.lib

@@ -1 @@
-// $Id: primdec.lib,v 1.8 1998-01-27 18:51:23 Singular Exp $
+// $Id: primdec.lib,v 1.9 1998-02-16 12:07:04 pfister Exp $
 ///////////////////////////////////////////////////////
 // primdec.lib
@@ -5 +5 @@
 // the ideas of Gianni,Trager,Zacharias
 // written by Gerhard Pfister
+//
+// algorithms for primary decomposition based on
+// the ideas of Shimoyama/Yokoyama
+// written by Wolfram Decker and Hans Schoenemann
 //////////////////////////////////////////////////////
 
@@ -12 +16 @@
   //computes the minimal associated primes of the ideal I
 
-  decomp (ideal I)
+  primdecGTZ (ideal I)
   // Computes a complete primary decomposition via Gianni,Trager,Zacharias
 
@@ -24 +28 @@
   //computes the radicals of the equidimensional parts of I
 
+  min_ass_prim_charsets (ideal I, int choose)
+  // minimal associated primes via the  characteristic set
+  // package written by Michael Messollen.
+  // The integer choose must be either 0 or 1.
+  // If choose=0, the given ordering of the variables is used.
+  // If choose=1, the system tries to find an "optimal ordering",
+  // which in some cases may considerably speed up the algorithm
+
+  // You may also may want to try one of the algorithms for
+  // minimal associated primes in the the library
+  // primdec.lib, written by Gerhard Pfister.
+  // These algorithms are variants of the algorithm
+  // of Gianni-Trager-Zacharias
+
+  primdecSY (ideal I, int choose)
+
+  // Computes a complete primary decomposition via
+  // a variant of the pseudoprimary approach of
+  // Shimoyama-Yokoyama.
+  // The integer choose must be either 0, 1, 2 or 3.
+  // If choose=0, min_ass_prim_charsets with the given
+  // ordering of the variables is used.
+  // If choose=1, min_ass_prim_charsets with the "optimized"
+  // ordering of the variables is used.
+  // If choose=2, minAssPrimes from primdec.lib is used
+  // If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
+
+LIB "general.lib";
+LIB "elim.lib";
+LIB "poly.lib";
 LIB "random.lib";
-LIB "elim.lib";
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -62 +95 @@
    if(ordstr(basering)[1,2]!="dp")
    {
-      execute "ring @Phelp=("+charstr(@P)+"),("+varstr(@P)+"),dp;";
+      execute "ring @Phelp=("+charstr(@P)+"),("+varstr(@P)+"),(C,dp);";
       ideal inew=std(imap(@P,id));
       ideal  @h=imap(@P,h);
@@ -175 +208 @@
    ideal fac=tsil[2];
    poly f=tsil[3];
-
+   
    ideal star=quotient(laedi,f);
-
    action=1;
 
@@ -190 +222 @@
       {
         g=1;
+        verg=laedi;
+        
          for(j=1;j<=size(fac);j++)
          {
@@ -198 +232 @@
          }
          verg=quotient(laedi,g);
+         
          if(specialIdealsEqual(verg,star)==1)
          {
@@ -211 +246 @@
       }
    }
-   l=star,fac,f;
+   l=star,fac,f;   
    return(l);
 }
-
 
 ////////////////////////////////////////////////////////////////////////////////
 proc testFactor(list act,poly p)
 {
+  poly keep=p;
+  
    int i;
    poly q=act[1][1]^act[2][1];
@@ -230 +266 @@
    {
       "ERROR IN FACTOR";
+      basering;
+      
       act;
+      keep;
+      pause;
+      
       p;
       q;
@@ -274 +315 @@
      return(@l);
   }
-  @l=factorize(p,2);
-  if(npars(basering)>0)
-  {
+//  @l=factorize(p,2);
+  @l=factorize(p);
+  // if(npars(basering)>0)
+  // {
      for(@j=1;@j<=size(@l[1]);@j++)
      {
@@ -320 +362 @@
         }
      }
-  }
+     // }
   return(@l);
 }
@@ -381 +423 @@
       }
       attrib(k2,"isSB",1);
-      if(size(reduce(k1,k2))==0)
+      if(size(reduce(k1,k2,1))==0)
       {
          return(1);
@@ -514 +556 @@
          m=m-1;
       }
-      //check whether i[m] =(c*var(n)+h)^e modulo prm for some
+      //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
@@ -524 +566 @@
          i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
 
-         if (reduce(i[m]-t^e,prm) !=0)
+         if (reduce(i[m]-t^e,prm,1) !=0)
          {
            return(ideal(0));
@@ -606 +648 @@
    {
       i++;
-      if((size(ser)>0)&&(size(reduce(ser,l[2*i-1]))==0))
+      if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==0))
       {
          l[2*i-1]=ideal(1);
@@ -736 +778 @@
    for(i=1;i<=size(l)/2;i++)
    {
-      if((size(l[2*i])==0)&&(specialIdealsEqual(keepprime[i],l[2*i-1])!=1))
+     attrib(l[2*i-1],"isSB",1);
+     
+     if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==0)&&(deg(l[2*i-1][1])>0))
+     {
+       "Achtung in split";
+       
+         l[2*i-1]=ideal(1);
+         l[2*i]=ideal(1);
+     }
+     if((size(l[2*i])==0)&&(specialIdealsEqual(keepprime[i],l[2*i-1])!=1))
       {
          keepprime[i]=std(keepprime[i]);
@@ -782 +833 @@
   def   @P = basering;
   int nva = nvars(basering);
-  int @k,@s,@n,@k1;
-  list primary,lres,act,@lh,@wh;
-  map phi,psi;
-  ideal jmap,helpprim,@qh,@qht;
+  int @k,@s,@n,@k1,zz;
+  list primary,lres,lres1,act,@lh,@wh;
+  map phi,psi,phi1,psi1;
+  ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1;
   intvec @vh,@hilb;
   string @ri;
@@ -796 +847 @@
      return(primary);
   }
-  j=interred(j);
+  
+  j=interred(j);  
   attrib(j,"isSB",1);
   if(vdim(j)==deg(j[1]))
-  {
-     if((size(ser)>0)&&(size(reduce(ser,j))==0))
-     {
-          primary[1]=ideal(1);
-          primary[2]=ideal(1);
-          return(primary);
-     }
+  {    
      act=factor(j[1]);
      for(@k=1;@k<=size(act[1]);@k++)
@@ -822 +868 @@
        @qh[1]=act[1][@k];
        primary[2*@k]=interred(@qh);
-     }
-        return(primary);
+       attrib( primary[2*@k-1],"isSB",1);
+       
+       if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],1))==0))
+       {
+          primary[2*@k-1]=ideal(1);
+          primary[2*@k]=ideal(1);          
+       }
+     }
+     return(primary);
   }
 
@@ -829 +882 @@
   {
      primary[1]=j;
-     if((size(ser)>0)&&(size(reduce(ser,j))==0))
+     if((size(ser)>0)&&(size(reduce(ser,j,1))==0))
      {
           primary[1]=ideal(1);
@@ -894 +947 @@
   }
   else
-  {
+  {  
      primary[1]=j;
      if((size(#)>0)&&(act[2][1]>1))
@@ -913 +966 @@
   if(size(#)==0)
   {
+    
      primary=splitPrimary(primary,ser,@wr,act);
+     
   }
 
@@ -922 +977 @@
 
   @k=0;
-  int zz;
   while(@k<(size(primary)/2))
   {
@@ -937 +991 @@
     }
   }
+    
   @k=0;
+  ideal keep;
   while(@k<(size(primary)/2))
   {
@@ -943 +999 @@
     if (size(primary[2*@k])==0)
     {
-//      "the univariat polynomials";
-//       int qwe=timer;
-//       system("finduni",primary[2*@k-1]);
 
        jmap=randomLast(100);
-//       timer-qwe;
-//       primary[2*@k-1];
-//       pause;
+       jmap1=maxideal(1);
+       jmap2=maxideal(1);
+       @qht=primary[2*@k-1];
+
        for(@n=2;@n<=size(primary[2*@k-1]);@n++)
        {
           if(deg(lead(primary[2*@k-1][@n]))==1)
           {
+             for(zz=1;zz<=nva;zz++)
+             {
+                if(lead(primary[2*@k-1][@n])/var(zz)!=0)
+                {
+                   jmap1[zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]
+                   +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]);
+                    jmap2[zz]=primary[2*@k-1][@n];
+                    @qht[@n]=var(zz);
+                     
+                }
+             }
              jmap[nva]=subst(jmap[nva],lead(primary[2*@k-1][@n]),0);
           }
        }
+       if(size(subst(jmap[nva],var(1),0)-var(nva))!=0)
+       {
+          jmap[nva]=subst(jmap[nva],var(1),0);
+       }
+       phi1=@P,jmap1;
        phi=@P,jmap;
-       jmap[nva]=-(jmap[nva]-2*var(nva));
+ 
+       for(@n=1;@n<=nva;@n++)
+       {
+          jmap[@n]=-(jmap[@n]-2*var(@n));
+       }
        psi=@P,jmap;
-       @qht=primary[2*@k-1];
+       psi1=@P,jmap2;
+
        @qh=phi(@qht);
+       
        if(npars(@P)>0)
        {
           @ri= "ring @Phelp ="
-                  +string(char(@P))+",("+varstr(@P)+","+parstr(@P)+",@t),dp;";
+                  +string(char(@P))+",
+                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
        }
        else
        {
           @ri= "ring @Phelp ="
-                  +string(char(@P))+",("+varstr(@P)+",@t),dp;";
+                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
        }
        execute(@ri);
@@ -979 +1056 @@
        @hilb=hilb(@qh1,1);
        @ri= "ring @Phelp1 ="
-                  +string(char(@P))+",("+varstr(@Phelp)+"),lp;";
+                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
        execute(@ri);
        ideal @qh=homog(imap(@P,@qh),@t);
@@ -989 +1066 @@
        kill @Phelp1;
        @qh=clearSB(@qh);
-       attrib(@qh,"isSB",1);
-
-       @lh=zero_decomp (@qh,psi(ser),@wr);
-
+       attrib(@qh,"isSB",1);       
+       ser1=phi1(ser);
+       @lh=zero_decomp (@qh,phi(ser1),@wr);
+//       @lh=zero_decomp (@qh,psi(ser),@wr);
+       
+       
        kill lres;
        list lres;
@@ -999 +1078 @@
           helpprim=@lh[2];
           lres[1]=primary[2*@k-1];
-          lres[2]=psi(helpprim);
-          if(size(reduce(lres[2],lres[1]))==0)
+          ser1=psi(helpprim);
+          lres[2]=psi1(ser1);
+          if(size(reduce(lres[2],lres[1],1))==0)
           {
              primary[2*@k]=primary[2*@k-1];
@@ -1014 +1094 @@
              {
                 @f=act[1][@n]^act[2][@n];
-                lres[2*@n-1]=interred(primary[2*@k-1]+psi(@f));
+                ser1=psi(@f);
+                lres[2*@n-1]=interred(primary[2*@k-1]+psi1(ser1));
                 helpprim=@lh[2*@n];
-                lres[2*@n]=psi(helpprim);
+                ser1=psi(helpprim);
+                lres[2*@n]=psi1(ser1);
              }
           }
           else
           {
-             lres=psi(@lh);
+             lres1=psi(@lh);
+             lres=psi1(lres1);
           }
        }
@@ -1027 +1110 @@
        {
           @ri= "ring @Phelp ="
-                  +string(char(@P))+",("+varstr(@P)+","+parstr(@P)+",@t),dp;";
+                  +string(char(@P))+",
+                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
        }
        else
        {
           @ri= "ring @Phelp ="
-                  +string(char(@P))+",("+varstr(@P)+",@t),dp;";
+                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
        }
        execute(@ri);
@@ -1069 +1153 @@
        }
        @ri= "ring @Phelp1 ="
-                  +string(char(@P))+",("+varstr(@Phelp)+"),lp;";
+                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
        execute(@ri);
        list @lr=imap(@Phelp,@lr);
@@ -1156 +1240 @@
         return(1)
      }
-     p=gcd(p,i[k]);
+     p=GCD(p,i[k]);
      if(deg(p)==0)
      {
@@ -1204 +1288 @@
    def bsr= basering;
    string @id=string(h);
-   execute "ring @r=0,("+@pa+","+varstr(bsr)+"),dp;";
+   execute "ring @r=0,("+@pa+","+varstr(bsr)+"),(C,dp);";
    execute "ideal @i="+@id+";";
    poly @p=lcmP(@i);
@@ -1250 +1334 @@
      if(deg(i[k])!=0)
      {
-        q=gcd(p,i[k]);
+        q=GCD(p,i[k]);
         if(deg(q)==0)
         {
@@ -1514 +1598 @@
     @jh=@jh+var(@n);
   }
-  quotring=quotring+"),lp;";
+  quotring=quotring+"),(C,lp);";
 
   return(quotring);
@@ -1577 +1661 @@
    #[1]=1;
    def @P=basering;
+   list qr=simplifyIdeal(i);
+   map phi=@P,qr[2];
+   i=qr[1];
+   
    execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
              +ordstr(basering)+");";
@@ -1593 +1681 @@
       setring @P;
       list @res=imap(gnir,tluser);
-      return(@res);
+      return(phi(@res));
    }
    list @res,empty;
@@ -1606 +1694 @@
          setring @P;
          list @res=maxideal(1);
-         return(@res);
+         return(phi(@res));
       }
       if(dim(@pr[1])>1)
       {
          setring @P;
-         kill gnir;
-         execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
+        // kill gnir;
+         execute "ring gnir1 = ("+charstr(basering)+"),
+                              ("+varstr(basering)+"),(C,lp);";
          ideal i=fetch(@P,i);
          list @pr=facstd(i);
-         ideal ser;
-      }
-   }
-   option(noredSB);
+        // ideal ser;
+         setring gnir;
+         @pr=fetch(gnir1,@pr);
+         kill gnir1;
+      }
+   }
+    option(noredSB);
    int j,k,odim,ndim,count;
    attrib(@pr[1],"isSB",1);
@@ -1662 +1754 @@
      for(j=1;j<=size(@pr);j++)
      {
-         @res[j]=decomp(@pr[j],2);
- //      @res[j]=decomp(@pr[j],2,@pr[j],ser);
- //      for(k=1;k<=size(@res[j]);k++)
- //      {
- //         ser=intersect1(ser,@res[j][k]);
- //      }
+//@pr[j];
+//pause;
+        @res[j]=decomp(@pr[j],2);
+//       @res[j]=decomp(@pr[j],2,@pr[j],ser);
+//       for(k=1;k<=size(@res[j]);k++)
+//       {
+//          ser=intersect(ser,@res[j][k]);
+//       }
      }
    }
@@ -1674 +1768 @@
    setring @P;
    list @res=imap(gnir,@res);
-   return(@res);
+   return(phi(@res));
 }
 example
@@ -1696 +1790 @@
   def P=basering;
 
-  execute "ring ir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
+  execute "ring ir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);";
   list l=fetch(P,li);
   list @erg;
@@ -1735 +1829 @@
         if(size(reduce(i1,i2,1))==0)
         {
-           if(size(reduce(@erg[i],@erg[k]))==0)
+           if(size(reduce(@erg[i],@erg[k],1))==0)
            {
               @erg[k]=ideal(1);
@@ -1743 +1837 @@
         if(size(reduce(i2,i1,1))==0)
         {
-           if(size(reduce(@erg[k],@erg[i]))==0)
+           if(size(reduce(@erg[k],@erg[i],1))==0)
            {
               break;
@@ -1780 +1874 @@
    return(radical(i,1));
 }
+
 ///////////////////////////////////////////////////////////////////////////////
-proc decomp (ideal i,list #)
+proc decomp(ideal i,list #)
 USAGE:   decomp(i); i ideal  (for primary decomposition)   (resp.
          decomp(i,1);        (for the minimal associated primes) )
@@ -1793 +1888 @@
   list primary,indep,ltras;
   intvec @vh,isat;
-  int @wr,@k,@n,@m,@n1,@n2,@n3,homo;
+  int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi;
   ideal peek=i;
   ideal ser,tras;
 
-  int @aa=timer;
-
-  homo=homog(i);
   if(size(#)>0)
   {
@@ -1807 +1899 @@
         if(size(#)>1)
         {
+          seri=1;
           peek=#[2];
           ser=#[3];
@@ -1813 +1906 @@
       else
       {
-         peek=#[1];
-         ser=#[2];
-      }
-  }
-
+        seri=1;
+        peek=#[1];
+        ser=#[2];
+      }
+  }
+  
+  homo=homog(i);
+  
   if(homo==1)
   {
-     ltras=mstd(i);
-     attrib(ltras[1],"isSB",1);
-     tras=ltras[1];
-     if(dim(tras)==0)
-     {
+    if(attrib(i,"isSB")!=1)
+    {
+      ltras=mstd(i);
+      attrib(ltras[1],"isSB",1);
+    }
+    else
+    {
+      ltras=i,i;
+    }
+    tras=ltras[1];
+    if(dim(tras)==0)
+    {
         primary[1]=ltras[2];
         primary[2]=maxideal(1);
@@ -1836 +1939 @@
         return(primary);
      }
-      intvec @hilb=hilb(tras,1);
+     intvec @hilb=hilb(tras,1);
+     intvec keephilb=@hilb;
   }
 
@@ -1853 +1957 @@
   //----------------------------------------------------------------
 
-  execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
-
+  execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);";
   option(redSB);
-  ideal ser=fetch(@P,ser);
-  ideal peek=std(fetch(@P,peek));
-  homo=homog(peek);
+
+  ideal ser=fetch(@P,ser);  
 
   if(homo==1)
   {
-     if(ordstr(@P)[1,2]!="lp")
-     {
-        ideal @j=std(fetch(@P,i),@hilb);
+     if((ordstr(@P)[1]!="(C,lp)")&&(ordstr(@P)[3]!="(C,lp)"))
+     {
+       ideal @j=std(fetch(@P,i),@hilb);
      }
      else
@@ -1876 +1978 @@
      ideal @j=std(fetch(@P,i));
   }
-
+  option(noredSB);
+  if(seri==1)
+  {
+    ideal peek=fetch(@P,peek);
+    attrib(peek,"isSB",1);
+  }
+  else
+  {
+    ideal peek=@j;
+  }
+  if(size(ser)==0)
+  {
+    ideal fried;
+    @n=size(@j);
+    for(@k=1;@k<=@n;@k++)
+    {
+      if(deg(lead(@j[@k]))==1)
+      {
+        fried[size(fried)+1]=@j[@k];
+        @j[@k]=0;
+      }
+    }
+    if(size(fried)>0)
+    {
+       @j=simplify(@j,2);
+       attrib(@j,"isSB",1);
+       list pr=decomp(@j);
+       for(@k=1;@k<=size(pr);@k++)
+       {
+         @j=pr[@k]+fried;
+         pr[@k]=@j;
+       }
+       setring @P;
+       return(fetch(gnir,pr));
+    }
+  }
+ 
   //----------------------------------------------------------------
   //j is the ring
@@ -1883 +2021 @@
   if (dim(@j)==-1)
   {
-     setring @P;
-     option(noredSB);
-     return(ideal(0));
+    setring @P;
+    return(ideal(0));
   }
 
@@ -1911 +2048 @@
      }
      setring @P;
-     option(noredSB);
      primary=fetch(gnir,gprimary);
 
      return(primary);
   }
-
+  
  //------------------------------------------------------------------
  //the zero-dimensional case
@@ -1923 +2059 @@
   if (dim(@j)==0)
   {
-      list gprimary= zero_decomp(@j,ser,@wr);
-
-      setring @P;
-      option(noredSB);
-      primary=fetch(gnir,gprimary);
-      if(size(ser)>0)
-      {
-         primary=cleanPrimary(primary);
-      }
-      return(primary);
-   }
+    option(redSB);
+    list gprimary= zero_decomp(@j,ser,@wr);
+    option(noredSB);
+    setring @P;
+    primary=fetch(gnir,gprimary);
+    if(size(ser)>0)
+    {
+      primary=cleanPrimary(primary);
+    }
+    return(primary);
+  }
 
   poly @gs,@gh,@p;
@@ -1941 +2077 @@
   int jdim=dim(@j);
   list fett;
-  int lauf,di;
+  int lauf,di,newtest;
   //------------------------------------------------------------------
   //search for a maximal independent set indep,i.e.
@@ -1970 +2106 @@
       indep=maxIndependSet(@j);
    }
-
+  
   ideal jkeep=@j;
 
@@ -1979 +2115 @@
   else
   {
-     execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),dp;";
-  }
-  ideal jwork=std(imap(gnir,@j));
+     execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),(C,dp);";
+  }
+
+  if(homo==1)
+  {
+    if((ordstr(@P)[3]=="d")||(ordstr(@P)[1]=="d")||(ordstr(@P)[1]=="w")
+       ||(ordstr(@P)[3]=="w"))
+    {
+      ideal jwork=imap(@P,tras);
+      attrib(jwork,"isSB",1);
+    }
+    else
+    {
+      ideal jwork=std(imap(gnir,@j),@hilb);
+    }
+    
+  }
+  else
+  {
+    ideal jwork=std(imap(gnir,@j));
+  }
+  list hquprimary;
   poly @p,@q;
-  ideal @h,fac;
+  ideal @h,fac,ser;
   di=dim(jwork);
+  keepdi=di;
+  
   setring gnir;
   for(@m=1;@m<=size(indep);@m++)
@@ -1990 +2147 @@
      isat=0;
      @n2=0;
+     option(redSB);
      if((indep[@m][1]==varstr(basering))&&(@m==1))
      //this is the good case, nothing to do, just to have the same notations
@@ -1998 +2156 @@
         ideal @j=fetch(gnir,@j);
         attrib(@j,"isSB",1);
+        ideal ser=fetch(gnir,ser);
+        
      }
      else
@@ -2013 +2173 @@
           ideal @j=std(phi(@j));
         }
-      }
+        ideal ser=phi(ser);
+        
+     }
+     option(noredSB);     
      if((deg(@j[1])==0)||(dim(@j)<jdim))
      {
@@ -2050 +2213 @@
     // @j considered in the quotientring
      ideal @j=imap(gnir1,@j);
-     ideal ser=imap(gnir,ser);
+     ideal ser=imap(gnir1,ser);
 
      kill gnir1;
@@ -2069 +2232 @@
         }
         //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
-
+        option(redSB);        
         list uprimary= zero_decomp(@j,ser,@wr);
-
+        option(noredSB);
      }
      else
      {
        list uprimary;
+       uprimary[1]=ideal(1);
        uprimary[2]=ideal(1);
-       uprimary[1]=ideal(1);
      }
 
@@ -2130 +2293 @@
      //here the intersection with the polynomialring
      //mentioned above is really computed
-
-    for(@n=@n3/2+1;@n<=@n2/2;@n++)
-     {
+     for(@n=@n3/2+1;@n<=@n2/2;@n++)
+     {       
         if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
         {
@@ -2147 +2309 @@
         }
      }
+    
      if(size(@h)>0)
      {
@@ -2155 +2318 @@
         @h=imap(gnir,@h);
         if(@wr!=1)
-//        if(@wr==0)
-        {
-           @q=minSat(jwork,@h)[2];
+        {
+          @q=minSat(jwork,@h)[2];   
         }
         else
@@ -2177 +2339 @@
         }
         jwork=std(jwork,@q);
-        if(dim(jwork)<di)
+        keepdi=dim(jwork);        
+        if(keepdi<di)
         {
            setring gnir;
@@ -2195 +2358 @@
   {
      @j=ideal(1);
+     quprimary[1]=ideal(1);
      quprimary[2]=ideal(1);
-     quprimary[1]=ideal(1);
-  }
-
+  }
+  if((size(quprimary)==0))
+  {
+    keepdi=di-1;
+  }  
   //---------------------------------------------------------------
   //notice that j=sat(j,gh) intersected with (j,gh^n)
@@ -2205 +2371 @@
   if((deg(@j[1])!=0)&&(@wr!=1))
   {
-     int uq=size(quprimary);
-     if(uq>0)
-     {
+     if(size(quprimary)>0)
+     {
+        setring @Phelp;
+        ser=imap(gnir,ser);
+        hquprimary=imap(gnir,quprimary);
         if(@wr==0)
         {
-           ideal htest=quprimary[1];
-
-           for (@n1=2;@n1<=size(quprimary)/2;@n1++)
+           ideal htest=hquprimary[1];
+           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
            {
-              htest=intersect(htest,quprimary[2*@n1-1]);
+              htest=intersect(htest,hquprimary[2*@n1-1]);
            }
         }
         else
         {
-           ideal htest=quprimary[2];
-
-           for (@n1=2;@n1<=size(quprimary)/2;@n1++)
+           ideal htest=hquprimary[2];
+
+           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
            {
-              htest=intersect(htest,quprimary[2*@n1]);
+              htest=intersect(htest,hquprimary[2*@n1]);
            }
         }
+
         if(size(ser)>0)
-        {
-           htest=intersect(htest,ser);
-        }
-        ser=std(htest);
-     }
-     //we are not ready yet
-     if (specialIdealsEqual(ser,peek)!=1)
-     {
+        { 
+           ser=intersect(htest,ser); 
+        }
+        else
+        {
+          ser=htest;
+        } 
+        setring gnir;
+        ser=imap(@Phelp,ser);
+     }
+     if(size(reduce(ser,peek,1))!=0)
+     {       
         for(@m=1;@m<=size(restindep);@m++)
         {
+         // if(restindep[@m][3]>=keepdi)
+         // { 
            isat=0;
            @n2=0;
+           option(redSB);
+           
            if(restindep[@m][1]==varstr(basering))
            //this is the good case, nothing to do, just to have the same notations
@@ -2256 +2432 @@
               if(homo==1)
               {
-                 ideal @j=std(phi(jkeep),@hilb);
+                 ideal @j=std(phi(jkeep),keephilb);
               }
               else
@@ -2262 +2438 @@
                 ideal @j=std(phi(jkeep));
               }
+              ideal ser=phi(ser);
            }
-
+           option(noredSB);
+           
            for (lauf=1;lauf<=size(@j);lauf++)
            {
@@ -2295 +2473 @@
            // @j considered in the quotientring
            ideal @j=imap(gnir1,@j);
-           ideal ser=imap(gnir,ser);
+           ideal ser=imap(gnir1,ser);
 
            kill gnir1;
@@ -2312 +2490 @@
            }
            //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
-
-            list uprimary= zero_decomp(@j,ser,@wr);
-
+           
+           option(redSB);
+           list uprimary= zero_decomp(@j,ser,@wr);
+           option(noredSB);
+           
+            
            //we need the intersection of the ideals in the list quprimary with the
            //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
@@ -2349 +2530 @@
            @n2=size(quprimary);
            @n3=@n2;
-
+           
            for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
            {
@@ -2362 +2543 @@
               }
            }
-
+           
+           
            //here the intersection with the polynomialring
            //mentioned above is really computed
@@ -2380 +2562 @@
                  }
                  quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
-             }
-             if(@wr==0)
-             {
-                ser=std(intersect(ser,quprimary[2*@n-1]));
-             }
-             else
-             {
-                ser=std(intersect(ser,quprimary[2*@n]));
-             }
+              }
            }
-        }
-        //we are not ready yet
-        if (specialIdealsEqual(ser,peek)!=1)
+           if(@n2>=@n3+2)
+           {
+              setring @Phelp;
+              ser=imap(gnir,ser);
+              hquprimary=imap(gnir,quprimary);
+              for(@n=@n3/2+1;@n<=@n2/2;@n++)
+              {
+                if(@wr==0)
+                {
+                   ser=intersect(ser,hquprimary[2*@n-1]);
+                }
+                else
+                {
+                   ser=intersect(ser,hquprimary[2*@n]);
+                }
+              }
+              setring gnir;
+              ser=imap(@Phelp,ser);
+           }
+           
+         // }
+        }             
+        if(size(reduce(ser,peek,1))!=0)
         {
            if(@wr>0)
@@ -2404 +2598 @@
            // 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++)
@@ -2412 +2605 @@
         }
      }
+     
    }
   //------------------------------------------------------------
@@ -2417 +2611 @@
   //the final result: primary
   //------------------------------------------------------------
+  
   setring @P;
   primary=imap(gnir,quprimary);
-
-  option(noredSB);
   return(primary);
 }
@@ -2438 +2631 @@
 
 ///////////////////////////////////////////////////////////////////////////////
+proc primdecGTZ(ideal i)
+{
+   return(decomp(i));
+}
+///////////////////////////////////////////////////////////////////////////////
+proc primdecSY(ideal i,int j)
+{
+   return(prim_dec(i,j));
+}
+///////////////////////////////////////////////////////////////////////////////
+
 proc radical(ideal i,list #)
 {
@@ -2444 +2648 @@
    if(size(#)>0)
    {
-     il=#[1];
+     il=#[1]; 
    }
    ideal re=1;
    option(redSB);
+   list qr=simplifyIdeal(i);
+   map phi=@P,qr[2];
+   i=qr[1];
+   
    list pr=facstd(i);
 
@@ -2456 +2664 @@
       {
          ideal @res=maxideal(1);
-         return(@res);
+         return(phi(@res));
       }
       if(dim(pr[1])>1)
       {
-         execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
+         execute "ring gnir = ("+charstr(basering)+"),
+                              ("+varstr(basering)+"),(C,lp);";
          ideal i=fetch(@P,i);
          list @pr=facstd(i);
@@ -2469 +2678 @@
    option(noredSB);
    int s=size(pr);
+
    if(s==1)
    {
-      return(radicalEHV(i,ideal(1),il));
+     i=radicalEHV(i,ideal(1),il);
+     return(phi(i));
    }
    intvec pos;
    pos[s]=0;
-
    if(il==1)
    {
@@ -2482 +2692 @@
      int odim=dim(pr[1]);
      int count=1;
-
+   
      for(j=2;j<=s;j++)
      {
@@ -2502 +2712 @@
      }
    }
-
    for(j=1;j<=s;j++)
    {
       if(pos[j]==0)
       {
-         re=intersect1(re,radicalEHV(pr[s+1-j],re,il));
-      }
-   }
-   return(re);
+         re=intersect(re,radicalEHV(pr[s+1-j],re,il));
+      }
+   }
+   return(phi(re));
 }
 
@@ -2516 +2725 @@
 {
    def R=basering;
-   execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+",@t),dp;";
+   execute "ring gnir = ("+charstr(basering)+"),
+                        ("+varstr(basering)+",@t),(C,dp);";
    ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))-1)*imap(R,j);
    ideal j=eliminate(i,var(nvars(basering)));
@@ -2523 +2733 @@
    return(phi(j));
 }
-
+ 
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -2544 +2754 @@
       return(ideal(0));
    }
-
+   
    def  @P = basering;
    list indep,allindep,restindep,fett,@mu;
@@ -2552 +2762 @@
    ideal rad=ideal(1);
    ideal te=ser;
+
    poly @p,@q;
    string @va,quotring;
@@ -2563 +2774 @@
    @j1=m[2];
    int jdim=dim(@j);
-   if(size(reduce(ser,@j))==0)
+   if(size(reduce(ser,@j,1))==0)
    {
       return(ser);
@@ -2595 +2806 @@
    {
       fac=factorize(@j[1],1);
-      poly @p=1;
+      @p=1;
       for(@k=1;@k<=size(fac);@k++)
       {
@@ -2603 +2814 @@
 
       return(ideal(@p));
-   }
+   }   
    //------------------------------------------------------------------
    //the case of a complete intersection
@@ -2619 +2830 @@
    if (jdim==0)
    {
-      @j1=system("finduni",@j);
+      @j1=system("finduni",@j); 
       for(@k=1;@k<=size(@j1);@k++)
       {
@@ -2634 +2845 @@
       return(@j);
    }
-
+  
    //------------------------------------------------------------------
    //search for a maximal independent set indep,i.e.
@@ -2643 +2854 @@
 
    indep=maxIndependSet(@j);
-
+  
    di=dim(@j);
 
@@ -2772 +2983 @@
       {
          fac=ideal(0);
-         for(lauf=1;lauf<=ncols(@h);lauf++)
+         for(lauf=1;lauf<=ncols(@h);lauf++)   
          {
             if(deg(@h[lauf])>0)
@@ -2786 +2997 @@
          }
 
-
+       
          @mu=mstd(quotient(@j+ideal(@q),rad));
          @j=@mu[1];
          attrib(@j,"isSB",1);
-
+       
       }
       if((deg(rad[1])>0)&&(deg(collectrad[1])>0))
       {
-int xyz=timer;
-"bei collecterad";
          rad=intersect(rad,collectrad);
-timer-xyz;
       }
       else
@@ -2808 +3016 @@
 
       te=simplify(reduce(te*rad,@j),2);
-
+ 
       if((dim(@j)<di)||(size(te)==0))
       {
@@ -2822 +3030 @@
    {
       return(rad);
-   }
-   ideal tes=radicalKL(@mu,rad,@wr);
-int sml=timer;
-"bei rad";
-   rad=intersect(rad,tes);
-timer-sml;
-  // rad=intersect(rad,radicalKL(@mu,ideal(1),@wr));
+   }   
+  // rad=intersect(rad,radicalKL(@mu,rad,@wr));
+   rad=intersect(rad,radicalKL(@mu,ideal(1),@wr));
 
 
@@ -2849 +3053 @@
       il=#[1];
    }
+   
    option(redSB);
    list m=mstd(i);
         I=m[2];
    option(noredSB);
-   if(size(reduce(re,m[1]))==0)
+   if(size(reduce(re,m[1],1))==0)
    {
       return(re);
    }
    int cod=nvars(basering)-dim(m[1]);
-   if(nvars(basering)<9)
+   if((nvars(basering)<=5)&&(size(m[2])<=5))
    {
       if(cod==size(m[2]))
       {
-         J=minor(jacob(I),cod);
-         return(quotient(I,J));
-      }
-
+        J=minor(jacob(I),cod);
+        return(quotient(I,J));
+      }
+   
       for(l=1;l<=cod;l++)
       {
@@ -2877 +3082 @@
          radI1=quotient(radI0,L);
 
-         if(size(reduce(radI1,m[1]))==0)
+         if(size(reduce(radI1,m[1],1))==0)
          {
             return(I);
@@ -2883 +3088 @@
          if(il==1)
          {
+      
             return(radI1);
          }
@@ -2892 +3098 @@
             return(radI1);
          }
-         return(intersect(radI1,radicalEHV(I2,re,il)));
+         return(intersect(radI1,radicalEHV(I2,re,il)));       
       }
    }
@@ -2901 +3107 @@
 {
   M=prune(M);  //to obtain a small embedding
-  return(quotient(M,freemodule(nrows(M))));
+  ideal ann=quotient1(M,freemodule(nrows(M)));
+  return(ann);
 }
 
@@ -2919 +3126 @@
   }
   return(ann);
-}
-
+} 
+ 
 //computes all equidimensional parts of the ideal i
 proc prepareAss(ideal i)
@@ -2937 +3144 @@
   {
      list re=mres(i,0);             //fehler in sres
-  }
+  }  
   for(e=cod;e<=nvars(basering);e++)
   {
      ann=AnnExt_R(e,re);
+
      if(nvars(basering)-dim(std(ann))==e)
      {
@@ -2947 +3155 @@
   }
   return(er);
-}
+}  
 
 //computes the annihilator of Ext^n(R/i,R) with given resolution re
@@ -2958 +3166 @@
      module k=res(f,2)[2];             //the kernel
      matrix g=transpose(re[n]);        //the image of Hom(_,R)
+
      ideal ann=quotient(g,k);           //the anihilator
   }
@@ -2964 +3173 @@
      ideal ann=Ann(transpose(re[n]));
   }
-  return(ann);
-}
-
+  return(ann);            
+}
+
+proc quotient1(module a,module b)
+{
+   int i;
+   ideal re=quotient(a,module(b[1]));
+   for(i=2;i<=size(b);i++)
+   {
+      re=intersect1(re,quotient(a,module(b[i])));
+   }
+   return(re);   
+}
+
+
+
+proc analyze(list pr)
+{ 
+   int ii,jj;
+   for(ii=1;ii<=size(pr)/2;ii++)
+   {
+      dim(std(pr[2*ii]));
+      idealsEqual(pr[2*ii-1],pr[2*ii]);
+      "===========================";
+   }
+
+   for(ii=size(pr)/2;ii>1;ii--)
+   {
+      for(jj=1;jj<ii;jj++)
+      {
+         if(size(reduce(pr[2*jj],std(pr[2*ii],1)))==0)
+         {
+            "eingebette Komponente";
+            jj;
+            ii;
+         }
+      }
+   }   
+}
+
+
+proc simplifyIdeal(ideal i)
+{
+  def r=basering;
+  
+  int j,k;
+  map phi;
+  poly p;
+  
+  ideal iwork=i;
+  ideal imap1=maxideal(1);
+  ideal imap2=maxideal(1);
+  
+
+  for(j=1;j<=nvars(basering);j++)
+  {
+    for(k=1;k<=size(i);k++)
+    {
+      if(deg(iwork[k]/var(j))==0)
+      {
+        p=-1/leadcoef(iwork[k]/var(j))*iwork[k];
+        imap1[j]=p+2*var(j);
+        phi=r,imap1;
+        iwork=phi(iwork);
+        iwork=subst(iwork,var(j),0);
+        iwork[k]=var(j);
+        imap1=maxideal(1);
+        imap2[j]=-p;        
+        break;
+      }
+    }
+  }
+  return(iwork,imap2);
+}
+
+        
+///////////////////////////////////////////////////////
+// ini_mod
+// input: a polynomial p
+// output: the initial term of p as needed
+// in the context of characteristic sets
+//////////////////////////////////////////////////////
+
+proc ini_mod(poly p)
+{
+  if (p==0)
+  {
+    return(0);
+  }
+  int n; matrix m;
+  for( n=nvars(basering); n>0; n=n-1)
+  {
+    m=coef(p,var(n));
+    if(m[1,1]!=1)
+    {
+      p=m[2,1];
+      break;
+    }
+  }
+  if(deg(p)==0)
+  {
+    p=0;
+  }
+  return(p);
+}
+///////////////////////////////////////////////////////
+// min_ass_prim_charsets
+// input: generators of an ideal PS and an integer cho
+// If cho=0, the given ordering of the variables is used.
+// Otherwise, the system tries to find an "optimal ordering",
+// which in some cases may considerably speed up the algorithm
+// output: the minimal associated primes of PS
+// algorithm: via characteriostic sets
+//////////////////////////////////////////////////////
+
+
+proc min_ass_prim_charsets (ideal PS, int cho)
+{
+  if((cho<0) and (cho>1))
+    {
+      "ERROR: <int> must be 0 or 1"
+      return();
+    }
+  if(system("version")>933)
+  {
+    option(notWarnSB);
+  }
+  if(cho==0)
+  {
+    return(min_ass_prim_charsets0(PS));
+  }
+  else
+  {
+     return(min_ass_prim_charsets1(PS));
+  }
+}
+///////////////////////////////////////////////////////
+// min_ass_prim_charsets0
+// input: generators of an ideal PS
+// output: the minimal associated primes of PS
+// algorithm: via characteristic sets
+// the given ordering of the variables is used
+//////////////////////////////////////////////////////
+
+
+proc min_ass_prim_charsets0 (ideal PS)
+{
+
+  matrix m=char_series(PS);  // We compute an irreducible
+                             // characteristic series
+  int i,j,k;
+  list PSI;
+  list PHI;  // the ideals given by the characteristic series
+  for(i=nrows(m);i>=1; i--)
+  {
+     PHI[i]=ideal(m[i,1..ncols(m)]);
+  }
+  // We compute the radical of each ideal in PHI
+  ideal I,JS,II;
+  int sizeJS, sizeII;
+  for(i=size(PHI);i>=1; i--)
+  {
+     I=0;
+     for(j=size(PHI[i]);j>0;j=j-1)
+     {
+       I=I+ini_mod(PHI[i][j]);
+     }
+     JS=std(PHI[i]);
+     sizeJS=size(JS);
+     for(j=size(I);j>0;j=j-1)
+     {
+       II=0;
+       sizeII=0;
+       k=0;
+       while(k<=sizeII)                  // successive saturation
+       {
+         option(returnSB);
+         II=quotient(JS,I[j]);
+         option(noreturnSB);
+//std
+//         II=std(II);
+         sizeII=size(II);
+         if(sizeII==sizeJS)
+         {
+            for(k=1;k<=sizeII;k++)
+            {
+              if(leadexp(II[k])!=leadexp(JS[k])) break;
+            }
+         }
+         JS=II;
+         sizeJS=sizeII;
+       }
+    }
+    PSI=insert(PSI,JS);
+  }
+  int sizePSI=size(PSI);
+  // We eliminate redundant ideals
+  for(i=1;i<sizePSI;i++)
+  {
+    for(j=i+1;j<=sizePSI;j++)
+    {
+      if(size(PSI[i])!=0)
+      {
+        if(size(PSI[j])!=0)
+        {
+          if(size(NF(PSI[i],PSI[j],1))==0)
+          {
+            PSI[j]=ideal(0);
+          }
+          else
+          {
+            if(size(NF(PSI[j],PSI[i],1))==0)
+            {
+              PSI[i]=ideal(0);
+            }
+          }
+        }
+      }
+    }
+  }
+  for(i=sizePSI;i>=1;i--)
+  {
+    if(size(PSI[i])==0)
+    {
+      PSI=delete(PSI,i);
+    }
+  }
+  return (PSI);
+}
+
+///////////////////////////////////////////////////////
+// min_ass_prim_charsets1
+// input: generators of an ideal PS
+// output: the minimal associated primes of PS
+// algorithm: via characteristic sets
+// input: generators of an ideal PS and an integer i
+// The system tries to find an "optimal ordering" of
+// the variables
+//////////////////////////////////////////////////////
+
+
+proc min_ass_prim_charsets1 (ideal PS)
+{
+  def oldring=basering;
+  string n=system("neworder",PS);
+  execute "ring r="+charstr(oldring)+",("+n+"),dp;";
+  ideal PS=imap(oldring,PS);
+  matrix m=char_series(PS);  // We compute an irreducible
+                             // characteristic series
+  int i,j,k;
+  ideal I;
+  list PSI;
+  list PHI;    // the ideals given by the characteristic series
+  list ITPHI;  // their initial terms
+  for(i=nrows(m);i>=1; i--)
+  {
+     PHI[i]=ideal(m[i,1..ncols(m)]);
+     I=0;
+     for(j=size(PHI[i]);j>0;j=j-1)
+     {
+       I=I,ini_mod(PHI[i][j]);
+     }
+      I=I[2..ncols(I)];
+      ITPHI[i]=I;
+  }
+  setring oldring;
+  matrix m=imap(r,m);
+  list PHI=imap(r,PHI);
+  list ITPHI=imap(r,ITPHI);
+  // We compute the radical of each ideal in PHI
+  ideal I,JS,II;
+  int sizeJS, sizeII;
+  for(i=size(PHI);i>=1; i--)
+  {
+     I=0;
+     for(j=size(PHI[i]);j>0;j=j-1)
+     {
+       I=I+ITPHI[i][j];
+     }
+     JS=std(PHI[i]);
+     sizeJS=size(JS);
+     for(j=size(I);j>0;j=j-1)
+     {
+       II=0;
+       sizeII=0;
+       k=0;
+       while(k<=sizeII)                  // successive iteration
+       {
+         option(returnSB);
+         II=quotient(JS,I[j]);
+         option(noreturnSB);
+//std
+//         II=std(II);
+         sizeII=size(II);
+         if(sizeII==sizeJS)
+         {
+            for(k=1;k<=sizeII;k++)
+            {
+              if(leadexp(II[k])!=leadexp(JS[k])) break;
+            }
+         }
+         JS=II;
+         sizeJS=sizeII;
+       }
+    }
+    PSI=insert(PSI,JS);
+  }
+  int sizePSI=size(PSI);
+  // We eliminate redundant ideals
+  for(i=1;i<sizePSI;i++)
+  {
+    for(j=i+1;j<=sizePSI;j++)
+    {
+      if(size(PSI[i])!=0)
+      {
+        if(size(PSI[j])!=0)
+        {
+          if(size(NF(PSI[i],PSI[j],1))==0)
+          {
+            PSI[j]=ideal(0);
+          }
+          else
+          {
+            if(size(NF(PSI[j],PSI[i],1))==0)
+            {
+              PSI[i]=ideal(0);
+            }
+          }
+        }
+      }
+    }
+  }
+  for(i=sizePSI;i>=1;i--)
+  {
+    if(size(PSI[i])==0)
+    {
+      PSI=delete(PSI,i);
+    }
+  }
+  return (PSI);
+}
+
+
+/////////////////////////////////////////////////////
+// proc prim_dec
+// input:  generators of an ideal I and an integer choose
+// If choose=0, min_ass_prim_charsets with the given
+// ordering of the variables is used.
+// If choose=1, min_ass_prim_charsets with the "optimized"
+// ordering of the variables is used.
+// If choose=2, minAssPrimes from primdec.lib is used
+// If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
+// output: a primary decomposition of I, i.e., a list
+// of pairs consisting of a standard basis of a primary component
+// of I and a standard basis of the corresponding associated prime.
+// To compute the minimal associated primes of a given ideal
+// min_ass_prim_l is called, i.e., the minimal associated primes
+// are computed via characteristic sets.
+// In the homogeneous case, the performance of the procedure
+// will be improved if I is already given by a minimal set of
+// generators. Apply minbase if necessary.
+//////////////////////////////////////////////////////////
+
+
+proc prim_dec(ideal I, int choose)
+{
+   if((choose<0) or (choose>3))
+   {
+     "ERROR: <int> must be 0 or 1 or 2 or 3";
+      return();
+   }
+   if(system("version")>933)
+   {
+      option(notWarnSB);
+    }
+  ideal H=1; // The intersection of the primary components
+  list U;    // the leaves of the decomposition tree, i.e.,
+             // pairs consisting of a primary component of I
+             // and the corresponding associated prime
+  list W;    // the non-leaf vertices in the decomposition tree.
+             // every entry has 6 components:
+                // 1- the vertex itself , i.e., a standard bais of the
+                //    given ideal I (type 1), or a standard basis of a
+                //    pseudo-primary component arising from
+                //    pseudo-primary decomposition (type 2), or a
+                //    standard basis of a remaining component arising from
+                //    pseudo-primary decomposition or extraction (type 3)
+                // 2- the type of the vertex as indicated above
+                // 3- the weighted_tree_depth of the vertex
+                // 4- the tester of the vertex
+                // 5- a standard basis of the associated prime
+                //    of a vertex of type 2, or 0 otherwise
+                // 6- a list of pairs consisting of a standard
+                //    basis of a minimal associated prime ideal
+                //    of the father of the vertex and the
+                //    irreducible factors of the "minimal
+                //    divisor" of the seperator or extractor
+                //    corresponding to the prime ideal
+                //    as computed by the procedure minsat,
+                //    if the vertex is of type 3, or
+                //    the empty list otherwise
+  ideal SI=std(I);
+  int ncolsSI=ncols(SI);
+  int ncolsH=1;
+  W[1]=list(I,1,0,poly(1),ideal(0),list()); // The root of the tree
+  int weighted_tree_depth;
+  int i,j;
+  int check;
+  list V;  // current vertex
+  list VV; // new vertex
+  list QQ;
+  list WI;
+  ideal Qi,SQ,SRest,fac;
+  poly tester;
+
+  while(1)
+  {
+    i=1;
+    while(1)
+    {
+      while(i<=size(W)) // find vertex V of smallest weighted tree-depth
+      {
+        if (W[i][3]<=weighted_tree_depth) break;
+        i++;
+      }
+      if (i<=size(W)) break;
+      i=1;
+      weighted_tree_depth++;
+    }
+    V=W[i];
+    W=delete(W,i); // delete V from W
+
+    // now proceed by type of vertex V
+
+    if (V[2]==2)  // extraction needed
+    {
+      SQ,SRest,fac=extraction(V[1],V[5]);
+                        // standard basis of primary component,
+                        // standard basis of remaining component,
+                        // irreducible factors of
+                        // the "minimal divisor" of the extractor
+                        // as computed by the procedure minsat,
+      check=0;
+      for(j=1;j<=ncolsH;j++)
+      {
+        if (NF(H[j],SQ,1)!=0) // Q is not redundant
+        {
+          check=1;
+          break;
+        }
+      }
+      if(check==1)             // Q is not redundant
+      {
+        QQ=list();
+        QQ[1]=list(SQ,V[5]);  // primary component, associated prime,
+                              // i.e., standard bases thereof
+        U=U+QQ;
+        H=intersect(H,SQ);
+        H=std(H);
+        ncolsH=ncols(H);
+        check=0;
+        if(ncolsH==ncolsSI)
+        {
+          for(j=1;j<=ncolsSI;j++)
+          {
+            if(leadexp(H[j])!=leadexp(SI[j]))
+            {
+              check=1;
+              break;
+            }
+          }
+        }
+        else
+        {
+          check=1;
+        }
+        if(check==0) // H==I => U is a primary decomposition
+        {
+          return(U);
+        }
+      }
+      if (SRest[1]!=1)        // the remaining component is not
+                              // the whole ring
+      {
+        if (rad_con(V[4],SRest)==0) // the new vertex is not the
+                                    // root of a redundant subtree
+        {
+          VV[1]=SRest;     // remaining component
+          VV[2]=3;         // pseudoprimdec_special
+          VV[3]=V[3]+1;    // weighted depth
+          VV[4]=V[4];      // the tester did not change
+          VV[5]=ideal(0);
+          VV[6]=list(list(V[5],fac));
+          W=insert(W,VV,size(W));
+        }
+      }
+    }
+    else
+    {
+      if (V[2]==3) // pseudo_prim_dec_special is needed
+      {
+        QQ,SRest=pseudo_prim_dec_special_charsets(V[1],V[6],choose);
+                         // QQ = quadruples:
+                         // standard basis of pseudo-primary component,
+                         // standard basis of corresponding prime,
+                         // seperator, irreducible factors of
+                         // the "minimal divisor" of the seperator
+                         // as computed by the procedure minsat,
+                         // SRest=standard basis of remaining component
+      }
+      else     // V is the root, pseudo_prim_dec is needed
+      {
+        QQ,SRest=pseudo_prim_dec_charsets(I,SI,choose);
+                         // QQ = quadruples:
+                         // standard basis of pseudo-primary component,
+                         // standard basis of corresponding prime,
+                         // seperator, irreducible factors of
+                         // the "minimal divisor" of the seperator
+                         // as computed by the procedure minsat,
+                         // SRest=standard basis of remaining component
+
+      }
+//check
+      for(i=size(QQ);i>=1;i--)
+      //for(i=1;i<=size(QQ);i++)
+      {
+        tester=QQ[i][3]*V[4];
+        Qi=QQ[i][2];
+        if(NF(tester,Qi,1)!=0)  // the new vertex is not the
+                                // root of a redundant subtree
+        {
+          VV[1]=QQ[i][1];
+          VV[2]=2;
+          VV[3]=V[3]+1;
+          VV[4]=tester;      // the new tester as computed above
+          VV[5]=Qi;          // QQ[i][2];
+          VV[6]=list();
+          W=insert(W,VV,size(W));
+        }
+      }
+      if (SRest[1]!=1)        // the remaining component is not
+                              // the whole ring
+      {
+        if (rad_con(V[4],SRest)==0) // the vertex is not the root
+                                    // of a redundant subtree
+        {
+          VV[1]=SRest;
+          VV[2]=3;
+          VV[3]=V[3]+2;
+          VV[4]=V[4];      // the tester did not change
+          VV[5]=ideal(0);
+          WI=list();
+          for(i=1;i<=size(QQ);i++)
+          {
+            WI=insert(WI,list(QQ[i][2],QQ[i][4]));
+          }
+          VV[6]=WI;
+          W=insert(W,VV,size(W));
+        }
+      }
+    }
+  }
+}
+
+//////////////////////////////////////////////////////////////////////////
+// proc pseudo_prim_dec_charsets
+// input: Generators of an arbitrary ideal I, a standard basis SI of I,
+// and an integer choo
+// If choo=0, min_ass_prim_charsets with the given
+// ordering of the variables is used.
+// If choo=1, min_ass_prim_charsets with the "optimized"
+// ordering of the variables is used.
+// If choo=2, minAssPrimes from primdec.lib is used
+// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
+// output: a pseudo primary decomposition of I, i.e., a list
+// of pseudo primary components together with a standard basis of the
+// remaining component. Each pseudo primary component is
+// represented by a quadrupel: A standard basis of the component,
+// a standard basis of the corresponding associated prime, the
+// seperator of the component, and the irreducible factors of the
+// "minimal divisor" of the seperator as computed by the procedure minsat,
+// calls  proc pseudo_prim_dec_i
+//////////////////////////////////////////////////////////////////////////
+
+
+proc pseudo_prim_dec_charsets (ideal I, ideal SI, int choo)
+{
+  list L;          // The list of minimal associated primes,
+                   // each one given by a standard basis
+  if((choo==0) or (choo==1))
+    {
+       L=min_ass_prim_charsets(I,choo);
+    }
+  else
+    {
+      if(choo==2)
+      {
+        L=minAssPrimes(I);
+      }
+      else
+      {
+        L=minAssPrimes(I,1);
+      }
+      for(int i=size(L);i>=1;i=i-1)
+        {
+          L[i]=std(L[i]);
+        }
+    }
+  return (pseudo_prim_dec_i(SI,L));
+}
+
+////////////////////////////////////////////////////////////////
+// proc pseudo_prim_dec_special_charsets
+// input: a standard basis of an ideal I whose radical is the
+// intersection of the radicals of ideals generated by one prime ideal
+// P_i together with one polynomial f_i, the list V6 must be the list of
+// pairs (standard basis of P_i, irreducible factors of f_i),
+// and an integer choo
+// If choo=0, min_ass_prim_charsets with the given
+// ordering of the variables is used.
+// If choo=1, min_ass_prim_charsets with the "optimized"
+// ordering of the variables is used.
+// If choo=2, minAssPrimes from primdec.lib is used
+// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
+// output: a pseudo primary decomposition of I, i.e., a list
+// of pseudo primary components together with a standard basis of the
+// remaining component. Each pseudo primary component is
+// represented by a quadrupel: A standard basis of the component,
+// a standard basis of the corresponding associated prime, the
+// seperator of the component, and the irreducible factors of the
+// "minimal divisor" of the seperator as computed by the procedure minsat,
+// calls  proc pseudo_prim_dec_i
+////////////////////////////////////////////////////////////////
+
+
+proc pseudo_prim_dec_special_charsets (ideal SI,list V6, int choo)
+{
+  int i,j,l;
+  list m;
+  list L;
+  int sizeL;
+  ideal P,SP; ideal fac;
+  int dimSP;
+  for(l=size(V6);l>=1;l--)   // creates a list of associated primes
+                             // of I, possibly redundant
+  {
+    P=V6[l][1];
+    fac=V6[l][2];
+    for(i=ncols(fac);i>=1;i--)
+    {
+      SP=P+fac[i];
+      SP=std(SP);
+      if(SP[1]!=1)
+      {
+        if((choo==0) or (choo==1))
+        {
+          m=min_ass_prim_charsets(SP,choo);  // a list of SB
+        }
+        else
+        {
+          if(choo==2)
+          {
+            m=minAssPrimes(SP);
+          }
+          else
+          {
+            m=minAssPrimes(SP,1);
+          }
+          for(j=size(m);j>=1;j=j-1)
+            {
+              m[j]=std(m[j]);
+            }
+        }
+        dimSP=dim(SP);
+        for(j=size(m);j>=1; j--)
+        {
+          if(dim(m[j])==dimSP)
+          {
+            L=insert(L,m[j],size(L));
+          }
+        }
+      }
+    }
+  }
+  sizeL=size(L);
+  for(i=1;i<sizeL;i++)     // get rid of redundant primes
+  {
+    for(j=i+1;j<=sizeL;j++)
+    {
+      if(size(L[i])!=0)
+      {
+        if(size(L[j])!=0)
+        {
+          if(size(NF(L[i],L[j],1))==0)
+          {
+            L[j]=ideal(0);
+          }
+          else
+          {
+            if(size(NF(L[j],L[i],1))==0)
+            {
+              L[i]=ideal(0);
+            }
+          }
+        }
+      }
+    }
+  }
+  for(i=sizeL;i>=1;i--)
+  {
+    if(size(L[i])==0)
+    {
+      L=delete(L,i);
+    }
+  }
+  return (pseudo_prim_dec_i(SI,L));
+}
+
+
+////////////////////////////////////////////////////////////////
+// proc pseudo_prim_dec_i
+// input: A standard basis of an arbitrary ideal I, and standard bases
+// of the minimal associated primes of I
+// output: a pseudo primary decomposition of I, i.e., a list
+// of pseudo primary components together with a standard basis of the
+// remaining component. Each pseudo primary component is
+// represented by a quadrupel: A standard basis of the component Q_i,
+// a standard basis of the corresponding associated prime P_i, the
+// seperator of the component, and the irreducible factors of the
+// "minimal divisor" of the seperator as computed by the procedure minsat,
+////////////////////////////////////////////////////////////////
+
+
+proc pseudo_prim_dec_i (ideal SI, list L)
+{
+  list Q;
+  if (size(L)==1)               // one minimal associated prime only
+                                // the ideal is already pseudo primary
+  {
+    Q=SI,L[1],1;
+    list QQ;
+    QQ[1]=Q;
+    return (QQ,ideal(1));
+  }
+
+  poly f0,f,g;
+  ideal fac;
+  int i,j,k,l;
+  ideal SQi;
+  ideal I'=SI;
+  list QP;
+  int sizeL=size(L);
+  for(i=1;i<=sizeL;i++)
+  {
+    fac=0;
+    for(j=1;j<=sizeL;j++)           // compute the seperator sep_i
+                                    // of the i-th component
+    {
+      if (i!=j)                       // search g not in L[i], but L[j]
+      {
+        for(k=1;k<=ncols(L[j]);k++)
+        {
+          if(NF(L[j][k],L[i],1)!=0)
+          {
+            break;
+          }
+        }
+        fac=fac+L[j][k];
+      }
+    }
+    // delete superfluous polynomials
+    fac=simplify(fac,8);
+    // saturation
+    SQi,f0,f,fac=minsat_ppd(SI,fac);
+    I'=I',f;
+    QP=SQi,L[i],f0,fac;
+             // the quadrupel:
+             // a standard basis of Q_i,
+             // a standard basis of P_i,
+             // sep_i,
+             // irreducible factors of
+             // the "minimal divisor" of the seperator
+             //  as computed by the procedure minsat,
+    Q[i]=QP;
+  }
+  I'=std(I');
+  return (Q, I');
+                   // I' = remaining component
+}
+
+
+////////////////////////////////////////////////////////////////
+// proc extraction
+// input: A standard basis of a pseudo primary ideal I, and a standard
+// basis of the unique minimal associated prime P of I
+// output: an extraction of I, i.e., a standard basis of the primary
+// component Q of I with associated prime P, a standard basis of the
+// remaining component, and the irreducible factors of the
+// "minimal divisor" of the extractor as computed by the procedure minsat
+////////////////////////////////////////////////////////////////
+
+
+proc extraction (ideal SI, ideal SP)
+{
+  list indsets=system("indsetall",SP,0);
+  poly f;
+  if(size(indsets)!=0)      //check, whether dim P != 0
+  {
+    intvec v;               // a maximal independent set of variables
+                            // modulo P
+    string U;               // the independent variables
+    string A;               // the dependent variables
+    int j,k;
+    int a;                  //  the size of A
+    int degf;
+    ideal g;
+    list polys;
+    int sizepolys;
+    list newpoly;
+    def R=basering;
+    //intvec hv=hilb(SI,1);
+    for (k=1;k<=size(indsets);k++)
+    {
+      v=indsets[k];
+      for (j=1;j<=nvars(R);j++)
+      {
+        if (v[j]==1)
+        {
+          U=U+varstr(j)+",";
+        }
+        else
+        {
+          A=A+varstr(j)+",";
+          a++;
+        }
+      }
+
+      U[size(U)]=")";           // we compute the extractor of I (w.r.t. U)
+      execute "ring RAU="+charstr(basering)+",("+A+U+",(dp("+string(a)+"),dp);";
+      ideal I=imap(R,SI);
+      //I=std(I,hv);            // the standard basis in (R[U])[A]
+      I=std(I);            // the standard basis in (R[U])[A]
+      A[size(A)]=")";
+      execute "ring Rloc=("+charstr(basering)+","+U+",("+A+",dp;";
+      ideal I=imap(RAU,I);
+      //"std in lokalisierung:"+newline,I;
+      ideal h;
+      for(j=ncols(I);j>=1;j--)
+      {
+        h[j]=leadcoef(I[j]);  // consider I in (R(U))[A]
+      }
+      setring R;
+      g=imap(Rloc,h);
+      kill RAU,Rloc;
+      U="";
+      A="";
+      a=0;
+      f=lcm(g);
+      newpoly[1]=f;
+      polys=polys+newpoly;
+      newpoly=list();
+    }
+    f=polys[1];
+    degf=deg(f);
+    sizepolys=size(polys);
+    for (k=2;k<=sizepolys;k++)
+    {
+      if (deg(polys[k])<degf)
+      {
+        f=polys[k];
+        degf=deg(f);
+      }
+    }
+  }
+  else
+  {
+    f=1;
+  }
+  poly f0,h0; ideal SQ; ideal fac;
+  if(f!=1)
+  {
+    SQ,f0,h0,fac=minsat(SI,f);
+    return(SQ,std(SI+h0),fac);
+             // the tripel
+             // a standard basis of Q,
+             // a standard basis of remaining component,
+             // irreducible factors of
+             // the "minimal divisor" of the extractor
+             // as computed by the procedure minsat
+  }
+  else
+  {
+    return(SI,ideal(1),ideal(1));
+  }
+}
+
+/////////////////////////////////////////////////////
+// proc minsat
+// input:  a standard basis of an ideal I and a polynomial p
+// output: a standard basis IS of the saturation of I w.r. to p,
+// the maximal squarefree factor f0 of p,
+// the "minimal divisor" f of f0 such that the saturation of
+// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
+// the irreducible factors of f
+//////////////////////////////////////////////////////////
+
+
+proc minsat(ideal SI, poly p)
+{
+  ideal fac=factorize(p,1);       //the irreducible factors of p
+  fac=sort(fac)[1];
+  int i,k;
+  poly f0=1;
+  for(i=ncols(fac);i>=1;i--)
+  {
+    f0=f0*fac[i];
+  }
+  poly f=1;
+  ideal iold;
+  list quotM;
+  quotM[1]=SI;
+  quotM[2]=fac;
+  quotM[3]=f0;
+  // we deal seperately with the first quotient;
+  // factors, which do not contribute to this one,
+  // are omitted
+  iold=quotM[1];
+  quotM=minquot(quotM);
+  fac=quotM[2];
+  if(quotM[3]==1)
+    {
+      return(quotM[1],f0,f,fac);
+    }
+  while(special_ideals_equal(iold,quotM[1])==0)
+    {
+      f=f*quotM[3];
+      iold=quotM[1];
+      quotM=minquot(quotM);
+    }
+  return(quotM[1],f0,f,fac);           // the quadrupel ((I:p),f0,f, irr. factors of f)
+}
+
+/////////////////////////////////////////////////////
+// proc minsat_ppd
+// input:  a standard basis of an ideal I and a polynomial p
+// output: a standard basis IS of the saturation of I w.r. to p,
+// the maximal squarefree factor f0 of p,
+// the "minimal divisor" f of f0 such that the saturation of
+// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
+// the irreducible factors of f
+//////////////////////////////////////////////////////////
+
+
+proc minsat_ppd(ideal SI, ideal fac)
+{
+  fac=sort(fac)[1];
+  int i,k;
+  poly f0=1;
+  for(i=ncols(fac);i>=1;i--)
+  {
+    f0=f0*fac[i];
+  }
+  poly f=1;
+  ideal iold;
+  list quotM;
+  quotM[1]=SI;
+  quotM[2]=fac;
+  quotM[3]=f0;
+  // we deal seperately with the first quotient;
+  // factors, which do not contribute to this one,
+  // are omitted
+  iold=quotM[1];
+  quotM=minquot(quotM);
+  fac=quotM[2];
+  if(quotM[3]==1)
+    {
+      return(quotM[1],f0,f,fac);
+    }
+  while(special_ideals_equal(iold,quotM[1])==0)
+  {
+    f=f*quotM[3];
+    iold=quotM[1];
+    quotM=minquot(quotM);
+    k++;
+  }
+  return(quotM[1],f0,f,fac);           // the quadrupel ((I:p),f0,f, irr. factors of f)
+}
+/////////////////////////////////////////////////////////////////
+// proc minquot
+// input: a list with 3 components: a standard basis
+// of an ideal I, a set of irreducible polynomials, and
+// there product f0
+// output: a standard basis of the ideal (I:f0), the irreducible
+// factors of the "minimal divisor" f of f0 with (I:f0) = (I:f),
+// the "minimal divisor" f
+/////////////////////////////////////////////////////////////////
+
+proc minquot(list tsil)
+{
+   int i,j,k,action;
+   ideal verg;
+   list l;
+   poly g;
+   ideal laedi=tsil[1];
+   ideal fac=tsil[2];
+   poly f=tsil[3];
+
+//std
+//   ideal star=quotient(laedi,f);
+//   star=std(star);
+   option(returnSB);
+   ideal star=quotient(laedi,f);
+   option(noreturnSB);
+   if(special_ideals_equal(laedi,star)==1)
+     {
+       return(laedi,ideal(1),1);
+     }
+   action=1;
+   while(action==1)
+   {
+      if(size(fac)==1)
+      {
+         action=0;
+         break;
+      }
+      for(i=1;i<=size(fac);i++)
+      {
+        g=1;
+         for(j=1;j<=size(fac);j++)
+         {
+            if(i!=j)
+            {
+               g=g*fac[j];
+            }
+         }
+//std
+//         verg=quotient(laedi,g);
+//         verg=std(verg);
+         option(returnSB);
+         verg=quotient(laedi,g);
+         option(noreturnSB);
+         if(special_ideals_equal(verg,star)==1)
+         {
+            f=g;
+            fac[i]=0;
+            fac=simplify(fac,2);
+            break;
+         }
+         if(i==size(fac))
+         {
+            action=0;
+         }
+      }
+   }
+   l=star,fac,f;
+   return(l);
+}
+/////////////////////////////////////////////////
+// proc special_ideals_equal
+// input: standard bases of ideal k1 and k2 such that
+// k1 is contained in k2, or k2 is contained ink1
+// output: 1, if k1 equals k2, 0 otherwise
+//////////////////////////////////////////////////
+
+proc special_ideals_equal( ideal k1, ideal k2)
+{
+   int j;
+   if(size(k1)==size(k2))
+   {
+      for(j=1;j<=size(k1);j++)
+      {
+         if(leadexp(k1[j])!=leadexp(k2[j]))
+         {
+            return(0);
+         }
+      }
+      return(1);
+   }
+   return(0);
+}
+        
+        
+