Changeset 654a23 in git

Timestamp:

Jul 23, 2015, 2:48:14 PM (8 years ago)

Author:

Stephan Oberfranz <oberfran@…>

Branches:

(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')

Children:

74fe3587102ed73b3858258975abbe6d63139abf

Parents:

8d7ca03d4444afafa6c2111db95d65216b93ad464132ee2e5b539d15463dcf3db09fbf9aa7c00877

Message:

update

Merge branch 'spielwiese' of github.com:Singular/Sources into spielwiese

Files:

• Singular/LIB/classify_aeq.lib (modified) (9 diffs)
• Singular/LIB/crypto.lib (modified) (2 diffs)
• Singular/LIB/gmssing.lib (modified) (1 diff)
• Singular/LIB/gradedModules.lib (modified) (88 diffs)
• Singular/LIB/grobcov.lib (modified) (24 diffs)
• Singular/LIB/grwalk.lib (modified) (10 diffs)
• Singular/LIB/modstd.lib (modified) (1 diff)
• Singular/LIB/modwalk.lib (modified) (1 diff, 1 prop)
• Singular/LIB/ncfactor.lib (modified) (1 diff)
• Singular/LIB/nctools.lib (modified) (1 diff)
• Singular/LIB/nfmodstd.lib (modified) (1 diff)
• Singular/LIB/primdec.lib (modified) (9 diffs)
• Singular/LIB/rwalk.lib (modified) (11 diffs, 1 prop)
• Singular/LIB/surf.lib (modified) (4 diffs)
• Singular/LIB/swalk.lib (modified) (1 prop)
• Singular/Makefile.am (modified) (1 diff)
• Singular/blackbox.cc (modified) (2 diffs)
• Singular/cntrlc.cc (modified) (4 diffs)
• Singular/dyn_modules/gfanlib/bbcone.cc (modified) (2 diffs)
• Singular/dyn_modules/singmathic/singmathic.cc (modified) (1 diff)
• Singular/dyn_modules/syzextra/mod_main.cc (modified) (3 diffs)
• Singular/extra.cc (modified) (1 diff)
• Singular/fevoices.cc (modified) (1 diff)
• Singular/ipassign.cc (modified) (1 diff)
• Singular/ipshell.cc (modified) (2 diffs)
• Singular/misc_ip.cc (modified) (2 diffs)
• Singular/omSingularConfig.h (deleted)
• Singular/scanner.cc (modified) (1 diff)
• Singular/scanner.ll (modified) (1 diff)
• Singular/sdb.cc (modified) (1 diff)
• Singular/test.cc (modified) (1 diff)
• Singular/tesths.cc (modified) (1 diff)
• Singular/walk.cc (modified) (171 diffs, 1 prop)
• Singular/walk.h (modified) (1 diff, 1 prop)
• Tst/Long/modulo_l.res.gz.uu (modified) (1 diff)
• Tst/Long/modulo_l.tst (modified) (1 diff)
• Tst/Long/nfmodstd_l.res.gz.uu (modified) (1 diff)
• Tst/Long/nfmodstd_l.tst (modified) (2 diffs)
• Tst/Manual.lst (modified) (1 diff)
• Tst/Manual/all.lst (modified) (1 diff)
• Tst/Manual/m.lst (modified) (1 diff)
• Tst/Manual/setglobalrings.res.gz.uu (deleted)
• Tst/Manual/setglobalrings.stat (deleted)
• Tst/Manual/setglobalrings.tst (deleted)
• Tst/New/grmodexamples.tst (modified) (2 diffs)
• Tst/Old/m8ex.res.gz.uu (modified) (1 diff)
• Tst/Old/m8ex.tst (modified) (1 diff)
• Tst/Old/preimage.res.gz.uu (modified) (2 diffs)
• Tst/Short/bug_newstruct2.tst (modified) (1 diff)
• Tst/Short/bug_tr678_680.tst (modified) (1 diff)
• Tst/Short/bug_tr704.tst (modified) (1 diff)
• Tst/Short/bug_tr709.tst (modified) (1 diff)
• Tst/Short/elimlib_s.res.gz.uu (modified) (1 diff)
• Tst/regress.lst (modified) (1 diff)
• configure.ac (modified) (1 diff)
• doc/NEWS.texi (modified) (2 diffs)
• doc/plural.doc (modified) (previous)
• doc/reference.doc (modified) (previous)
• doc/singular.doc (modified) (previous)
• doc/start.doc (modified) (previous)
• factory/NTLconvert.cc (modified) (1 diff)
• kernel/GBEngine/kspoly.cc (modified) (1 diff)
• kernel/GBEngine/kstd2.cc (modified) (1 diff)
• kernel/GBEngine/syz0.cc (modified) (4 diffs)
• kernel/preimage.cc (modified) (1 diff)
• libpolys/coeffs/longrat.h (modified) (2 diffs)
• libpolys/coeffs/modulop.cc (modified) (1 diff)
• libpolys/polys/ext_fields/transext.cc (modified) (27 diffs)
• libpolys/polys/monomials/p_polys.cc (modified) (6 diffs)
• libpolys/polys/monomials/ring.cc (modified) (4 diffs)
• libpolys/polys/monomials/ring.h (modified) (1 diff)
• libpolys/polys/simpleideals.h (modified) (1 diff)
• libpolys/reporter/reporter.h (modified) (2 diffs)
• omalloc/omDefaultConfig.h (modified) (1 diff)

Singular/LIB/classify_aeq.lib

@@ -7 +7 @@
 AUTHORS: Faira Kanwal Janjua     fairakanwaljanjua@gmail.com
          Gerhard Pfister         pfister@mathematik.uni-kl.de
+         Khawar Mehmood          khawar1073@gmail.com                         NEU
 
 OVERVIEW: A library for classifying the simple singularities
@@ -22 +23 @@
   [4] Bruce, J.W.,Gaffney, T.J.: Simple singularities of mappings (C, 0) ->(C2,0).
   J. London Math. Soc. (2) 26 (1982), 465-474.
-
+  [5] Ishikawa,G; Janeczko,S.: The Complex Symplectic Moduli Spaces of Unimodal Parametric Plane Curve  NEU Singularities. Insitute of Mathematics of the Polish Academy of Sciences,Preprint 664(2006)
 PROCEDURES:
           sagbiAlg(G);    Compute the Sagbi-basis of the Algebra.
@@ -662 +663 @@
 
 ////////////////////////////////////////////////////////////////////////////////
-proc sagbiAlg(ideal I)
+proc sagbiAlg(ideal I,list #)
 "USAGE":  sagbiAlg(I);  I ideal
 RETURN: An ideal.The sagbi bases of I.
@@ -676 +677 @@
 
     if(size(I)==0){return(I);}
-    int b=ConductorBound(I);
-
-  // int b=200;
+    if(size(#)==0)
+    {
+       int b=ConductorBound(I);
+    // int b=200;
   //   b=correctBound(I,b);
+    }
+    else
+    {
+       int b=#[1];
+    }
    ideal S=interReduceSagbi(I,b) ;
   // b=correctBound(S,b);
@@ -699 +706 @@
    return(S);
 }
-
 example
 {
@@ -709 +715 @@
 sagbiAlg(I);
 }
-
+//////////////////////////////////////////////////////////////////////////////// NEU
+proc reducedSagbiAlg(ideal I,list #)
+{
+
+   I=sagbiAlg(I,#);
+   intvec L=semiGroup(I)[3];
+   if(size(#)==0)
+   {
+      int b=findConductor(L);
+   }
+   else
+   {
+      int b=#[1];
+   }
+   int i;
+   poly q;
+   for(i=1;i<=size(I);i++)
+   {
+      q=I[i]-lead(I[i]);
+      q=sagbiNF(q,I,b);
+      I[i]=lead(I[i])+q;
+   }
+   return(I);
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+ring R=0,t,ds;
+ideal I=t4+2t9,t9+t10+19/18t11-3t12+t13-t14;
+reducedSagbiAlg(I);
+}
+////////////////////////////////////////////////////////////////////////////////   NEU
+proc classifyAEQunimodal(ideal I)
+"USAGE":  classifyAEQunimodal(I);  I ideal generated by 2 polynomials
+RETURN: An ideal.Ideal is one of the singularity in the list of Ishikawa and Jenczko [5]
+EXAMPLE: example classifyAEQunimodal;  shows an example
+"
+{
+    def R=basering;
+    ring @S=0,t,ds;
+    ideal I=fetch(R,I);
+    ideal J;
+    poly g;
+    if(size(I)>=3){ERROR("not a plane curve");}
+    I=simplify(I,2);   //deletes zero’s i I
+    I=simplify(I,1);   //creates monic generators in I
+    if(ord(I[1])>ord(I[2])){poly q=I[2];I[2]=I[1];I[1]=q;}
+    if((ord(I[1])>=6)||(ord(I[1])<=3)){return("not in the unimodal list");}
+    //compute estimate of the term with the modulus
+    int c=ord(sagbiNF(I[2],ideal(I[1]),10));
+    if(c==10)
+    {
+       if(ord(I[1])!=4){return("not in the unimodal list");}
+       c=ord(I[2][2])+2;
+    }
+    else
+    {
+       c=0;
+       intvec v=ord(I[1]),ord(I[2]);
+       if(v==intvec(5,6)){c=14;}
+       if(v==intvec(5,7)){c=18;}
+       if(v==intvec(5,8)){c=22;}
+       if(v==intvec(4,9)){c=19;}
+       if(v==intvec(4,11)){c=25;}
+       if(c==0){return("not in the unimodal list");}
+    }
+    while(size(I[1])>1)
+    {
+      I=jet(subst(I,t,t-number(1)/number(ord(I[1]))*leadcoef(I[1][2])*t^(ord(I[1][2])-ord(I[1])+1)),c);
+    }
+    ideal G=I;
+    G[2]=sagbiNF(G[2],ideal(G[1]),c);
+    ideal M=sagbiMod(diff(G,t),G);
+    list K=semiMod(M,G);
+
+    if(K[1]==intvec(4,9))
+    {
+       if(K[4]==intvec(3,8)){J=t4,t9;}
+       if(K[4]==intvec(3,8,22)){J=t4,t9+t19;}
+       if(K[4]==intvec(3,8,18)){J=t4,t9+t15;}
+       if(K[4]==intvec(3,8,14)){J=t4,t9+t11;}
+       if(K[4]==intvec(3,8,13))
+       {
+          G=reducedSagbiAlg(G,15);
+          if(ord(G[2][4])==14)
+          {
+             //kill the term t14 by some transformation
+             G=subst(G,t,t-leadcoef(G[2][4])/9*t^6);
+             G=jet(G,15);
+             G[1]=sagbiNF(G[1],ideal(G[2]),15);
+             //arrange the first element to be t4
+             while(size(G[1])>1)
+             {
+                G=subst(G,t,t-1/4*(G[1]-lead(G[1]))/t^3);
+                G=jet(G,15);
+             }
+          }
+          G[2]=sagbiNF(G[2],ideal(G[1]),15);
+          //arrange the coefficient of t10 to become 1
+          number m=leadcoef(G[2]-lead(G[2]));
+          G[2]=m^9*subst(G[2],t,1/m*t);
+          J=G;
+       }
+       if(K[4]==intvec(3,8,13,18))
+       {
+          G=reducedSagbiAlg(G,11);
+          number m=leadcoef(G[2]-lead(G[2]));
+          G[2]=m^9*subst(G[2],t,1/m*t);
+          J=G;
+       }
+    }
+    if(K[1]==intvec(4,11))
+    {
+       if(K[4]==intvec(3,10)){J=t4,t11;}
+       if(K[4]==intvec(3,10,28)){J=t4,t11+t25;}
+       if(K[4]==intvec(3,10,24)){J=t4,t11+t21;}
+       if(K[4]==intvec(3,10,20)){J=t4,t11+t17;}
+       if(K[4]==intvec(3,10,16))
+       {
+          G=reducedSagbiAlg(G,14);
+          number m=leadcoef(G[2]-lead(G[2]));
+          number l=leadcoef(G[2][3]);
+          //lambda^2=l^2/m^3
+          J=G;
+       }
+       if(K[4]==intvec(3,10,17))
+       {
+          G=reducedSagbiAlg(G,21);
+          if(ord(G[2][4])==18)
+          {
+             //kill the term t18 by some transformation
+             G=subst(G,t,t-leadcoef(G[2][4])/11*t^8);
+             G=jet(G,21);
+             G[1]=sagbiNF(G[1],ideal(G[2]),21);
+             //arrange the first element to be t4
+             while(size(G[1])>1)
+             {
+                G=subst(G,t,t-1/4*(G[1]-lead(G[1]))/t^3);
+                G=jet(G,21);
+             }
+          }
+          G[2]=sagbiNF(G[2],ideal(G[1]),21);
+          //arrange the coefficient of t14 to become 1
+          number m=leadcoef(G[2]-lead(G[2]));
+          number l=leadcoef(G[2][4]);
+          //lambda^2=l^3/m^10
+          J=G;
+
+
+       }
+       if(K[4]==intvec(3,10,17,24))
+       {
+          G=reducedSagbiAlg(G,18);
+          //arrange the coefficient of t14 to become 1
+          number m=leadcoef(G[2]-lead(G[2]));
+          G[2]=t11+t14+leadcoef(G[2][3])/m^2*t17;
+          J=G;
+       }
+    }
+    if((size(K[1])==3)&&(K[1][1]==4)&&(K[1][2]==10))
+    {
+       int l=(K[1][3]-19) div 2;
+       G=reducedSagbiAlg(G,2*l+12);
+       number m=leadcoef(G[2]-lead(G[2]));
+       number s=leadcoef(G[2][3]);
+       //lambda^(2l-1)=s^(2l-1)/m^(2l+1)
+       J=G;
+    }
+    if(K[1]==intvec(5,6))
+    {
+       if(K[4]==intvec(4,5)){J=t5,t6;}
+       if(K[4]==intvec(4,5,18)){J=t5,t6+t14;}
+       if(K[4]==intvec(4,5,13)){J=t5,t6+t9;}
+       if(K[4]==intvec(4,5,12))
+       {
+          G=reducedSagbiAlg(G,9);
+          if(ord(G[2][2])==7)
+          {
+             //kill the term t7 by some transformation
+             G=subst(G,t,t-leadcoef(G[2][2])/6*t^2);
+             G=jet(G,10);
+             G[1]=sagbiNF(G[1],ideal(G[2]),9);
+             //arrange the first element to be t4
+             while(size(G[1])>1)
+             {
+                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
+                G=jet(G,9);
+             }
+          }
+          G[2]=sagbiNF(G[2],ideal(G[1]),9);
+          //arrange the coefficient of t8 to become 1
+          number m=leadcoef(G[2]-lead(G[2]));
+          number l=leadcoef(G[2][3]);
+          //lambda^2=l^2/m^3
+          J=G;
+
+       }
+    }
+    if(K[1]==intvec(5,7))
+    {
+       if(K[4]==intvec(4,6)){J=t5,t7;}
+       if(K[4]==intvec(4,6,22)){J=t5,t7+t18;}
+       if(K[4]==intvec(4,6,17)){J=t5,t7+t13;}
+       if(K[4]==intvec(4,6,12))
+       {
+          G=reducedSagbiAlg(G,11);
+          if(ord(G[2][3])==9)
+          {
+             //kill the term t9 by some transformation
+             G=subst(G,t,t-leadcoef(G[2][3])/7*t^3);
+             G=jet(G,11);
+             G[1]=sagbiNF(G[1],ideal(G[2]),11);
+             //arrange the first element to be t4
+             while(size(G[1])>1)
+             {
+                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
+                G=jet(G,11);
+             }
+          }
+          G[2]=sagbiNF(G[2],ideal(G[1]),11);
+          //arrange the coefficient of t8 to become 1
+          number m=leadcoef(G[2]-lead(G[2]));
+          G[2]=m^7*subst(G[2],t,1/m*t);
+          J=G;
+       }
+       if(K[4]==intvec(4,6,15))
+       {
+          G=reducedSagbiAlg(G,14);
+          if(ord(G[2][2])==9)
+          {
+             //kill the term t9 by some transformation
+             G=subst(G,t,t-leadcoef(G[2][2])/7*t^3);
+             G=jet(G,14);
+             G[1]=sagbiNF(G[1],ideal(G[2]),14);
+             //arrange the first element to be t4
+             while(size(G[1])>1)
+             {
+                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
+                G=jet(G,14);
+             }
+          }
+          G[2]=sagbiNF(G[2],ideal(G[1]),14);
+          //arrange the coefficient of t11 to become 1
+          number m=leadcoef(G[2]-lead(G[2]));
+          number l=leadcoef(G[2][3]);
+          //lambda^2=l^2/m^3
+          J=G;
+       }
+
+    }
+    if(K[1]==intvec(5,8))
+    {
+       if(K[4]==intvec(4,7)){J=t5,t8;}
+       if(K[4]==intvec(4,7,26)){J=t5,t8+t22;}
+       if(K[4]==intvec(4,7,21)){J=t5,t8+t17;}
+       if(K[4]==intvec(4,7,13))
+       {
+          G=reducedSagbiAlg(G,12);
+          if(ord(G[2][3])==11)
+          {
+             //kill the term t11 by some transformation
+             G=subst(G,t,t-leadcoef(G[2][3])/8*t^4);
+             G=jet(G,12);
+             G[1]=sagbiNF(G[1],ideal(G[2]),12);
+             //arrange the first element to be t4
+             while(size(G[1])>1)
+             {
+                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
+                G=jet(G,12);
+             }
+          }
+          G[2]=sagbiNF(G[2],ideal(G[1]),12);
+          //arrange the coefficient of t9 to become 1
+          number m=leadcoef(G[2]-lead(G[2]));
+          G[2]=m^8*subst(G[2],t,1/m*t);
+          J=G;
+       }
+
+       if(K[4]==intvec(4,7,16))
+       {
+          G=reducedSagbiAlg(G,14);
+          if(ord(G[2][2])==11)
+          {
+             //kill the term t11 by some transformation
+             G=subst(G,t,t-leadcoef(G[2][2])/8*t^4);
+             G=jet(G,14);
+             G[1]=sagbiNF(G[1],ideal(G[2]),14);
+             //arrange the first element to be t4
+             while(size(G[1])>1)
+             {
+                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
+                G=jet(G,14);
+             }
+          }
+          G[2]=sagbiNF(G[2],ideal(G[1]),14);
+          //arrange the coefficient of t12 to become 1
+          number m=leadcoef(G[2]-lead(G[2]));
+          number l=leadcoef(G[2][3]);
+          //lambda^2=l^2/m^3
+          J=G;
+
+       }
+       if(K[4]==intvec(4,7,18))
+       {
+          G=reducedSagbiAlg(G,17);
+          if(ord(G[2][2])==11)
+          {
+             //kill the term t11 by some transformation
+             G=subst(G,t,t-leadcoef(G[2][2])/8*t^4);
+             G=jet(G,17);
+             G[1]=sagbiNF(G[1],ideal(G[2]),17);
+             //arrange the first element to be t4
+             while(size(G[1])>1)
+             {
+                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
+                G=jet(G,17);
+             }
+          }
+          G[2]=sagbiNF(G[2],ideal(G[1]),17);
+          //arrange the coefficient of t12 to become 1
+          number m=leadcoef(G[2]-lead(G[2]));
+          number l=leadcoef(G[2][3]);
+          //lambda^2=l^2/m^3
+          J=G;
+       }
+    }
+    setring R;
+    ideal J=fetch(@S,J);
+    if(size(J)==0)
+    {
+       return("not in the unimodal list");
+    }
+    return(J);
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+ring R=0,t,ds;
+ideal I=t4,9t9+18t10+38t11-216t12+144t13-288t14;
+classifyAEQunimodal(I);
+I=t4,9t9+18t10+40t11-216t12+144t13-288t14;
+classifyAEQunimodal(I);
+I=t4,t11+t12+3t14+2t15+7t16+7t17;
+classifyAEQunimodal(I);
+I=t4,t11+t14+25/22t17+3t18+4t21;
+classifyAEQunimodal(I);
+I=t5,t6+2t7+t8+3t9;
+classifyAEQunimodal(I);
+I=t5,t7+3t8+3t9+5t10;
+classifyAEQunimodal(I);
+I=t5,t7+3t11+3t12+5t13;
+classifyAEQunimodal(I);
+I=t5,t8+3t9+5t10+2t11+3t12+5t13;
+classifyAEQunimodal(I);
+I=t5,t8+5t11+3t12+7t13+5t14;
+classifyAEQunimodal(I);
+I=t5,t8+5t11+7t13+5t14+7t15+2t16+8t17;
+classifyAEQunimodal(I);
+I=subst(I,t,t+t2);
+classifyAEQunimodal(I);
+I=t4+2t5+3t6+5t7+t8,9t9+18t10+40t11-216t12+144t13-288t14;
+classifyAEQunimodal(I);
+}
+////////////////////////////////////////////////////////////////////////////////   NEU
+proc computeModulus(poly p)
+"USAGE":  computeModulus(p);  p monic poly with 3 or 4 monomials
+RETURN: A polynomial with first and second coefficient 1
+EXAMPLE: computeModulus;  shows an example
+ASSUME: the basering has one vearable and one parameter
+"
+{
+   def R=basering;
+   int a1=ord(p);
+   int a2=ord(p-lead(p));
+   number m=leadcoef(p-lead(p));
+
+   poly q=par(1)^(a2-a1)-1/m;
+   ring S=(0,a),t,ds;
+   number m=fetch(R,m);
+   minpoly=par(1)^(a2-a1)-1/m;
+   poly p=fetch(R,p);
+   p=1/par(1)^a1*subst(p,var(1),par(1)*var(1));
+   setring R;
+   p=imap(S,p);
+   return(list(p,q));
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+ring R=(0,a),t,ds;
+poly p=t8-395/16t14+4931/32t17;
+computeModulus(p);
+p=t8+3t12-395/16t14;
+computeModulus(p);
+p=t8-395/16t14+4931/32t17;
+computeModulus(p);
+
+}
+
+////////////////////////////////////////////////////////////////////////////////   NEU
+static proc n_thRoot(poly p,int n, int b)
+{
+   //computes the n-th root of 1+p up to order b
+   //assumes that p(0)=0
+   poly s=1;
+   poly q=jet(p,b);
+   if(q==0){return(s);}
+   int i;
+   for(i=1;i<=b;i++)
+   {
+      s=s+bino(n,i)*q;
+      q=jet(q*p,b);
+   }
+   return(jet(s,b));
+}
+////////////////////////////////////////////////////////////////////////////////   NEU
+static proc bino(number n, int i)
+{
+//computes the i-th binomial coefficient of 1/n
+   if(i==0){return(1);}
+   return(bino(n,i-1)*(1/n-i+1)/i);
+}
 ////////////////////////////////////////////////////////////////////////////////
 proc sagbiMod(ideal I,ideal G)
@@ -1283 +1710 @@
     execute
     ("ring T=("+charstr(br)+",x(1),z(1..n)),(x(2),y(1..m)),dp;");
-    execute
-    ("ring R=("+charstr(br)+"),(x(1..2),y(1..m),z(1..n)),(lp(m+2),dp(n));");
+     execute
+    ("ring R=("+charstr(br)+"),(x(1..2),y(1..m),z(1..n)),(lp(2),dp(m),dp(n));");
+
      map phi = br,x(1);
      ideal P = phi(G1);
@@ -1305 +1733 @@
      M[size(M)+1]=check-x(2);
      check=check*keep;
+
      option(redSB);
      M=std(M);
      int j,s;
+
      for(i=1;i<=size(M);i++)
      {
@@ -1342 +1772 @@
          f=sagbiNFMODO(f,G,I,b);
     }
-    return(f+sagbiNFMOD(p-lead(p),G,I,b));
+    return(lead(f)+sagbiNFMOD(p-lead(p),G,I,b));
 }
 ////////////////////////////////////////////////////////////////////////////////

Singular/LIB/crypto.lib

@@ -69 +69 @@
   merkle_hellman_transformation(list knapsack, int key, int mod1                generates a hard knapsack for the  Merkle-Hellman Kryptosystem for a given easy knapsack , a multiplicator key and a modulus mod1
   merkle_hellman_encryption(list knapsack, list message)                    encrypts a message with the Merkle-Hellman Kryptosystem, using a hard knapsack and a message encoded as binary list
-  merkle_hellman_decryption(list knapsack, int key, int mod1, int message)          decrypts a message with the multiplicative Merkle-Hellman Kryptosystem, using the hard knapsack, the key, the modulus mod1 and the message encoded as integer   
+  merkle_hellman_decryption(list knapsack, int key, int mod1, int message)          decrypts a message with the multiplicative Merkle-Hellman Kryptosystem, using the hard knapsack, the key, the modulus mod1 and the message encoded as integer
   super_increasing_knapsack(int ksize)                            Creates the smallest super-increasing knapsack of given size ksize
   h_increasing_knapsack(int ksize, int h)                            Creates the smallest h-increasing knapsack of given size ksize and h
@@ -1272 +1272 @@
    if((k mod 2)==0)
    {
-      resu=ellipticMult(N,a,b,P,k/2);
+      resu=ellipticMult(N,a,b,P,k div 2);
       return(ellipticAdd(N,a,b,resu,resu));
    }

Singular/LIB/gmssing.lib

@@ -110 +110 @@
 
   M=transpose(simplify(M,2));
-  I=M[1];
+  I=ideal(M[1]);
   attrib(I,"isSB",1);
   M=M[2..ncols(M)];

Singular/LIB/gradedModules.lib

@@ -31 +31 @@
     grsum(M,N)      direct sum of two graded modules coker(M) + coker(N)
     grpower(M,p)    direct p-th power of graded module coker(M)
-    grtranspose(M)  un-ordered graded transpose of map M        
+    grtranspose(M)  un-ordered graded transpose of map M
     grgens(M)       try to compute submodule generators of coker(M)
     grpres(F)       presentation of submodule generated by columns of F
@@ -51 +51 @@
     mappingcone(M,N)   mapping cone?
     grlifting3(A,B)    RND! chain lifting? probably wrong one
-    mappingcone3(A,B)  mapping cone3? 
+    mappingcone3(A,B)  mapping cone3?
     grrange(M)         get the row-weightings
     grneg(A)           graded object given by -A
     matrixpres(a)      matrix presentation of direct sum of Omega^{a[i]}(i)
 ";
-    
+
 //    grisequal(A,B)  check whether A is exactly eqal to B? TODO: isomorphic!
 
@@ -86 +86 @@
 
   v = -v; // NOTE: due to Mathematical meanings of Singular data
-  
+
 
   int lst = v[1];
-  int cnt = 1; 
+  int cnt = 1;
 
   string p = R;
@@ -95 +95 @@
 
   int k, d;
-  for (k = 2; k <= n; k++ ) 
+  for (k = 2; k <= n; k++ )
   {
     d = v[k];
     if( d == lst ) { cnt = cnt + 1; }
-    else 
+    else
     {
       if (cnt > 1){ p = p + "^" + string(cnt); }
@@ -120 +120 @@
   def E = grtwist(2, 0);
   def v = grrange(E); // grdeg(E);
-  grsumstr("", v );  
+  grsumstr("", v );
 }
 
@@ -180 +180 @@
   {
     int n = size(N); ASSUME(0, n > 0);
-    
+
     string msg1 = "";
     if( size(R) >= 2 )
@@ -186 +186 @@
       msg1 = msg1 + "(let R:="+R+")";
       R = "R"; // !!!
-    }    
+    }
     msg1 = msg1 + ": " ;
 
-    
+
 
     list arr; arr[n] = list();
     int exact = (1==1);
-    
+
     int i = 1;
-    
+
     ASSUME(1, grtest(N[i]));
-    
+
     string dst = grsumstr(R, grrange(N[i]));
     string src = grsumstr(R, grdeg(N[i]));
-    
+
     arr[i] = list(dst,  src);
 
     i = i + 1;
-    
+
     while( i <= n )
     {
       ASSUME(1, grtest(N[i]));
-      
+
       dst = grsumstr(R, grrange(N[i]));
-      
+
       if( exact && (src != dst) )
       {
@@ -218 +218 @@
 
       src = grsumstr(R, grdeg(N[i]));
-      
+
       arr[i] = list(dst,  src);
 
@@ -236 +236 @@
         msg = msg + newline + arr[i][1] + " <-- "+o+"_" + string(i) + " --";
       };
-      
+
       msg =  msg + newline + arr[n][2];
-      msg = msg + ", given by maps: "; 
+      msg = msg + ", given by maps: ";
     } else
     {
 //      print(arr);
 
-      msg = msg + "-object collection";      
+      msg = msg + "-object collection";
       o = "o";
-      
+
 //      for( i = 1; i <= n; i++ )
 //      {
 //        msg = msg + newline + arr[i][1] + " <-- "+o+"_" + string(i) + " -- " + arr[i][2];
-//      };      
-      msg = msg + ", given by the following maps (named here as "+o+"_[1 .. "+string(n)+"]): "; 
+//      };
+      msg = msg + ", given by the following maps (named here as "+o+"_[1 .. "+string(n)+"]): ";
     }
 
@@ -257 +257 @@
 
     for( i = 1; i <= size(N); i++ )
-    { 
+    {
       print( o+"_" + string(i) + " :" );
-      grview( N[i] ); 
+      grview( N[i] );
     };
 
@@ -268 +268 @@
 
   ASSUME(1, grtest(N) );
-  
+
   msg = msg + " homomorphism";
   if( size(R) >= 2 )
@@ -280 +280 @@
   intvec gr = grrange(N); // grading weights?
   string dst = grsumstr(R, gr);
-  
+
   intvec G = grdeg(N);
   string src = grsumstr(R, G);
-  
+
   if( ncols(N) == 0 )
   {
     src = "0";
-  } 
+  }
 
   lst = msg;
 
-  if( (size(lst) + size(dst) + size(src) + 4) > 80 ) 
-  { 
+  if( (size(lst) + size(dst) + size(src) + 4) > 80 )
+  {
     if( (size(lst) + size(dst)) > 80 ) { msg = msg + newline; lst = ""; }
 
@@ -312 +312 @@
 //  lst = lst + ", given by ";
 
-  
+
   int nc = ncols(N); int nr = nrows(N);
 
@@ -323 +323 @@
   {
     ASSUME(0, nc > 0);
-    
+
     matrix M = module(N);
-    
+
     int r,c;
     int d = 1; // number of extra cols/rows for extra info around the central degree(N) block in D
@@ -361 +361 @@
       }
 
-    } else 
-    { 
+    } else
+    {
       msg = msg + "a matrix";
     }
 
-    print(msg + ", with degrees: " );    
-    draw(D, d); // print it nicely!    
+    print(msg + ", with degrees: " );
+    draw(D, d); // print it nicely!
   }
 }
@@ -462 +462 @@
 "
 {
-  ASSUME(1, grtest(M) ); 
-  
+  ASSUME(1, grtest(M) );
+
   if ( typeof(attrib(M, "degHomog")) == "intvec" )
   {
@@ -469 +469 @@
 
     if( size(t) == 0 ){ return (t); } // ZERO!
-    
+
     ASSUME(2, ncols(M) == size(t) );
     return (t);
@@ -629 +629 @@
   "Input degrees: "; grview(I);
 
-  def RR = grres(I, 0, 1); list L = RR; 
+  def RR = grres(I, 0, 1); list L = RR;
 
   " = Non-minimal betti numbers: "; print(betti(L, 0), "betti");
@@ -830 +830 @@
 
   ring r=32003,(x,y,z),dp;
-  
+
   grview( grtwists ( intvec(-4, 1, 6 )) );
 
@@ -843 +843 @@
 "
 {
-  ASSUME(0, a > 0);  
+  ASSUME(0, a > 0);
   def Z = grtwists( intvec(d:a) ); // will set the rank as well
 //  ASSUME(2, grisequal(Z, grpower( grshift(grzero(), d), a ) )); // optional check
@@ -931 +931 @@
 
   def SS = grobj(S, c, dc);
-  
+
   ASSUME(0, size( grrange(SS) ) == (size(a) + size(b)) );
   ASSUME(0, size( grdeg(SS) ) == (size(da) + size(db)) );
   ASSUME(0, ncols( SS ) == size( grdeg(SS) ) );
   ASSUME(0, nrows( SS ) == size( grrange(SS) ) );
-  
+
   return(SS);
 }
@@ -984 +984 @@
 {
   ASSUME(1, grtest(M) );
-           
-           
+
+
 
   intvec a = grrange(M);
@@ -993 +993 @@
   {
     "!! Warning: shifting '0 <- 0' leaves it as it unchanged!";
-    return (M);    
-  }
-  
+    return (M);
+  }
+
   a = a - intvec(d:size(a));
   t = t - intvec(d:size(t));
@@ -1056 +1056 @@
 {
   ASSUME(0, size(w) >= nrows(A) );
-  
+
   module M = module(A);
 
@@ -1063 +1063 @@
 
   intvec @ww = 0:0;
-  
+
   if( size(#) > 0 )
   {
     ASSUME(0, typeof(#[1]) == "intvec" );
-    
+
     @ww = intvec( #[1] );
 
@@ -1077 +1077 @@
       }
     }
-    
+
     ASSUME(0, size(@ww) == ncols(M) );
   }
@@ -1091 +1091 @@
         M = freemodule(0);
       }
-      
+
       attrib( M, "rank", size(w) );
       attrib( M, "isHomog", w );
-      
+
 //      ASSUME(0, /* PROBLEM WITH ZERO COLUMNS / THEIR DEGREES! */ (ncols(M) == 0) );
     }
   }
-  
+
 //  type(@ww);  type(M);
-  
+
   ASSUME(0, size(@ww) == ncols(M) ); // ?!
-  
+
   attrib(M, "degHomog", @ww);
 
@@ -1132 +1132 @@
   grview( grobj( freemodule(0), intvec(1,2,3) ) );
 
-  matrix z1[0][3]; grview( grobj( z1, 0:0, intvec(1,2,3) ) ); 
+  matrix z1[0][3]; grview( grobj( z1, 0:0, intvec(1,2,3) ) );
   grview( grobj( freemodule(0), 0:0, intvec(1,2,3) ) );
 
   matrix z0[0][0]; grview( grobj( z0, 0:0 ) );
   grview( grobj( freemodule(0), 0:0 ) );
-  
-  
+
+
 
 }
@@ -1145 +1145 @@
 "USAGE:  grtest(M[,b]), anyting M, optionally int b
 RETURN:  1 if M is a valid graded object, 0 otherwise
-PURPOSE: validate a graded object. Print an invalid object message if b is not given 
+PURPOSE: validate a graded object. Print an invalid object message if b is not given
 NOTE: M should be an ideal or module or matrix, with weighting attribute
    'isHomog' and optionally total graded degrees attribute 'degHomog'.
@@ -1169 +1169 @@
   if ( nrows(N) != size(gr) )
   {
-    if(b) { "   ? grtest: Input has wrong number of rows!"; };    
+    if(b) { "   ? grtest: Input has wrong number of rows!"; };
     return (0);
   };
@@ -1177 +1177 @@
     return(1);
   }
-  
+
 //  if( attrib(N, "rank") != size(gr) ){ return (0); } // wrong rank :(
 
@@ -1189 +1189 @@
       return (1);
     }
-    
+
     if ( ncols(N) != size(T) )
     {
-      if(b) { "   ? grtest: Input has wrong number of cols!"; };    
+      if(b) { "   ? grtest: Input has wrong number of cols!"; };
       return (0);
     };
-    
+
     int k = nvars(basering) + 1; // index of mod. column in the leadexp
 
@@ -1259 +1259 @@
 {
   ASSUME(0, d >= 0 );
-  
+
   if( d == 0 ) { return (A); }
-  
+
   if( ncols(A) == 0 )
   {
     matrix B[nrows(A) + d][0];
-    return (B);    
-  }
-  
+    return (B);
+  }
+
   module T; T[d] = 0;
   T = T, module(transpose(A));
@@ -1281 +1281 @@
   matrix m[0][2];
   type( align(m, 3) );
-}  
+}
 
 
@@ -1302 +1302 @@
   module A = grobj( module([x+y, x, 0, 0], [0, x+y, y, 0]), intvec(0,0,0,1) );
   grview(A);
-  
+
   module B = grgroebner(A);
   grview(B);
@@ -1324 +1324 @@
   ring r=32003,(x,y,z),dp;
 
-  module A = grgroebner( grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ) ); 
+  module A = grgroebner( grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ) );
   grview(A);
-  
+
   grview(grsyz(A));
-  
-  module X = grgroebner( grobj( module([x]), intvec(2) ) ); 
+
+  module X = grgroebner( grobj( module([x]), intvec(2) ) );
   grview(X);
 
   // syzygy module should be zero!
   grview(grsyz(X));
-  
-  
+
+
 }
 
@@ -1351 +1351 @@
   intvec a = grdeg(A);
   intvec b = grrange(B);
-  
+
   ASSUME(0, (size(a) == size(b)) && (a == b));  // == for intvec :(
 
@@ -1361 +1361 @@
   ring r=32003,(x,y,z),dp;
 
-  module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); 
+  module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) );
   grview(A);
-  
+
   A = grgroebner(A);
   grview(A);
-  
+
   module B = grsyz(A);
   grview(B);
   print(B);
-  
+
   module D = grprod( A, B );
   grview(D);
@@ -1384 +1384 @@
 "USAGE:  grres(M, l[, b]), graded object M, int l, int b
 RETURN:  graded resolution = list of graded objects
-PURPOSE: compute graded resolution of M (of length l) and minimise it if b was given                             
+PURPOSE: compute graded resolution of M (of length l) and minimise it if b was given
 EXAMPLE: example grres; shows an example
 "
@@ -1392 +1392 @@
 
   intvec v = grrange(A);
-  
+
   int b = (size(#) > 0);
   if(b) { list r = res(A, l, #[1]); } else { list r = res(A, l); }
 
   l = size(r);
-  
+
   int i;
-  
+
   for ( i = 1; i <= l; i++ )
   {
@@ -1411 +1411 @@
     r[i] = grobj(r[i], v); v = grdeg(r[i]);
   }
-  
+
   i = i-1;
 
@@ -1421 +1421 @@
   ring r=32003,(x,y,z),dp;
 
-  module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); 
+  module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) );
   grview(A);
-  
+
   module B = grgroebner(A);
   grview(B);
@@ -1438 +1438 @@
 RETURN:  graded object
 PURPOSE: graded transpose of M
-NOTE:    no reordering is performend by this procedure   
+NOTE:    no reordering is performend by this procedure
 EXAMPLE: example grtranspose; shows an example
 "
@@ -1457 +1457 @@
 
   module K = grsyz( L ); grview(K);
-  
+
 
   // Corner cases: 0 <- 0!
@@ -1463 +1463 @@
   grview( grtranspose( Z ) );
 
-  
+
   // Corner cases: * <- 0
   matrix M1[3][0];
-  
+
   module Z1 = grobj( M1, intvec(-1, 0, 1) ); grview(Z1);
   grview( grtranspose( Z1 ) );
-  
-  
+
+
   // Corner cases: 0 <- *
   matrix M2[0][3];
@@ -1476 +1476 @@
   module Z2 = grobj( M2, 0:0, intvec(-1, 0, 1) ); grview(Z2);
   grview( grtranspose( Z2 ) );
-  
+
 }
 
@@ -1482 +1482 @@
 proc grgens(def M)
 "USAGE:   grgens(M), graded object M (map)
-RETURN:  graded object  
+RETURN:  graded object
 PURPOSE: try compute graded generators of coker(M) and return them as columns
          of a graded map.
@@ -1492 +1492 @@
 
   module N = grtranspose( grsyz( grtranspose(M) ) );
-  
+
 //  ASSUME(3, grisequal( grgroebner(M), grgroebner( grpres( N ) ) ) ); // FIXME: not always true!?
-  
+
   return ( N );
 }
@@ -1505 +1505 @@
 
   module N = grgens(M);
-  
+
   grview( N ); print(N); // fine == M
 
 
-  module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); 
+  module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) );
 
   A = grgroebner(A); grview(A);
@@ -1518 +1518 @@
 
   grview( grgens( grzero() ) );
-  
+
 }
 
@@ -1534 +1534 @@
 
 //  ASSUME(3, grisequal( M, grgens( N ) ) );
-  
+
   return ( N );
 }
@@ -1545 +1545 @@
   grview(A);
 
-  "graded transpose: "; def B = grtranspose(A); grview( B ); print(B); 
+  "graded transpose: "; def B = grtranspose(A); grview( B ); print(B);
 
   "... syzygy: "; def C = grsyz(B); grview(C);
@@ -1553 +1553 @@
   "... and back to presentation: "; def E = grsyz( D ); grview(E); print(E);
 
-  def F = grgens( E ); grview(F); print(F); 
+  def F = grgens( E ); grview(F); print(F);
   def G = grpres( F ); grview(G); print(G);
 
 
   def M = grtwists( intvec(-2, 0, 4, 4) ); grview(M);
-  
+
   def N = grgens(M); grview( N ); print(N);
-  
-  def L = grpres( N ); grview( L ); print(L);  
+
+  def L = grpres( N ); grview( L ); print(L);
 }
 
@@ -1652 +1652 @@
   module N;
 
-  if(i>n) 
+  if(i>n)
     { // no middle part
       if(a[1]>0)
@@ -1659 +1659 @@
 
           if(a[n+1]>0)
-            { N=grsum(N,grtwist(a[n+1],-1));}        
-        }  
+            { N=grsum(N,grtwist(a[n+1],-1));}
+        }
       else
         { N=grtwist(a[n+1],-1);}
-      
-      return (N); // grorder(N)); 
+
+      return (N); // grorder(N));
     }
-  else // i <= n: middle part is present, a_i != 0 
+  else // i <= n: middle part is present, a_i != 0
     { // a = a1  ... |  i:2, a_2 ..... i: n, a_n | .... i: n+1a_(n+1)
       j = i - 1;
@@ -1697 +1697 @@
 
 
-      return ((N)); //      return (grorder(N)); 
+      return ((N)); //      return (grorder(N));
     }
 }
@@ -1704 +1704 @@
   ring r = 32003,(x(0..4)),dp;
 
-  def N1 = KeneshlouMatrixPresentation(intvec(2,0,0,0,0)); 
+  def N1 = KeneshlouMatrixPresentation(intvec(2,0,0,0,0));
   grview(N1);
 
@@ -1738 +1738 @@
   ASSUME(1, grtest(B));
   ASSUME(0, grrange(A)==grrange(B));
- 
+
   intvec v = grrange(A);
   intvec w=grdeg(A),grdeg(B);
@@ -1746 +1746 @@
 { "EXAMPLE:"; echo = 2;
   ring r;
-  
+
   module R=grobj(module([x,y,z]),intvec(0:3));
   grview(R);
@@ -1771 +1771 @@
 //  matrix U;
   matrix L =lift(A,B/*,U*/);  //  module(B*U) = module(matrix(A)*L)
-  
+
   return(grobj(L, grdeg(A), grdeg(B)));
 }
@@ -1782 +1782 @@
 
 
-  module D=grobj(module([y,0,z],[x2+y2,z,0]),intvec(0,1,0)); 
+  module D=grobj(module([y,0,z],[x2+y2,z,0]),intvec(0,1,0));
   grview(D);
 
@@ -1788 +1788 @@
   grview(G);
 
-  ASSUME(0, grisequal( grprod(D, G), P) );  
+  ASSUME(0, grisequal( grprod(D, G), P) );
 }
 
@@ -1807 +1807 @@
 
   module Z = grobj(freemodule(0),intvec(0:0),intvec(0:0));
-  
+
   grrange(Z);
   grdeg(Z);
-  
+
   grview(Z);
 
@@ -1838 +1838 @@
        grtranspose( grprod( N,  alpha1 ) ) )
        ) ); // alpha0!
-    
+
 }
 example
@@ -1899 +1899 @@
 
   ASSUME(0, t >= 2);
-  
+
   list P;
 
   "t: ", t;
-  
+
   P[1]= grrndmap( rM[1], rN[1] ); // alpha1
 
-  if(t==2){return(P[1]);}  
-    
+  if(t==2){return(P[1]);}
+
   for(k=2; k<=t; k++)
   {
-    P[k] = grlift( grprod(P[k-1],rM[k]), rN[k] ); 
-     grview(P[k]);  
-    
-  }
-      
+    P[k] = grlift( grprod(P[k-1],rM[k]), rN[k] );
+     grview(P[k]);
+
+  }
+
   return(P);
-    
+
 }
 example
@@ -1923 +1923 @@
   ring r=32003,(x,y,z),dp;
 
-  module P=grobj(module([xy,0,xz]),intvec(0,1,0)); 
+  module P=grobj(module([xy,0,xz]),intvec(0,1,0));
   grview(P);
 
-  module D=grobj(module([y,0,z],[x2+y2,z,0]),intvec(0,1,0)); 
+  module D=grobj(module([y,0,z],[x2+y2,z,0]),intvec(0,1,0));
   grview(D);
 
@@ -1937 +1937 @@
   module D=grobj(module([y,0,z],[x2+y2,z,0], [z3, xy, xy2]),intvec(0,1,0));
   D = grgroebner(D);
-  grview( grres(D, 0)); 
+  grview( grres(D, 0));
 
   def G=grlifting(D, D);
@@ -1951 +1951 @@
   def M = KeneshlouMatrixPresentation(intvec(0,0,1,0));
   grview( grres(M, 0) );
-  
+
 //   module N = grshift(kos[3], 1); // psi, Syz_2(K(1))
   def N = KeneshlouMatrixPresentation(intvec(0,1,0,0));
@@ -1966 +1966 @@
 proc mappingcone(M,N)
 "USAGE: mappingcone(M,N), M,N graded objects
-RETURN: chain complex (as a list) 
+RETURN: chain complex (as a list)
 PURPOSE: construct a free resolution of the cokernel of a random map between Img(M), and Img(N).
 EXAMPLE: example mappingcone; shows an example
@@ -2026 +2026 @@
 
 // correct
-proc grrndmap(def S, def D, list #) 
+proc grrndmap(def S, def D, list #)
 "USAGE: grrndmap(S,D), graded objects S and D
 RETURN: graded object
@@ -2092 +2092 @@
 // 0<---M<----F0<----F1<----F2<----F3<----
 //              |p1   |p2
-// 
+//
 // 0<---N<----G0<----G1<----G2<----G3<----
 //                B(g1)      g2     g3
-// 
+//
 proc grlifting2(A,B)
 "USAGE: grlifting2(A,B), graded objects A and B (matrices defining maps)
-RETURN: map of chain complexes (as a list) 
+RETURN: map of chain complexes (as a list)
 PURPOSE: construct a map of chain complexes between free resolution of
 M=coker(A) and N=coker(B).
@@ -2127 +2127 @@
   P[1]=grrndmap2(A,B);
 
-  // A(or B)=0 
+  // A(or B)=0
   if(t==1){return(P[1])};
- 
+
   for(k=2;k<=t;k++)
   {
@@ -2161 +2161 @@
 def T=grlifting2(DD,PP); T;
 
-// def Z=grlifting2(P,D); Z; // WRONG!!!  
-             
+// def Z=grlifting2(P,D); Z; // WRONG!!!
+
 }
 
@@ -2171 +2171 @@
 proc mappingcone2(A,B)
 "USAGE: mappingcone2(A,B), graded objects A and B (matrices defining maps)
-RETURN: chain complex (as a list) 
+RETURN: chain complex (as a list)
 PURPOSE: construct the free resolution of a cokernel of a random map between M=coker(A), and N=coker(B)
 EXAMPLE: example mappingcone2;
@@ -2234 +2234 @@
 */
 
-proc grlifting3(A,B) 
+proc grlifting3(A,B)
 "TODO: grlifting4 was newer and had more documentation than this proc, but was removed... Please verify and update!
 "
@@ -2247 +2247 @@
   list rN = grres(B,0,1);
   print( betti(rN), "betti");
-  
+
   int i,j,k;
 
@@ -2260 +2260 @@
   }
   int t=min(i,j);
-  
+
   list P;
 
   "t: ", t;
 //  grview(rM[t]);  grview(rN[t]);
-  
+
   P[t]= grrndmap2(rM[t],rN[t]);
   grview(P[t]);
@@ -2312 +2312 @@
 //def I=KeneshlouMatrixPresentation(intvec(2,3,0,6,2));
 //def J=KeneshlouMatrixPresentation(intvec(4,0,1,2,1));
-//def N=grlifting3(I,J); grview(N); 
+//def N=grlifting3(I,J); grview(N);
 }
 
@@ -2318 +2318 @@
 "USAGE: grneg(A), graded object A
 RETURN: graded object
-PURPOSE: graded map defined by -A 
+PURPOSE: graded map defined by -A
 EXAMPLE: example grneg; shows an example
 "
@@ -2341 +2341 @@
 proc mappingcone3(A,B)
 "USAGE: mappingcone3(A,B), graded objects A and B (matrices defining maps)
-RETURN: chain complex (as a list) 
+RETURN: chain complex (as a list)
 PURPOSE: construct a free resolution of the cokernel of a random map between M=coker(A), and N=coker(B)
 EXAMPLE: example mappingcone3; shows an example
@@ -2376 +2376 @@
 
     T[i]=grtranspose(D);
-    
+
     kill A, B, C, D, v, w, r, s, zero;
   }
@@ -2394 +2394 @@
 def F=grlifting3(A,T); grview(F);
 
-// BUG in the proc  
+// BUG in the proc
 def G=mappingcone3(A,T); grview(G);
 
@@ -2403 +2403 @@
 ideal P=groebner(flatten(U[2]));
 resolution L=mres(P,0);
-print(betti(L),"betti");    
+print(betti(L),"betti");
 */
 
@@ -2421 +2421 @@
 def I=KeneshlouMatrixPresentation(intvec(2,3,0,6,2));
 def J=KeneshlouMatrixPresentation(intvec(4,0,1,2,1));
-// def N=grlifting3(I,J); 
+// def N=grlifting3(I,J);
 // 2nd module does not lie in the first:
 // def NN=mappingcone3(I,J); // ????????
@@ -2458 +2458 @@
   module N;
 
-  if(i>n) 
+  if(i>n)
     { // no middle part
       if(a[1]>0)
@@ -2465 +2465 @@
 
           if(a[n+1]>0)
-            { N=grsum(N,grtwist(a[n+1],0));}         
-        }  
+            { N=grsum(N,grtwist(a[n+1],0));}
+        }
       else
         { N=grtwist(a[n+1],0);}
-      
-      return (N); // grorder(N)); 
+
+      return (N); // grorder(N));
     }
 
-else // i <= n: middle part is present, a_i != 0 
+else // i <= n: middle part is present, a_i != 0
     { // a = a1  ... |  i:2, a_2 ..... i: n, a_n | .... i: n+1a_(n+1)
       module I = maxideal(1);
@@ -2505 +2505 @@
 
 
-      return ((N)); //      return (grorder(N)); 
+      return ((N)); //      return (grorder(N));
     }
 }
@@ -2518 +2518 @@
 grview(S);
 
-def N1 = matrixpres(intvec(2,0,0,0,0)); 
+def N1 = matrixpres(intvec(2,0,0,0,0));
 grview(N1);
 

Singular/LIB/grobcov.lib

@@ -1 +1 @@
 //
-version="version grobcov.lib 4.0.1.2 Jan_2015 "; // $Id$
+version="version grobcov.lib 4.0.2.0 Jul_2015 "; // $Id$
+           // version L;  July_2015;
 category="General purpose";
 info="
@@ -54 +55 @@
             Groebner Cover, and new theoretical developments have been done.
 
+            The actual version also includes a routine (ConsLevels)
+            for computing the canonical form of a constructible set, given as a
+            union of locally closed sets. It is used in the new version for the
+            computation of loci and envelops. It determines the canonical locally closed
+            level sets of a constructuble. They will be described in a forthcoming paper:
+
+             J.M. Brunat, A. Montes,
+            \"Computing the canonical representation of constructible sets\".
+            Submited to Mathematics in Computer Science. July 2015.
+
             A new set of routines (locus, locusdg, locusto) has been included to
             compute loci of points. The routines are used in the Dynamic
@@ -69 +80 @@
              \''Envelops in Dynamic Geometry using the Groebner cover\''.
 
-            The actual version also includes two routines (AddCons and AddconsP)
-            for computing the canonical form of a constructible set, given as a
-            union of locally closed sets. They are used in the new version for the
-            computation of loci and envelops. It determines the canonical locally closed
-            level sets of a constructuble. They will be described in a forthcoming paper:
-
-            A. Montes, J.M. Brunat,
-            \"Canonical representations of constructible sets\".
-
-            This version was finished on 31/01/2015
+
+            This version was finished on 31/07/2015
 
 NOTATIONS: All given and determined polynomials and ideals are in the
@@ -85 +88 @@
 @*         grobcov, cgsdr,
 @*         generate the global rings
-@*         Grobcov::@R   (Q[a][x]),
-@*         Grobcov::@P   (Q[a]),
-@*         Grobcov::@RP  (Q[x,a])
+@*         @R   (Q[a][x]),
+@*         @P   (Q[a]),
+@*         @RP  (Q[x,a])
 @*         that are used inside and killed before the output.
-@*         If you want to use some internal routine you must
-@*         create before the above rings by calling setglobalrings();
-@*         because some of the internal routines use these rings.
-@*         The call to the basic routines grobcov, cgsdr will
-@*         kill these rings.
 
 PROCEDURES:
@@ -109 +107 @@
              (the own routine of 2010 that is no more used).
              Now, KSW algorithm is used.
-
-setglobalrings();  Generates the global rings @R, @P and @PR that
-              are respectively the rings Q[a][x], Q[a], Q[x,a].  (a=parameters,
-             x=variables) It is called inside each of the fundamental
-             routines of the library: grobcov, cgsdr, locus, locusdg and
-             killed before the output.
 
 pdivi(f,F); Performs a pseudodivision of a parametric polynomial
@@ -138 +130 @@
              the bases are computed, and one can obtain the full
              representation using extend.
+
+ConsLevels(L); Given a list L of locally closed sets, it returns the canonical levels
+             of the constructible set of the union of them, as well as the levels
+             of the complement. It is described in the paper
+
+             J.M. Brunat, A. Montes,
+            \"Computing the canonical representation of constructible sets\".
+            Submited to Mathematics in Computer Science. July 2015.
 
 locus(G);    Special routine for determining the geometrical locus of points
@@ -178 +178 @@
              \''Envelops in Dynamic Geometry using the Gr\"obner cover\''.
 
-
 envelopdg(ev); Is a special routine to determine the 'Relevant' components
              of the envelop of a family of curves to be used in Dynamic Geometry.
              It must be called to the output of envelop(F,C).
 
-locusto(L); Transforms the output of locus, locusdg, envelop and  envelopdg
+locusto(L); Transforms the output of locus, locusdg, envelop and envelopdg
              into a string that can be reed from different computational systems.
 
-AddCons(L); Uses the routine AddConsP. Given a set of locally closed sets as
-             difference of of varieties (does not need to be in C-representation)
-             it returns the canonical P-representation of the constructible set
-             formed by the union of them. The result is formed by a set of embedded,
-             disjoint, locally closed sets (levels).
-
-AddConsP(L);  Given a set of locally closed P-components, it returns the
-             canonical P-representation of the constructible set
-             formed by the union of them, consisting of a set of embedded,
-             disjoint locally closed sets (levels).
-             The routines AddCons and AddConsP and the canonical structure
-             of the constructible sets will be described in a forthcoming paper.
-
-             A. Montes, J.M. Brunat,
-             \"Canonical representations of constructible sets\".
 
 SEE ALSO: compregb_lib
@@ -228 +212 @@
 // Uses KSW algorithm for cgsdr
 // Final data: 21-11-2013
-// Release K: (public)
-// Updated locus: 8-1-2015
-// Updated AddConsP and AddCons: 8-1-2015
-// Reformed many routines, examples and helps: 8-1-2015
+// Release L: (public)
+// New routine ConsLevels: 10-7-2015
+// Updated locus: 10-7-2015 (uses Conslevels)
 // New routines for computing the envelop of a family of curves: 22-1-2015
-// Final data: 22-1-2015
+// Final data: 22-7-2015
 
 //*************Auxiliary routines**************
@@ -297 +280 @@
   if (size(l)==1 and i==1){return(L);}
   // L=l[1];
-  if(i==1)
-  {
-    for(j=2;j<=size(l);j++)
-    {
-      L[j-1]=l[j];
-    }
-  }
-  else
+  if(i>1)
   {
     for(j=1;j<=i-1;j++)
     {
-      L[j]=l[j];
-    }
-    for(j=i+1;j<=size(l);j++)
-    {
-      L[j-1]=l[j];
-    }
+      L[size(L)+1]=l[j];
+    }
+  }
+  for(j=i+1;j<=size(l);j++)
+  {
+    L[size(L)+1]=l[j];
   }
   return(L);
@@ -773 +749 @@
 //}
 
-proc setglobalrings()
-"USAGE:   setglobalrings();
-          No arguments
-RETURN:   After its call the rings @R=Q[a][x], @P=Q[a], @RP=Q[x,a] are
-          defined as global variables.  (a=parameters, x=variables)
-NOTE: It is called internally by many basic routines of the
-          library grobcov, cgsdr, extend, pdivi, pnormalf, locus, locusdg,
-          envelop, envelopdg, and killed before the output.
-          The user does not need to call it, except when it is interested
-          in using some internal routine of the library that
-          uses these rings.
-          The basering R, must be of the form Q[a][x], (a=parameters,
-          x=variables), and should be defined previously.
-KEYWORDS: ring, rings
-EXAMPLE:  setglobalrings; shows an example"
+static proc setglobalrings()
+// "USAGE:   setglobalrings();
+//           No arguments
+// RETURN: After its call the rings Grobcov::@R=Q[a][x], Grobcov::@P=Q[a],
+//           Grobcov::@RP=Q[x,a] are defined as global variables.
+//           (a=parameters, x=variables)
+// NOTE: It is called internally by many basic routines of the
+//           library grobcov, cgsdr, extend, pdivi, pnormalf, locus, locusdg,
+//           envelop, envelopdg, and killed before the output.
+//           The user does not need to call it, except when it is interested
+//           in using some internal routine of the library that
+//           uses these rings.
+//           The basering R, must be of the form Q[a][x], (a=parameters,
+//           x=variables), and should be defined previously.
+// KEYWORDS: ring, rings
+// EXAMPLE:  setglobalrings; shows an example"
 {
   if (defined(@P))
@@ -810 +787 @@
   setring(RR);
 };
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring R=(0,a,b),(x,y,z),dp;
-  setglobalrings();
-  R;
-  Grobcov::@R;
-  Grobcov::@P;
-  Grobcov::@RP;
-  ringlist(Grobcov::@R);
-  ringlist(Grobcov::@P);
-  ringlist(Grobcov::@RP);
-}
+// example
+// {
+//   "EXAMPLE:"; echo = 2;
+//   ring R=(0,a,b),(x,y,z),dp;
+//   setglobalrings();
+//   " ";"R=";R;
+//   " ";"Grobcov::@R=";Grobcov::@R;
+//   " ";"Grobcov::@P=";Grobcov::@P;
+//   " ";"Grobcov::@RP=";Grobcov::@RP;
+//  " "; "ringlist(Grobcov::@R)=";  ringlist(Grobcov::@R);
+//  " "; "ringlist(Grobcov::@P)=";  ringlist(Grobcov::@P);
+//  " "; "ringlist(Grobcov::@RP)=";  ringlist(Grobcov::@RP);
+// }
 
 // cld : clears denominators of an ideal and normalizes to content 1
@@ -1515 +1492 @@
     }
   }
+  //"T_abans="; prep;
   if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));}
+  //"T_Prep="; prep;
+  //def Lout=CompleteA(prep,prep);
+  //"T_Lout="; Lout;
   return(prep);
 }
@@ -3934 +3915 @@
 
 //********************* End KapurSunWang *************************
+
+//********************* Begin ConsLevels ***************************
+
+static proc zeroone(int n)
+{
+  list L; list L2;
+  intvec e; intvec e2; intvec e3;
+  int j;
+  if(n==1)
+  {
+    e[1]=0;
+    L[1]=e;
+    e[1]=1;
+    L[2]=e;
+    return(L);
+  }
+  if(n>1)
+  {
+    L=zeroone(n-1);
+    for(j=1;j<=size(L);j++)
+    {
+      e2=L[j];
+      e3=e2;
+      e3[size(e3)+1]=0;
+      L2[size(L2)+1]=e3;
+      e3=e2;
+      e3[size(e3)+1]=1;
+      L2[size(L2)+1]=e3;
+    }
+  }
+  return(L2);
+}
+
+// Auxiliary routine
+// subsets: the list of subsets of (1,..n)
+static proc subsets(int n)
+{
+  list L; list L1;
+  int i; int j;
+  L=zeroone(n);
+  intvec e; intvec e1;
+  for(i=1;i<=size(L);i++)
+  {
+    e=L[i];
+    kill e1; intvec e1;
+    for(j=1;j<=n;j++)
+    {
+      if(e[n+1-j]==1)
+      {
+        if(e1==intvec(0)){e1[1]=j;}
+        else{e1[size(e1)+1]=j};
+      }
+    }
+    L1[i]=e1;
+  }
+  return(L1);
+}
+
+// Input a list A of locally closed sets in C-rep
+// Output a list B of a simplified list of A
+static proc SimplifyUnion(list A)
+{
+  int i; int j;
+  list L=A;
+  int n=size(L);
+  if(n<2){return(A);}
+  for(i=1;i<=size(L);i++)
+  {
+    for(j=1;j<=size(L);j++)
+    {
+      if(i != j)
+      {
+        if(equalideals(L[i][2],L[j][1])==1)
+        {
+          L[i][2]=L[j][2];
+        }
+      }
+    }
+  }
+  ideal T=ideal(1);
+  intvec v;
+  for(i=1;i<=size(L);i++)
+  {
+    if(equalideals(L[i][2],ideal(1)))
+    {
+      v[size(v)+1]=i;
+      T=intersect(T,L[i][1]);
+    }
+  }
+  if(size(v)>0)
+  {
+    for(i=1; i<=size(v);i++)
+    {
+      j=v[size(v)+1-i];
+      L=elimfromlist(L, j);
+    }
+  }
+  if(equalideals(T,ideal(1))==0){L[size(L)+1]=list(std(T),ideal(1))};
+  //string("T_n=",n," new n0",size(L));
+  return(L);
+}
+
+// Input: list(A)
+//          A is a list of locally closed sets in Crep. A=[[P1,Q1],[P2,Q2],..,[Pr,Qr]]
+//          whose union is a constructible set from
+// Output list [Lev,C]:
+//          where Lev is the Crep of the canonical first level of A, and
+//          C is the complement of the first level Lev wrt to the closure of A. The elements are given in Crep,
+static proc FirstLevel(list A)
+{
+  int n=size(A);
+  list T=zeroone(n);
+  ideal P; ideal Q;
+  list Cb;  ideal Cc=ideal(1);
+  int i; int j;
+  intvec t;
+  ideal Z=ideal(1);
+  list C;
+  for(i=1;i<=size(A);i++)
+  {
+    Z=intersect(Z,A[i][1]);
+  }
+  for(i=2; i<=size(T);i++)
+  {
+    t=T[i];
+    P=ideal(1); Q=ideal(0);
+    for(j=1;j<=n;j++)
+    {
+      if(t[n+1-j]==1)
+      {
+        Q=Q+A[j][2];
+      }
+      else
+      {
+        P=intersect(P,A[j][1]);
+      }
+    }
+    //"T_Q="; Q; "T_P="; P;
+    Cb=Crep(Q,P);
+    //"T_Cb="; Cb;
+    if(Cb[1][1]<>1)
+    {
+      C[size(C)+1]=Cb;
+      Cc=intersect(Cc,Cb[1]);
+    }
+  }
+  list Lev=list(Z,Cc);                //Crep(Z,Cc);
+  if(size(C)>1){C=SimplifyUnion(C);}
+  return(list(Lev,C));
+}
+
+// Input: list(A)
+//          A is a list of locally closed sets in Crep. A=[[P1,Q1],[P2,Q2],..,[Pr,Qr]]
+//          whose union is a constructible set from which the algorithm computes its canonical form.
+// Output:
+//     list [L,C] where
+//          where L is the list of canonical levels of A,
+//          and C is the list of canonical levels of the complement of A wrt to the closure of A.
+//          All locally closed sets are given in Crep.
+//          L=[[1,[p1,p2],[3,[p3,p4],..,[2r-1,[p_{2r-1},p_2r]]]] is the list of levels of A in Crep.
+//          C=[[2,p2,p3],[4,[p4,p5],..,[2s,[p_{2s},p_{2s+1}]]]  is the list of levels of C,
+//                                              the complement of S wrt the closure of A, in Crep.
+proc ConsLevels(list A)
+"USAGE:   ConsLevels(A);
+          A is a list of locally closed sets in Crep. A=[[P1,Q1],[P2,Q2],..,[Pr,Qr]]
+          whose union is a constructible set from which the algorithm computes its
+          canonical form.
+RETURN: A list with two elements: [L,C]
+          where L is the list of canonical levels of A, and C is the list of canonical
+          levels of the
+          complement of A wrt to the closure of A.
+          All locally closed sets are given in Crep.
+          L=[[1,[p1,p2],[3,[p3,p4],..,[2r-1,[p_{2r-1},p_2r]]]]
+          C=[[2,p2,p3],[4,[p4,p5],..,[2s,[p_{2s},p_{2s+1}]]]
+NOTE: The basering R, must be of the form Q[a]
+KEYWORDS: locally closed set, constructible set
+EXAMPLE:  ConsLevels; shows an example"
+{
+  list L; list C; int i;
+  list B; list T;
+  for(i=1; i<=size(A);i++)
+  {
+    T=Crep(A[i][1],A[i][2]);
+    B[size(B)+1]=T;
+  }
+  int level=0;
+  list K;
+  while(size(B)>0)
+  {
+    level++;
+    K=FirstLevel(B);
+    if(level mod 2==1){L[size(L)+1]=list(level,K[1]);}
+    else{C[size(C)+1]=list(level,K[1]);}
+    //"T_L="; L;
+    //"T_C="; C;
+    B=K[2];
+    //"T_size(B)="; size(B);
+  }
+  return(list(L,C));
+}
+example
+{ "EXAMPLE:"; echo = 2;
+// Example of ConsLevels with nice geometrical interpretetion and 2 levels
+
+if (defined(R)){kill R;}
+if (defined(@P)){kill @P; kill @R; kill @RP;}
+
+  ring R=0,(x,y,z),lp;
+  short=0;
+  ideal P1=x*(x^2+y^2+z^2-1);
+  ideal Q1=z,x^2+y^2-1;
+  ideal P2=y,x^2+z^2-1;
+  ideal Q2=z*(z+1),y,x*(x+1);
+
+  list Cr1=Crep(P1,Q1);
+  list Cr2=Crep(P2,Q2);
+
+  list L=list(Cr1,Cr2);
+  L;
+  // ConsLevels(L)=
+  ConsLevels(L);
+
+//----------------------------------------------------------------------
+// Begin Problem S93
+// Automatic theorem proving
+// Generalized Steiner-Lehmus Theorem
+// A.Montes, T.Recio
+
+// Bisector of A(-1,0) = Bisector of B(1,0) varying C(a,b)
+
+if(defined(R1)){kill R1;}
+ring R1=(0,a,b),(x1,y1,x2,y2,p,r),dp;
+ideal S93=(a+1)^2+b^2-(p+1)^2,
+                    (a-1)^2+b^2-(1-r)^2,
+                    a*y1-b*x1-y1+b,
+                    a*y2-b*x2+y2-b,
+                    -2*y1+b*x1-(a+p)*y1+b,
+                    2*y2+b*x2-(a+r)*y2-b,
+                    (x1+1)^2+y1^2-(x2-1)^2-y2^2;
+
+short=0;
+def GC93=grobcov(S93,"nonnull",ideal(b),"rep",1);
+//"grobcov(S93,'nonnull',ideal(b),"rep",1)="; GC93;
+
+list L0;
+for(int i=1;i<=size(GC93);i++)
+{
+   L0[size(L0)+1]=GC93[i][3];
+}
+
+def GC93ab=grobcov(S93,"nonnull",ideal(a*b),"rep",1);
+
+ring RR=0,(a,b),lp;
+
+list L1;
+L1=imap(R1,L0);
+// L1=Total elements of the grobcov(S93,'nonnull',ideal(b),'rep',1) to be added=
+L1;
+
+// Adding all the elements of grobcov(S93,'nonnull',ideal(b),'rep',1)=
+ConsLevels(L1);
+}
+
+//**************************** End ConsLevels ******************
 
 //******************** Begin locus ******************************
@@ -4522 +4767 @@
   locus(grobcov(S));
   kill R;
-  "********************************************";
+  //********************************************
 
   ring R=(0,x,y),(x1,x2),dp;
@@ -4611 +4856 @@
   locusdg(locus(grobcov(S96)));
   kill R;
-  "********************************************";
+  //********************************************
   ring R=(0,a,b),(x4,x3,x2,x1),dp;
   ideal S=(x1-3)^2+(x2-1)^2-9,
@@ -4622 +4867 @@
   locusdg(locus(grobcov(S)));
   kill R;
-  "********************************************";
+  //********************************************
 
   ring R=(0,x,y),(x1,x2),dp;
@@ -4633 +4878 @@
 }
 
-// locusto: Transforms the output of locus to a string that
-//      can be read by different computational systems.
+// locusto: Transforms the output of locus, locusdg, envelop and envelopdg
+//             into a string that can be reed from different computational systems.
 // input:
 //     list L: The output of locus
@@ -4717 +4962 @@
   locusto(locusdg(locus(grobcov(S))));
   kill R;
-  "********************************************";
+  //********************************************
 
   // 1. Take a fixed line l: x1-y1=0  and consider
@@ -4737 +4982 @@
   locusto(envelopdg(envelop(F,C)));
   kill R;
-  "********************************************";
+  //********************************************
 
   // Steiner Deltoid
@@ -4804 +5049 @@
   return(d);
 }
-
-// // locusdgto: Transforms the output of locusdg to a string that
-// //      can be read by different computational systems.
-// // input:
-// //     list L: The output of locus
-// // output:
-// //     string s: The output of locus converted to a string readable by other programs
-// //                   Outputs only the relevant dynamical geometry components.
-// //                   Without unnecessary parenthesis
-// proc locusdgto(list LL)
-// "USAGE: locusdgto(L);
-//           The argument must be the output of locusdg of a parametrical ideal
-//           It transforms the output into a string in standard form
-//           readable in many languages (Geogebra).
-// RETURN: The locusdg in string standard form
-// NOTE: It can only be called after computing the locusdg(grobcov(F)) of the
-//           parametrical ideal.
-//           The basering R, must be of the form Q[a,b,..][x,y,..].
-// KEYWORDS: geometrical locus, locus, loci.
-// EXAMPLE:  locusdgto; shows an example"
-// {
-//   int i; int j; int k; int short0=short; int ojo;
-//   int te=0;
-//   short=0;
-//   if(size(LL)==0){ojo=1; list L;}
-//   else
-//   {
-//     def L=LL;
-//   }
-//   string s="["; string sf="]"; string st=s+sf;
-//   if(size(L)==0){return(st);}
-//   ideal p;
-//   ideal q;
-//   for(i=1;i<=size(L);i++)
-//   {
-//     if(L[i][3]=="Relevant")
-//     {
-//       s=string(s,"[[");
-//       for (j=1;j<=size(L[i][1]);j++)
-//       {
-//         s=string(s,L[i][1][j],",");
-//       }
-//       s[size(s)]="]";
-//       s=string(s,",[");
-//       for(j=1;j<=size(L[i][2]);j++)
-//       {
-//         s=string(s,"[");
-//         for(k=1;k<=size(L[i][2][j]);k++)
-//         {
-//           s=string(s,L[i][2][j][k],",");
-//         }
-//         s[size(s)]="]";
-//         s=string(s,",");
-//       }
-//       s[size(s)]="]";
-//       s=string(s,"]");
-//       s[size(s)]="]";
-//       s=string(s,",");
-//     }
-//   }
-//   if(s=="["){s="[]";}
-//   else{s[size(s)]="]";}
-//   short=short0;
-//   return(s);
-// }
-// example
-// {"EXAMPLE:"; echo = 2;
-//   ring R=(0,a,b),(x,y),dp;
-//   short=0;
-//   ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
-//   "System="; S96; " ";
-//   "locusdgto(locusdg(grobcov(S96)))=";
-//   locusdgto(locusdg(grobcov(S96)));
-// }
 
 static proc norspec(ideal F)
@@ -4898 +5069 @@
   exportto(Top,@RP);     // global ring K[x,a] with product order
   setring(RR);
-
 }
 
@@ -5059 +5229 @@
 }
 
-//********************* End locus ****************************
-
-//********************* Begin AddCons **********************
-
-// Input: L1,L2: lists of components with common top
-// Output L: list of the union of L1 and L2.
-static proc Add2ComWithCommonTop(list L1, list L2)
-{
-  int i; int j; ideal pij; list L; list Lp; list PR; int k;
-  for(i=1;i<=size(L1[2]);i++)
-  {
-    for(j=1;j<=size(L2[2]);j++)
-    {
-      pij=std(L1[2][i]+L2[2][j]);
-      PR=minGTZ(pij);
-      for(k=1;k<=size(PR);k++)
-      {
-        Lp[size(Lp)+1]=PR[k];
-      }
-    }
-  }
-  for(i=1; i<=size(Lp);i++)
-  {
-    for(j=i+1;j<=size(Lp);j++)
-    {
-      if(idcontains(Lp[i],Lp[j])) {Lp=delete(Lp,j);}
-    }
-    for(j=1;j<i;j++)
-    {
-      if(idcontains(Lp[i],Lp[j])){Lp=delete(Lp,j); j=j-1; i=i-1;}
-    }
-  }
-  L[1]=L1[1];
-  L[2]=Lp;
-  return(L);
-}
-
-// Input: L list od P-rep of components to be added. L[i]=[p_i,[p_{i1},...p_{ir_i}]]
-// Output: lists A,B,L
-//       where no top in the lists are repeated
-//       A: list of components with higher tops
-//       B: list of components with lower tops
-//       L1: A,B
-static proc SepareAB(list L)
-{
-  int i;  int j; list Ln=L;
-  for(i=1;i<=size(Ln);i++)
-  {
-    for(j=i+1;j<=size(Ln);j++)
-    {
-      if (equalideals(Ln[j][1],Ln[i][1]))
-      {
-        Ln[i]=Add2ComWithCommonTop(Ln[i],Ln[j]);
-        Ln=delete(Ln,j);
-        j=j-1;
-      }
-    }
-  }
-  list A; list B; int clas; list T; list L1;
-  for(i=1;i<=size(Ln);i++)
-  {
-    j=1;
-    clas=0;
-    while((clas==0) and  (j<=size(Ln)))
-    {
-      if(j!=i)
-      {
-        if(idcontains(Ln[i][1],Ln[j][1]) ) {B[size(B)+1]=Ln[i]; clas=1;}
-      }
-      j++;
-    }
-    if(clas==0) {A[size(A)+1]=Ln[i];}
-  }
-  L1=A; for(j=1;j<=size(B);j++){L1[size(L1)+1]=B[j];}
-  T[1]=A; T[2]=B; T[3]=L1;
-  return(T);
-}
-
-// Input:
-//  A1: list of high set of P-reps to be completed by the remaining P-reps
-//  L1: the list A1, completed with the list of lower P-reps.
-// Output:
-//  A: list A1 completed with all possible parts of the remaining parts of L1 with the
-//      condition of building locally closed sets.
-static proc CompleteA(list A1,list L1)
-{
-  int modif; int i; int j; int k; int l;
-  ideal pij; ideal pk; ideal pijkl; ideal pkl;
-  int n; list nl; int te; int clas; list vvv; list AAA;
-  list Lp; int m; ideal Pij;
-  list A0;
-  modif=1;
-  list A=A1;
-  while(modif==1)
-  {
-      modif=0;
-      A0=A;
-      for(i=1;i<=size(A);i++)
-      {
-          for(j=1;j<=size(A[i][2]); j++)
-          {
-              pij=A[i][2][j];
-             for(k=1;k<=size(L1);k++)
-             {
-                 if(k!=i)
-                 {
-                     pk=L1[k][1];
-                     if(idcontains(pij,pk)==1)
-                     {
-                         te=0;
-                         kill nl;
-                         list nl; l=1;
-                         while((te==0) and (l<=size(L1[k][2])))
-                         {
-                              pkl=L1[k][2][l];
-                              if((equalideals(pij,pkl)==1) or (idcontains(pij,pkl)==1)) {te=1;}
-                              l++;
-                         }   // end while ((te=0) and (l>...
-                         //"T_te="; te; pause();
-                         if(te==0)
-                         {
-                           modif=1;
-                           kill Pij; ideal Pij=1;
-                           for(l=1; l<=size(L1[k][2]);l++)
-                           {
-                             pkl=L1[k][2][l];
-                             pijkl=std(pij+pkl);
-                             Pij=intersect(Pij,pijkl);
-                           }
-                           kill Lp; list Lp;
-                           Lp=minGTZ(Pij);
-                           for(m=1;m<=size(Lp);m++)
-                            {
-                               nl[size(nl)+1]=Lp[m];
-                            }
-                            A[i][2]=delete(A[i][2], j);
-                            for(n=1;n<=size(nl);n++)
-                            {
-                              A[i][2][size(A[i][2])+1]=nl[n];
-                            }
-                          } // end if(te==0)
-                     } // end if(idcontains(pij,pk)==1)
-                 }  // end if (k!=i)
-             } // end for k
-         }  // end for j
-         kill vvv; list vvv;
-         if(modif==1)
-         // Select the maximal ideals of the set A[I][2][j]
-         {
-             kill nl; list nl;
-             nl=selectminideals(A[i][2]);
-             kill AAA; list AAA;
-             for(j=1;j<=size(nl);j++)
-             {
-               AAA[size(AAA)+1]=A[i][2][nl[j]];
-             }
-             A[i][2]=AAA;
-         } // end if(modif=1)
-      }  // end for i
-      modif=1-equallistsAall(A,A0);
-  } // end while(modif==1)
-  return(A);
-}
-
-// Input:
-//   A: list of the top P-reps of one level
-//   B: list of remaining lower P-reps that have not yeen be possible to add
-// Output:
-//   Bn: list B where the elements that are completely included in A are removed and the parts that are
-//         included in A also.
-static proc ReduceB(list A,list B)
-{
-     int i; int j; list nl; list Bn; int te; int k; int elim;
-     ideal pC; ideal pD; ideal pCD; ideal pBC; list nB; int j0;
-     for(i=1;i<=size(B);i++)
-     {
-         j=1; te=0;
-         while((te==0) and (j<=size(A)))
-         {
-             if(idcontains(B[i][1],A[j][1])){te=1; j0=j;}
-             else{j++;}
-         }
-         pD=B[i][2][1];
-         for(k=2;k<=size(B[i][2]);k++){pD=intersect(pD,B[i][2][k]);}
-         pC=A[j0][2][1];
-         for(k=2;k<=size(A[j0][2]);k++) {pC=intersect(pC,A[j0][2][k]);}
-         pCD=std(pD+pC);
-         pBC=std(B[i][1]+pC);
-         elim=0;
-         if(idcontains(pBC,pCD)){elim=1;}   // B=delfromlist(B,i);i=i-1;
-         else
-         {
-              nB=Prep(pBC,pCD);
-              if(equalideals(nB[1][1],ideal(1))==0)
-              {
-                  for(k=1;k<=size(nB);k++){Bn[size(Bn)+1]=nB[k];}
-              }
-         }
-    }   // end for i
-    return(Bn);
-}
-
-// AddConsP: given a set of components of locally closed sets in P-representation, it builds the
-//       canonical P-representation of the corresponding constructible set, of its union,
-//       including levels it they are.
-proc AddConsP(list L)
-"USAGE:   AddConsP(L)
-      Input L: list of components in P-rep to be added
-      [  [[p_1,[p_11,..,p_1,r1]],..[p_k,[p_k1,..,p_kr_k]]  ]
-RETURN:
-     list of lists of levels of the different locally closed sets of
-     the canonical P-rep of the constructible.
-     [  [level_1,[ [Comp_11,..Comp_1r_1] ] ], .. ,
-        [level_s,[ [Comp_s1,..Comp_sr_1] ]
-     ]
-     where level_i=i,   Comp_ij=[ p_i,[p_i1,..,p_it_i] ] is a prime component.
-NOTE:     Operates in a ring R=Q[u_1,..,u_m]
-KEYWORDS: Constructible set, Canoncial form
-EXAMPLE:  AddConsP; shows an example"
-{
-  list LL; list A; list B; list L1; list T; int level=0; list h;
-  LL=L; int i;
-  while(size(LL)!=0)
-  {
-    level++;
-    L1=SepareAB(LL);
-    A=L1[1]; B=L1[2]; LL=L1[3];
-    A=CompleteA(A,LL);
-    for(i=1;i<=size(A);i++)
-    {
-      LL[i]=A[i];
-    }
-    h[1]=level; h[2]=A;
-    T[size(T)+1]=h;
-    LL=ReduceB(A,B);
-  }
-  return(T);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  if (defined(Grobcov::@P)){kill Grobcov::@P; kill Grobcov::@R; kill Grobcov::@RP;}
-  ring R=0,(x,y,z),lp;
-  short=0;
-
-  ideal P1=x;
-  ideal Q1=x,y;
-  ideal P2=y;
-  ideal Q2=y,z;
-
-  list L=list(Prep(P1,Q1)[1],Prep(P2,Q2)[1]);
-  L;
-  AddConsP(L);
-}
-
-// AddCons:  given a set of  locally closed sets by pairs of ideal, it builds the
-//       canonical P-representation of the corresponding constructible set, of its union,
-//       including levels it they are.
-//       It makes a call to AddConsP after transforming the input.
-// Input list of lists of pairs of ideals representing locally colsed sets:
-//     L=  [ [E1,N1], .. , [E_s,N_s] ]
-// Output: The canonical frepresentation of the constructible set union of the V(E_i) \ V(N_i)
-//     T=[  [level_1,[ [p_1,[p_11,..,p_1s_1]],.., [p_k,[p_k1,..,p_ks_k]] ]],, .. , [level_r,[..       ]]  ]
-proc AddCons(list L)
-"USAGE:   AddCons(L)
-      Calls internally AddConsP and allows a different input.
-      Input L: list of pairs of of ideals [ [P_1,Q_1], .., [Pr,Qr] ]
-      representing a set of locally closed setsV(P_i) \ V(Q_i)
-      to be added.
-RETURN:
-      list of lists of levels of the different locally closed sets of
-      the canonical P-rep of the constructible.
-      [  [level_1,[ [Comp_11,..Comp_1r_1] ] ], .. ,
-         [level_s,[ [Comp_s1,..Comp_sr_1] ]
-      ]
-      where level_i=i,   Comp_ij=[ p_i,[p_i1,..,p_it_i] ] is a prime component.
-NOTE: Operates in a ring R=Q[u_1,..,u_m]
-KEYWORDS: Constructible set, Canoncial form
-EXAMPLE: AddCons; shows an example"
-{
-  int i; list H; list P; int j;
-  for(i=1;i<=size(L);i++)
-  {
-    P=Prep(L[i][1],L[i][2]);
-    for(j=1;j<=size(P);j++)
-    {
-      H[size(H)+1]=P[j];
-    }
-  }
-  return(AddConsP(H));
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  if (defined(Grobcov::@P)){kill Grobcov::@P; kill Grobcov::@R; kill Grobcov::@RP;}
-  ring R=0,(x,y,z),lp;
-  short=0;
-
-  ideal P1=x;
-  ideal Q1=x,y;
-  ideal P2=y;
-  ideal Q2=y,z;
-
-  list L=list(list(P1,Q1), list(P2,Q2) );
-  L;
-
-  AddCons(L);
-}
-
 // AddLocus: auxilliary routine for locus0 that computes the components of the constructible:
 // Input:  the list of locally closed sets to be added, each with its type as third argument
@@ -5397 +5258 @@
     }
   }
-  L3=AddConsP(L1);
+  L3=LocusConsLevels(L1);
   list L4; int level;
   ideal p1; ideal pp1; int t; int k; int k0; string typ; list l4;
@@ -5426 +5287 @@
 }
 
-//********************* End AddCons **********************
-;
+// Input L: list of components in P-rep to be added
+//         [  [[p_1,[p_11,..,p_1,r1]],..[p_k,[p_k1,..,p_kr_k]]  ]
+// Output:
+//          list of lists of levels of the different locally closed sets of
+//          the canonical P-rep of the constructible.
+//          [  [level_1,[ [Comp_11,..Comp_1r_1] ] ], .. ,
+//             [level_s,[ [Comp_s1,..Comp_sr_1] ]
+//          ]
+//          where level_i=i,   Comp_ij=[ p_i,[p_i1,..,p_it_i] ] is a prime component.
+// LocusConsLevels: given a set of components of locally closed sets in P-representation, it builds the
+//       canonical P-representation of the corresponding constructible set of its union,
+//       including levels it they are.
+static proc LocusConsLevels(list L)
+{
+  list Lc; list Sc;
+  int i;
+  for(i=1;i<=size(L);i++)
+  {
+    Sc=PtoCrep(list(L[i]));
+    Lc[size(Lc)+1]=Sc;
+  }
+  list S=ConsLevels(Lc)[1];
+  list Sout;
+  list Lev;
+  for(i=1;i<=size(S);i++)
+  {
+    Lev=list(i,Prep(S[i][2][1],S[i][2][2]));
+    Sout[size(Sout)+1]=Lev;
+  }
+  return(Sout);
+}
+
+//********************* End locus ****************************
+

Singular/LIB/grwalk.lib

@@ -255 +255 @@
 }
 
-proc gwalk(ideal Go, list #)
-"SYNTAX: gwalk(ideal i);
-         gwalk(ideal i, intvec v, intvec w);
+proc gwalk(ideal Go, int reduction,int printout, list #)
+"SYNTAX: gwalk(ideal i, int reduction, int printout);
+         gwalk(ideal i, int reduction, int printout, intvec v, intvec w);
 TYPE:    ideal
 PURPOSE: compute the standard basis of the ideal, calculated via
@@ -274 +274 @@
    string ord_str =   OSCTW[2];
    intvec curr_weight   =   OSCTW[3]; /* original weight vector */
-   intvec target_weight =   OSCTW[4]; /* terget weight vector */
+   intvec target_weight =   OSCTW[4]; /* target weight vector */
    kill OSCTW;
    option(redSB);
@@ -284 +284 @@
    //print("//** help ring = " + string(basering));
    ideal G = fetch(xR, Go);
-   G = system("Mwalk", G, curr_weight, target_weight,basering);
+   G = system("Mwalk", G, curr_weight, target_weight,basering,reduction,printout);
 
    setring xR;
@@ -299 +299 @@
   //** compute a Groebner basis of I w.r.t. lp.
   ring r = 32003,(z,y,x), lp;
-  ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
-  gwalk(I);
+  ideal I = zy2+yx2+yx+3,
+            z3x+y3+zyx-yx2-yx-3,
+            z2yx3-y5+z2yx2+y3x2+y2x3+y3x+y2x2+3z2x+3y2+3yx,
+            zyx5+y6-y4x2-y3x3+2zyx4-y4x-y3x2+zyx3-3z2yx+3zx3-3y3-3y2x+3zx2,
+            yx7-y7+y5x2+y4x3+3yx6+y5x+y4x2+3yx5-6zyx3+yx4+3x5+3y4+3y3x-6zyx2+6x4+3x3-9zx;
+  gwalk(I,0,1);
 }
 
@@ -346 +350 @@
 }
 
-proc fwalk(ideal Go, list #)
-"SYNTAX: fwalk(ideal i);
-         fwalk(ideal i, intvec v, intvec w);
+proc fwalk(ideal Go, int reduction, int printout, list #)
+"SYNTAX: fwalk(ideal i,int reductioin);
+         fwalk(ideal i, int reduction intvec v, intvec w);
 TYPE:    ideal
 PURPOSE: compute the standard basis of the ideal w.r.t. the
@@ -372 +376 @@
 
    ideal G = fetch(xR, Go);
-   G = system("Mfwalk", G, curr_weight, target_weight);
+   G = system("Mfwalk", G, curr_weight, target_weight, reduction, printout);
 
    setring xR;
@@ -387 +391 @@
     ring r = 32003,(z,y,x), lp;
     ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
-    fwalk(I);
+    int reduction = 1;
+    int printout = 1;
+    fwalk(I,reduction,printout);
 }
 
@@ -437 +443 @@
 }
 
-proc pwalk(ideal Go, int n1, int n2, list #)
-"SYNTAX: pwalk(int d, ideal i, int n1, int n2);
-         pwalk(int d, ideal i, int n1, int n2, intvec v, intvec w);
+proc pwalk(ideal Go, int n1, int n2, int reduction, int printout, list #)
+"SYNTAX: pwalk(int d, ideal i, int n1, int n2, int reduction, int printout);
+         pwalk(int d, ideal i, int n1, int n2, int reduction, int printout, intvec v, intvec w);
 TYPE:    ideal
 PURPOSE: compute the standard basis of the ideal, calculated via
@@ -477 +483 @@
   ideal G = fetch(xR, Go);
 
-  G = system("Mpwalk", G, n1, n2, curr_weight, target_weight,nP);
-
+  G = system("Mpwalk",G,n1,n2,curr_weight,target_weight,nP,reduction,printout);
+ 
   setring xR;
-  //kill Go;
+  //kill Go; //unused
 
   keepring basering;
@@ -492 +498 @@
     ring r = 32003,(z,y,x), lp;
     ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
-    //I = std(I);
-    //ring rr = 32003,(z,y,x),lp;
-    //ideal I = fetch(r,I);
-    pwalk(I,2,2);
+    int reduction = 1;
+    int printout = 2;
+    pwalk(I,2,2,reduction,printout);
 }
 

Singular/LIB/modstd.lib

@@ -283 +283 @@
  * (same as size(reduce(I_reduce, G_reduce))).
  * Uses parallelization. */
-static proc reduce_parallel(ideal I_reduce, ideal G_reduce)
+static proc reduce_parallel(def I_reduce, def G_reduce)
 {
     exportto(Modstd, I_reduce);

Singular/LIB/modwalk.lib

@@ -16 +16 @@
 
 PROCEDURES:
-modWalk(ideal,int,int[,int,int,int,int]);        standard basis conversion of I using modular methods (chinese remainder)
+
+modWalk(I,#);                   standard basis conversion of I by Groebner Walk using modular methods
+modrWalk(I,radius,pertdeg,#);   standard basis conversion of I by Random Walk using modular methods
+modfWalk(I,#);                  standard basis conversion of I by Fractal Walk using modular methods
+modfrWalk(I,radius,#);          standard basis conversion of I by Random Fractal Walk using modular methods
 ";
 
-LIB "poly.lib";
-LIB "ring.lib";
-LIB "parallel.lib";
 LIB "rwalk.lib";
 LIB "grwalk.lib";
-LIB "modstd.lib";
-
-
-////////////////////////////////////////////////////////////////////////////////
-
-static proc modpWalk(def II, int p, int variant, list #)
-"USAGE:  modpWalk(I,p,#); I ideal, p integer, variant integer
-ASSUME:  If size(#) > 0, then
-           #[1] is an intvec describing the current weight vector
-           #[2] is an intvec describing the target weight vector
-RETURN:  ideal - a standard basis of I mod p, integer - p
-NOTE:    The procedure computes a standard basis of the ideal I modulo p and
-         fetches the result to the basering.
-EXAMPLE: example modpWalk; shows an example
-"
-{
-  option(redSB);
-  int k,nvar@r;
-  def R0 = basering;
-  string ordstr_R0 = ordstr(R0);
-  list rl = ringlist(R0);
-  int sizerl = size(rl);
-  int neg = 1 - attrib(R0,"global");
-  if(typeof(II) == "ideal")
-  {
-    ideal I = II;
+LIB "modular.lib";
+
+proc modWalk(ideal I, list #)
+"USAGE:   modWalk(I, [, v, w]); I ideal, v intvec, w intvec
+RETURN:   a standard basis of I
+NOTE:     The procedure computes a standard basis of I (over the rational
+          numbers) by using modular methods.
+SEE ALSO: modular
+EXAMPLE:  example modWalk; shows an example"
+{
+    /* read optional parameter */
+    if (size(#) > 0) {
+        if (size(#) == 1) {
+            intvec w = #[1];
+        }
+        if (size(#) == 2) {
+            intvec v = #[1];
+            intvec w = #[2];
+        }
+        if (size(#) > 2 || typeof(#[1]) != "intvec") {
+            ERROR("wrong optional parameter");
+        }
+    }
+
+    /* save options */
+    intvec opt = option(get);
+    option(redSB);
+
+    /* set additional parameters */
+    int reduction = 1;
+    int printout = 0;
+
+    /* call modular() */
+    if (size(#) > 0) {
+        I = modular("gwalk", list(I,reduction,printout,#));
+    }
+    else {
+        I = modular("gwalk", list(I,reduction,printout));
+    }
+
+    /* return the result */
+    attrib(I, "isSB", 1);
+    option(set, opt);
+    return(I);
+}
+example
+{
+    "EXAMPLE:";
+    echo = 2;
+    ring R1 = 0, (x,y,z,t), dp;
+    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
+    I = std(I);
+    ring R2 = 0, (x,y,z,t), lp;
+    ideal I = fetch(R1, I);
+    ideal J = modWalk(I);
+    J;
+}
+
+proc modrWalk(ideal I, int radius, int pertdeg, list #)
+"USAGE:   modrWalk(I, radius, pertdeg[, v, w]);
+          I ideal, radius int, pertdeg int, v intvec, w intvec
+RETURN:   a standard basis of I
+NOTE:     The procedure computes a standard basis of I (over the rational
+          numbers) by using modular methods.
+SEE ALSO: modular
+EXAMPLE:  example modrWalk; shows an example"
+{
+    /* read optional parameter */
+    if (size(#) > 0) {
+        if (size(#) == 1) {
+            intvec w = #[1];
+        }
+        if (size(#) == 2) {
+            intvec v = #[1];
+            intvec w = #[2];
+        }
+        if (size(#) > 2 || typeof(#[1]) != "intvec") {
+            ERROR("wrong optional parameter");
+        }
+    }
+
+    /* save options */
+    intvec opt = option(get);
+    option(redSB);
+
+    /* set additional parameters */
+    int reduction = 1;
+    int printout = 0;
+
+    /* call modular() */
+    if (size(#) > 0) {
+        I = modular("rwalk", list(I,radius,pertdeg,reduction,printout,#));
+    }
+    else {
+        I = modular("rwalk", list(I,radius,pertdeg,reduction,printout));
+    }
+
+    /* return the result */
+    attrib(I, "isSB", 1);
+    option(set, opt);
+    return(I);
+}
+example
+{
+    "EXAMPLE:";
+    echo = 2;
+    ring R1 = 0, (x,y,z,t), dp;
+    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
+    I = std(I);
+    ring R2 = 0, (x,y,z,t), lp;
+    ideal I = fetch(R1, I);
     int radius = 2;
-    int pert_deg = 2;
-  }
-  if(typeof(II) == "list" && typeof(II[1]) == "ideal")
-  {
-    ideal I = II[1];
-    if(size(II) == 2)
-    {
-      int radius = II[2];
-      int pert_deg = 2;
-    }
-    if(size(II) == 3)
-    {
-      int radius = II[2];
-      int pert_deg = II[3];
-    }
-  }
-  rl[1] = p;
-  int h = homog(I);
-  def @r = ring(rl);
-  setring @r;
-  ideal i = fetch(R0,I);
-  string order;
-  if(system("nblocks") <= 2)
-  {
-    if(find(ordstr_R0, "M") + find(ordstr_R0, "lp") + find(ordstr_R0, "rp") <= 0)
-    {
-      order = "simple";
-    }
-  }
-
-//-------------------------  make i homogeneous  -----------------------------
-  if(!mixedTest() && !h)
-  {
-    if(!((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg))
-    {
-      if(!((order == "simple") || (sizerl > 4)))
-      {
-        list rl@r = ringlist(@r);
-        nvar@r = nvars(@r);
-        intvec w;
-        for(k = 1; k <= nvar@r; k++)
-        {
-          w[k] = deg(var(k));
-        }
-        w[nvar@r + 1] = 1;
-        rl@r[2][nvar@r + 1] = "homvar";
-        rl@r[3][2][2] = w;
-        def HomR = ring(rl@r);
-        setring HomR;
-        ideal i = imap(@r, i);
-        i = homog(i, homvar);
-      }
-    }
-  }
-
-//-------------------------  compute a standard basis mod p  -----------------------------
-
-  if(variant == 1)
-  {
-    if(size(#)>0)
-    {
-      i = rwalk(i,radius,pert_deg,#);
-     // rwalk(i,radius,pert_deg,#); std(i);
-    }
-    else
-    {
-      i = rwalk(i,radius,pert_deg);
-    }
-  }
-  if(variant == 2)
-  {
-    if(size(#) == 2)
-    {
-      i = gwalk(i,#);
-    }
-    else
-    {
-      i = gwalk(i);
-    }
-  }
-  if(variant == 3)
-  {
-    if(size(#) == 2)
-    {
-      i = frandwalk(i,radius,#);
-    }
-    else
-    {
-      i = frandwalk(i,radius);
-    }
-  }
-  if(variant == 4)
-  {
-    if(size(#) == 2)
-    {
-      i=fwalk(i,#);
-    }
-    else
-    {
-      i=fwalk(i);
-    }
-  }
-  if(variant == 5)
-  {
-    if(size(#) == 2)
-    {
-     i=prwalk(i,radius,pert_deg,pert_deg,#);
-    }
-    else
-    {
-      i=prwalk(i,radius,pert_deg,pert_deg);
-    }
-  }
-  if(variant == 6)
-  {
-    if(size(#) == 2)
-    {
-      i=pwalk(i,pert_deg,pert_deg,#);
-    }
-    else
-    {
-      i=pwalk(i,pert_deg,pert_deg);
-    }
-  }
-
-  if(!mixedTest() && !h)
-  {
-    if(!((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg))
-    {
-      if(!((order == "simple") || (sizerl > 4)))
-      {
-        i = subst(i, homvar, 1);
-        i = simplify(i, 34);
-        setring @r;
-        i = imap(HomR, i);
-        i = interred(i);
-        kill HomR;
-      }
-    }
-  }
-  setring R0;
-  return(list(fetch(@r,i),p));
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  option(redSB);
-
-  int p = 181;
-  intvec a = 2,1,3,4;
-  intvec b = 1,9,1,1;
-  ring ra = 0,x(1..4),(a(a),lp);
-  ideal I = std(cyclic(4));
-  ring rb = 0,x(1..4),(a(b),lp);
-  ideal I = imap(ra,I);
-  modpWalk(I,p,1,a,b);
-  std(I);
-}
-
-////////////////////////////////////////////////////////////////////////////////
-
-proc modWalk(def II, int variant, list #)
-"USAGE:  modWalk(II); II ideal or list(ideal,int)
-ASSUME:  If variant =
-@*       - 1 the Random Walk algorithm with radius II[2] is applied
-           to II[1] if II = list(ideal, int). It is applied to II with radius 2
-           if II is an ideal.
-@*       - 2, the Groebner Walk algorithm is applied to II[1] or to II, respectively.
-@*       - 3, the Fractal Walk algorithm with random element is applied to II[1] or II,
-           respectively.
-@*       - 4, the Fractal Walk algorithm is applied to II[1] or II, respectively.
-@*       - 5, the Perturbation Walk algorithm with random element is applied to II[1]
-           or II, respectively, with radius II[3] and perturbation degree II[2].
-@*       - 6, the Perturbation Walk algorithm is applied to II[1] or II, respectively,
-           with perturbation degree II[3].
-         If size(#) > 0, then # contains either 1, 2 or 4 integers such that
-@*       - #[1] is the number of available processors for the computation,
-@*       - #[2] is an optional parameter for the exactness of the computation,
-                if #[2] = 1, the procedure computes a standard basis for sure,
-@*       - #[3] is the number of primes until the first lifting,
-@*       - #[4] is the constant number of primes between two liftings until
-           the computation stops.
-RETURN:  a standard basis of I if no warning appears.
-NOTE:    The procedure converts a standard basis of I (over the rational
-         numbers) from the ordering \"a(v),lp\", "dp\" or \"Dp\" to the ordering
-         \"(a(w),lp\" or \"a(1,0,...,0),lp\" by using modular methods.
-         By default the procedure computes a standard basis of I for sure, but
-         if the optional parameter #[2] = 0, it computes a standard basis of I
-         with high probability.
-EXAMPLE: example modWalk; shows an example
-"
-{
-  int TT = timer;
-  int RT = rtimer;
-  int i,j,pTest,sizeTest,weighted,n1;
-  bigint N;
-
-  def R0 = basering;
-  list rl = ringlist(R0);
-  if((npars(R0) > 0) || (rl[1] > 0))
-  {
-    ERROR("Characteristic of basering should be zero, basering should have no parameters.");
-  }
-
-  if(typeof(II) == "ideal")
-  {
-    ideal I = II;
-    kill II;
-    list II;
-    II[1] = I;
-    II[2] = 2;
-    II[3] = 2;
-  }
-  else
-  {
-    if(typeof(II) == "list" && typeof(II[1]) == "ideal")
-    {
-      ideal I = II[1];
-      if(size(II) == 1)
-      {
-        II[2] = 2;
-        II[3] = 2;
-      }
-      if(size(II) == 2)
-      {
-        II[3] = 2;
-      }
-
-    }
-    else
-    {
-      ERROR("Unexpected type of input.");
-    }
-  }
-
-//--------------------  Initialize optional parameters  ------------------------
-  n1 = system("--cpus");
-  if(size(#) == 0)
-  {
-    int exactness = 1;
-    int n2 = 10;
-    int n3 = 10;
-  }
-  else
-  {
-    if(size(#) == 1)
-    {
-      if(typeof(#[1]) == "int")
-      {
-        if(#[1] < n1)
-        {
-          n1 = #[1];
-        }
-        int exactness = 1;
-        if(n1 >= 10)
-        {
-          int n2 = n1 + 1;
-          int n3 = n1;
-        }
-        else
-        {
-          int n2 = 10;
-          int n3 = 10;
-        }
-      }
-      else
-      {
-        ERROR("Unexpected type of input.");
-      }
-    }
-    if(size(#) == 2)
-    {
-      if(typeof(#[1]) == "int" && typeof(#[2]) == "int")
-      {
-        if(#[1] < n1)
-        {
-          n1 = #[1];
-        }
-        int exactness = #[2];
-        if(n1 >= 10)
-        {
-          int n2 = n1 + 1;
-          int n3 = n1;
-        }
-        else
-        {
-          int n2 = 10;
-          int n3 = 10;
-        }
-      }
-      else
-      {
-        if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec")
-        {
-          intvec curr_weight = #[1];
-          intvec target_weight = #[2];
-          weighted = 1;
-          int n2 = 10;
-          int n3 = 10;
-          int exactness = 1;
-        }
-        else
-        {
-          ERROR("Unexpected type of input.");
-        }
-      }
-    }
-    if(size(#) == 3)
-    {
-      if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int")
-      {
-        intvec curr_weight = #[1];
-        intvec target_weight = #[2];
-        weighted = 1;
-        n1 = #[3];
-        int n2 = 10;
-        int n3 = 10;
-        int exactness = 1;
-      }
-      else
-      {
-        ERROR("Unexpected type of input.");
-      }
-    }
-    if(size(#) == 4)
-    {
-      if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int"
-          && typeof(#[4]) == "int")
-      {
-        intvec curr_weight = #[1];
-        intvec target_weight = #[2];
-        weighted = 1;
-        if(#[1] < n1)
-        {
-          n1 = #[3];
-        }
-        int exactness = #[4];
-        if(n1 >= 10)
-        {
-          int n2 = n1 + 1;
-          int n3 = n1;
-        }
-        else
-        {
-          int n2 = 10;
-          int n3 = 10;
-        }
-      }
-      else
-      {
-        if(typeof(#[1]) == "int" && typeof(#[2]) == "int" && typeof(#[3]) == "int" && typeof(#[4]) == "int")
-        {
-          if(#[1] < n1)
-          {
-            n1 = #[1];
-          }
-          int exactness = #[2];
-          if(n1 >= #[3])
-          {
-            int n2 = n1 + 1;
-          }
-          else
-          {
-            int n2 = #[3];
-          }
-          if(n1 >= #[4])
-          {
-            int n3 = n1;
-          }
-          else
-          {
-            int n3 = #[4];
-          }
-        }
-        else
-        {
-          ERROR("Unexpected type of input.");
-        }
-      }
-    }
-    if(size(#) == 6)
-    {
-      if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int" && typeof(#[4]) == "int" && typeof(#[5]) == "int" && typeof(#[6]) == "int")
-      {
-        intvec curr_weight = #[1];
-        intvec target_weight = #[2];
-        weighted = 1;
-        if(#[3] < n1)
-        {
-          n1 = #[3];
-        }
-        int exactness = #[4];
-        if(n1 >= #[5])
-        {
-          int n2 = n1 + 1;
-        }
-        else
-        {
-          int n2 = #[5];
-        }
-        if(n1 >= #[6])
-        {
-          int n3 = n1;
-        }
-        else
-        {
-          int n3 = #[6];
-        }
-      }
-      else
-      {
-        ERROR("Expected list(intvec,intvec,int,int,int,int) as optional parameter list.");
-      }
-    }
-    if(size(#) == 1 || size(#) == 5 || size(#) > 6)
-    {
-      ERROR("Expected 0,2,3,4 or 5 optional arguments.");
-    }
-  }
-  if(printlevel >= 10)
-  {
-  "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)+", exactness = "+string(exactness);
-  }
-
-//-------------------------  Save current options  -----------------------------
-  intvec opt = option(get);
-  //option(redSB);
-
-//--------------------  Initialize the list of primes  -------------------------
-  int tt = timer;
-  int rt = rtimer;
-  int en = 2134567879;
-  int an = 1000000000;
-  intvec L = primeList(I,n2);
-  if(n2 > 4)
-  {
-  //  L[5] = prime(random(an,en));
-  }
-  if(printlevel >= 10)
-  {
-    "CPU-time for primeList: "+string(timer-tt)+" seconds.";
-    "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
-  }
-  int h = homog(I);
-  list P,T1,T2,LL,Arguments,PP;
-  ideal J,K,H;
-
-//-------------------  parallelized Groebner Walk in positive characteristic  --------------------
-
-  if(weighted)
-  {
-    for(i=1; i<=size(L); i++)
-    {
-      Arguments[i] = list(II,L[i],variant,list(curr_weight,target_weight));
-    }
-  }
-  else
-  {
-    for(i=1; i<=size(L); i++)
-    {
-      Arguments[i] = list(II,L[i],variant);
-    }
-  }
-  P = parallelWaitAll("modpWalk",Arguments);
-  for(i=1; i<=size(P); i++)
-  {
-    T1[i] = P[i][1];
-    T2[i] = bigint(P[i][2]);
-  }
-
-  while(1)
-  {
-    LL = deleteUnluckyPrimes(T1,T2,h);
-    T1 = LL[1];
-    T2 = LL[2];
-//-------------------  Now all leading ideals are the same  --------------------
-//-------------------  Lift results to basering via farey  ---------------------
-
-    tt = timer; rt = rtimer;
-    N = T2[1];
-    for(i=2; i<=size(T2); i++)
-    {
-      N = N*T2[i];
-    }
-    H = chinrem(T1,T2);
-    J = parallelFarey(H,N,n1);
-    //J=farey(H,N);
-    if(printlevel >= 10)
-    {
-      "CPU-time for lifting-process is "+string(timer - tt)+" seconds.";
-      "Real-time for lifting-process is "+string(rtimer - rt)+" seconds.";
-    }
-
-//----------------  Test if we already have a standard basis of I --------------
-
-    tt = timer; rt = rtimer;
-    pTest = pTestSB(I,J,L,variant);
-    //pTest = primeTestSB(I,J,L,variant);
-    if(printlevel >= 10)
-    {
-      "CPU-time for pTest is "+string(timer - tt)+" seconds.";
-      "Real-time for pTest is "+string(rtimer - rt)+" seconds.";
-    }
-    if(pTest)
-    {
-      if(printlevel >= 10)
-      {
-        "CPU-time for computation without final tests is "+string(timer - TT)+" seconds.";
-        "Real-time for computation without final tests is "+string(rtimer - RT)+" seconds.";
-      }
-      attrib(J,"isSB",1);
-      if(exactness == 0)
-      {
-        option(set, opt);
-        return(J);
-      }
-      else
-      {
-        tt = timer;
-        rt = rtimer;
-        sizeTest = 1 - isIdealIncluded(I,J,n1);
-        if(printlevel >= 10)
-        {
-          "CPU-time for checking if I subset <G> is "+string(timer - tt)+" seconds.";
-          "Real-time for checking if I subset <G> is "+string(rtimer - rt)+" seconds.";
-        }
-        if(sizeTest == 0)
-        {
-          tt = timer;
-          rt = rtimer;
-          K = std(J);
-          if(printlevel >= 10)
-          {
-            "CPU-time for last std-computation is "+string(timer - tt)+" seconds.";
-            "Real-time for last std-computation is "+string(rtimer - rt)+" seconds.";
-          }
-          if(size(reduce(K,J)) == 0)
-          {
-            option(set, opt);
-            return(J);
-          }
-        }
-      }
-    }
-//--------------  We do not already have a standard basis of I, therefore do the main computation for more primes  --------------
-
-    T1 = H;
-    T2 = N;
-    j = size(L)+1;
-    tt = timer; rt = rtimer;
-    L = primeList(I,n3,L,n1);
-    if(printlevel >= 10)
-    {
-      "CPU-time for primeList: "+string(timer-tt)+" seconds.";
-      "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
-    }
-    Arguments = list();
-    PP = list();
-    if(weighted)
-    {
-      for(i=j; i<=size(L); i++)
-      {
-        //Arguments[i-j+1] = list(II,L[i],variant,list(curr_weight,target_weight));
-        Arguments[size(Arguments)+1] = list(II,L[i],variant,list(curr_weight,target_weight));
-      }
-    }
-    else
-    {
-      for(i=j; i<=size(L); i++)
-      {
-        //Arguments[i-j+1] = list(II,L[i],variant);
-        Arguments[size(Arguments)+1] = list(II,L[i],variant);
-      }
-    }
-    PP = parallelWaitAll("modpWalk",Arguments);
-    if(printlevel >= 10)
-    {
-      "parallel modpWalk";
-    }
-    for(i=1; i<=size(PP); i++)
-    {
-      //P[size(P) + 1] = PP[i];
-      T1[size(T1) + 1] = PP[i][1];
-      T2[size(T2) + 1] = bigint(PP[i][2]);
-    }
-  }
-  if(printlevel >= 10)
-  {
-    "CPU-time for computation with final tests is "+string(timer - TT)+" seconds.";
-    "Real-time for computation with final tests is "+string(rtimer - RT)+" seconds.";
-  }
-}
-
-example
-{
-  "EXAMPLE:";
-  echo = 2;
-  ring R=0,(x,y,z),lp;
-  ideal I=-x+y2z-z,xz+1,x2+y2-1;
-  // I is a standard basis in dp
-  ideal J = modWalk(I,1);
-  J;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-static proc isIdealIncluded(ideal I, ideal J, int n1)
-"USAGE:  isIdealIncluded(I,J,int n1); I ideal, J ideal, n1 integer
-"
-{
-  if(n1 > 1)
-  {
-    int k;
-    list args,results;
-    for(k=1; k<=size(I); k++)
-    {
-      args[k] = list(ideal(I[k]),J,1);
-    }
-    results = parallelWaitAll("reduce",args);
-    for(k=1; k<=size(results); k++)
-    {
-      if(results[k] == 0)
-      {
-        return(1);
-      }
-    }
-    return(0);
-  }
-  else
-  {
-    if(reduce(I,J,1) == 0)
-    {
-      return(1);
-    }
-    else
-    {
-      return(0);
-    }
-  }
-}
-
-////////////////////////////////////////////////////////////////////////////////
-static proc parallelChinrem(list T1, list T2, int n1)
-"USAGE:  parallelChinrem(T1,T2); T1 list of ideals, T2 list of primes, n1 integer"
-{
-  int i,j,k;
-
-  ideal H,J;
-
-  list arguments_chinrem,results_chinrem;
-  for(i=1; i<=size(T1); i++)
-  {
-    J = ideal(T1[i]);
-    attrib(J,"isSB",1);
-    arguments_chinrem[size(arguments_chinrem)+1] = list(list(J),T2);
-  }
-  results_chinrem = parallelWaitAll("chinrem",arguments_chinrem);
-    for(j=1; j <= size(results_chinrem); j++)
-    {
-      J = results_chinrem[j];
-      attrib(J,"isSB",1);
-      if(isIdealIncluded(J,H,n1) == 0)
-      {
-        if(H == 0)
-        {
-          H = J;
-        }
-        else
-        {
-          H = H,J;
-        }
-      }
-    }
-  return(H);
-}
-
-////////////////////////////////////////////////////////////////////////////////
-static proc parallelFarey(ideal H, bigint N, int n1)
-"USAGE:  parallelFarey(H,N,n1); H ideal, N bigint, n1 integer
-"
-{
-  int i,j;
-  int ii = 1;
-  list arguments_farey,results_farey;
-  for(i=1; i<=size(H); i++)
-  {
-    for(j=1; j<=size(H[i]); j++)
-    {
-      arguments_farey[size(arguments_farey)+1] = list(H[i][j],N);
-    }
-  }
-  results_farey = parallelWaitAll("farey",arguments_farey);
-  ideal J,K;
-  poly f_farey;
-  while(ii<=size(results_farey))
-  {
-    for(i=1; i<=size(H); i++)
-    {
-      f_farey = 0;
-      for(j=1; j<=size(H[i]); j++)
-      {
-        f_farey = f_farey + results_farey[ii][1];
-        ii++;
-      }
-      K = ideal(f_farey);
-      attrib(K,"isSB",1);
-      attrib(J,"isSB",1);
-      if(isIdealIncluded(K,J,n1) == 0)
-      {
-        if(J == 0)
-        {
-          J = K;
-        }
-        else
-        {
-          J = J,K;
-        }
-      }
-    }
-  }
-  return(J);
-}
-//////////////////////////////////////////////////////////////////////////////////
-static proc primeTestSB(def II, ideal J, list L, int variant, list #)
-"USAGE:  primeTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int
-RETURN:  1 (resp. 0) if for a randomly chosen prime p that is not in L
-         J mod p is (resp. is not) a standard basis of I mod p
-EXAMPLE: example primeTestSB; shows an example
-"
-{
-if(typeof(II) == "ideal")
-  {
-  ideal I = II;
-  int radius = 2;
-  }
-if(typeof(II) == "list")
-  {
-  ideal I = II[1];
-  int radius = II[2];
-  }
-
-int i,j,k,p;
-def R = basering;
-list r = ringlist(R);
-
-while(!j)
-  {
-  j = 1;
-  p = prime(random(1000000000,2134567879));
-  for(i = 1; i <= size(L); i++)
-    {
-    if(p == L[i])
-      {
-      j = 0;
-      break;
-      }
-    }
-  if(j)
-    {
-    for(i = 1; i <= ncols(I); i++)
-      {
-      for(k = 2; k <= size(I[i]); k++)
-        {
-        if((denominator(leadcoef(I[i][k])) mod p) == 0)
-          {
-          j = 0;
-          break;
-          }
-        }
-      if(!j)
-        {
-        break;
-        }
-      }
-    }
-  if(j)
-    {
-    if(!primeTest(I,p))
-      {
-      j = 0;
-      }
-    }
-  }
-r[1] = p;
-def @R = ring(r);
-setring @R;
-ideal I = imap(R,I);
-ideal J = imap(R,J);
-attrib(J,"isSB",1);
-
-int t = timer;
-j = 1;
-if(isIncluded(I,J) == 0)
-  {
-  j = 0;
-  }
-if(printlevel >= 11)
-  {
-  "isIncluded(I,J) takes "+string(timer - t)+" seconds";
-  "j = "+string(j);
-  }
-t = timer;
-if(j)
-  {
-  if(size(#) > 0)
-    {
-    ideal K = modpWalk(I,p,variant,#)[1];
-    }
-  else
-    {
-    ideal K = modpWalk(I,p,variant)[1];
-    }
-  t = timer;
-  if(isIncluded(J,K) == 0)
-    {
-    j = 0;
-    }
-  if(printlevel >= 11)
-    {
-    "isIncluded(K,J) takes "+string(timer - t)+" seconds";
-    "j = "+string(j);
-    }
-  }
-setring R;
-
-return(j);
-}
-example
-{ "EXAMPLE:"; echo = 2;
-   intvec L = 2,3,5;
-   ring r = 0,(x,y,z),lp;
-   ideal I = x+1,x+y+1;
-   ideal J = x+1,y;
-   primeTestSB(I,I,L,1);
-   primeTestSB(I,J,L,1);
-}
-
-////////////////////////////////////////////////////////////////////////////////
-static proc pTestSB(ideal I, ideal J, list L, int variant, list #)
-"USAGE:  pTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int
-RETURN:  1 (resp. 0) if for a randomly chosen prime p that is not in L
-         J mod p is (resp. is not) a standard basis of I mod p
-EXAMPLE: example pTestSB; shows an example
-"
-{
-   int i,j,k,p;
-   def R = basering;
-   list r = ringlist(R);
-
-   while(!j)
-   {
-      j = 1;
-      p = prime(random(1000000000,2134567879));
-      for(i = 1; i <= size(L); i++)
-      {
-         if(p == L[i]) { j = 0; break; }
-      }
-      if(j)
-      {
-         for(i = 1; i <= ncols(I); i++)
-         {
-            for(k = 2; k <= size(I[i]); k++)
-            {
-               if((denominator(leadcoef(I[i][k])) mod p) == 0) { j = 0; break; }
-            }
-            if(!j){ break; }
-         }
-      }
-      if(j)
-      {
-         if(!primeTest(I,p)) { j = 0; }
-      }
-   }
-   r[1] = p;
-   def @R = ring(r);
-   setring @R;
-   ideal I = imap(R,I);
-   ideal J = imap(R,J);
-   attrib(J,"isSB",1);
-
-   int t = timer;
-   j = 1;
-   if(isIncluded(I,J) == 0) { j = 0; }
-
-   if(printlevel >= 11)
-   {
-      "isIncluded(I,J) takes "+string(timer - t)+" seconds";
-      "j = "+string(j);
-   }
-
-   t = timer;
-   if(j)
-   {
-      if(size(#) > 0)
-      {
-         ideal K = modpStd(I,p,variant,#[1])[1];
-      }
-      else
-      {
-         ideal K = groebner(I);
-      }
-      t = timer;
-      if(isIncluded(J,K) == 0) { j = 0; }
-
-      if(printlevel >= 11)
-      {
-         "isIncluded(J,K) takes "+string(timer - t)+" seconds";
-         "j = "+string(j);
-      }
-   }
-   setring R;
-   return(j);
-}
-example
-{ "EXAMPLE:"; echo = 2;
-   intvec L = 2,3,5;
-   ring r = 0,(x,y,z),dp;
-   ideal I = x+1,x+y+1;
-   ideal J = x+1,y;
-   pTestSB(I,I,L,2);
-   pTestSB(I,J,L,2);
-}
-////////////////////////////////////////////////////////////////////////////////
-static proc mixedTest()
-"USAGE:  mixedTest();
-RETURN:  1 if ordering of basering is mixed, 0 else
-EXAMPLE: example mixedTest(); shows an example
-"
-{
-   int i,p,m;
-   for(i = 1; i <= nvars(basering); i++)
-   {
-      if(var(i) > 1)
-      {
-         p++;
-      }
-      else
-      {
-         m++;
-      }
-   }
-   if((p > 0) && (m > 0)) { return(1); }
-   return(0);
-}
-example
-{ "EXAMPLE:"; echo = 2;
-   ring R1 = 0,(x,y,z),dp;
-   mixedTest();
-   ring R2 = 31,(x(1..4),y(1..3)),(ds(4),lp(3));
-   mixedTest();
-   ring R3 = 181,x(1..9),(dp(5),lp(4));
-   mixedTest();
-}
+    int pertdeg = 3;
+    ideal J = modrWalk(I,radius,pertdeg);
+    J;
+}
+
+proc modfWalk(ideal I, list #)
+"USAGE:   modfWalk(I, [, v, w]); I ideal, v intvec, w intvec
+RETURN:   a standard basis of I
+NOTE:     The procedure computes a standard basis of I (over the rational
+          numbers) by using modular methods.
+SEE ALSO: modular
+EXAMPLE:  example modfWalk; shows an example"
+{
+    /* read optional parameter */
+    if (size(#) > 0) {
+        if (size(#) == 1) {
+            intvec w = #[1];
+        }
+        if (size(#) == 2) {
+            intvec v = #[1];
+            intvec w = #[2];
+        }
+        if (size(#) > 2 || typeof(#[1]) != "intvec") {
+            ERROR("wrong optional parameter");
+        }
+    }
+
+    /* save options */
+    intvec opt = option(get);
+    option(redSB);
+
+    /* set additional parameters */
+    int reduction = 1;
+    int printout = 0;
+
+    /* call modular() */
+    if (size(#) > 0) {
+        I = modular("fwalk", list(I,reduction,printout,#));
+    }
+    else {
+        I = modular("fwalk", list(I,reduction,printout));
+    }
+
+    /* return the result */
+    attrib(I, "isSB", 1);
+    option(set, opt);
+    return(I);
+}
+example
+{
+    "EXAMPLE:";
+    echo = 2;
+    ring R1 = 0, (x,y,z,t), dp;
+    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
+    I = std(I);
+    ring R2 = 0, (x,y,z,t), lp;
+    ideal I = fetch(R1, I);
+    ideal J = modfWalk(I);
+    J;
+}
+
+proc modfrWalk(ideal I, int radius, list #)
+"USAGE:   modfrWalk(I, radius [, v, w]); I ideal, radius int, v intvec, w intvec
+RETURN:   a standard basis of I
+NOTE:     The procedure computes a standard basis of I (over the rational
+          numbers) by using modular methods.
+SEE ALSO: modular
+EXAMPLE:  example modfrWalk; shows an example"
+{
+    /* read optional parameter */
+    if (size(#) > 0) {
+        if (size(#) == 1) {
+            intvec w = #[1];
+        }
+        if (size(#) == 2) {
+            intvec v = #[1];
+            intvec w = #[2];
+        }
+        if (size(#) > 2 || typeof(#[1]) != "intvec") {
+            ERROR("wrong optional parameter");
+        }
+    }
+
+    /* save options */
+    intvec opt = option(get);
+    option(redSB);
+
+    /* set additional parameters */
+    int reduction = 1;
+    int printout = 0;
+
+    /* call modular() */
+    if (size(#) > 0) {
+        I = modular("frandwalk", list(I,radius,reduction,printout,#));
+    }
+    else {
+        I = modular("frandwalk", list(I,radius,reduction,printout));
+    }
+
+    /* return the result */
+    attrib(I, "isSB", 1);
+    option(set, opt);
+    return(I);
+}
+example
+{
+    "EXAMPLE:";
+    echo = 2;
+    ring R1 = 0, (x,y,z,t), dp;
+    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
+    I = std(I);
+    ring R2 = 0, (x,y,z,t), lp;
+    ideal I = fetch(R1, I);
+    int radius = 2;
+    ideal J = modfrWalk(I,radius);
+    J;
+}

Singular/LIB/ncfactor.lib

@@ -6498 +6498 @@
     def Firstweyl = nc_algebra(1,1);
     setring Firstweyl;
-    ideal M = 0:nvars(@r);
+    ideal M = ideal(0:nvars(@r));
     M[varnumX]=var(2);
     M[varnumD]=var(1);

Singular/LIB/nctools.lib

@@ -1597 +1597 @@
   matrix U = UpOneMatrix(N);
   if (size(ideal(C-U)) <> 0) { return(notW); } // lt(xy)<>lt(yx)
-  ideal I = D;
+  ideal I = ideal(D);
   if (size(I) <> n) { return(notW); } // not n entries<>0
   I = simplify(I,4+2);

Singular/LIB/nfmodstd.lib

@@ -524 +524 @@
     I;
 }
+
+//////////////////////////////////////////////////////////////////////////////
+
+/*
+Benchmark Problems from
+
+Boku, Decker, Fieker, Steenpass: Groebner Bases over Algebraic Number
+Fields.
+
+// 1
+ring R = (0,a), (x,y,z), dp;
+minpoly = (a^2+1);
+poly f1 = (a+8)*x^2*y^2+5*x*y^3+(-a+3)*x^3*z
+          +x^2*y*z;
+poly f2 = x^5+2*y^3*z^2+13*y^2*z^3+5*y*z^4;
+poly f3 = 8*x^3+(a+12)*y^3+x*z^2+3;
+poly f4 = (-a+7)*x^2*y^4+y^3*z^3+18*y^3*z^2;
+ideal I1 = f1,f2,f3,f4;
+
+// 2
+ring R = (0,a), (x,y,z), dp;
+minpoly = (a^5+a^2+2);
+poly f1 = 2*x*y^4*z^2+(a-1)*x^2*y^3*z
+          +(2*a)*x*y*z^2+7*y^3+(7*a+1);
+poly f2 = 2*x^2*y^4*z+(a)*x^2*y*z^2-x*y^2*z^2
+          +(2*a+3)*x^2*y*z-12*x+(12*a)*y;
+poly f3 = (2*a)*y^5*z+x^2*y^2*z-x*y^3*z
+          +(-a)*x*y^3+y^4+2*y^2*z;
+poly f4 = (3*a)*x*y^4*z^3+(a+1)*x^2*y^2*z
+          -x*y^3*z+4*y^3*z^2+(3*a)*x*y*z^3
+          +4*z^2-x+(a)*y;
+ideal I2 = f1,f2,f3,f4;
+
+// 3a
+ring R = (0,a), (v,w,x,y,z), dp;
+minpoly = (a^7-7*a+3);
+poly f1 = (a)*v+(a-1)*w+x+(a+2)*y+z;
+poly f2 = v*w+(a-1)*w*x+(a+2)*v*y+x*y+(a)*y*z;
+poly f3 = (a)*v*w*x+(a+5)*w*x*y+(a)*v*w*z
+          +(a+2)*v*y*z+(a)*x*y*z;
+poly f4 = (a-11)*v*w*x*y+(a+5)*v*w*x*z
+          +(a)*v*w*y*z+(a)*v*x*y*z
+          +(a)*w*x*y*z;
+poly f5 = (a+3)*v*w*x*y*z+(a+23);
+ideal I3a = f1,f2,f3,f4,f5;
+
+// 3b
+ring R = (0,a), (u,v,w,x,y,z), dp;
+minpoly = (a^7-7*a+3);
+poly f1 = (a)*u+(a+2)*v+w+x+y+z;
+poly f2 = u*v+v*w+w*x+x*y+(a+3)*u*z+y*z;
+poly f3 = u*v*w+v*w*x+(a+1)*w*x*y+u*v*z+u*y*z
+          +x*y*z;
+poly f4 = (a-1)*u*v*w*x+v*w*x*y+u*v*w*z
+          +u*v*y*z+u*x*y*z+w*x*y*z;
+poly f5 = u*v*w*x*y+(a+1)*u*v*w*x*z+u*v*w*y*z
+          +u*v*x*y*z+u*w*x*y*z+v*w*x*y*z;
+poly f6 = u*v*w*x*y*z+(-a+2);
+ideal I3b = f1,f2,f3,f4,f5,f6;
+
+// 4
+ring R = (0,a), (w,x,y,z), dp;
+minpoly = (a^6+a^5+a^4+a^3+a^2+a+1);
+poly f1 = (a+5)*w^3*x^2*y+(a-3)*w^2*x^3*y
+          +(a+7)*w*x^2*y^2;
+poly f2 = (a)*w^5+(a+3)*w*x^2*y^2
+          +(a^2+11)*x^2*y^2*z;
+poly f3 = (a+7)*w^3+12*x^3+4*w*x*y+(a)*z^3;
+poly f4 = 3*w^3+(a-4)*x^3+x*y^2;
+ideal I4 = f1,f2,f3,f4;
+
+// 5
+ring R = (0,a), (w,x,y,z), dp;
+minpoly = (a^12-5*a^11+24*a^10-115*a^9+551*a^8
+          -2640*a^7+12649*a^6-2640*a^5+551*a^4
+          -115*a^3+24*a^2-5*a+1);
+poly f1 = (2*a+3)*w*x^4*y^2+(a+1)*w^2*x^3*y*z
+          +2*w*x*y^2*z^3+(7*a-1)*x^3*z^4;
+poly f2 = 2*w^2*x^4*y+w^2*x*y^2*z^2
+          +(-a)*w*x^2*y^2*z^2
+          +(a+11)*w^2*x*y*z^3-12*w*z^6
+          +12*x*z^6;
+poly f3 = 2*x^5*y+w^2*x^2*y*z-w*x^3*y*z
+          -w*x^3*z^2+(a)*x^4*z^2+2*x^2*y*z^3;
+poly f4 = 3*w*x^4*y^3+w^2*x^2*y*z^3
+          -w*x^3*y*z^3+(a+4)*x^3*y^2*z^3
+          +3*w*x*y^3*z^3+(4*a)*y^2*z^6-w*z^7
+          +x*z^7;
+ideal I5 = f1,f2,f3,f4;
+
+// 6
+ring R = (0,a), (u,v,w,x,y,z), dp;
+minpoly = (a^2+5*a+1);
+poly f1 = u+v+w+x+y+z+(a);
+poly f2 = u*v+v*w+w*x+x*y+y*z+(a)*u+(a)*z;
+poly f3 = u*v*w+v*w*x+w*x*y+x*y*z+(a)*u*v
+          +(a)*u*z+(a)*y*z;
+poly f4 = u*v*w*x+v*w*x*y+w*x*y*z+(a)*u*v*w
+          +(a)*u*v*z+(a)*u*y*z+(a)*x*y*z;
+poly f5 = u*v*w*x*y+v*w*x*y*z+(a)*u*v*w*x
+          +(a)*u*v*w*z+(a)*u*v*y*z+(a)*u*x*y*z
+          +(a)*w*x*y*z;
+poly f6 = u*v*w*x*y*z+(a)*u*v*w*x*y
+          +(a)*u*v*w*x*z+(a)*u*v*w*y*z
+          +(a)*u*v*x*y*z+(a)*u*w*x*y*z
+          +(a)*v*w*x*y*z;
+poly f7 = (a)*u*v*w*x*y*z-1;
+ideal I6 = f1,f2,f3,f4,f5,f6,f7;
+
+// 7
+ring R = (0,a), (w,x,y,z), dp;
+minpoly = (a^8-16*a^7+19*a^6-a^5-5*a^4+13*a^3
+          -9*a^2+13*a+17);
+poly f1 = (-a^2-1)*x^2*y+2*w*x*z-2*w
+          +(a^2+1)*y;
+poly f2 = (a^3-a-3)*w^3*y+4*w*x^2*y+4*w^2*x*z
+          +2*x^3*z+(a)*w^2-10*x^2+4*w*y-10*x*z
+          +(2*a^2+a);
+poly f3 = (a^2+a+11)*x*y*z+w*z^2-w-2*y;
+poly f4 = -w*y^3+4*x*y^2*z+4*w*y*z^2+2*x*z^3
+          +(2*a^3+a^2)*w*y+4*y^2-10*x*z-10*z^2
+          +(3*a^2+5);
+ideal I7 = f1,f2,f3,f4;
+
+// 8
+ring R = (0,a), (t,u,v,w,x,y,z), dp;
+minpoly = (a^7+10*a^5+5*a^3+10*a+1);
+poly f1 = v*x+w*y-x*z-w-y;
+poly f2 = v*w-u*x+x*y-w*z+v+x+z;
+poly f3 = t*w-w^2+x^2-t;
+poly f4 = (-a)*v^2-u*y+y^2-v*z-z^2+u;
+poly f5 = t*v+v*w+(-a^2-a-5)*x*y-t*z+w*z+v+x+z
+          +(a+1);
+poly f6 = t*u+u*w+(-a-11)*v*x-t*y+w*y-x*z-t-u
+          +w+y;
+poly f7 = w^2*y^3-w*x*y^3+x^2*y^3+w^2*y^2*z
+          -w*x*y^2*z+x^2*y^2*z+w^2*y*z^2
+          -w*x*y*z^2+x^2*y*z^2+w^2*z^3-w*x*z^3
+          +x^2*z^3;
+poly f8 = t^2*u^3+t^2*u^2*v+t^2*u*v^2+t^2*v^3
+          -t*u^3*x-t*u^2*v*x-t*u*v^2*x-t*v^3*x
+          +u^3*x^2+u^2*v*x^2+u*v^2*x^2
+          +v^3*x^2;
+ideal I8 = f1,f2,f3,f4,f5,f6,f7,f8;
+*/

Singular/LIB/primdec.lib

@@ -874 +874 @@
   int uytrewq;
   int nva = nvars(basering);
-  int @k,@s,@n,@k1,zz;
+  int @k,@s,@n,@k1,@zz;
   list primary,lres0,lres1,act,@lh,@wh;
   map phi,psi,phi1,psi1;
@@ -1055 +1055 @@
   {
     poly randp;
-    for(zz=1;zz<nvars(basering);zz++)
+    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);
+              +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(@zz);
     }
     randp=randp+var(nvars(basering));
@@ -1068 +1068 @@
     if (size(primary[2*@k])==0)
     {
-      for(zz=1;zz<size(primary[2*@k-1])-1;zz++)
+      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]))
+        if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][@zz]))
         {
           primary[2*@k]=primary[2*@k-1];
@@ -1099 +1099 @@
         if(deg(lead(primary[2*@k-1][@n]))==1)
         {
-          for(zz=1;zz<=nva;zz++)
+          for(@zz=1;@zz<=nva;@zz++)
           {
-            if(lead(primary[2*@k-1][@n])/var(zz)!=0)
+            if(lead(primary[2*@k-1][@n])/var(@zz)!=0)
             {
-              jmap1[zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]
+              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);
+              jmap2[@zz]=primary[2*@k-1][@n];
+              @qht[@n]=var(@zz);
             }
           }
@@ -3413 +3413 @@
 
   ideal jkeep=@j;
-  if(ordstr(@P)[1]=="w")
+  if((ordstr(@P)[1]=="w")&&(size(ringlist(@P)[3])==2)) // weighted ordering
   {
     list gnir_l=ringlist(gnir);
@@ -9023 +9023 @@
   int uytrewq;
   int nva = nvars(basering);
-  int @k,@s,@n,@k1,zz;
+  int @k,@s,@n,@k1,@zz;
   list primary,lres0,lres1,act,@lh,@wh;
   map phi,psi,phi1,psi1;
@@ -9206 +9206 @@
   {
     poly randp;
-    for(zz=1;zz<nvars(basering);zz++)
+    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);
+              +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(@zz);
     }
     randp=randp+var(nvars(basering));
@@ -9219 +9219 @@
     if (size(primary[2*@k])==0)
     {
-      for(zz=1;zz<size(primary[2*@k-1])-1;zz++)
+      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]))
+        if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][@zz]))
         {
           primary[2*@k]=primary[2*@k-1];
@@ -9250 +9250 @@
         if(deg(lead(primary[2*@k-1][@n]))==1)
         {
-          for(zz=1;zz<=nva;zz++)
+          for(@zz=1;@zz<=nva;@zz++)
           {
-            if(lead(primary[2*@k-1][@n])/var(zz)!=0)
+            if(lead(primary[2*@k-1][@n])/var(@zz)!=0)
             {
-              jmap1[zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]
+              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);
+              jmap2[@zz]=primary[2*@k-1][@n];
+              @qht[@n]=var(@zz);
             }
           }

Singular/LIB/rwalk.lib

@@ -10 +10 @@
 rwalk(ideal,int,int[,intvec,intvec]);   standard basis of ideal via Random Walk algorithm
 rwalk(ideal,int[,intvec,intvec]);       standard basis of ideal via Random Perturbation Walk algorithm
-frwalk(ideal,int[,intvec,intvec]);      standard basis of ideal via Random Fractal Walk algorithm
+frandwalk(ideal,int[,intvec,intvec]);      standard basis of ideal via Random Fractal Walk algorithm
 ";
 
@@ -141 +141 @@
  * Random Walk  *
  ****************/
-proc rwalk(ideal Go, int radius, int pert_deg, list #)
+proc rwalk(ideal Go, int radius, int pert_deg, int reduction, int printout, list #)
 "SYNTAX: rwalk(ideal i, int radius);
          if size(#)>0 then rwalk(ideal i, int radius, intvec v, intvec w);
+         intermediate Groebner bases are not reduced if reduction = 0
 TYPE:    ideal
 PURPOSE: compute the standard basis of the ideal, calculated via
@@ -178 +179 @@
 
 ideal G = fetch(xR, Go);
-G = system("Mrwalk", G, curr_weight, target_weight, radius, pert_deg, basering);
+G = system("Mrwalk", G, curr_weight, target_weight, radius, pert_deg, reduction, printout);
 
 setring xR;
@@ -196 +197 @@
   int radius = 1;
   int perturb_deg = 2;
-  rwalk(I,radius,perturb_deg);
+  int reduction = 0;
+  int printout = 1;
+  rwalk(I,radius,perturb_deg,reduction,printout);
 }
 
@@ -202 +205 @@
  * Perturbation Walk with random element *
  *****************************************/
-proc prwalk(ideal Go, int radius, int o_pert_deg, int t_pert_deg, list #)
+proc prwalk(ideal Go, int radius, int o_pert_deg, int t_pert_deg, int reduction, int printout, list #)
 "SYNTAX: rwalk(ideal i, int radius);
          if size(#)>0 then rwalk(ideal i, int radius, intvec v, intvec w);
@@ -227 +230 @@
   OSCTW= OrderStringalp_NP("al", #);
   }
+int nP = OSCTW[1];
 string ord_str = OSCTW[2];
 intvec curr_weight = OSCTW[3]; // original weight vector
@@ -238 +242 @@
 
 ideal G = fetch(xR, Go);
-G = system("Mprwalk", G, curr_weight, target_weight, radius, o_pert_deg, t_pert_deg, basering);
+G = system("Mprwalk", G, curr_weight, target_weight, radius, o_pert_deg, t_pert_deg,
+           nP, reduction, printout);
 
 setring xR;
@@ -257 +262 @@
   int o_perturb_deg = 2;
   int t_perturb_deg = 2;
-  prwalk(I,radius,o_perturb_deg,t_perturb_deg);
+  int reduction = 0;
+  int printout = 2;
+  prwalk(I,radius,o_perturb_deg,t_perturb_deg,reduction,printout);
 }
 
@@ -263 +270 @@
  * Fractal Walk with random element *
  ************************************/
-proc frandwalk(ideal Go, int radius, list #)
-"SYNTAX: frwalk(ideal i, int radius);
-         frwalk(ideal i, int radius, intvec v, intvec w);
+proc frandwalk(ideal Go, int radius, int reduction, int printout, list #)
+"SYNTAX: frwalk(ideal i, int radius, int reduction, int printout);
+         frwalk(ideal i, int radius, int reduction, int printout, intvec v, intvec w);
 TYPE:    ideal
 PURPOSE: compute the standard basis of the ideal, calculated via
@@ -299 +306 @@
    ideal G = fetch(xR, Go);
    int pert_deg = 2;
-   G = system("Mfrwalk", G, curr_weight, target_weight, radius);
+
+   G = system("Mfrwalk", G, curr_weight, target_weight, radius, reduction, printout);
 
    setring xR;
@@ -314 +322 @@
     ring r = 0,(z,y,x), lp;
     ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
-    frandwalk(I,2);
-}
+    int reduction = 0;
+    frandwalk(I,2,0,1);
+}

Singular/LIB/surf.lib

@@ -274 +274 @@
 
   string surf_call; i = 0;
-  
+
   if (isWindows())
   {
@@ -303 +303 @@
   else
   {
-    surf_call = "surfer";    
+    surf_call = "surfer";
     surf_call = surf_call + " " + l + " >/dev/null 2>&1";
     "Close window to exit from `surfer`.";
     i = system("sh", surf_call);
-    
-    if ( (i != 0) && isMacOSX() ) 
+
+    if ( (i != 0) && isMacOSX() )
     {
       "*!* Sorry: calling `surfer` failed ['"+surf_call+"']" + newline
       + " (The shell returned the error code " + string(i) + "." + newline
-      + "But since we are on Mac OS X, let us try to open SURFER.app instead..." + newline 
+      + "But since we are on Mac OS X, let us try to open SURFER.app instead..." + newline
       + "Appropriate SURFER.app is available for instance at http://www.mathematik.uni-kl.de/~motsak/files/SURFER.dmg";
-      
+
       // fallback, will only work if SURFER.app is available (e.g. in /Applications)
       // get SURFER.app e.g. from http://www.mathematik.uni-kl.de/~motsak/files/SURFER.dmg
       // note that the newer (Java-based) variant of Surfer may not support command line usage yet :(
-      
+
       surf_call = "open -a SURFER -W --args -t -s";
       surf_call = surf_call + " " + l + " >/dev/null 2>&1";
@@ -324 +324 @@
       i = system("sh", surf_call);
     }
-    
-    
-    
+
+
+
   }
   system("sh", "/bin/rm " + l);
@@ -381 +381 @@
 {
   string s = system("uname");
-  
+
   for (int i = 1; i <= size(s)-2; i = i + 1)
   {

Singular/Makefile.am

@@ -98 +98 @@
    mmalloc.h \
    mod_lib.h \
-   omSingularConfig.h \
    links/ndbm.h \
    newstruct.h \

Singular/blackbox.cc

@@ -54 +54 @@
 BOOLEAN blackbox_default_serialize(blackbox */*b*/, void */*d*/, si_link /*f*/)
 {
+  WerrorS("blackbox_serialize is not implemented");
   return TRUE;
 }
@@ -59 +60 @@
 BOOLEAN blackbox_default_deserialize(blackbox **/*b*/, void **/*d*/, si_link /*f*/)
 {
+  WerrorS("blackbox_deserialize is not implemented");
   return TRUE;
 }

Singular/cntrlc.cc

@@ -5 +5 @@
 * ABSTRACT - interupt handling
 */
-#include <kernel/mod2.h>
-
-/* includes */
-#ifdef DecAlpha_OSF1
-#define _XOPEN_SOURCE_EXTENDED
-#endif /* MP3-Y2 0.022UF */
-
-#include <omalloc/omalloc.h>
-
-#include <reporter/si_signals.h>
-#include <Singular/fevoices.h>
-
-#include <Singular/tok.h>
-#include <Singular/ipshell.h>
-void sig_chld_hdl(int sig); /*#include <Singular/links/ssiLink.h>*/
-#include <Singular/cntrlc.h>
-#include <Singular/feOpt.h>
-#include <Singular/misc_ip.h>
-#include <Singular/links/silink.h>
-#include <Singular/links/ssiLink.h>
-
 #include <stdio.h>
 #include <stddef.h>
@@ -33 +12 @@
 #include <sys/types.h>
 #include <sys/wait.h>
+
+#include <kernel/mod2.h>
+
+#include <omalloc/omalloc.h>
+
+#include <reporter/si_signals.h>
+#include <Singular/fevoices.h>
+
+#include <Singular/tok.h>
+#include <Singular/ipshell.h>
+#include <Singular/cntrlc.h>
+#include <Singular/feOpt.h>
+#include <Singular/misc_ip.h>
+#include <Singular/links/silink.h>
+#include <Singular/links/ssiLink.h>
 
 /* undef, if you don't want GDB to come up on error */
@@ -352 +346 @@
       if (feOptValue(FE_OPT_EMACS) == NULL)
       {
-        fputs("abort after this command(a), abort immediately(r), print backtrace(b), continue(c) or quit Singular(q) ?",stderr);fflush(stderr);
+        fputs("abort after this command(a), abort immediately(r), print backtrace(b), continue(c) or quit Singular(q) ?",stderr);
+        fflush(stderr);fflush(stdin);
         c = fgetc(stdin);
       }
@@ -371 +366 @@
                   fputs("** Warning: Singular should be restarted as soon as possible **\n",stderr);
                   fflush(stderr);
+                  extern void my_yy_flush();
+                  my_yy_flush();
+                  currentVoice=feInitStdin(NULL);
                   longjmp(si_start_jmpbuf,1);
                 }

Singular/dyn_modules/gfanlib/bbcone.cc

@@ -1739 +1739 @@
 }
 
-
+#if 0
+BOOLEAN bbcone_serialize(blackbox *b, void *d, si_link f)
+{
+  ssiInfo *dd = (ssiInfo *)f->data;
+  sleftv l;
+  memset(&l,0,sizeof(l));
+  l.rtyp=STRING_CMD;
+  l.data=(void*)"cone";
+  f->m->Write(f, &l);
+  gfan::ZCone *Z=((gfan::ZCone*) d;
+  /* AMBIENT_DIM */ fprintf(dd->f_write("%d ",Z->ambientDimension());
+  /* FACETS or INEQUALITIES */ fprintf(dd->f_write("%d ",Z->areFacetsKnown());
+  gfan::ZMatrix i=Z->getInequalities();
+....
+  /* LINEAR_SPAN or EQUATIONS */ fprintf(dd->f_write("%d ",Z->areImpliedEquationsKnown());
+  gfan::ZMatrix e=Z->getEquations();
+....
+  /* RAYS */
+  gfan::ZMatrix r=Z->extremeRays();
+....
+  /* LINEALITY_SPACE */
+  gfan::ZMatrix l=Z->generatorsOfLinealitySpace();
+....
+  return FALSE;
+}
+BOOLEAN bbcone_deserialize(blackbox **b, void **d, si_link f)
+{
+  ssiInfo *dd = (ssiInfo *)f->data;
+  gfan::ZCone *Z;
+  /* AMBIENT_DIM */ = s_readint(dd->f_read);
+  /* areFacetsKnown: */ = s_readint(dd->f_read);
+  if (areFacetsKnown)
+  ....FACETS
+  else
+  ....INEQUALITIES
+  /* areImpliedEquationsKnown*/ = s_readint(dd->f_read);
+  if(areImpliedEquationsKnown)
+  ....EQUATIONS
+  else
+  ...LINEAR_SPAN
+  ...RAYS
+  ...LINEALITY_SPACE
+  *d=Z;
+  return FALSE;
+}
+#endif
 void bbcone_setup(SModulFunctions* p)
 {
@@ -1753 +1798 @@
   b->blackbox_Assign=bbcone_Assign;
   b->blackbox_Op2=bbcone_Op2;
+  //b->blackbox_serialize=bbcone_serialize;
+  //b->blackbox_deserialize=bbcone_deserialize;
   p->iiAddCproc("","coneViaInequalities",FALSE,coneViaNormals);
   p->iiAddCproc("","coneViaPoints",FALSE,coneViaRays);

Singular/dyn_modules/singmathic/singmathic.cc

@@ -468 +468 @@
   result->rtyp=NONE;
 
-  if (arg == NULL || arg->next != NULL || 
+  if (arg == NULL || arg->next != NULL ||
   ((arg->Typ() != IDEAL_CMD) &&(arg->Typ() != MODUL_CMD))){
     WerrorS("Syntax: mathicgb(<ideal>/<module>)");

Singular/dyn_modules/syzextra/mod_main.cc

@@ -82 +82 @@
     for (int j=0; j<l; j++)
       if (id->m[j] != NULL && p_GetComp(id->m[j], r) > 0)
-        return TRUE;    
+        return TRUE;
 
     return FALSE; // rank: 1, only zero or no entries? can be an ideal OR module... BUT in the use-case should better be an ideal!
@@ -91 +91 @@
 
 
-   
+
 
 static inline void NoReturn(leftv& res)
@@ -1550 +1550 @@
   const int iLimit = r->typ[pos].data.is.limit;
   const ideal F = r->typ[pos].data.is.F;
-  
+
   ideal FF = id_Copy(F, r);
 

Singular/extra.cc

@@ -3794 +3794 @@
     }
     else
-#endif    
+#endif
 /*==================== Error =================*/
       Werror( "(extended) system(\"%s\",...) %s", sys_cmd, feNotImplemented );

Singular/fevoices.cc

@@ -17 +17 @@
 #include <Singular/subexpr.h>
 #include <Singular/ipshell.h>
+#include <Singular/sdb.h>
 
 #include <stdlib.h>

Singular/ipassign.cc

@@ -844 +844 @@
   if (res->data!=NULL) idDelete((ideal*)&res->data);
   matrix m=(matrix)a->CopyD(MATRIX_CMD);
+  if (TEST_V_ALLWARN)
+    if (MATROWS(m)>1)
+      Warn("assign matrix with %d rows to an ideal in >>%s<<",MATROWS(m),my_yylinebuf);
   IDELEMS((ideal)m)=MATROWS(m)*MATCOLS(m);
   ((ideal)m)->rank=1;

Singular/ipshell.cc

@@ -65 +65 @@
 #include <Singular/subexpr.h>
 #include <Singular/fevoices.h>
+#include <Singular/sdb.h>
 
 #include <math.h>
@@ -1029 +1030 @@
 void iiDebug()
 {
+#ifdef HAVE_SDB
+  sdb_flags=1;
+#endif
   Print("\n-- break point in %s --\n",VoiceName());
   if (iiDebugMarker) VoiceBackTrack();

Singular/misc_ip.cc

@@ -716 +716 @@
   } while (v!=NULL);
 
-#ifdef OM_SINGULAR_CONFIG_H
    // set global variable to show memory usage
   extern int om_sing_opt_show_mem;
   if (BVERBOSE(V_SHOW_MEM)) om_sing_opt_show_mem = 1;
   else om_sing_opt_show_mem = 0;
-#endif
 
   return FALSE;
@@ -1095 +1093 @@
 #endif   // HAVE_SIMPLEIPC
     fe_reset_input_mode();
+    monitor(NULL,0);
 #ifdef PAGE_TEST
     mmEndStat();

Singular/scanner.cc

@@ -2330 +2330 @@
   //yy_flush_buffer((YY_BUFFER_STATE)oldb);
 }
+
+void my_yy_flush() { YY_FLUSH_BUFFER;BEGIN(0); }
+

Singular/scanner.ll

@@ -380 +380 @@
   //yy_flush_buffer((YY_BUFFER_STATE)oldb);
 }
+
+void my_yy_flush() { YY_FLUSH_BUFFER;BEGIN(0); }

Singular/sdb.cc

@@ -244 +244 @@
           "p <var> - display type and value of the variable <var>\n"
           "q <flags> - quit debugger, set debugger flags(0,1,2)\n"
+          "   0: stop debug, 1:continue, 2: throw an error, return to toplevel\n"
           "Q - quit Singular\n");
           int i;

Singular/test.cc

@@ -194 +194 @@
 #include <Singular/links/ndbm.h>
 #include <Singular/newstruct.h>
-#include <Singular/omSingularConfig.h>
 #include <Singular/pcv.h>
 #include <Singular/links/pipeLink.h>

Singular/tesths.cc

@@ -146 +146 @@
   {
     if (feOptValue(FE_OPT_SORT)) On(SW_USE_NTL_SORT);
-#ifdef HAVE_SDB
-    sdb_flags = 0;
-#endif
     dup2(1,2);
     /* alternative:

Singular/walk.cc

@@ -16 +16 @@
 
 //#define TEST_OVERFLOW
-//#define CHECK_IDEAL
-//#define CHECK_IDEAL_MWALK
+
+#define CHECK_IDEAL_MWALK //to print intermediate results
 
 //#define NEXT_VECTORS_CC
-//#define PRINT_VECTORS //to print vectors (sigma, tau, omega)
+#define PRINT_VECTORS //to print weight vectors
 
 #define INVEPS_SMALL_IN_FRACTAL  //to choose the small invers of epsilon
@@ -27 +27 @@
 
 #define FIRST_STEP_FRACTAL // to define the first step of the fractal
-//#define MSTDCC_FRACTAL // apply Buchberger alg to compute a red GB, if
-//                          tau doesn't stay in the correct cone
+#define MSTDCC_FRACTAL // apply Buchberger alg to compute a red GB, if tau doesn't stay in the correct cone
 
 //#define TIME_TEST // print the used time of each subroutine
@@ -42 +41 @@
 #include <Singular/ipshell.h>
 #include <Singular/ipconv.h>
+#include <coeffs/ffields.h>
 #include <coeffs/coeffs.h>
 #include <Singular/subexpr.h>
+#include <polys/templates/p_Procs.h>
 
 #include <polys/monomials/maps.h>
@@ -120 +121 @@
  ************************************/
 // unused
-#if 0
+/*
 static void initSSpecialCC (ideal F, ideal Q, ideal P,kStrategy strat)
 {
@@ -268 +269 @@
 #endif
 }
-#endif
+*/
 
 /*****************
@@ -277 +278 @@
   int j;
   kStrategy strat = new skStrategy;
-
-//  if (TEST_OPT_PROT)
-//  {
-//    writeTime("start InterRed:");
-//    mflush();
-//  }
-  //strat->syzComp     = 0;
+/*
+  if (TEST_OPT_PROT)
+  {
+    writeTime("start InterRed:");
+    mflush();
+  }
+  strat->syzComp     = 0;
+*/
   strat->kHEdgeFound = (currRing->ppNoether) != NULL;
   strat->kNoether=pCopy((currRing->ppNoether));
@@ -346 +348 @@
     strat->fromQ = NULL;
   }
-//  if (TEST_OPT_PROT)
-//  {
-//    writeTime("end Interred:");
-//    mflush();
-//  }
+/*
+  if (TEST_OPT_PROT)
+  {
+    writeTime("end Interred:");
+    mflush();
+  }
+*/
   ideal shdl=strat->Shdl;
   idSkipZeroes(shdl);
@@ -358 +362 @@
 }
 
-//unused
-#if 0
+#ifdef TIME_TEST
 static void TimeString(clock_t tinput, clock_t tostd, clock_t tif,clock_t tstd,
                        clock_t tlf,clock_t tred, clock_t tnw, int step)
@@ -397 +400 @@
         ((((double) xtextra)/1000000)/totm)*100);
 }
-#endif
-
-//unused
-#if 0
+
 static void TimeStringFractal(clock_t tinput, clock_t tostd, clock_t tif,clock_t tstd,
                        clock_t textra, clock_t tlf,clock_t tred, clock_t tnw)
@@ -442 +442 @@
 }
 #endif
-
+/*
 #if defined(CHECK_IDEAL_MWALK) || defined(ENDWALKS)
 static void headidString(ideal L, char* st)
@@ -496 +496 @@
 }
 #endif
-
+*/
 
 static void ivString(intvec* iv, const char* ch)
@@ -510 +510 @@
 }
 
-//unused
-#if 0
+#ifdef PRINT_VECTORS
 static void MivString(intvec* iva, intvec* ivb, intvec* ivc)
 {
@@ -558 +557 @@
   }
   return p0;
+}
+
+/*****************************************************************************
+ * compute the gcd of the entries of the vectors curr_weight and diff_weight *
+ *****************************************************************************/
+static int simplify_gcd(intvec* curr_weight, intvec* diff_weight)
+{
+  int j;
+  int nRing = currRing->N;
+  int gcd_tmp = (*curr_weight)[0];
+  for (j=1; j<nRing; j++)
+  {
+    gcd_tmp = gcd(gcd_tmp, (*curr_weight)[j]);
+    if(gcd_tmp == 1)
+    {
+      break;
+    }
+  }
+  if(gcd_tmp != 1)
+  {
+    for (j=0; j<nRing; j++)
+    {
+    gcd_tmp = gcd(gcd_tmp, (*diff_weight)[j]);
+    if(gcd_tmp == 1)
+      {
+        break;
+      }
+    }
+  }
+  return gcd_tmp;
 }
 
@@ -774 +803 @@
   for(i=nG-1; i>=0; i--)
   {
-    mi = MpolyInitialForm(G->m[i], iv);
-    gi = G->m[i];
-
+    mi = pHead(MpolyInitialForm(G->m[i], iv));
+    //Print("\n **// test_w_in_ConeCC: lm(initial)= %s \n",pString(mi));
+    gi = pHead(G->m[i]);
+    //Print("\n **// test_w_in_ConeCC: lm(ideal)= %s \n",pString(gi));
     if(mi == NULL)
     {
@@ -953 +983 @@
 }
 
-/*****************************************************************************
-* create a weight matrix order as intvec of an extra weight vector (a(iv),lp)*
-******************************************************************************/
+/*********************************************************************************
+* create a weight matrix order as intvec of an extra weight vector (a(iv),M(iw)) *
+**********************************************************************************/
 intvec* MivMatrixOrderRefine(intvec* iv, intvec* iw)
 {
-  assume(iv->length() == iw->length());
-  int i, nR = iv->length();
-
+  assume((iv->length())*(iv->length()) == iw->length());
+  int i,j, nR = iv->length();
+  
   intvec* ivm = new intvec(nR*nR);
 
@@ -966 +996 @@
   {
     (*ivm)[i] = (*iv)[i];
-    (*ivm)[i+nR] = (*iw)[i];
-  }
-  for(i=2; i<nR; i++)
-  {
-    (*ivm)[i*nR+i-2] = 1;
+  }
+  for(i=1; i<nR; i++)
+  {
+    for(j=0; j<nR; j++)
+    {
+      (*ivm)[j+i*nR] = (*iw)[j+i*nR];
+    }
   }
   return ivm;
@@ -1005 +1037 @@
  * print the max total degree and the max coefficient of G                   *
  *****************************************************************************/
-#if 0
+/*
 static void checkComplexity(ideal G, char* cG)
 {
@@ -1046 +1078 @@
   PrintLn();
 }
-#endif
+*/
 
 /*****************************************************************************
@@ -1068 +1100 @@
   intvec* v_null =  new intvec(nV);
 
-
   // Check that the perturbed degree is valid
   if(pdeg > nV || pdeg <= 0)
@@ -1082 +1113 @@
   }
   mpz_t *pert_vector = (mpz_t*)omAlloc(nV*sizeof(mpz_t));
-  //mpz_t *pert_vector1 = (mpz_t*)omAlloc(nV*sizeof(mpz_t));
+  mpz_t *pert_vector1 = (mpz_t*)omAlloc(nV*sizeof(mpz_t));
 
   for(i=0; i<nV; i++)
   {
     mpz_init_set_si(pert_vector[i], (*ivtarget)[i]);
-   // mpz_init_set_si(pert_vector1[i], (*ivtarget)[i]);
+    mpz_init_set_si(pert_vector1[i], (*ivtarget)[i]);
   }
   // Calculate max1 = Max(A2)+Max(A3)+...+Max(Apdeg),
@@ -1167 +1198 @@
     }
   }
+
+  // 2147483647 is max. integer representation in SINGULAR
+  mpz_t sing_int;
+  mpz_init_set_ui(sing_int,  2147483647);
+
+  mpz_t check_int;
+  mpz_init_set_ui(check_int,  100000);
+
   mpz_t ztemp;
   mpz_init(ztemp);
@@ -1186 +1225 @@
   }
 
-  intvec *pert_vector1= new intvec(nV);
-  j = 0;
   for(i=0; i<nV; i++)
   {
-    (* pert_vector1)[i] = mpz_get_si(pert_vector[i]);
-    (* pert_vector1)[i] = 0.1*(* pert_vector1)[i];
-    (* pert_vector1)[i] = floor((* pert_vector1)[i] + 0.5);
-    if((* pert_vector1)[i] == 0)
-    {
-      j++;
-    }
-  }
-  if(j > nV - 1)
-  {
-    // Print("\n//  MPertVectors: geaenderter vector gleich Null! \n");
-    delete pert_vector1;
-    goto CHECK_OVERFLOW;
-  }
-
-// check that the perturbed weight vector lies in the Groebner cone
-  if(test_w_in_ConeCC(G,pert_vector1) != 0)
-  {
-    // Print("\n//  MPertVectors: geaenderter vector liegt in Groebnerkegel! \n");
+    if(mpz_cmp(pert_vector[i], check_int)>=0)
+    {
+      for(j=0; j<nV; j++)
+      {
+        mpz_fdiv_q_ui(pert_vector1[j], pert_vector[j], 100);
+      }
+    }
+  }
+
+  intvec* result = new intvec(nV);
+
+  int ntrue=0;
+
+  for(i=0; i<nV; i++)
+  {
+    (*result)[i] = mpz_get_si(pert_vector1[i]);
+    if(mpz_cmp(pert_vector1[i], sing_int)>=0)
+    {
+      ntrue++;
+    }
+  }
+  if(ntrue > 0 || test_w_in_ConeCC(G,result)==0)
+  {
+    ntrue=0;
     for(i=0; i<nV; i++)
     {
-      mpz_set_si(pert_vector[i], (*pert_vector1)[i]);
-    }
-  }
-  else
-  {
-    //Print("\n// MpertVectors: geaenderter vector liegt nicht in Groebnerkegel! \n");
-  }
-  delete pert_vector1;
-
-  CHECK_OVERFLOW:
-  intvec* result = new intvec(nV);
-
-  /* 2147483647 is max. integer representation in SINGULAR */
-  mpz_t sing_int;
-  mpz_init_set_ui(sing_int,  2147483647);
-
-  int ntrue=0;
-  for(i=0; i<nV; i++)
-  {
-    (*result)[i] = mpz_get_si(pert_vector[i]);
-    if(mpz_cmp(pert_vector[i], sing_int)>=0)
-    {
-      ntrue++;
-      if(Overflow_Error == FALSE)
-      {
-        Overflow_Error = TRUE;
-        PrintS("\n// ** OVERFLOW in \"MPertvectors\": ");
-        mpz_out_str( stdout, 10, pert_vector[i]);
-        PrintS(" is greater than 2147483647 (max. integer representation)");
-        Print("\n//  So vector[%d] := %d is wrong!!", i+1, (*result)[i]);
-      }
-    }
-  }
-
-  if(Overflow_Error == TRUE)
-  {
-    ivString(result, "pert_vector");
-    Print("\n// %d element(s) of it is overflow!!", ntrue);
+      (*result)[i] = mpz_get_si(pert_vector[i]);
+      if(mpz_cmp(pert_vector[i], sing_int)>=0)
+      {
+        ntrue++;
+        if(Overflow_Error == FALSE)
+        {
+          Overflow_Error = TRUE;
+          PrintS("\n// ** OVERFLOW in \"MPertvectors\": ");
+          mpz_out_str( stdout, 10, pert_vector[i]);
+          PrintS(" is greater than 2147483647 (max. integer representation)");
+          Print("\n//  So vector[%d] := %d is wrong!!", i+1, (*result)[i]);
+        }
+      }
+    }
+
+    if(Overflow_Error == TRUE)
+    {
+      ivString(result, "pert_vector");
+      Print("\n// %d element(s) of it is overflow!!", ntrue);
+    }
   }
 
   mpz_clear(ztemp);
   mpz_clear(sing_int);
+  mpz_clear(check_int);
   omFree(pert_vector);
-  //omFree(pert_vector1);
+  omFree(pert_vector1);
   mpz_clear(tot_deg);
   mpz_clear(maxdeg);
@@ -1456 +1481 @@
 
 //unused
-#if 0
+/*
 static intvec* MatrixOrderdp(int nV)
 {
@@ -1472 +1497 @@
   return(ivM);
 }
-#endif
+*/
 
 intvec* MivUnit(int nV)
@@ -1549 +1574 @@
     mpz_cdiv_q_ui(inveps, inveps, nV);
   }
-  //PrintS("\n// choose the \"small\" inverse epsilon!");
+  // choose the small inverse epsilon
 #endif
 
@@ -1583 +1608 @@
 
     for(j=0; j<nV; j++)
-      {
+    {
       mpz_init_set(pert_vector[i*nV+j],ivtemp[j]);
-      }
-  }
-
-  /* 2147483647 is max. integer representation in SINGULAR */
+    }
+  }
+
+  // 2147483647 is max. integer representation in SINGULAR
   mpz_t sing_int;
   mpz_init_set_ui(sing_int,  2147483647);
@@ -1611 +1636 @@
     mpz_divexact(pert_vector[i], pert_vector[i], ztmp);
     (* result)[i] = mpz_get_si(pert_vector[i]);
-  }
-
-  j = 0;
-  for(i=0; i<nV; i++)
-  {
-    (* result1)[i] = mpz_get_si(pert_vector[i]);
-    (* result1)[i] = 0.1*(* result1)[i];
-    (* result1)[i] = floor((* result1)[i] + 0.5);
-    if((* result1)[i] == 0)
-    {
-      j++;
-    }
-  }
-  if(j > nV - 1)
-  {
-    // Print("\n//  MfPertwalk: geaenderter vector gleich Null! \n");
-    delete result1;
-    goto CHECK_OVERFLOW;
-  }
-
-// check that the perturbed weight vector lies in the Groebner cone
-  if(test_w_in_ConeCC(G,result1) != 0)
-  {
-    // Print("\n//  MfPertwalk: geaenderter vector liegt in Groebnerkegel! \n");
-    delete result;
-    result = result1;
-    for(i=0; i<nV; i++)
-    {
-      mpz_set_si(pert_vector[i], (*result1)[i]);
-    }
-  }
-  else
-  {
-    delete result1;
-    // Print("\n// Mfpertwalk: geaenderter vector liegt nicht in Groebnerkegel! \n");
   }
 
@@ -1685 +1675 @@
     while(p!=NULL)
     {
-      p_Setm(p,currRing); pIter(p);
+      p_Setm(p,currRing);
+      pIter(p);
     }
   }
@@ -1768 +1759 @@
 
 //unused
-#if 0
+/*
 static void checkidealCC(ideal G, char* Ch)
 {
@@ -1794 +1785 @@
   PrintLn();
 }
-#endif
+*/
 
 //unused
-#if 0
+/*
 static void HeadidString(ideal L, char* st)
 {
@@ -1809 +1800 @@
   Print(" %s;\n", pString(pHead(L->m[nL])));
 }
-#endif
-
+
+*/
 static inline int MivComp(intvec* iva, intvec* ivb)
 {
@@ -1859 +1850 @@
 }
 
+
+/**************************************************************
+ * Look for the position of the smallest absolut value in vec *
+ **************************************************************/
+static int MivAbsMaxArg(intvec* vec)
+{
+  int k = MivAbsMax(vec);
+  int i=0;
+  while(1)
+  {
+    if((*vec)[i] == k || (*vec)[i] == -k)
+    {
+      break;
+    }
+    i++;
+  }
+  return i;
+}
+
+
 /**********************************************************************
  * Compute a next weight vector between curr_weight and target_weight *
  * with respect to an ideal <G>.                                      *
 **********************************************************************/
+/*
 static intvec* MwalkNextWeightCC(intvec* curr_weight, intvec* target_weight,
                                  ideal G)
@@ -1873 +1885 @@
 
   int nRing = currRing->N;
-  int checkRed, j, kkk, nG = IDELEMS(G);
+  int checkRed, j, nG = IDELEMS(G);
   intvec* ivtemp;
 
@@ -1911 +1923 @@
   mpz_init(dcw);
 
-  //int tn0, tn1, tz1, ncmp, gcd_tmp, ntmp;
   int gcd_tmp;
   intvec* diff_weight = MivSub(target_weight, curr_weight);
@@ -1917 +1928 @@
   intvec* diff_weight1 = MivSub(target_weight, curr_weight);
   poly g;
-  //poly g, gw;
+
   for (j=0; j<nG; j++)
   {
@@ -1934 +1945 @@
         mpz_set_si(MwWd, MivDotProduct(ivtemp, curr_weight));
         mpz_sub(s_zaehler, deg_w0_p1, MwWd);
-
         if(mpz_cmp(s_zaehler, t_null) != 0)
         {
           mpz_set_si(MwWd, MivDotProduct(ivtemp, diff_weight));
           mpz_sub(s_nenner, MwWd, deg_d0_p1);
-
           // check for 0 < s <= 1
           if( (mpz_cmp(s_zaehler,t_null) > 0 &&
@@ -1979 +1988 @@
     }
   }
-//Print("\n// Alloc Size = %d \n", nRing*sizeof(mpz_t));
+  //Print("\n// Alloc Size = %d \n", nRing*sizeof(mpz_t));
   mpz_t *vec=(mpz_t*)omAlloc(nRing*sizeof(mpz_t));
 
 
-  // there is no 0<t<1 and define the next weight vector that is equal to the current weight vector
+  // there is no 0<t<1 and define the next weight vector that is equal
+  // to the current weight vector
   if(mpz_cmp(t_nenner, t_null) == 0)
   {
-    #ifndef SING_NDEBUG
-    Print("\n//MwalkNextWeightCC: t_nenner ist Null!");
-    #endif
+#ifndef SING_NDEBUG
+    Print("\n//MwalkNextWeightCC: t_nenner=0\n");
+#endif
     delete diff_weight;
     diff_weight = ivCopy(curr_weight);//take memory
@@ -2054 +2064 @@
 #endif
 
-  //  BOOLEAN isdwpos;
-
-  // construct a new weight vector
+// construct a new weight vector and check whether vec[j] is overflow,
+// i.e. vec[j] > 2^31.
+// If vec[j] doesn't overflow, define a weight vector. Otherwise,
+// report that overflow appears. In the second case, test whether the
+// the correctness of the new vector plays an important role
+
   for (j=0; j<nRing; j++)
   {
@@ -2100 +2113 @@
     }
   }
-
+  // reduce the vector with the gcd
+  if(mpz_cmp_si(ggt,1) != 0)
+  {
+    for (j=0; j<nRing; j++)
+    {
+      mpz_divexact(vec[j], vec[j], ggt);
+    }
+  }
 #ifdef  NEXT_VECTORS_CC
   PrintS("\n// gcd of elements of the vector: ");
@@ -2106 +2126 @@
 #endif
 
-/**********************************************************************
-* construct a new weight vector and check whether vec[j] is overflow, *
-* i.e. vec[j] > 2^31.                                                 *
-* If vec[j] doesn't overflow, define a weight vector. Otherwise,      *
-* report that overflow appears. In the second case, test whether the  *
-* the correctness of the new vector plays an important role           *
-**********************************************************************/
-  kkk=0;
   for(j=0; j<nRing; j++)
-    {
+  {
     if(mpz_cmp(vec[j], sing_int_half) >= 0)
-      {
+    {
       goto REDUCTION;
-      }
-    }
+    }
+  }
   checkRed = 1;
   for (j=0; j<nRing; j++)
@@ -2129 +2141 @@
 
   REDUCTION:
+  checkRed = 1;
   for (j=0; j<nRing; j++)
   {
-    (*diff_weight)[j] = mpz_get_si(vec[j]);
-  }
-  while(MivAbsMax(diff_weight) >= 5)
-  {
-    for (j=0; j<nRing; j++)
-    {
-      if(mpz_cmp_si(ggt,1)==0)
-      {
-        (*diff_weight1)[j] = floor(0.1*(*diff_weight)[j] + 0.5);
-        // Print("\n//  vector[%d] = %d \n",j+1, (*diff_weight1)[j]);
-      }
-      else
-      {
-        mpz_divexact(vec[j], vec[j], ggt);
-        (*diff_weight1)[j] = floor(0.1*(*diff_weight)[j] + 0.5);
-        // Print("\n//  vector[%d] = %d \n",j+1, (*diff_weight1)[j]);
-      }
-/*
-      if((*diff_weight1)[j] == 0)
-      {
-        kkk = kkk + 1;
-      }
-*/
-    }
-
-
-/*
-    if(kkk > nRing - 1)
-    {
-      // diff_weight was reduced to zero
-      // Print("\n // MwalkNextWeightCC: geaenderter Vector gleich Null! \n");
-      goto TEST_OVERFLOW;
-    }
-*/
-
-    if(test_w_in_ConeCC(G,diff_weight1) != 0)
-    {
-      Print("\n// MwalkNextWeightCC: geaenderter vector liegt in Groebnerkegel! \n");
-      for (j=0; j<nRing; j++)
-      {
-        (*diff_weight)[j] = (*diff_weight1)[j];
-      }
-      if(MivAbsMax(diff_weight) < 5)
-      {
-        checkRed = 1;
-        goto SIMPLIFY_GCD;
-      }
-    }
-    else
-    {
-      // Print("\n// MwalkNextWeightCC: geaenderter vector liegt nicht in Groebnerkegel! \n");
-      break;
-    }
+    (*diff_weight1)[j] = mpz_get_si(vec[j]);
+  }
+  while(test_w_in_ConeCC(G,diff_weight1))
+  {
+    for(j=0; j<nRing; j++)
+    {
+      (*diff_weight)[j] = (*diff_weight1)[j];
+      mpz_set_si(vec[j], (*diff_weight)[j]);      
+    }
+    for(j=0; j<nRing; j++)
+    {
+      (*diff_weight1)[j] = floor(0.1*(*diff_weight)[j] + 0.5);
+    }
+  }
+  if(MivAbsMax(diff_weight)>10000)
+  {
+    for(j=0; j<nRing; j++)
+    {
+      (*diff_weight1)[j] = (*diff_weight)[j];
+    }
+    j = 0;
+    while(test_w_in_ConeCC(G,diff_weight1))
+    {
+      (*diff_weight)[j] = (*diff_weight1)[j];
+      mpz_set_si(vec[j], (*diff_weight)[j]);
+      j = MivAbsMaxArg(diff_weight1);
+      (*diff_weight1)[j] = floor(0.1*(*diff_weight1)[j] + 0.5);
+    }
+    goto SIMPLIFY_GCD;
   }
 
@@ -2222 +2211 @@
    mpz_clear(t_null);
 
-
-
   if(Overflow_Error == FALSE)
   {
     Overflow_Error = nError;
   }
- rComplete(currRing);
-  for(kkk=0; kkk<IDELEMS(G);kkk++)
-  {
-    poly p=G->m[kkk];
+  rComplete(currRing);
+  for(j=0; j<IDELEMS(G); j++)
+  {
+    poly p=G->m[j];
     while(p!=NULL)
     {
@@ -2240 +2227 @@
 return diff_weight;
 }
+*/
+/**********************************************************************
+ * Compute a next weight vector between curr_weight and target_weight *
+ * with respect to an ideal <G>.                                      *
+**********************************************************************/
+static intvec* MwalkNextWeightCC(intvec* curr_weight, intvec* target_weight,
+                                 ideal G)
+{
+  BOOLEAN nError = Overflow_Error;
+  Overflow_Error = FALSE;
+
+  assume(currRing != NULL && curr_weight != NULL &&
+         target_weight != NULL && G != NULL);
+
+  int nRing = currRing->N;
+  int checkRed, j, nG = IDELEMS(G);
+  intvec* ivtemp;
+
+  mpz_t t_zaehler, t_nenner;
+  mpz_init(t_zaehler);
+  mpz_init(t_nenner);
+
+  mpz_t s_zaehler, s_nenner, temp, MwWd;
+  mpz_init(s_zaehler);
+  mpz_init(s_nenner);
+  mpz_init(temp);
+  mpz_init(MwWd);
+
+  mpz_t sing_int;
+  mpz_init(sing_int);
+  mpz_set_si(sing_int,  2147483647);
+
+  mpz_t sing_int_half;
+  mpz_init(sing_int_half);
+  mpz_set_si(sing_int_half,  3*(1073741824/2));
+
+  mpz_t deg_w0_p1, deg_d0_p1;
+  mpz_init(deg_w0_p1);
+  mpz_init(deg_d0_p1);
+
+  mpz_t sztn, sntz;
+  mpz_init(sztn);
+  mpz_init(sntz);
+
+  mpz_t t_null;
+  mpz_init(t_null);
+
+  mpz_t ggt;
+  mpz_init(ggt);
+
+  mpz_t dcw;
+  mpz_init(dcw);
+
+  int gcd_tmp;
+  //intvec* diff_weight = MivSub(target_weight, curr_weight);
+
+  intvec* diff_weight1 = new intvec(nRing); //MivSub(target_weight, curr_weight);
+  poly g;
+
+  // reduce the size of the entries of the current weight vector
+  if(TEST_OPT_REDSB)
+  {
+    for (j=0; j<nRing; j++)
+    {
+      (*diff_weight1)[j] = (*curr_weight)[j];
+    }
+    while(MivAbsMax(diff_weight1)>10000 && test_w_in_ConeCC(G,diff_weight1)==1)
+    {
+      for(j=0; j<nRing; j++)
+      {
+        (*curr_weight)[j] = (*diff_weight1)[j];  
+      }
+      for(j=0; j<nRing; j++)
+      {
+        (*diff_weight1)[j] = floor(0.1*(*diff_weight1)[j] + 0.5);
+      }
+    }
+
+    if(MivAbsMax(curr_weight)>100000)
+    {
+      for(j=0; j<nRing; j++)
+      {
+        (*diff_weight1)[j] = (*curr_weight)[j];
+      }
+      j = 0;
+      while(test_w_in_ConeCC(G,diff_weight1)==1 && MivAbsMax(diff_weight1)>1000)
+      {
+        (*curr_weight)[j] = (*diff_weight1)[j];
+        j = MivAbsMaxArg(diff_weight1);
+        (*diff_weight1)[j] = floor(0.1*(*diff_weight1)[j] + 0.5);
+      }
+    }
+
+  }
+  intvec* diff_weight = MivSub(target_weight, curr_weight);
+
+  // compute a suitable next weight vector
+  for (j=0; j<nG; j++)
+  {
+    g = G->m[j];
+    if (g != NULL)
+    {
+      ivtemp = MExpPol(g);
+      mpz_set_si(deg_w0_p1, MivDotProduct(ivtemp, curr_weight));
+      mpz_set_si(deg_d0_p1, MivDotProduct(ivtemp, diff_weight));
+      delete ivtemp;
+
+      pIter(g);
+      while (g != NULL)
+      {
+        ivtemp = MExpPol(g);
+        mpz_set_si(MwWd, MivDotProduct(ivtemp, curr_weight));
+        mpz_sub(s_zaehler, deg_w0_p1, MwWd);
+        if(mpz_cmp(s_zaehler, t_null) != 0)
+        {
+          mpz_set_si(MwWd, MivDotProduct(ivtemp, diff_weight));
+          mpz_sub(s_nenner, MwWd, deg_d0_p1);
+          // check for 0 < s <= 1
+          if( (mpz_cmp(s_zaehler,t_null) > 0 &&
+               mpz_cmp(s_nenner, s_zaehler)>=0) ||
+              (mpz_cmp(s_zaehler, t_null) < 0 &&
+               mpz_cmp(s_nenner, s_zaehler)<=0))
+          {
+            // make both positive
+            if (mpz_cmp(s_zaehler, t_null) < 0)
+            {
+              mpz_neg(s_zaehler, s_zaehler);
+              mpz_neg(s_nenner, s_nenner);
+            }
+
+            //compute a simple fraction of s
+            cancel(s_zaehler, s_nenner);
+
+            if(mpz_cmp(t_nenner, t_null) != 0)
+            {
+              mpz_mul(sztn, s_zaehler, t_nenner);
+              mpz_mul(sntz, s_nenner, t_zaehler);
+
+              if(mpz_cmp(sztn,sntz) < 0)
+              {
+                mpz_add(t_nenner, t_null, s_nenner);
+                mpz_add(t_zaehler,t_null, s_zaehler);
+              }
+            }
+            else
+            {
+              mpz_add(t_nenner, t_null, s_nenner);
+              mpz_add(t_zaehler,t_null, s_zaehler);
+            }
+          }
+        }
+        pIter(g);
+        delete ivtemp;
+      }
+    }
+  }
+  //Print("\n// Alloc Size = %d \n", nRing*sizeof(mpz_t));
+  mpz_t *vec=(mpz_t*)omAlloc(nRing*sizeof(mpz_t));
+
+
+  // there is no 0<t<1 and define the next weight vector that is equal
+  // to the current weight vector
+  if(mpz_cmp(t_nenner, t_null) == 0)
+  {
+#ifndef SING_NDEBUG
+    Print("\n//MwalkNextWeightCC: t_nenner=0\n");
+#endif
+    delete diff_weight;
+    diff_weight = ivCopy(curr_weight);//take memory
+    goto FINISH;
+  }
+
+  // define the target vector as the next weight vector, if t = 1
+  if(mpz_cmp_si(t_nenner, 1)==0 && mpz_cmp_si(t_zaehler,1)==0)
+  {
+    delete diff_weight;
+    diff_weight = ivCopy(target_weight); //this takes memory
+    goto FINISH;
+  }
+
+   //checkRed = 0;
+
+  SIMPLIFY_GCD:
+
+  // simplify the vectors curr_weight and diff_weight (C-int)
+  gcd_tmp = (*curr_weight)[0];
+
+  for (j=1; j<nRing; j++)
+  {
+    gcd_tmp = gcd(gcd_tmp, (*curr_weight)[j]);
+    if(gcd_tmp == 1)
+    {
+      break;
+    }
+  }
+  if(gcd_tmp != 1)
+  {
+    for (j=0; j<nRing; j++)
+    {
+      gcd_tmp = gcd(gcd_tmp, (*diff_weight)[j]);
+      if(gcd_tmp == 1)
+      {
+        break;
+      }
+    }
+  }
+  if(gcd_tmp != 1)
+  {
+    for (j=0; j<nRing; j++)
+    {
+      (*curr_weight)[j] =  (*curr_weight)[j]/gcd_tmp;
+      (*diff_weight)[j] =  (*diff_weight)[j]/gcd_tmp;
+    }
+  }
+  if(checkRed > 0)
+  {
+    for (j=0; j<nRing; j++)
+    {
+      mpz_set_si(vec[j], (*diff_weight)[j]);
+    }
+    goto TEST_OVERFLOW;
+  }
+
+#ifdef  NEXT_VECTORS_CC
+  Print("\n// gcd of the weight vectors (current and target) = %d", gcd_tmp);
+  ivString(curr_weight, "new cw");
+  ivString(diff_weight, "new dw");
+
+  PrintS("\n// t_zaehler: ");  mpz_out_str( stdout, 10, t_zaehler);
+  PrintS(", t_nenner: ");  mpz_out_str( stdout, 10, t_nenner);
+#endif
+
+// construct a new weight vector and check whether vec[j] is overflow, i.e. vec[j] > 2^31.
+// If vec[j] doesn't overflow, define a weight vector. Otherwise, report that overflow
+// appears. In the second case, test whether the the correctness of the new vector plays
+// an important role
+
+  for (j=0; j<nRing; j++)
+  {
+    mpz_set_si(dcw, (*curr_weight)[j]);
+    mpz_mul(s_nenner, t_nenner, dcw);
+
+    if( (*diff_weight)[j]>0)
+    {
+      mpz_mul_ui(s_zaehler, t_zaehler, (*diff_weight)[j]);
+    }
+    else
+    {
+      mpz_mul_ui(s_zaehler, t_zaehler, -(*diff_weight)[j]);
+      mpz_neg(s_zaehler, s_zaehler);
+    }
+    mpz_add(sntz, s_nenner, s_zaehler);
+    mpz_init_set(vec[j], sntz);
+
+#ifdef NEXT_VECTORS_CC
+    Print("\n//   j = %d ==> ", j);
+    PrintS("(");
+    mpz_out_str( stdout, 10, t_nenner);
+    Print(" * %d)", (*curr_weight)[j]);
+    Print(" + ("); mpz_out_str( stdout, 10, t_zaehler);
+    Print(" * %d) =  ",  (*diff_weight)[j]);
+    mpz_out_str( stdout, 10, s_nenner);
+    PrintS(" + ");
+    mpz_out_str( stdout, 10, s_zaehler);
+    PrintS(" = "); mpz_out_str( stdout, 10, sntz);
+    Print(" ==> vector[%d]: ", j); mpz_out_str(stdout, 10, vec[j]);
+#endif
+
+    if(j==0)
+    {
+      mpz_set(ggt, sntz);
+    }
+    else
+    {
+      if(mpz_cmp_si(ggt,1) != 0)
+      {
+        mpz_gcd(ggt, ggt, sntz);
+      }
+    }
+  }
+  // reduce the vector with the gcd
+  if(mpz_cmp_si(ggt,1) != 0)
+  {
+    for (j=0; j<nRing; j++)
+    {
+      mpz_divexact(vec[j], vec[j], ggt);
+    }
+  }
+#ifdef  NEXT_VECTORS_CC
+  PrintS("\n// gcd of elements of the vector: ");
+  mpz_out_str( stdout, 10, ggt);
+#endif
+
+  for (j=0; j<nRing; j++)
+  {
+    (*diff_weight)[j] = mpz_get_si(vec[j]);
+  }
+
+ TEST_OVERFLOW:
+
+  for (j=0; j<nRing; j++)
+  {
+    if(mpz_cmp(vec[j], sing_int)>=0)
+    {
+      if(Overflow_Error == FALSE)
+      {
+        Overflow_Error = TRUE;
+        PrintS("\n// ** OVERFLOW in \"MwalkNextWeightCC\": ");
+        mpz_out_str( stdout, 10, vec[j]);
+        PrintS(" is greater than 2147483647 (max. integer representation)\n");
+        Print("//  So vector[%d] := %d is wrong!!\n",j+1, vec[j]);// vec[j] is mpz_t
+      }
+    }
+  }
+
+ FINISH:
+   delete diff_weight1;
+   mpz_clear(t_zaehler);
+   mpz_clear(t_nenner);
+   mpz_clear(s_zaehler);
+   mpz_clear(s_nenner);
+   mpz_clear(sntz);
+   mpz_clear(sztn);
+   mpz_clear(temp);
+   mpz_clear(MwWd);
+   mpz_clear(deg_w0_p1);
+   mpz_clear(deg_d0_p1);
+   mpz_clear(ggt);
+   omFree(vec);
+   mpz_clear(sing_int_half);
+   mpz_clear(sing_int);
+   mpz_clear(dcw);
+   mpz_clear(t_null);
+
+  if(Overflow_Error == FALSE)
+  {
+    Overflow_Error = nError;
+  }
+  rComplete(currRing);
+  for(j=0; j<IDELEMS(G); j++)
+  {
+    poly p=G->m[j];
+    while(p!=NULL)
+    {
+      p_Setm(p,currRing);
+      pIter(p);
+    }
+  }
+return diff_weight;
+}
+
 
 /**********************************************************************
@@ -2271 +2609 @@
 }
 
-/**************************************************************
+/********************************************************************
  * define and execute a new ring which order is (a(vb),a(va),lp,C)  *
- * ************************************************************/
-#if 0
-// unused
+ * ******************************************************************/
 static void VMrHomogeneous(intvec* va, intvec* vb)
 {
@@ -2353 +2689 @@
   rChangeCurrRing(r);
 }
-#endif
+
 
 /**************************************************************
@@ -2428 +2764 @@
 }
 
-static ring VMrDefault1(intvec* va)
-{
-
-  if ((currRing->ppNoether)!=NULL)
-  {
-    pDelete(&(currRing->ppNoether));
-  }
-  if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) ||
-      ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data))))
-  {
-    sLastPrinted.CleanUp();
-  }
-
-  ring r = (ring) omAlloc0Bin(sip_sring_bin);
-  int i, nv = currRing->N;
-
-  r->cf  = currRing->cf;
-  r->N   = currRing->N;
-
-  int nb = 4;
-
-  //names
-  char* Q; // In order to avoid the corrupted memory, do not change.
-  r->names = (char **) omAlloc0(nv * sizeof(char_ptr));
-  for(i=0; i<nv; i++)
-  {
-    Q = currRing->names[i];
-    r->names[i]  = omStrDup(Q);
-  }
-
-  /*weights: entries for 3 blocks: NULL Made:???*/
-  r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
-  r->wvhdl[0] = (int*) omAlloc(nv*sizeof(int));
-  for(i=0; i<nv; i++)
-    r->wvhdl[0][i] = (*va)[i];
-
-  /* order: a,lp,C,0 */
-  r->order = (int *) omAlloc(nb * sizeof(int *));
-  r->block0 = (int *)omAlloc0(nb * sizeof(int *));
-  r->block1 = (int *)omAlloc0(nb * sizeof(int *));
-
-  // ringorder a for the first block: var 1..nv
-  r->order[0]  = ringorder_a;
-  r->block0[0] = 1;
-  r->block1[0] = nv;
-
-  // ringorder lp for the second block: var 1..nv
-  r->order[1]  = ringorder_lp;
-  r->block0[1] = 1;
-  r->block1[1] = nv;
-
-  // ringorder C for the third block
-  // it is very important within "idLift",
-  // especially, by ring syz_ring=rCurrRingAssure_SyzComp();
-  // therefore, nb  must be (nBlocks(currRing)  + 1)
-  r->order[2]  = ringorder_C;
-
-  // the last block: everything is 0
-  r->order[3]  = 0;
-
-  // polynomial ring
-  r->OrdSgn    = 1;
-
-  // complete ring intializations
-
-  rComplete(r);
-
-  //rChangeCurrRing(r);
-  return r;
-}
-
 /****************************************************************
  * define and execute a new ring with ordering (a(va),Wp(vb),C) *
  * **************************************************************/
-
 static ring VMrRefine(intvec* va, intvec* vb)
 {
@@ -2576 +2840 @@
 
   // complete ring intializations
-
+  
   rComplete(r);
 
@@ -2806 +3070 @@
 }
 
-
-/* define a ring with parameters und change to it */
-/* DefRingPar and DefRingParlp corrupt still memory */
+/***************************************************
+* define a ring with parameters und change to it   *
+* DefRingPar and DefRingParlp corrupt still memory *
+****************************************************/
 static void DefRingPar(intvec* va)
 {
@@ -2956 +3221 @@
 }
 
-//unused
-/**************************************************************
- * check wheather one or more components of a vector are zero *
- **************************************************************/
-#if 0
+/*************************************************************
+ * check whether one or more components of a vector are zero *
+ *************************************************************/
 static int isNolVector(intvec* hilb)
 {
@@ -2973 +3236 @@
   return 0;
 }
-#endif
+
+/*************************************************************
+ * check whether one or more components of a vector are <= 0 *
+ *************************************************************/
+static int isNegNolVector(intvec* hilb)
+{
+  int i;
+  for(i=hilb->length()-1; i>=0; i--)
+  {
+    if((* hilb)[i]<=0)
+    {
+      return 1;
+    }
+  }
+  return 0;
+}
+
+/**************************************************************************
+* Gomega is the initial ideal of G w. r. t. the current weight vector     *
+* curr_weight. Check whether curr_weight lies on a border of the Groebner *
+* cone, i. e. check whether a monomial is divisible by a leading monomial *
+***************************************************************************/
+static ideal middleOfCone(ideal G, ideal Gomega)
+{
+  BOOLEAN middle = FALSE;
+  int i,j,N = IDELEMS(Gomega);
+  poly p,lm,factor1,factor2;
+
+  ideal Go = idCopy(G);
+  
+  // check whether leading monomials of G and Gomega coincide
+  // and return NULL if not
+  for(i=0; i<N; i++)
+  {
+    if(!pIsConstant(pSub(pCopy(Gomega->m[i]),pCopy(pHead(G->m[i])))))
+    {
+      idDelete(&Go);
+      return NULL; 
+    }
+  }
+  for(i=0; i<N; i++)
+  {
+    for(j=0; j<N; j++)
+    {
+      if(i!=j)
+      {
+        p = pCopy(Gomega->m[i]);
+        lm = pCopy(Gomega->m[j]);
+        pIter(p);
+        while(p!=NULL)
+        {
+          if(pDivisibleBy(lm,p))
+          {
+            if(middle == FALSE)
+            {
+              middle = TRUE;
+            }
+            factor1 = singclap_pdivide(pHead(p),lm,currRing);
+            factor2 = pMult(pCopy(factor1),pCopy(Go->m[j]));
+            pDelete(&factor1);
+            Go->m[i] = pAdd((Go->m[i]),pNeg(pCopy(factor2)));
+            pDelete(&factor2);
+          }
+          pIter(p);
+        }
+        pDelete(&lm);
+        pDelete(&p);
+      }
+    }
+  }
+
+  if(middle == TRUE)
+  {
+    return Go;
+  }
+  idDelete(&Go);
+  return NULL; 
+}
 
 /******************************  Februar 2002  ****************************
@@ -3104 +3444 @@
     {
       Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
+/*
       idElements(Gomega, "Gw");
       headidString(Gomega, "Gw");
+*/
     }
 #endif
@@ -3128 +3470 @@
     else
     {