Changeset 47362dc in git


Ignore:
Timestamp:
Jul 23, 2015, 2:50:29 PM (8 years ago)
Author:
Stephan Oberfranz <oberfran@…>
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
53535c6cbef445ffa8464eb62c4971b17a05b3da
Parents:
38301d00726a0671f6dbc3f82923b5c9f6964c5d4132ee2e5b539d15463dcf3db09fbf9aa7c00877
Message:
update
Merge branch 'spielwiese' of github.com:Singular/Sources into spielwiese
Files:
3 deleted
42 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/classify_aeq.lib

    r38301d r47362dc  
    77AUTHORS: Faira Kanwal Janjua     fairakanwaljanjua@gmail.com
    88         Gerhard Pfister         pfister@mathematik.uni-kl.de
     9         Khawar Mehmood          khawar1073@gmail.com                         NEU
    910
    1011OVERVIEW: A library for classifying the simple singularities
     
    2223  [4] Bruce, J.W.,Gaffney, T.J.: Simple singularities of mappings (C, 0) ->(C2,0).
    2324  J. London Math. Soc. (2) 26 (1982), 465-474.
    24 
     25  [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)
    2526PROCEDURES:
    2627          sagbiAlg(G);    Compute the Sagbi-basis of the Algebra.
     
    662663
    663664////////////////////////////////////////////////////////////////////////////////
    664 proc sagbiAlg(ideal I)
     665proc sagbiAlg(ideal I,list #)
    665666"USAGE":  sagbiAlg(I);  I ideal
    666667RETURN: An ideal.The sagbi bases of I.
     
    676677
    677678    if(size(I)==0){return(I);}
    678     int b=ConductorBound(I);
    679 
    680   // int b=200;
     679    if(size(#)==0)
     680    {
     681       int b=ConductorBound(I);
     682    // int b=200;
    681683  //   b=correctBound(I,b);
     684    }
     685    else
     686    {
     687       int b=#[1];
     688    }
    682689   ideal S=interReduceSagbi(I,b) ;
    683690  // b=correctBound(S,b);
     
    699706   return(S);
    700707}
    701 
    702708example
    703709{
     
    709715sagbiAlg(I);
    710716}
    711 
     717//////////////////////////////////////////////////////////////////////////////// NEU
     718proc reducedSagbiAlg(ideal I,list #)
     719{
     720
     721   I=sagbiAlg(I,#);
     722   intvec L=semiGroup(I)[3];
     723   if(size(#)==0)
     724   {
     725      int b=findConductor(L);
     726   }
     727   else
     728   {
     729      int b=#[1];
     730   }
     731   int i;
     732   poly q;
     733   for(i=1;i<=size(I);i++)
     734   {
     735      q=I[i]-lead(I[i]);
     736      q=sagbiNF(q,I,b);
     737      I[i]=lead(I[i])+q;
     738   }
     739   return(I);
     740}
     741example
     742{
     743  "EXAMPLE:"; echo=2;
     744ring R=0,t,ds;
     745ideal I=t4+2t9,t9+t10+19/18t11-3t12+t13-t14;
     746reducedSagbiAlg(I);
     747}
     748////////////////////////////////////////////////////////////////////////////////   NEU
     749proc classifyAEQunimodal(ideal I)
     750"USAGE":  classifyAEQunimodal(I);  I ideal generated by 2 polynomials
     751RETURN: An ideal.Ideal is one of the singularity in the list of Ishikawa and Jenczko [5]
     752EXAMPLE: example classifyAEQunimodal;  shows an example
     753"
     754{
     755    def R=basering;
     756    ring @S=0,t,ds;
     757    ideal I=fetch(R,I);
     758    ideal J;
     759    poly g;
     760    if(size(I)>=3){ERROR("not a plane curve");}
     761    I=simplify(I,2);   //deletes zero’s i I
     762    I=simplify(I,1);   //creates monic generators in I
     763    if(ord(I[1])>ord(I[2])){poly q=I[2];I[2]=I[1];I[1]=q;}
     764    if((ord(I[1])>=6)||(ord(I[1])<=3)){return("not in the unimodal list");}
     765    //compute estimate of the term with the modulus
     766    int c=ord(sagbiNF(I[2],ideal(I[1]),10));
     767    if(c==10)
     768    {
     769       if(ord(I[1])!=4){return("not in the unimodal list");}
     770       c=ord(I[2][2])+2;
     771    }
     772    else
     773    {
     774       c=0;
     775       intvec v=ord(I[1]),ord(I[2]);
     776       if(v==intvec(5,6)){c=14;}
     777       if(v==intvec(5,7)){c=18;}
     778       if(v==intvec(5,8)){c=22;}
     779       if(v==intvec(4,9)){c=19;}
     780       if(v==intvec(4,11)){c=25;}
     781       if(c==0){return("not in the unimodal list");}
     782    }
     783    while(size(I[1])>1)
     784    {
     785      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);
     786    }
     787    ideal G=I;
     788    G[2]=sagbiNF(G[2],ideal(G[1]),c);
     789    ideal M=sagbiMod(diff(G,t),G);
     790    list K=semiMod(M,G);
     791
     792    if(K[1]==intvec(4,9))
     793    {
     794       if(K[4]==intvec(3,8)){J=t4,t9;}
     795       if(K[4]==intvec(3,8,22)){J=t4,t9+t19;}
     796       if(K[4]==intvec(3,8,18)){J=t4,t9+t15;}
     797       if(K[4]==intvec(3,8,14)){J=t4,t9+t11;}
     798       if(K[4]==intvec(3,8,13))
     799       {
     800          G=reducedSagbiAlg(G,15);
     801          if(ord(G[2][4])==14)
     802          {
     803             //kill the term t14 by some transformation
     804             G=subst(G,t,t-leadcoef(G[2][4])/9*t^6);
     805             G=jet(G,15);
     806             G[1]=sagbiNF(G[1],ideal(G[2]),15);
     807             //arrange the first element to be t4
     808             while(size(G[1])>1)
     809             {
     810                G=subst(G,t,t-1/4*(G[1]-lead(G[1]))/t^3);
     811                G=jet(G,15);
     812             }
     813          }
     814          G[2]=sagbiNF(G[2],ideal(G[1]),15);
     815          //arrange the coefficient of t10 to become 1
     816          number m=leadcoef(G[2]-lead(G[2]));
     817          G[2]=m^9*subst(G[2],t,1/m*t);
     818          J=G;
     819       }
     820       if(K[4]==intvec(3,8,13,18))
     821       {
     822          G=reducedSagbiAlg(G,11);
     823          number m=leadcoef(G[2]-lead(G[2]));
     824          G[2]=m^9*subst(G[2],t,1/m*t);
     825          J=G;
     826       }
     827    }
     828    if(K[1]==intvec(4,11))
     829    {
     830       if(K[4]==intvec(3,10)){J=t4,t11;}
     831       if(K[4]==intvec(3,10,28)){J=t4,t11+t25;}
     832       if(K[4]==intvec(3,10,24)){J=t4,t11+t21;}
     833       if(K[4]==intvec(3,10,20)){J=t4,t11+t17;}
     834       if(K[4]==intvec(3,10,16))
     835       {
     836          G=reducedSagbiAlg(G,14);
     837          number m=leadcoef(G[2]-lead(G[2]));
     838          number l=leadcoef(G[2][3]);
     839          //lambda^2=l^2/m^3
     840          J=G;
     841       }
     842       if(K[4]==intvec(3,10,17))
     843       {
     844          G=reducedSagbiAlg(G,21);
     845          if(ord(G[2][4])==18)
     846          {
     847             //kill the term t18 by some transformation
     848             G=subst(G,t,t-leadcoef(G[2][4])/11*t^8);
     849             G=jet(G,21);
     850             G[1]=sagbiNF(G[1],ideal(G[2]),21);
     851             //arrange the first element to be t4
     852             while(size(G[1])>1)
     853             {
     854                G=subst(G,t,t-1/4*(G[1]-lead(G[1]))/t^3);
     855                G=jet(G,21);
     856             }
     857          }
     858          G[2]=sagbiNF(G[2],ideal(G[1]),21);
     859          //arrange the coefficient of t14 to become 1
     860          number m=leadcoef(G[2]-lead(G[2]));
     861          number l=leadcoef(G[2][4]);
     862          //lambda^2=l^3/m^10
     863          J=G;
     864
     865
     866       }
     867       if(K[4]==intvec(3,10,17,24))
     868       {
     869          G=reducedSagbiAlg(G,18);
     870          //arrange the coefficient of t14 to become 1
     871          number m=leadcoef(G[2]-lead(G[2]));
     872          G[2]=t11+t14+leadcoef(G[2][3])/m^2*t17;
     873          J=G;
     874       }
     875    }
     876    if((size(K[1])==3)&&(K[1][1]==4)&&(K[1][2]==10))
     877    {
     878       int l=(K[1][3]-19) div 2;
     879       G=reducedSagbiAlg(G,2*l+12);
     880       number m=leadcoef(G[2]-lead(G[2]));
     881       number s=leadcoef(G[2][3]);
     882       //lambda^(2l-1)=s^(2l-1)/m^(2l+1)
     883       J=G;
     884    }
     885    if(K[1]==intvec(5,6))
     886    {
     887       if(K[4]==intvec(4,5)){J=t5,t6;}
     888       if(K[4]==intvec(4,5,18)){J=t5,t6+t14;}
     889       if(K[4]==intvec(4,5,13)){J=t5,t6+t9;}
     890       if(K[4]==intvec(4,5,12))
     891       {
     892          G=reducedSagbiAlg(G,9);
     893          if(ord(G[2][2])==7)
     894          {
     895             //kill the term t7 by some transformation
     896             G=subst(G,t,t-leadcoef(G[2][2])/6*t^2);
     897             G=jet(G,10);
     898             G[1]=sagbiNF(G[1],ideal(G[2]),9);
     899             //arrange the first element to be t4
     900             while(size(G[1])>1)
     901             {
     902                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
     903                G=jet(G,9);
     904             }
     905          }
     906          G[2]=sagbiNF(G[2],ideal(G[1]),9);
     907          //arrange the coefficient of t8 to become 1
     908          number m=leadcoef(G[2]-lead(G[2]));
     909          number l=leadcoef(G[2][3]);
     910          //lambda^2=l^2/m^3
     911          J=G;
     912
     913       }
     914    }
     915    if(K[1]==intvec(5,7))
     916    {
     917       if(K[4]==intvec(4,6)){J=t5,t7;}
     918       if(K[4]==intvec(4,6,22)){J=t5,t7+t18;}
     919       if(K[4]==intvec(4,6,17)){J=t5,t7+t13;}
     920       if(K[4]==intvec(4,6,12))
     921       {
     922          G=reducedSagbiAlg(G,11);
     923          if(ord(G[2][3])==9)
     924          {
     925             //kill the term t9 by some transformation
     926             G=subst(G,t,t-leadcoef(G[2][3])/7*t^3);
     927             G=jet(G,11);
     928             G[1]=sagbiNF(G[1],ideal(G[2]),11);
     929             //arrange the first element to be t4
     930             while(size(G[1])>1)
     931             {
     932                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
     933                G=jet(G,11);
     934             }
     935          }
     936          G[2]=sagbiNF(G[2],ideal(G[1]),11);
     937          //arrange the coefficient of t8 to become 1
     938          number m=leadcoef(G[2]-lead(G[2]));
     939          G[2]=m^7*subst(G[2],t,1/m*t);
     940          J=G;
     941       }
     942       if(K[4]==intvec(4,6,15))
     943       {
     944          G=reducedSagbiAlg(G,14);
     945          if(ord(G[2][2])==9)
     946          {
     947             //kill the term t9 by some transformation
     948             G=subst(G,t,t-leadcoef(G[2][2])/7*t^3);
     949             G=jet(G,14);
     950             G[1]=sagbiNF(G[1],ideal(G[2]),14);
     951             //arrange the first element to be t4
     952             while(size(G[1])>1)
     953             {
     954                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
     955                G=jet(G,14);
     956             }
     957          }
     958          G[2]=sagbiNF(G[2],ideal(G[1]),14);
     959          //arrange the coefficient of t11 to become 1
     960          number m=leadcoef(G[2]-lead(G[2]));
     961          number l=leadcoef(G[2][3]);
     962          //lambda^2=l^2/m^3
     963          J=G;
     964       }
     965
     966    }
     967    if(K[1]==intvec(5,8))
     968    {
     969       if(K[4]==intvec(4,7)){J=t5,t8;}
     970       if(K[4]==intvec(4,7,26)){J=t5,t8+t22;}
     971       if(K[4]==intvec(4,7,21)){J=t5,t8+t17;}
     972       if(K[4]==intvec(4,7,13))
     973       {
     974          G=reducedSagbiAlg(G,12);
     975          if(ord(G[2][3])==11)
     976          {
     977             //kill the term t11 by some transformation
     978             G=subst(G,t,t-leadcoef(G[2][3])/8*t^4);
     979             G=jet(G,12);
     980             G[1]=sagbiNF(G[1],ideal(G[2]),12);
     981             //arrange the first element to be t4
     982             while(size(G[1])>1)
     983             {
     984                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
     985                G=jet(G,12);
     986             }
     987          }
     988          G[2]=sagbiNF(G[2],ideal(G[1]),12);
     989          //arrange the coefficient of t9 to become 1
     990          number m=leadcoef(G[2]-lead(G[2]));
     991          G[2]=m^8*subst(G[2],t,1/m*t);
     992          J=G;
     993       }
     994
     995       if(K[4]==intvec(4,7,16))
     996       {
     997          G=reducedSagbiAlg(G,14);
     998          if(ord(G[2][2])==11)
     999          {
     1000             //kill the term t11 by some transformation
     1001             G=subst(G,t,t-leadcoef(G[2][2])/8*t^4);
     1002             G=jet(G,14);
     1003             G[1]=sagbiNF(G[1],ideal(G[2]),14);
     1004             //arrange the first element to be t4
     1005             while(size(G[1])>1)
     1006             {
     1007                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
     1008                G=jet(G,14);
     1009             }
     1010          }
     1011          G[2]=sagbiNF(G[2],ideal(G[1]),14);
     1012          //arrange the coefficient of t12 to become 1
     1013          number m=leadcoef(G[2]-lead(G[2]));
     1014          number l=leadcoef(G[2][3]);
     1015          //lambda^2=l^2/m^3
     1016          J=G;
     1017
     1018       }
     1019       if(K[4]==intvec(4,7,18))
     1020       {
     1021          G=reducedSagbiAlg(G,17);
     1022          if(ord(G[2][2])==11)
     1023          {
     1024             //kill the term t11 by some transformation
     1025             G=subst(G,t,t-leadcoef(G[2][2])/8*t^4);
     1026             G=jet(G,17);
     1027             G[1]=sagbiNF(G[1],ideal(G[2]),17);
     1028             //arrange the first element to be t4
     1029             while(size(G[1])>1)
     1030             {
     1031                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
     1032                G=jet(G,17);
     1033             }
     1034          }
     1035          G[2]=sagbiNF(G[2],ideal(G[1]),17);
     1036          //arrange the coefficient of t12 to become 1
     1037          number m=leadcoef(G[2]-lead(G[2]));
     1038          number l=leadcoef(G[2][3]);
     1039          //lambda^2=l^2/m^3
     1040          J=G;
     1041       }
     1042    }
     1043    setring R;
     1044    ideal J=fetch(@S,J);
     1045    if(size(J)==0)
     1046    {
     1047       return("not in the unimodal list");
     1048    }
     1049    return(J);
     1050}
     1051example
     1052{
     1053  "EXAMPLE:"; echo=2;
     1054ring R=0,t,ds;
     1055ideal I=t4,9t9+18t10+38t11-216t12+144t13-288t14;
     1056classifyAEQunimodal(I);
     1057I=t4,9t9+18t10+40t11-216t12+144t13-288t14;
     1058classifyAEQunimodal(I);
     1059I=t4,t11+t12+3t14+2t15+7t16+7t17;
     1060classifyAEQunimodal(I);
     1061I=t4,t11+t14+25/22t17+3t18+4t21;
     1062classifyAEQunimodal(I);
     1063I=t5,t6+2t7+t8+3t9;
     1064classifyAEQunimodal(I);
     1065I=t5,t7+3t8+3t9+5t10;
     1066classifyAEQunimodal(I);
     1067I=t5,t7+3t11+3t12+5t13;
     1068classifyAEQunimodal(I);
     1069I=t5,t8+3t9+5t10+2t11+3t12+5t13;
     1070classifyAEQunimodal(I);
     1071I=t5,t8+5t11+3t12+7t13+5t14;
     1072classifyAEQunimodal(I);
     1073I=t5,t8+5t11+7t13+5t14+7t15+2t16+8t17;
     1074classifyAEQunimodal(I);
     1075I=subst(I,t,t+t2);
     1076classifyAEQunimodal(I);
     1077I=t4+2t5+3t6+5t7+t8,9t9+18t10+40t11-216t12+144t13-288t14;
     1078classifyAEQunimodal(I);
     1079}
     1080////////////////////////////////////////////////////////////////////////////////   NEU
     1081proc computeModulus(poly p)
     1082"USAGE":  computeModulus(p);  p monic poly with 3 or 4 monomials
     1083RETURN: A polynomial with first and second coefficient 1
     1084EXAMPLE: computeModulus;  shows an example
     1085ASSUME: the basering has one vearable and one parameter
     1086"
     1087{
     1088   def R=basering;
     1089   int a1=ord(p);
     1090   int a2=ord(p-lead(p));
     1091   number m=leadcoef(p-lead(p));
     1092
     1093   poly q=par(1)^(a2-a1)-1/m;
     1094   ring S=(0,a),t,ds;
     1095   number m=fetch(R,m);
     1096   minpoly=par(1)^(a2-a1)-1/m;
     1097   poly p=fetch(R,p);
     1098   p=1/par(1)^a1*subst(p,var(1),par(1)*var(1));
     1099   setring R;
     1100   p=imap(S,p);
     1101   return(list(p,q));
     1102}
     1103example
     1104{
     1105  "EXAMPLE:"; echo=2;
     1106ring R=(0,a),t,ds;
     1107poly p=t8-395/16t14+4931/32t17;
     1108computeModulus(p);
     1109p=t8+3t12-395/16t14;
     1110computeModulus(p);
     1111p=t8-395/16t14+4931/32t17;
     1112computeModulus(p);
     1113
     1114}
     1115
     1116////////////////////////////////////////////////////////////////////////////////   NEU
     1117static proc n_thRoot(poly p,int n, int b)
     1118{
     1119   //computes the n-th root of 1+p up to order b
     1120   //assumes that p(0)=0
     1121   poly s=1;
     1122   poly q=jet(p,b);
     1123   if(q==0){return(s);}
     1124   int i;
     1125   for(i=1;i<=b;i++)
     1126   {
     1127      s=s+bino(n,i)*q;
     1128      q=jet(q*p,b);
     1129   }
     1130   return(jet(s,b));
     1131}
     1132////////////////////////////////////////////////////////////////////////////////   NEU
     1133static proc bino(number n, int i)
     1134{
     1135//computes the i-th binomial coefficient of 1/n
     1136   if(i==0){return(1);}
     1137   return(bino(n,i-1)*(1/n-i+1)/i);
     1138}
    7121139////////////////////////////////////////////////////////////////////////////////
    7131140proc sagbiMod(ideal I,ideal G)
     
    12831710    execute
    12841711    ("ring T=("+charstr(br)+",x(1),z(1..n)),(x(2),y(1..m)),dp;");
    1285     execute
    1286     ("ring R=("+charstr(br)+"),(x(1..2),y(1..m),z(1..n)),(lp(m+2),dp(n));");
     1712     execute
     1713    ("ring R=("+charstr(br)+"),(x(1..2),y(1..m),z(1..n)),(lp(2),dp(m),dp(n));");
     1714
    12871715     map phi = br,x(1);
    12881716     ideal P = phi(G1);
     
    13051733     M[size(M)+1]=check-x(2);
    13061734     check=check*keep;
     1735
    13071736     option(redSB);
    13081737     M=std(M);
    13091738     int j,s;
     1739
    13101740     for(i=1;i<=size(M);i++)
    13111741     {
     
    13421772         f=sagbiNFMODO(f,G,I,b);
    13431773    }
    1344     return(f+sagbiNFMOD(p-lead(p),G,I,b));
     1774    return(lead(f)+sagbiNFMOD(p-lead(p),G,I,b));
    13451775}
    13461776////////////////////////////////////////////////////////////////////////////////
  • Singular/LIB/gmssing.lib

    r38301d r47362dc  
    110110
    111111  M=transpose(simplify(M,2));
    112   I=M[1];
     112  I=ideal(M[1]);
    113113  attrib(I,"isSB",1);
    114114  M=M[2..ncols(M)];
  • Singular/LIB/grobcov.lib

    r38301d r47362dc  
    11//
    2 version="version grobcov.lib 4.0.1.2 Jan_2015 "; // $Id$
     2version="version grobcov.lib 4.0.2.0 Jul_2015 "; // $Id$
     3           // version L;  July_2015;
    34category="General purpose";
    45info="
     
    5455            Groebner Cover, and new theoretical developments have been done.
    5556
     57            The actual version also includes a routine (ConsLevels)
     58            for computing the canonical form of a constructible set, given as a
     59            union of locally closed sets. It is used in the new version for the
     60            computation of loci and envelops. It determines the canonical locally closed
     61            level sets of a constructuble. They will be described in a forthcoming paper:
     62
     63             J.M. Brunat, A. Montes,
     64            \"Computing the canonical representation of constructible sets\".
     65            Submited to Mathematics in Computer Science. July 2015.
     66
    5667            A new set of routines (locus, locusdg, locusto) has been included to
    5768            compute loci of points. The routines are used in the Dynamic
     
    6980             \''Envelops in Dynamic Geometry using the Groebner cover\''.
    7081
    71             The actual version also includes two routines (AddCons and AddconsP)
    72             for computing the canonical form of a constructible set, given as a
    73             union of locally closed sets. They are used in the new version for the
    74             computation of loci and envelops. It determines the canonical locally closed
    75             level sets of a constructuble. They will be described in a forthcoming paper:
    76 
    77             A. Montes, J.M. Brunat,
    78             \"Canonical representations of constructible sets\".
    79 
    80             This version was finished on 31/01/2015
     82
     83            This version was finished on 31/07/2015
    8184
    8285NOTATIONS: All given and determined polynomials and ideals are in the
     
    8588@*         grobcov, cgsdr,
    8689@*         generate the global rings
    87 @*         Grobcov::@R   (Q[a][x]),
    88 @*         Grobcov::@P   (Q[a]),
    89 @*         Grobcov::@RP  (Q[x,a])
     90@*         @R   (Q[a][x]),
     91@*         @P   (Q[a]),
     92@*         @RP  (Q[x,a])
    9093@*         that are used inside and killed before the output.
    91 @*         If you want to use some internal routine you must
    92 @*         create before the above rings by calling setglobalrings();
    93 @*         because some of the internal routines use these rings.
    94 @*         The call to the basic routines grobcov, cgsdr will
    95 @*         kill these rings.
    9694
    9795PROCEDURES:
     
    109107             (the own routine of 2010 that is no more used).
    110108             Now, KSW algorithm is used.
    111 
    112 setglobalrings();  Generates the global rings @R, @P and @PR that
    113               are respectively the rings Q[a][x], Q[a], Q[x,a].  (a=parameters,
    114              x=variables) It is called inside each of the fundamental
    115              routines of the library: grobcov, cgsdr, locus, locusdg and
    116              killed before the output.
    117109
    118110pdivi(f,F); Performs a pseudodivision of a parametric polynomial
     
    138130             the bases are computed, and one can obtain the full
    139131             representation using extend.
     132
     133ConsLevels(L); Given a list L of locally closed sets, it returns the canonical levels
     134             of the constructible set of the union of them, as well as the levels
     135             of the complement. It is described in the paper
     136
     137             J.M. Brunat, A. Montes,
     138            \"Computing the canonical representation of constructible sets\".
     139            Submited to Mathematics in Computer Science. July 2015.
    140140
    141141locus(G);    Special routine for determining the geometrical locus of points
     
    178178             \''Envelops in Dynamic Geometry using the Gr\"obner cover\''.
    179179
    180 
    181180envelopdg(ev); Is a special routine to determine the 'Relevant' components
    182181             of the envelop of a family of curves to be used in Dynamic Geometry.
    183182             It must be called to the output of envelop(F,C).
    184183
    185 locusto(L); Transforms the output of locus, locusdg, envelop and  envelopdg
     184locusto(L); Transforms the output of locus, locusdg, envelop and envelopdg
    186185             into a string that can be reed from different computational systems.
    187186
    188 AddCons(L); Uses the routine AddConsP. Given a set of locally closed sets as
    189              difference of of varieties (does not need to be in C-representation)
    190              it returns the canonical P-representation of the constructible set
    191              formed by the union of them. The result is formed by a set of embedded,
    192              disjoint, locally closed sets (levels).
    193 
    194 AddConsP(L);  Given a set of locally closed P-components, it returns the
    195              canonical P-representation of the constructible set
    196              formed by the union of them, consisting of a set of embedded,
    197              disjoint locally closed sets (levels).
    198              The routines AddCons and AddConsP and the canonical structure
    199              of the constructible sets will be described in a forthcoming paper.
    200 
    201              A. Montes, J.M. Brunat,
    202              \"Canonical representations of constructible sets\".
    203187
    204188SEE ALSO: compregb_lib
     
    228212// Uses KSW algorithm for cgsdr
    229213// Final data: 21-11-2013
    230 // Release K: (public)
    231 // Updated locus: 8-1-2015
    232 // Updated AddConsP and AddCons: 8-1-2015
    233 // Reformed many routines, examples and helps: 8-1-2015
     214// Release L: (public)
     215// New routine ConsLevels: 10-7-2015
     216// Updated locus: 10-7-2015 (uses Conslevels)
    234217// New routines for computing the envelop of a family of curves: 22-1-2015
    235 // Final data: 22-1-2015
     218// Final data: 22-7-2015
    236219
    237220//*************Auxiliary routines**************
     
    297280  if (size(l)==1 and i==1){return(L);}
    298281  // L=l[1];
    299   if(i==1)
    300   {
    301     for(j=2;j<=size(l);j++)
    302     {
    303       L[j-1]=l[j];
    304     }
    305   }
    306   else
     282  if(i>1)
    307283  {
    308284    for(j=1;j<=i-1;j++)
    309285    {
    310       L[j]=l[j];
    311     }
    312     for(j=i+1;j<=size(l);j++)
    313     {
    314       L[j-1]=l[j];
    315     }
     286      L[size(L)+1]=l[j];
     287    }
     288  }
     289  for(j=i+1;j<=size(l);j++)
     290  {
     291    L[size(L)+1]=l[j];
    316292  }
    317293  return(L);
     
    773749//}
    774750
    775 proc setglobalrings()
    776 "USAGE:   setglobalrings();
    777           No arguments
    778 RETURN:   After its call the rings @R=Q[a][x], @P=Q[a], @RP=Q[x,a] are
    779           defined as global variables.  (a=parameters, x=variables)
    780 NOTE: It is called internally by many basic routines of the
    781           library grobcov, cgsdr, extend, pdivi, pnormalf, locus, locusdg,
    782           envelop, envelopdg, and killed before the output.
    783           The user does not need to call it, except when it is interested
    784           in using some internal routine of the library that
    785           uses these rings.
    786           The basering R, must be of the form Q[a][x], (a=parameters,
    787           x=variables), and should be defined previously.
    788 KEYWORDS: ring, rings
    789 EXAMPLE:  setglobalrings; shows an example"
     751static proc setglobalrings()
     752// "USAGE:   setglobalrings();
     753//           No arguments
     754// RETURN: After its call the rings Grobcov::@R=Q[a][x], Grobcov::@P=Q[a],
     755//           Grobcov::@RP=Q[x,a] are defined as global variables.
     756//           (a=parameters, x=variables)
     757// NOTE: It is called internally by many basic routines of the
     758//           library grobcov, cgsdr, extend, pdivi, pnormalf, locus, locusdg,
     759//           envelop, envelopdg, and killed before the output.
     760//           The user does not need to call it, except when it is interested
     761//           in using some internal routine of the library that
     762//           uses these rings.
     763//           The basering R, must be of the form Q[a][x], (a=parameters,
     764//           x=variables), and should be defined previously.
     765// KEYWORDS: ring, rings
     766// EXAMPLE:  setglobalrings; shows an example"
    790767{
    791768  if (defined(@P))
     
    810787  setring(RR);
    811788};
    812 example
    813 {
    814   "EXAMPLE:"; echo = 2;
    815   ring R=(0,a,b),(x,y,z),dp;
    816   setglobalrings();
    817   R;
    818   Grobcov::@R;
    819   Grobcov::@P;
    820   Grobcov::@RP;
    821   ringlist(Grobcov::@R);
    822   ringlist(Grobcov::@P);
    823   ringlist(Grobcov::@RP);
    824 }
     789// example
     790// {
     791//   "EXAMPLE:"; echo = 2;
     792//   ring R=(0,a,b),(x,y,z),dp;
     793//   setglobalrings();
     794//   " ";"R=";R;
     795//   " ";"Grobcov::@R=";Grobcov::@R;
     796//   " ";"Grobcov::@P=";Grobcov::@P;
     797//   " ";"Grobcov::@RP=";Grobcov::@RP;
     798//  " "; "ringlist(Grobcov::@R)=";  ringlist(Grobcov::@R);
     799//  " "; "ringlist(Grobcov::@P)=";  ringlist(Grobcov::@P);
     800//  " "; "ringlist(Grobcov::@RP)=";  ringlist(Grobcov::@RP);
     801// }
    825802
    826803// cld : clears denominators of an ideal and normalizes to content 1
     
    15151492    }
    15161493  }
     1494  //"T_abans="; prep;
    15171495  if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));}
     1496  //"T_Prep="; prep;
     1497  //def Lout=CompleteA(prep,prep);
     1498  //"T_Lout="; Lout;
    15181499  return(prep);
    15191500}
     
    39343915
    39353916//********************* End KapurSunWang *************************
     3917
     3918//********************* Begin ConsLevels ***************************
     3919
     3920static proc zeroone(int n)
     3921{
     3922  list L; list L2;
     3923  intvec e; intvec e2; intvec e3;
     3924  int j;
     3925  if(n==1)
     3926  {
     3927    e[1]=0;
     3928    L[1]=e;
     3929    e[1]=1;
     3930    L[2]=e;
     3931    return(L);
     3932  }
     3933  if(n>1)
     3934  {
     3935    L=zeroone(n-1);
     3936    for(j=1;j<=size(L);j++)
     3937    {
     3938      e2=L[j];
     3939      e3=e2;
     3940      e3[size(e3)+1]=0;
     3941      L2[size(L2)+1]=e3;
     3942      e3=e2;
     3943      e3[size(e3)+1]=1;
     3944      L2[size(L2)+1]=e3;
     3945    }
     3946  }
     3947  return(L2);
     3948}
     3949
     3950// Auxiliary routine
     3951// subsets: the list of subsets of (1,..n)
     3952static proc subsets(int n)
     3953{
     3954  list L; list L1;
     3955  int i; int j;
     3956  L=zeroone(n);
     3957  intvec e; intvec e1;
     3958  for(i=1;i<=size(L);i++)
     3959  {
     3960    e=L[i];
     3961    kill e1; intvec e1;
     3962    for(j=1;j<=n;j++)
     3963    {
     3964      if(e[n+1-j]==1)
     3965      {
     3966        if(e1==intvec(0)){e1[1]=j;}
     3967        else{e1[size(e1)+1]=j};
     3968      }
     3969    }
     3970    L1[i]=e1;
     3971  }
     3972  return(L1);
     3973}
     3974
     3975// Input a list A of locally closed sets in C-rep
     3976// Output a list B of a simplified list of A
     3977static proc SimplifyUnion(list A)
     3978{
     3979  int i; int j;
     3980  list L=A;
     3981  int n=size(L);
     3982  if(n<2){return(A);}
     3983  for(i=1;i<=size(L);i++)
     3984  {
     3985    for(j=1;j<=size(L);j++)
     3986    {
     3987      if(i != j)
     3988      {
     3989        if(equalideals(L[i][2],L[j][1])==1)
     3990        {
     3991          L[i][2]=L[j][2];
     3992        }
     3993      }
     3994    }
     3995  }
     3996  ideal T=ideal(1);
     3997  intvec v;
     3998  for(i=1;i<=size(L);i++)
     3999  {
     4000    if(equalideals(L[i][2],ideal(1)))
     4001    {
     4002      v[size(v)+1]=i;
     4003      T=intersect(T,L[i][1]);
     4004    }
     4005  }
     4006  if(size(v)>0)
     4007  {
     4008    for(i=1; i<=size(v);i++)
     4009    {
     4010      j=v[size(v)+1-i];
     4011      L=elimfromlist(L, j);
     4012    }
     4013  }
     4014  if(equalideals(T,ideal(1))==0){L[size(L)+1]=list(std(T),ideal(1))};
     4015  //string("T_n=",n," new n0",size(L));
     4016  return(L);
     4017}
     4018
     4019// Input: list(A)
     4020//          A is a list of locally closed sets in Crep. A=[[P1,Q1],[P2,Q2],..,[Pr,Qr]]
     4021//          whose union is a constructible set from
     4022// Output list [Lev,C]:
     4023//          where Lev is the Crep of the canonical first level of A, and
     4024//          C is the complement of the first level Lev wrt to the closure of A. The elements are given in Crep,
     4025static proc FirstLevel(list A)
     4026{
     4027  int n=size(A);
     4028  list T=zeroone(n);
     4029  ideal P; ideal Q;
     4030  list Cb;  ideal Cc=ideal(1);
     4031  int i; int j;
     4032  intvec t;
     4033  ideal Z=ideal(1);
     4034  list C;
     4035  for(i=1;i<=size(A);i++)
     4036  {
     4037    Z=intersect(Z,A[i][1]);
     4038  }
     4039  for(i=2; i<=size(T);i++)
     4040  {
     4041    t=T[i];
     4042    P=ideal(1); Q=ideal(0);
     4043    for(j=1;j<=n;j++)
     4044    {
     4045      if(t[n+1-j]==1)
     4046      {
     4047        Q=Q+A[j][2];
     4048      }
     4049      else
     4050      {
     4051        P=intersect(P,A[j][1]);
     4052      }
     4053    }
     4054    //"T_Q="; Q; "T_P="; P;
     4055    Cb=Crep(Q,P);
     4056    //"T_Cb="; Cb;
     4057    if(Cb[1][1]<>1)
     4058    {
     4059      C[size(C)+1]=Cb;
     4060      Cc=intersect(Cc,Cb[1]);
     4061    }
     4062  }
     4063  list Lev=list(Z,Cc);                //Crep(Z,Cc);
     4064  if(size(C)>1){C=SimplifyUnion(C);}
     4065  return(list(Lev,C));
     4066}
     4067
     4068// Input: list(A)
     4069//          A is a list of locally closed sets in Crep. A=[[P1,Q1],[P2,Q2],..,[Pr,Qr]]
     4070//          whose union is a constructible set from which the algorithm computes its canonical form.
     4071// Output:
     4072//     list [L,C] where
     4073//          where L is the list of canonical levels of A,
     4074//          and C is the list of canonical levels of the complement of A wrt to the closure of A.
     4075//          All locally closed sets are given in Crep.
     4076//          L=[[1,[p1,p2],[3,[p3,p4],..,[2r-1,[p_{2r-1},p_2r]]]] is the list of levels of A in Crep.
     4077//          C=[[2,p2,p3],[4,[p4,p5],..,[2s,[p_{2s},p_{2s+1}]]]  is the list of levels of C,
     4078//                                              the complement of S wrt the closure of A, in Crep.
     4079proc ConsLevels(list A)
     4080"USAGE:   ConsLevels(A);
     4081          A is a list of locally closed sets in Crep. A=[[P1,Q1],[P2,Q2],..,[Pr,Qr]]
     4082          whose union is a constructible set from which the algorithm computes its
     4083          canonical form.
     4084RETURN: A list with two elements: [L,C]
     4085          where L is the list of canonical levels of A, and C is the list of canonical
     4086          levels of the
     4087          complement of A wrt to the closure of A.
     4088          All locally closed sets are given in Crep.
     4089          L=[[1,[p1,p2],[3,[p3,p4],..,[2r-1,[p_{2r-1},p_2r]]]]
     4090          C=[[2,p2,p3],[4,[p4,p5],..,[2s,[p_{2s},p_{2s+1}]]]
     4091NOTE: The basering R, must be of the form Q[a]
     4092KEYWORDS: locally closed set, constructible set
     4093EXAMPLE:  ConsLevels; shows an example"
     4094{
     4095  list L; list C; int i;
     4096  list B; list T;
     4097  for(i=1; i<=size(A);i++)
     4098  {
     4099    T=Crep(A[i][1],A[i][2]);
     4100    B[size(B)+1]=T;
     4101  }
     4102  int level=0;
     4103  list K;
     4104  while(size(B)>0)
     4105  {
     4106    level++;
     4107    K=FirstLevel(B);
     4108    if(level mod 2==1){L[size(L)+1]=list(level,K[1]);}
     4109    else{C[size(C)+1]=list(level,K[1]);}
     4110    //"T_L="; L;
     4111    //"T_C="; C;
     4112    B=K[2];
     4113    //"T_size(B)="; size(B);
     4114  }
     4115  return(list(L,C));
     4116}
     4117example
     4118{ "EXAMPLE:"; echo = 2;
     4119// Example of ConsLevels with nice geometrical interpretetion and 2 levels
     4120
     4121if (defined(R)){kill R;}
     4122if (defined(@P)){kill @P; kill @R; kill @RP;}
     4123
     4124  ring R=0,(x,y,z),lp;
     4125  short=0;
     4126  ideal P1=x*(x^2+y^2+z^2-1);
     4127  ideal Q1=z,x^2+y^2-1;
     4128  ideal P2=y,x^2+z^2-1;
     4129  ideal Q2=z*(z+1),y,x*(x+1);
     4130
     4131  list Cr1=Crep(P1,Q1);
     4132  list Cr2=Crep(P2,Q2);
     4133
     4134  list L=list(Cr1,Cr2);
     4135  L;
     4136  // ConsLevels(L)=
     4137  ConsLevels(L);
     4138
     4139//----------------------------------------------------------------------
     4140// Begin Problem S93
     4141// Automatic theorem proving
     4142// Generalized Steiner-Lehmus Theorem
     4143// A.Montes, T.Recio
     4144
     4145// Bisector of A(-1,0) = Bisector of B(1,0) varying C(a,b)
     4146
     4147if(defined(R1)){kill R1;}
     4148ring R1=(0,a,b),(x1,y1,x2,y2,p,r),dp;
     4149ideal S93=(a+1)^2+b^2-(p+1)^2,
     4150                    (a-1)^2+b^2-(1-r)^2,
     4151                    a*y1-b*x1-y1+b,
     4152                    a*y2-b*x2+y2-b,
     4153                    -2*y1+b*x1-(a+p)*y1+b,
     4154                    2*y2+b*x2-(a+r)*y2-b,
     4155                    (x1+1)^2+y1^2-(x2-1)^2-y2^2;
     4156
     4157short=0;
     4158def GC93=grobcov(S93,"nonnull",ideal(b),"rep",1);
     4159//"grobcov(S93,'nonnull',ideal(b),"rep",1)="; GC93;
     4160
     4161list L0;
     4162for(int i=1;i<=size(GC93);i++)
     4163{
     4164   L0[size(L0)+1]=GC93[i][3];
     4165}
     4166
     4167def GC93ab=grobcov(S93,"nonnull",ideal(a*b),"rep",1);
     4168
     4169ring RR=0,(a,b),lp;
     4170
     4171list L1;
     4172L1=imap(R1,L0);
     4173// L1=Total elements of the grobcov(S93,'nonnull',ideal(b),'rep',1) to be added=
     4174L1;
     4175
     4176// Adding all the elements of grobcov(S93,'nonnull',ideal(b),'rep',1)=
     4177ConsLevels(L1);
     4178}
     4179
     4180//**************************** End ConsLevels ******************
    39364181
    39374182//******************** Begin locus ******************************
     
    45224767  locus(grobcov(S));
    45234768  kill R;
    4524   "********************************************";
     4769  //********************************************
    45254770
    45264771  ring R=(0,x,y),(x1,x2),dp;
     
    46114856  locusdg(locus(grobcov(S96)));
    46124857  kill R;
    4613   "********************************************";
     4858  //********************************************
    46144859  ring R=(0,a,b),(x4,x3,x2,x1),dp;
    46154860  ideal S=(x1-3)^2+(x2-1)^2-9,
     
    46224867  locusdg(locus(grobcov(S)));
    46234868  kill R;
    4624   "********************************************";
     4869  //********************************************
    46254870
    46264871  ring R=(0,x,y),(x1,x2),dp;
     
    46334878}
    46344879
    4635 // locusto: Transforms the output of locus to a string that
    4636 //      can be read by different computational systems.
     4880// locusto: Transforms the output of locus, locusdg, envelop and envelopdg
     4881//             into a string that can be reed from different computational systems.
    46374882// input:
    46384883//     list L: The output of locus
     
    47174962  locusto(locusdg(locus(grobcov(S))));
    47184963  kill R;
    4719   "********************************************";
     4964  //********************************************
    47204965
    47214966  // 1. Take a fixed line l: x1-y1=0  and consider
     
    47374982  locusto(envelopdg(envelop(F,C)));
    47384983  kill R;
    4739   "********************************************";
     4984  //********************************************
    47404985
    47414986  // Steiner Deltoid
     
    48045049  return(d);
    48055050}
    4806 
    4807 // // locusdgto: Transforms the output of locusdg to a string that
    4808 // //      can be read by different computational systems.
    4809 // // input:
    4810 // //     list L: The output of locus
    4811 // // output:
    4812 // //     string s: The output of locus converted to a string readable by other programs
    4813 // //                   Outputs only the relevant dynamical geometry components.
    4814 // //                   Without unnecessary parenthesis
    4815 // proc locusdgto(list LL)
    4816 // "USAGE: locusdgto(L);
    4817 //           The argument must be the output of locusdg of a parametrical ideal
    4818 //           It transforms the output into a string in standard form
    4819 //           readable in many languages (Geogebra).
    4820 // RETURN: The locusdg in string standard form
    4821 // NOTE: It can only be called after computing the locusdg(grobcov(F)) of the
    4822 //           parametrical ideal.
    4823 //           The basering R, must be of the form Q[a,b,..][x,y,..].
    4824 // KEYWORDS: geometrical locus, locus, loci.
    4825 // EXAMPLE:  locusdgto; shows an example"
    4826 // {
    4827 //   int i; int j; int k; int short0=short; int ojo;
    4828 //   int te=0;
    4829 //   short=0;
    4830 //   if(size(LL)==0){ojo=1; list L;}
    4831 //   else
    4832 //   {
    4833 //     def L=LL;
    4834 //   }
    4835 //   string s="["; string sf="]"; string st=s+sf;
    4836 //   if(size(L)==0){return(st);}
    4837 //   ideal p;
    4838 //   ideal q;
    4839 //   for(i=1;i<=size(L);i++)
    4840 //   {
    4841 //     if(L[i][3]=="Relevant")
    4842 //     {
    4843 //       s=string(s,"[[");
    4844 //       for (j=1;j<=size(L[i][1]);j++)
    4845 //       {
    4846 //         s=string(s,L[i][1][j],",");
    4847 //       }
    4848 //       s[size(s)]="]";
    4849 //       s=string(s,",[");
    4850 //       for(j=1;j<=size(L[i][2]);j++)
    4851 //       {
    4852 //         s=string(s,"[");
    4853 //         for(k=1;k<=size(L[i][2][j]);k++)
    4854 //         {
    4855 //           s=string(s,L[i][2][j][k],",");
    4856 //         }
    4857 //         s[size(s)]="]";
    4858 //         s=string(s,",");
    4859 //       }
    4860 //       s[size(s)]="]";
    4861 //       s=string(s,"]");
    4862 //       s[size(s)]="]";
    4863 //       s=string(s,",");
    4864 //     }
    4865 //   }
    4866 //   if(s=="["){s="[]";}
    4867 //   else{s[size(s)]="]";}
    4868 //   short=short0;
    4869 //   return(s);
    4870 // }
    4871 // example
    4872 // {"EXAMPLE:"; echo = 2;
    4873 //   ring R=(0,a,b),(x,y),dp;
    4874 //   short=0;
    4875 //   ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
    4876 //   "System="; S96; " ";
    4877 //   "locusdgto(locusdg(grobcov(S96)))=";
    4878 //   locusdgto(locusdg(grobcov(S96)));
    4879 // }
    48805051
    48815052static proc norspec(ideal F)
     
    48985069  exportto(Top,@RP);     // global ring K[x,a] with product order
    48995070  setring(RR);
    4900 
    49015071}
    49025072
     
    50595229}
    50605230
    5061 //********************* End locus ****************************
    5062 
    5063 //********************* Begin AddCons **********************
    5064 
    5065 // Input: L1,L2: lists of components with common top
    5066 // Output L: list of the union of L1 and L2.
    5067 static proc Add2ComWithCommonTop(list L1, list L2)
    5068 {
    5069   int i; int j; ideal pij; list L; list Lp; list PR; int k;
    5070   for(i=1;i<=size(L1[2]);i++)
    5071   {
    5072     for(j=1;j<=size(L2[2]);j++)
    5073     {
    5074       pij=std(L1[2][i]+L2[2][j]);
    5075       PR=minGTZ(pij);
    5076       for(k=1;k<=size(PR);k++)
    5077       {
    5078         Lp[size(Lp)+1]=PR[k];
    5079       }
    5080     }
    5081   }
    5082   for(i=1; i<=size(Lp);i++)
    5083   {
    5084     for(j=i+1;j<=size(Lp);j++)
    5085     {
    5086       if(idcontains(Lp[i],Lp[j])) {Lp=delete(Lp,j);}
    5087     }
    5088     for(j=1;j<i;j++)
    5089     {
    5090       if(idcontains(Lp[i],Lp[j])){Lp=delete(Lp,j); j=j-1; i=i-1;}
    5091     }
    5092   }
    5093   L[1]=L1[1];
    5094   L[2]=Lp;
    5095   return(L);
    5096 }
    5097 
    5098 // Input: L list od P-rep of components to be added. L[i]=[p_i,[p_{i1},...p_{ir_i}]]
    5099 // Output: lists A,B,L
    5100 //       where no top in the lists are repeated
    5101 //       A: list of components with higher tops
    5102 //       B: list of components with lower tops
    5103 //       L1: A,B
    5104 static proc SepareAB(list L)
    5105 {
    5106   int i;  int j; list Ln=L;
    5107   for(i=1;i<=size(Ln);i++)
    5108   {
    5109     for(j=i+1;j<=size(Ln);j++)
    5110     {
    5111       if (equalideals(Ln[j][1],Ln[i][1]))
    5112       {
    5113         Ln[i]=Add2ComWithCommonTop(Ln[i],Ln[j]);
    5114         Ln=delete(Ln,j);
    5115         j=j-1;
    5116       }
    5117     }
    5118   }
    5119   list A; list B; int clas; list T; list L1;
    5120   for(i=1;i<=size(Ln);i++)
    5121   {
    5122     j=1;
    5123     clas=0;
    5124     while((clas==0) and  (j<=size(Ln)))
    5125     {
    5126       if(j!=i)
    5127       {
    5128         if(idcontains(Ln[i][1],Ln[j][1]) ) {B[size(B)+1]=Ln[i]; clas=1;}
    5129       }
    5130       j++;
    5131     }
    5132     if(clas==0) {A[size(A)+1]=Ln[i];}
    5133   }
    5134   L1=A; for(j=1;j<=size(B);j++){L1[size(L1)+1]=B[j];}
    5135   T[1]=A; T[2]=B; T[3]=L1;
    5136   return(T);
    5137 }
    5138 
    5139 // Input:
    5140 //  A1: list of high set of P-reps to be completed by the remaining P-reps
    5141 //  L1: the list A1, completed with the list of lower P-reps.
    5142 // Output:
    5143 //  A: list A1 completed with all possible parts of the remaining parts of L1 with the
    5144 //      condition of building locally closed sets.
    5145 static proc CompleteA(list A1,list L1)
    5146 {
    5147   int modif; int i; int j; int k; int l;
    5148   ideal pij; ideal pk; ideal pijkl; ideal pkl;
    5149   int n; list nl; int te; int clas; list vvv; list AAA;
    5150   list Lp; int m; ideal Pij;
    5151   list A0;
    5152   modif=1;
    5153   list A=A1;
    5154   while(modif==1)
    5155   {
    5156       modif=0;
    5157       A0=A;
    5158       for(i=1;i<=size(A);i++)
    5159       {
    5160           for(j=1;j<=size(A[i][2]); j++)
    5161           {
    5162               pij=A[i][2][j];
    5163              for(k=1;k<=size(L1);k++)
    5164              {
    5165                  if(k!=i)
    5166                  {
    5167                      pk=L1[k][1];
    5168                      if(idcontains(pij,pk)==1)
    5169                      {
    5170                          te=0;
    5171                          kill nl;
    5172                          list nl; l=1;
    5173                          while((te==0) and (l<=size(L1[k][2])))
    5174                          {
    5175                               pkl=L1[k][2][l];
    5176                               if((equalideals(pij,pkl)==1) or (idcontains(pij,pkl)==1)) {te=1;}
    5177                               l++;
    5178                          }   // end while ((te=0) and (l>...
    5179                          //"T_te="; te; pause();
    5180                          if(te==0)
    5181                          {
    5182                            modif=1;
    5183                            kill Pij; ideal Pij=1;
    5184                            for(l=1; l<=size(L1[k][2]);l++)
    5185                            {
    5186                              pkl=L1[k][2][l];
    5187                              pijkl=std(pij+pkl);
    5188                              Pij=intersect(Pij,pijkl);
    5189                            }
    5190                            kill Lp; list Lp;
    5191                            Lp=minGTZ(Pij);
    5192                            for(m=1;m<=size(Lp);m++)
    5193                             {
    5194                                nl[size(nl)+1]=Lp[m];
    5195                             }
    5196                             A[i][2]=delete(A[i][2], j);
    5197                             for(n=1;n<=size(nl);n++)
    5198                             {
    5199                               A[i][2][size(A[i][2])+1]=nl[n];
    5200                             }
    5201                           } // end if(te==0)
    5202                      } // end if(idcontains(pij,pk)==1)
    5203                  }  // end if (k!=i)
    5204              } // end for k
    5205          }  // end for j
    5206          kill vvv; list vvv;
    5207          if(modif==1)
    5208          // Select the maximal ideals of the set A[I][2][j]
    5209          {
    5210              kill nl; list nl;
    5211              nl=selectminideals(A[i][2]);
    5212              kill AAA; list AAA;
    5213              for(j=1;j<=size(nl);j++)
    5214              {
    5215                AAA[size(AAA)+1]=A[i][2][nl[j]];
    5216              }
    5217              A[i][2]=AAA;
    5218          } // end if(modif=1)
    5219       }  // end for i
    5220       modif=1-equallistsAall(A,A0);
    5221   } // end while(modif==1)
    5222   return(A);
    5223 }
    5224 
    5225 // Input:
    5226 //   A: list of the top P-reps of one level
    5227 //   B: list of remaining lower P-reps that have not yeen be possible to add
    5228 // Output:
    5229 //   Bn: list B where the elements that are completely included in A are removed and the parts that are
    5230 //         included in A also.
    5231 static proc ReduceB(list A,list B)
    5232 {
    5233      int i; int j; list nl; list Bn; int te; int k; int elim;
    5234      ideal pC; ideal pD; ideal pCD; ideal pBC; list nB; int j0;
    5235      for(i=1;i<=size(B);i++)
    5236      {
    5237          j=1; te=0;
    5238          while((te==0) and (j<=size(A)))
    5239          {
    5240              if(idcontains(B[i][1],A[j][1])){te=1; j0=j;}
    5241              else{j++;}
    5242          }
    5243          pD=B[i][2][1];
    5244          for(k=2;k<=size(B[i][2]);k++){pD=intersect(pD,B[i][2][k]);}
    5245          pC=A[j0][2][1];
    5246          for(k=2;k<=size(A[j0][2]);k++) {pC=intersect(pC,A[j0][2][k]);}
    5247          pCD=std(pD+pC);
    5248          pBC=std(B[i][1]+pC);
    5249          elim=0;
    5250          if(idcontains(pBC,pCD)){elim=1;}   // B=delfromlist(B,i);i=i-1;
    5251          else
    5252          {
    5253               nB=Prep(pBC,pCD);
    5254               if(equalideals(nB[1][1],ideal(1))==0)
    5255               {
    5256                   for(k=1;k<=size(nB);k++){Bn[size(Bn)+1]=nB[k];}
    5257               }
    5258          }
    5259     }   // end for i
    5260     return(Bn);
    5261 }
    5262 
    5263 // AddConsP: given a set of components of locally closed sets in P-representation, it builds the
    5264 //       canonical P-representation of the corresponding constructible set, of its union,
    5265 //       including levels it they are.
    5266 proc AddConsP(list L)
    5267 "USAGE:   AddConsP(L)
    5268       Input L: list of components in P-rep to be added
    5269       [  [[p_1,[p_11,..,p_1,r1]],..[p_k,[p_k1,..,p_kr_k]]  ]
    5270 RETURN:
    5271      list of lists of levels of the different locally closed sets of
    5272      the canonical P-rep of the constructible.
    5273      [  [level_1,[ [Comp_11,..Comp_1r_1] ] ], .. ,
    5274         [level_s,[ [Comp_s1,..Comp_sr_1] ]
    5275      ]
    5276      where level_i=i,   Comp_ij=[ p_i,[p_i1,..,p_it_i] ] is a prime component.
    5277 NOTE:     Operates in a ring R=Q[u_1,..,u_m]
    5278 KEYWORDS: Constructible set, Canoncial form
    5279 EXAMPLE:  AddConsP; shows an example"
    5280 {
    5281   list LL; list A; list B; list L1; list T; int level=0; list h;
    5282   LL=L; int i;
    5283   while(size(LL)!=0)
    5284   {
    5285     level++;
    5286     L1=SepareAB(LL);
    5287     A=L1[1]; B=L1[2]; LL=L1[3];
    5288     A=CompleteA(A,LL);
    5289     for(i=1;i<=size(A);i++)
    5290     {
    5291       LL[i]=A[i];
    5292     }
    5293     h[1]=level; h[2]=A;
    5294     T[size(T)+1]=h;
    5295     LL=ReduceB(A,B);
    5296   }
    5297   return(T);
    5298 }
    5299 example
    5300 {
    5301   "EXAMPLE:"; echo = 2;
    5302   if (defined(Grobcov::@P)){kill Grobcov::@P; kill Grobcov::@R; kill Grobcov::@RP;}
    5303   ring R=0,(x,y,z),lp;
    5304   short=0;
    5305 
    5306   ideal P1=x;
    5307   ideal Q1=x,y;
    5308   ideal P2=y;
    5309   ideal Q2=y,z;
    5310 
    5311   list L=list(Prep(P1,Q1)[1],Prep(P2,Q2)[1]);
    5312   L;
    5313   AddConsP(L);
    5314 }
    5315 
    5316 // AddCons:  given a set of  locally closed sets by pairs of ideal, it builds the
    5317 //       canonical P-representation of the corresponding constructible set, of its union,
    5318 //       including levels it they are.
    5319 //       It makes a call to AddConsP after transforming the input.
    5320 // Input list of lists of pairs of ideals representing locally colsed sets:
    5321 //     L=  [ [E1,N1], .. , [E_s,N_s] ]
    5322 // Output: The canonical frepresentation of the constructible set union of the V(E_i) \ V(N_i)
    5323 //     T=[  [level_1,[ [p_1,[p_11,..,p_1s_1]],.., [p_k,[p_k1,..,p_ks_k]] ]],, .. , [level_r,[..       ]]  ]
    5324 proc AddCons(list L)
    5325 "USAGE:   AddCons(L)
    5326       Calls internally AddConsP and allows a different input.
    5327       Input L: list of pairs of of ideals [ [P_1,Q_1], .., [Pr,Qr] ]
    5328       representing a set of locally closed setsV(P_i) \ V(Q_i)
    5329       to be added.
    5330 RETURN:
    5331       list of lists of levels of the different locally closed sets of
    5332       the canonical P-rep of the constructible.
    5333       [  [level_1,[ [Comp_11,..Comp_1r_1] ] ], .. ,
    5334          [level_s,[ [Comp_s1,..Comp_sr_1] ]
    5335       ]
    5336       where level_i=i,   Comp_ij=[ p_i,[p_i1,..,p_it_i] ] is a prime component.
    5337 NOTE: Operates in a ring R=Q[u_1,..,u_m]
    5338 KEYWORDS: Constructible set, Canoncial form
    5339 EXAMPLE: AddCons; shows an example"
    5340 {
    5341   int i; list H; list P; int j;
    5342   for(i=1;i<=size(L);i++)
    5343   {
    5344     P=Prep(L[i][1],L[i][2]);
    5345     for(j=1;j<=size(P);j++)
    5346     {
    5347       H[size(H)+1]=P[j];
    5348     }
    5349   }
    5350   return(AddConsP(H));
    5351 }
    5352 example
    5353 {
    5354   "EXAMPLE:"; echo = 2;
    5355   if (defined(Grobcov::@P)){kill Grobcov::@P; kill Grobcov::@R; kill Grobcov::@RP;}
    5356   ring R=0,(x,y,z),lp;
    5357   short=0;
    5358 
    5359   ideal P1=x;
    5360   ideal Q1=x,y;
    5361   ideal P2=y;
    5362   ideal Q2=y,z;
    5363 
    5364   list L=list(list(P1,Q1), list(P2,Q2) );
    5365   L;
    5366 
    5367   AddCons(L);
    5368 }
    5369 
    53705231// AddLocus: auxilliary routine for locus0 that computes the components of the constructible:
    53715232// Input:  the list of locally closed sets to be added, each with its type as third argument
     
    53975258    }
    53985259  }
    5399   L3=AddConsP(L1);
     5260  L3=LocusConsLevels(L1);
    54005261  list L4; int level;
    54015262  ideal p1; ideal pp1; int t; int k; int k0; string typ; list l4;
     
    54265287}
    54275288
    5428 //********************* End AddCons **********************
    5429 ;
     5289// Input L: list of components in P-rep to be added
     5290//         [  [[p_1,[p_11,..,p_1,r1]],..[p_k,[p_k1,..,p_kr_k]]  ]
     5291// Output:
     5292//          list of lists of levels of the different locally closed sets of
     5293//          the canonical P-rep of the constructible.
     5294//          [  [level_1,[ [Comp_11,..Comp_1r_1] ] ], .. ,
     5295//             [level_s,[ [Comp_s1,..Comp_sr_1] ]
     5296//          ]
     5297//          where level_i=i,   Comp_ij=[ p_i,[p_i1,..,p_it_i] ] is a prime component.
     5298// LocusConsLevels: given a set of components of locally closed sets in P-representation, it builds the
     5299//       canonical P-representation of the corresponding constructible set of its union,
     5300//       including levels it they are.
     5301static proc LocusConsLevels(list L)
     5302{
     5303  list Lc; list Sc;
     5304  int i;
     5305  for(i=1;i<=size(L);i++)
     5306  {
     5307    Sc=PtoCrep(list(L[i]));
     5308    Lc[size(Lc)+1]=Sc;
     5309  }
     5310  list S=ConsLevels(Lc)[1];
     5311  list Sout;
     5312  list Lev;
     5313  for(i=1;i<=size(S);i++)
     5314  {
     5315    Lev=list(i,Prep(S[i][2][1],S[i][2][2]));
     5316    Sout[size(Sout)+1]=Lev;
     5317  }
     5318  return(Sout);
     5319}
     5320
     5321//********************* End locus ****************************
     5322
  • Singular/LIB/grwalk.lib

    r4132ee r47362dc  
    255255}
    256256
    257 proc gwalk(ideal Go, list #)
    258 "SYNTAX: gwalk(ideal i);
    259          gwalk(ideal i, intvec v, intvec w);
     257proc gwalk(ideal Go, int reduction,int printout, list #)
     258"SYNTAX: gwalk(ideal i, int reduction, int printout);
     259         gwalk(ideal i, int reduction, int printout, intvec v, intvec w);
    260260TYPE:    ideal
    261261PURPOSE: compute the standard basis of the ideal, calculated via
     
    274274   string ord_str =   OSCTW[2];
    275275   intvec curr_weight   =   OSCTW[3]; /* original weight vector */
    276    intvec target_weight =   OSCTW[4]; /* terget weight vector */
     276   intvec target_weight =   OSCTW[4]; /* target weight vector */
    277277   kill OSCTW;
    278278   option(redSB);
     
    284284   //print("//** help ring = " + string(basering));
    285285   ideal G = fetch(xR, Go);
    286    G = system("Mwalk", G, curr_weight, target_weight,basering);
     286   G = system("Mwalk", G, curr_weight, target_weight,basering,reduction,printout);
    287287
    288288   setring xR;
     
    299299  //** compute a Groebner basis of I w.r.t. lp.
    300300  ring r = 32003,(z,y,x), lp;
    301   ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
    302   gwalk(I);
     301  ideal I = zy2+yx2+yx+3,
     302            z3x+y3+zyx-yx2-yx-3,
     303            z2yx3-y5+z2yx2+y3x2+y2x3+y3x+y2x2+3z2x+3y2+3yx,
     304            zyx5+y6-y4x2-y3x3+2zyx4-y4x-y3x2+zyx3-3z2yx+3zx3-3y3-3y2x+3zx2,
     305            yx7-y7+y5x2+y4x3+3yx6+y5x+y4x2+3yx5-6zyx3+yx4+3x5+3y4+3y3x-6zyx2+6x4+3x3-9zx;
     306  gwalk(I,0,1);
    303307}
    304308
     
    346350}
    347351
    348 proc fwalk(ideal Go, list #)
    349 "SYNTAX: fwalk(ideal i);
    350          fwalk(ideal i, intvec v, intvec w);
     352proc fwalk(ideal Go, int reduction, int printout, list #)
     353"SYNTAX: fwalk(ideal i,int reductioin);
     354         fwalk(ideal i, int reduction intvec v, intvec w);
    351355TYPE:    ideal
    352356PURPOSE: compute the standard basis of the ideal w.r.t. the
     
    372376
    373377   ideal G = fetch(xR, Go);
    374    G = system("Mfwalk", G, curr_weight, target_weight);
     378   G = system("Mfwalk", G, curr_weight, target_weight, reduction, printout);
    375379
    376380   setring xR;
     
    387391    ring r = 32003,(z,y,x), lp;
    388392    ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
    389     fwalk(I);
     393    int reduction = 1;
     394    int printout = 1;
     395    fwalk(I,reduction,printout);
    390396}
    391397
     
    437443}
    438444
    439 proc pwalk(ideal Go, int n1, int n2, list #)
    440 "SYNTAX: pwalk(int d, ideal i, int n1, int n2);
    441          pwalk(int d, ideal i, int n1, int n2, intvec v, intvec w);
     445proc pwalk(ideal Go, int n1, int n2, int reduction, int printout, list #)
     446"SYNTAX: pwalk(int d, ideal i, int n1, int n2, int reduction, int printout);
     447         pwalk(int d, ideal i, int n1, int n2, int reduction, int printout, intvec v, intvec w);
    442448TYPE:    ideal
    443449PURPOSE: compute the standard basis of the ideal, calculated via
     
    477483  ideal G = fetch(xR, Go);
    478484
    479   G = system("Mpwalk", G, n1, n2, curr_weight, target_weight,nP);
    480 
     485  G = system("Mpwalk",G,n1,n2,curr_weight,target_weight,nP,reduction,printout);
     486 
    481487  setring xR;
    482   //kill Go;
     488  //kill Go; //unused
    483489
    484490  keepring basering;
     
    492498    ring r = 32003,(z,y,x), lp;
    493499    ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
    494     //I = std(I);
    495     //ring rr = 32003,(z,y,x),lp;
    496     //ideal I = fetch(r,I);
    497     pwalk(I,2,2);
     500    int reduction = 1;
     501    int printout = 2;
     502    pwalk(I,2,2,reduction,printout);
    498503}
    499504
  • Singular/LIB/modstd.lib

    r38301d r47362dc  
    283283 * (same as size(reduce(I_reduce, G_reduce))).
    284284 * Uses parallelization. */
    285 static proc reduce_parallel(ideal I_reduce, ideal G_reduce)
     285static proc reduce_parallel(def I_reduce, def G_reduce)
    286286{
    287287    exportto(Modstd, I_reduce);
  • Singular/LIB/modwalk.lib

    • Property mode changed from 100644 to 100755
    r4132ee r47362dc  
    1616
    1717PROCEDURES:
    18 modWalk(ideal,int,int[,int,int,int,int]);        standard basis conversion of I using modular methods (chinese remainder)
     18
     19modWalk(I,#);                   standard basis conversion of I by Groebner Walk using modular methods
     20modrWalk(I,radius,pertdeg,#);   standard basis conversion of I by Random Walk using modular methods
     21modfWalk(I,#);                  standard basis conversion of I by Fractal Walk using modular methods
     22modfrWalk(I,radius,#);          standard basis conversion of I by Random Fractal Walk using modular methods
    1923";
    2024
    21 LIB "poly.lib";
    22 LIB "ring.lib";
    23 LIB "parallel.lib";
    2425LIB "rwalk.lib";
    2526LIB "grwalk.lib";
    26 LIB "modstd.lib";
    27 
    28 
    29 ////////////////////////////////////////////////////////////////////////////////
    30 
    31 static proc modpWalk(def II, int p, int variant, list #)
    32 "USAGE:  modpWalk(I,p,#); I ideal, p integer, variant integer
    33 ASSUME:  If size(#) > 0, then
    34            #[1] is an intvec describing the current weight vector
    35            #[2] is an intvec describing the target weight vector
    36 RETURN:  ideal - a standard basis of I mod p, integer - p
    37 NOTE:    The procedure computes a standard basis of the ideal I modulo p and
    38          fetches the result to the basering.
    39 EXAMPLE: example modpWalk; shows an example
    40 "
    41 {
    42   option(redSB);
    43   int k,nvar@r;
    44   def R0 = basering;
    45   string ordstr_R0 = ordstr(R0);
    46   list rl = ringlist(R0);
    47   int sizerl = size(rl);
    48   int neg = 1 - attrib(R0,"global");
    49   if(typeof(II) == "ideal")
    50   {
    51     ideal I = II;
     27LIB "modular.lib";
     28
     29proc modWalk(ideal I, list #)
     30"USAGE:   modWalk(I, [, v, w]); I ideal, v intvec, w intvec
     31RETURN:   a standard basis of I
     32NOTE:     The procedure computes a standard basis of I (over the rational
     33          numbers) by using modular methods.
     34SEE ALSO: modular
     35EXAMPLE:  example modWalk; shows an example"
     36{
     37    /* read optional parameter */
     38    if (size(#) > 0) {
     39        if (size(#) == 1) {
     40            intvec w = #[1];
     41        }
     42        if (size(#) == 2) {
     43            intvec v = #[1];
     44            intvec w = #[2];
     45        }
     46        if (size(#) > 2 || typeof(#[1]) != "intvec") {
     47            ERROR("wrong optional parameter");
     48        }
     49    }
     50
     51    /* save options */
     52    intvec opt = option(get);
     53    option(redSB);
     54
     55    /* set additional parameters */
     56    int reduction = 1;
     57    int printout = 0;
     58
     59    /* call modular() */
     60    if (size(#) > 0) {
     61        I = modular("gwalk", list(I,reduction,printout,#));
     62    }
     63    else {
     64        I = modular("gwalk", list(I,reduction,printout));
     65    }
     66
     67    /* return the result */
     68    attrib(I, "isSB", 1);
     69    option(set, opt);
     70    return(I);
     71}
     72example
     73{
     74    "EXAMPLE:";
     75    echo = 2;
     76    ring R1 = 0, (x,y,z,t), dp;
     77    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
     78    I = std(I);
     79    ring R2 = 0, (x,y,z,t), lp;
     80    ideal I = fetch(R1, I);
     81    ideal J = modWalk(I);
     82    J;
     83}
     84
     85proc modrWalk(ideal I, int radius, int pertdeg, list #)
     86"USAGE:   modrWalk(I, radius, pertdeg[, v, w]);
     87          I ideal, radius int, pertdeg int, v intvec, w intvec
     88RETURN:   a standard basis of I
     89NOTE:     The procedure computes a standard basis of I (over the rational
     90          numbers) by using modular methods.
     91SEE ALSO: modular
     92EXAMPLE:  example modrWalk; shows an example"
     93{
     94    /* read optional parameter */
     95    if (size(#) > 0) {
     96        if (size(#) == 1) {
     97            intvec w = #[1];
     98        }
     99        if (size(#) == 2) {
     100            intvec v = #[1];
     101            intvec w = #[2];
     102        }
     103        if (size(#) > 2 || typeof(#[1]) != "intvec") {
     104            ERROR("wrong optional parameter");
     105        }
     106    }
     107
     108    /* save options */
     109    intvec opt = option(get);
     110    option(redSB);
     111
     112    /* set additional parameters */
     113    int reduction = 1;
     114    int printout = 0;
     115
     116    /* call modular() */
     117    if (size(#) > 0) {
     118        I = modular("rwalk", list(I,radius,pertdeg,reduction,printout,#));
     119    }
     120    else {
     121        I = modular("rwalk", list(I,radius,pertdeg,reduction,printout));
     122    }
     123
     124    /* return the result */
     125    attrib(I, "isSB", 1);
     126    option(set, opt);
     127    return(I);
     128}
     129example
     130{
     131    "EXAMPLE:";
     132    echo = 2;
     133    ring R1 = 0, (x,y,z,t), dp;
     134    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
     135    I = std(I);
     136    ring R2 = 0, (x,y,z,t), lp;
     137    ideal I = fetch(R1, I);
    52138    int radius = 2;
    53     int pert_deg = 2;
    54   }
    55   if(typeof(II) == "list" && typeof(II[1]) == "ideal")
    56   {
    57     ideal I = II[1];
    58     if(size(II) == 2)
    59     {
    60       int radius = II[2];
    61       int pert_deg = 2;
    62     }
    63     if(size(II) == 3)
    64     {
    65       int radius = II[2];
    66       int pert_deg = II[3];
    67     }
    68   }
    69   rl[1] = p;
    70   int h = homog(I);
    71   def @r = ring(rl);
    72   setring @r;
    73   ideal i = fetch(R0,I);
    74   string order;
    75   if(system("nblocks") <= 2)
    76   {
    77     if(find(ordstr_R0, "M") + find(ordstr_R0, "lp") + find(ordstr_R0, "rp") <= 0)
    78     {
    79       order = "simple";
    80     }
    81   }
    82 
    83 //-------------------------  make i homogeneous  -----------------------------
    84   if(!mixedTest() && !h)
    85   {
    86     if(!((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg))
    87     {
    88       if(!((order == "simple") || (sizerl > 4)))
    89       {
    90         list rl@r = ringlist(@r);
    91         nvar@r = nvars(@r);
    92         intvec w;
    93         for(k = 1; k <= nvar@r; k++)
    94         {
    95           w[k] = deg(var(k));
    96         }
    97         w[nvar@r + 1] = 1;
    98         rl@r[2][nvar@r + 1] = "homvar";
    99         rl@r[3][2][2] = w;
    100         def HomR = ring(rl@r);
    101         setring HomR;
    102         ideal i = imap(@r, i);
    103         i = homog(i, homvar);
    104       }
    105     }
    106   }
    107 
    108 //-------------------------  compute a standard basis mod p  -----------------------------
    109 
    110   if(variant == 1)
    111   {
    112     if(size(#)>0)
    113     {
    114       i = rwalk(i,radius,pert_deg,#);
    115      // rwalk(i,radius,pert_deg,#); std(i);
    116     }
    117     else
    118     {
    119       i = rwalk(i,radius,pert_deg);
    120     }
    121   }
    122   if(variant == 2)
    123   {
    124     if(size(#) == 2)
    125     {
    126       i = gwalk(i,#);
    127     }
    128     else
    129     {
    130       i = gwalk(i);
    131     }
    132   }
    133   if(variant == 3)
    134   {
    135     if(size(#) == 2)
    136     {
    137       i = frandwalk(i,radius,#);
    138     }
    139     else
    140     {
    141       i = frandwalk(i,radius);
    142     }
    143   }
    144   if(variant == 4)
    145   {
    146     if(size(#) == 2)
    147     {
    148       i=fwalk(i,#);
    149     }
    150     else
    151     {
    152       i=fwalk(i);
    153     }
    154   }
    155   if(variant == 5)
    156   {
    157     if(size(#) == 2)
    158     {
    159      i=prwalk(i,radius,pert_deg,pert_deg,#);
    160     }
    161     else
    162     {
    163       i=prwalk(i,radius,pert_deg,pert_deg);
    164     }
    165   }
    166   if(variant == 6)
    167   {
    168     if(size(#) == 2)
    169     {
    170       i=pwalk(i,pert_deg,pert_deg,#);
    171     }
    172     else
    173     {
    174       i=pwalk(i,pert_deg,pert_deg);
    175     }
    176   }
    177 
    178   if(!mixedTest() && !h)
    179   {
    180     if(!((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg))
    181     {
    182       if(!((order == "simple") || (sizerl > 4)))
    183       {
    184         i = subst(i, homvar, 1);
    185         i = simplify(i, 34);
    186         setring @r;
    187         i = imap(HomR, i);
    188         i = interred(i);
    189         kill HomR;
    190       }
    191     }
    192   }
    193   setring R0;
    194   return(list(fetch(@r,i),p));
    195 }
    196 example
    197 {
    198   "EXAMPLE:"; echo = 2;
    199   option(redSB);
    200 
    201   int p = 181;
    202   intvec a = 2,1,3,4;
    203   intvec b = 1,9,1,1;
    204   ring ra = 0,x(1..4),(a(a),lp);
    205   ideal I = std(cyclic(4));
    206   ring rb = 0,x(1..4),(a(b),lp);
    207   ideal I = imap(ra,I);
    208   modpWalk(I,p,1,a,b);
    209   std(I);
    210 }
    211 
    212 ////////////////////////////////////////////////////////////////////////////////
    213 
    214 proc modWalk(def II, int variant, list #)
    215 "USAGE:  modWalk(II); II ideal or list(ideal,int)
    216 ASSUME:  If variant =
    217 @*       - 1 the Random Walk algorithm with radius II[2] is applied
    218            to II[1] if II = list(ideal, int). It is applied to II with radius 2
    219            if II is an ideal.
    220 @*       - 2, the Groebner Walk algorithm is applied to II[1] or to II, respectively.
    221 @*       - 3, the Fractal Walk algorithm with random element is applied to II[1] or II,
    222            respectively.
    223 @*       - 4, the Fractal Walk algorithm is applied to II[1] or II, respectively.
    224 @*       - 5, the Perturbation Walk algorithm with random element is applied to II[1]
    225            or II, respectively, with radius II[3] and perturbation degree II[2].
    226 @*       - 6, the Perturbation Walk algorithm is applied to II[1] or II, respectively,
    227            with perturbation degree II[3].
    228          If size(#) > 0, then # contains either 1, 2 or 4 integers such that
    229 @*       - #[1] is the number of available processors for the computation,
    230 @*       - #[2] is an optional parameter for the exactness of the computation,
    231                 if #[2] = 1, the procedure computes a standard basis for sure,
    232 @*       - #[3] is the number of primes until the first lifting,
    233 @*       - #[4] is the constant number of primes between two liftings until
    234            the computation stops.
    235 RETURN:  a standard basis of I if no warning appears.
    236 NOTE:    The procedure converts a standard basis of I (over the rational
    237          numbers) from the ordering \"a(v),lp\", "dp\" or \"Dp\" to the ordering
    238          \"(a(w),lp\" or \"a(1,0,...,0),lp\" by using modular methods.
    239          By default the procedure computes a standard basis of I for sure, but
    240          if the optional parameter #[2] = 0, it computes a standard basis of I
    241          with high probability.
    242 EXAMPLE: example modWalk; shows an example
    243 "
    244 {
    245   int TT = timer;
    246   int RT = rtimer;
    247   int i,j,pTest,sizeTest,weighted,n1;
    248   bigint N;
    249 
    250   def R0 = basering;
    251   list rl = ringlist(R0);
    252   if((npars(R0) > 0) || (rl[1] > 0))
    253   {
    254     ERROR("Characteristic of basering should be zero, basering should have no parameters.");
    255   }
    256 
    257   if(typeof(II) == "ideal")
    258   {
    259     ideal I = II;
    260     kill II;
    261     list II;
    262     II[1] = I;
    263     II[2] = 2;
    264     II[3] = 2;
    265   }
    266   else
    267   {
    268     if(typeof(II) == "list" && typeof(II[1]) == "ideal")
    269     {
    270       ideal I = II[1];
    271       if(size(II) == 1)
    272       {
    273         II[2] = 2;
    274         II[3] = 2;
    275       }
    276       if(size(II) == 2)
    277       {
    278         II[3] = 2;
    279       }
    280 
    281     }
    282     else
    283     {
    284       ERROR("Unexpected type of input.");
    285     }
    286   }
    287 
    288 //--------------------  Initialize optional parameters  ------------------------
    289   n1 = system("--cpus");
    290   if(size(#) == 0)
    291   {
    292     int exactness = 1;
    293     int n2 = 10;
    294     int n3 = 10;
    295   }
    296   else
    297   {
    298     if(size(#) == 1)
    299     {
    300       if(typeof(#[1]) == "int")
    301       {
    302         if(#[1] < n1)
    303         {
    304           n1 = #[1];
    305         }
    306         int exactness = 1;
    307         if(n1 >= 10)
    308         {
    309           int n2 = n1 + 1;
    310           int n3 = n1;
    311         }
    312         else
    313         {
    314           int n2 = 10;
    315           int n3 = 10;
    316         }
    317       }
    318       else
    319       {
    320         ERROR("Unexpected type of input.");
    321       }
    322     }
    323     if(size(#) == 2)
    324     {
    325       if(typeof(#[1]) == "int" && typeof(#[2]) == "int")
    326       {
    327         if(#[1] < n1)
    328         {
    329           n1 = #[1];
    330         }
    331         int exactness = #[2];
    332         if(n1 >= 10)
    333         {
    334           int n2 = n1 + 1;
    335           int n3 = n1;
    336         }
    337         else
    338         {
    339           int n2 = 10;
    340           int n3 = 10;
    341         }
    342       }
    343       else
    344       {
    345         if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec")
    346         {
    347           intvec curr_weight = #[1];
    348           intvec target_weight = #[2];
    349           weighted = 1;
    350           int n2 = 10;
    351           int n3 = 10;
    352           int exactness = 1;
    353         }
    354         else
    355         {
    356           ERROR("Unexpected type of input.");
    357         }
    358       }
    359     }
    360     if(size(#) == 3)
    361     {
    362       if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int")
    363       {
    364         intvec curr_weight = #[1];
    365         intvec target_weight = #[2];
    366         weighted = 1;
    367         n1 = #[3];
    368         int n2 = 10;
    369         int n3 = 10;
    370         int exactness = 1;
    371       }
    372       else
    373       {
    374         ERROR("Unexpected type of input.");
    375       }
    376     }
    377     if(size(#) == 4)
    378     {
    379       if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int"
    380           && typeof(#[4]) == "int")
    381       {
    382         intvec curr_weight = #[1];
    383         intvec target_weight = #[2];
    384         weighted = 1;
    385         if(#[1] < n1)
    386         {
    387           n1 = #[3];
    388         }
    389         int exactness = #[4];
    390         if(n1 >= 10)
    391         {
    392           int n2 = n1 + 1;
    393           int n3 = n1;
    394         }
    395         else
    396         {
    397           int n2 = 10;
    398           int n3 = 10;
    399         }
    400       }
    401       else
    402       {
    403         if(typeof(#[1]) == "int" && typeof(#[2]) == "int" && typeof(#[3]) == "int" && typeof(#[4]) == "int")
    404         {
    405           if(#[1] < n1)
    406           {
    407             n1 = #[1];
    408           }
    409           int exactness = #[2];
    410           if(n1 >= #[3])
    411           {
    412             int n2 = n1 + 1;
    413           }
    414           else
    415           {
    416             int n2 = #[3];
    417           }
    418           if(n1 >= #[4])
    419           {
    420             int n3 = n1;
    421           }
    422           else
    423           {
    424             int n3 = #[4];
    425           }
    426         }
    427         else
    428         {
    429           ERROR("Unexpected type of input.");
    430         }
    431       }
    432     }
    433     if(size(#) == 6)
    434     {
    435       if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int" && typeof(#[4]) == "int" && typeof(#[5]) == "int" && typeof(#[6]) == "int")
    436       {
    437         intvec curr_weight = #[1];
    438         intvec target_weight = #[2];
    439         weighted = 1;
    440         if(#[3] < n1)
    441         {
    442           n1 = #[3];
    443         }
    444         int exactness = #[4];
    445         if(n1 >= #[5])
    446         {
    447           int n2 = n1 + 1;
    448         }
    449         else
    450         {
    451           int n2 = #[5];
    452         }
    453         if(n1 >= #[6])
    454         {
    455           int n3 = n1;
    456         }
    457         else
    458         {
    459           int n3 = #[6];
    460         }
    461       }
    462       else
    463       {
    464         ERROR("Expected list(intvec,intvec,int,int,int,int) as optional parameter list.");
    465       }
    466     }
    467     if(size(#) == 1 || size(#) == 5 || size(#) > 6)
    468     {
    469       ERROR("Expected 0,2,3,4 or 5 optional arguments.");
    470     }
    471   }
    472   if(printlevel >= 10)
    473   {
    474   "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)+", exactness = "+string(exactness);
    475   }
    476 
    477 //-------------------------  Save current options  -----------------------------
    478   intvec opt = option(get);
    479   //option(redSB);
    480 
    481 //--------------------  Initialize the list of primes  -------------------------
    482   int tt = timer;
    483   int rt = rtimer;
    484   int en = 2134567879;
    485   int an = 1000000000;
    486   intvec L = primeList(I,n2);
    487   if(n2 > 4)
    488   {
    489   //  L[5] = prime(random(an,en));
    490   }
    491   if(printlevel >= 10)
    492   {
    493     "CPU-time for primeList: "+string(timer-tt)+" seconds.";
    494     "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
    495   }
    496   int h = homog(I);
    497   list P,T1,T2,LL,Arguments,PP;
    498   ideal J,K,H;
    499 
    500 //-------------------  parallelized Groebner Walk in positive characteristic  --------------------
    501 
    502   if(weighted)
    503   {
    504     for(i=1; i<=size(L); i++)
    505     {
    506       Arguments[i] = list(II,L[i],variant,list(curr_weight,target_weight));
    507     }
    508   }
    509   else
    510   {
    511     for(i=1; i<=size(L); i++)
    512     {
    513       Arguments[i] = list(II,L[i],variant);
    514     }
    515   }
    516   P = parallelWaitAll("modpWalk",Arguments);
    517   for(i=1; i<=size(P); i++)
    518   {
    519     T1[i] = P[i][1];
    520     T2[i] = bigint(P[i][2]);
    521   }
    522 
    523   while(1)
    524   {
    525     LL = deleteUnluckyPrimes(T1,T2,h);
    526     T1 = LL[1];
    527     T2 = LL[2];
    528 //-------------------  Now all leading ideals are the same  --------------------
    529 //-------------------  Lift results to basering via farey  ---------------------
    530 
    531     tt = timer; rt = rtimer;
    532     N = T2[1];
    533     for(i=2; i<=size(T2); i++)
    534     {
    535       N = N*T2[i];
    536     }
    537     H = chinrem(T1,T2);
    538     J = parallelFarey(H,N,n1);
    539     //J=farey(H,N);
    540     if(printlevel >= 10)
    541     {
    542       "CPU-time for lifting-process is "+string(timer - tt)+" seconds.";
    543       "Real-time for lifting-process is "+string(rtimer - rt)+" seconds.";
    544     }
    545 
    546 //----------------  Test if we already have a standard basis of I --------------
    547 
    548     tt = timer; rt = rtimer;
    549     pTest = pTestSB(I,J,L,variant);
    550     //pTest = primeTestSB(I,J,L,variant);
    551     if(printlevel >= 10)
    552     {
    553       "CPU-time for pTest is "+string(timer - tt)+" seconds.";
    554       "Real-time for pTest is "+string(rtimer - rt)+" seconds.";
    555     }
    556     if(pTest)
    557     {
    558       if(printlevel >= 10)
    559       {
    560         "CPU-time for computation without final tests is "+string(timer - TT)+" seconds.";
    561         "Real-time for computation without final tests is "+string(rtimer - RT)+" seconds.";
    562       }
    563       attrib(J,"isSB",1);
    564       if(exactness == 0)
    565       {
    566         option(set, opt);
    567         return(J);
    568       }
    569       else
    570       {
    571         tt = timer;
    572         rt = rtimer;
    573         sizeTest = 1 - isIdealIncluded(I,J,n1);
    574         if(printlevel >= 10)
    575         {
    576           "CPU-time for checking if I subset <G> is "+string(timer - tt)+" seconds.";
    577           "Real-time for checking if I subset <G> is "+string(rtimer - rt)+" seconds.";
    578         }
    579         if(sizeTest == 0)
    580         {
    581           tt = timer;
    582           rt = rtimer;
    583           K = std(J);
    584           if(printlevel >= 10)
    585           {
    586             "CPU-time for last std-computation is "+string(timer - tt)+" seconds.";
    587             "Real-time for last std-computation is "+string(rtimer - rt)+" seconds.";
    588           }
    589           if(size(reduce(K,J)) == 0)
    590           {
    591             option(set, opt);
    592             return(J);
    593           }
    594         }
    595       }
    596     }
    597 //--------------  We do not already have a standard basis of I, therefore do the main computation for more primes  --------------
    598 
    599     T1 = H;
    600     T2 = N;
    601     j = size(L)+1;
    602     tt = timer; rt = rtimer;
    603     L = primeList(I,n3,L,n1);
    604     if(printlevel >= 10)
    605     {
    606       "CPU-time for primeList: "+string(timer-tt)+" seconds.";
    607       "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
    608     }
    609     Arguments = list();
    610     PP = list();
    611     if(weighted)
    612     {
    613       for(i=j; i<=size(L); i++)
    614       {
    615         //Arguments[i-j+1] = list(II,L[i],variant,list(curr_weight,target_weight));
    616         Arguments[size(Arguments)+1] = list(II,L[i],variant,list(curr_weight,target_weight));
    617       }
    618     }
    619     else
    620     {
    621       for(i=j; i<=size(L); i++)
    622       {
    623         //Arguments[i-j+1] = list(II,L[i],variant);
    624         Arguments[size(Arguments)+1] = list(II,L[i],variant);
    625       }
    626     }
    627     PP = parallelWaitAll("modpWalk",Arguments);
    628     if(printlevel >= 10)
    629     {
    630       "parallel modpWalk";
    631     }
    632     for(i=1; i<=size(PP); i++)
    633     {
    634       //P[size(P) + 1] = PP[i];
    635       T1[size(T1) + 1] = PP[i][1];
    636       T2[size(T2) + 1] = bigint(PP[i][2]);
    637     }
    638   }
    639   if(printlevel >= 10)
    640   {
    641     "CPU-time for computation with final tests is "+string(timer - TT)+" seconds.";
    642     "Real-time for computation with final tests is "+string(rtimer - RT)+" seconds.";
    643   }
    644 }
    645 
    646 example
    647 {
    648   "EXAMPLE:";
    649   echo = 2;
    650   ring R=0,(x,y,z),lp;
    651   ideal I=-x+y2z-z,xz+1,x2+y2-1;
    652   // I is a standard basis in dp
    653   ideal J = modWalk(I,1);
    654   J;
    655 }
    656 
    657 ////////////////////////////////////////////////////////////////////////////////
    658 static proc isIdealIncluded(ideal I, ideal J, int n1)
    659 "USAGE:  isIdealIncluded(I,J,int n1); I ideal, J ideal, n1 integer
    660 "
    661 {
    662   if(n1 > 1)
    663   {
    664     int k;
    665     list args,results;
    666     for(k=1; k<=size(I); k++)
    667     {
    668       args[k] = list(ideal(I[k]),J,1);
    669     }
    670     results = parallelWaitAll("reduce",args);
    671     for(k=1; k<=size(results); k++)
    672     {
    673       if(results[k] == 0)
    674       {
    675         return(1);
    676       }
    677     }
    678     return(0);
    679   }
    680   else
    681   {
    682     if(reduce(I,J,1) == 0)
    683     {
    684       return(1);
    685     }
    686     else
    687     {
    688       return(0);
    689     }
    690   }
    691 }
    692 
    693 ////////////////////////////////////////////////////////////////////////////////
    694 static proc parallelChinrem(list T1, list T2, int n1)
    695 "USAGE:  parallelChinrem(T1,T2); T1 list of ideals, T2 list of primes, n1 integer"
    696 {
    697   int i,j,k;
    698 
    699   ideal H,J;
    700 
    701   list arguments_chinrem,results_chinrem;
    702   for(i=1; i<=size(T1); i++)
    703   {
    704     J = ideal(T1[i]);
    705     attrib(J,"isSB",1);
    706     arguments_chinrem[size(arguments_chinrem)+1] = list(list(J),T2);
    707   }
    708   results_chinrem = parallelWaitAll("chinrem",arguments_chinrem);
    709     for(j=1; j <= size(results_chinrem); j++)
    710     {
    711       J = results_chinrem[j];
    712       attrib(J,"isSB",1);
    713       if(isIdealIncluded(J,H,n1) == 0)
    714       {
    715         if(H == 0)
    716         {
    717           H = J;
    718         }
    719         else
    720         {
    721           H = H,J;
    722         }
    723       }
    724     }
    725   return(H);
    726 }
    727 
    728 ////////////////////////////////////////////////////////////////////////////////
    729 static proc parallelFarey(ideal H, bigint N, int n1)
    730 "USAGE:  parallelFarey(H,N,n1); H ideal, N bigint, n1 integer
    731 "
    732 {
    733   int i,j;
    734   int ii = 1;
    735   list arguments_farey,results_farey;
    736   for(i=1; i<=size(H); i++)
    737   {
    738     for(j=1; j<=size(H[i]); j++)
    739     {
    740       arguments_farey[size(arguments_farey)+1] = list(H[i][j],N);
    741     }
    742   }
    743   results_farey = parallelWaitAll("farey",arguments_farey);
    744   ideal J,K;
    745   poly f_farey;
    746   while(ii<=size(results_farey))
    747   {
    748     for(i=1; i<=size(H); i++)
    749     {
    750       f_farey = 0;
    751       for(j=1; j<=size(H[i]); j++)
    752       {
    753         f_farey = f_farey + results_farey[ii][1];
    754         ii++;
    755       }
    756       K = ideal(f_farey);
    757       attrib(K,"isSB",1);
    758       attrib(J,"isSB",1);
    759       if(isIdealIncluded(K,J,n1) == 0)
    760       {
    761         if(J == 0)
    762         {
    763           J = K;
    764         }
    765         else
    766         {
    767           J = J,K;
    768         }
    769       }
    770     }
    771   }
    772   return(J);
    773 }
    774 //////////////////////////////////////////////////////////////////////////////////
    775 static proc primeTestSB(def II, ideal J, list L, int variant, list #)
    776 "USAGE:  primeTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int
    777 RETURN:  1 (resp. 0) if for a randomly chosen prime p that is not in L
    778          J mod p is (resp. is not) a standard basis of I mod p
    779 EXAMPLE: example primeTestSB; shows an example
    780 "
    781 {
    782 if(typeof(II) == "ideal")
    783   {
    784   ideal I = II;
    785   int radius = 2;
    786   }
    787 if(typeof(II) == "list")
    788   {
    789   ideal I = II[1];
    790   int radius = II[2];
    791   }
    792 
    793 int i,j,k,p;
    794 def R = basering;
    795 list r = ringlist(R);
    796 
    797 while(!j)
    798   {
    799   j = 1;
    800   p = prime(random(1000000000,2134567879));
    801   for(i = 1; i <= size(L); i++)
    802     {
    803     if(p == L[i])
    804       {
    805       j = 0;
    806       break;
    807       }
    808     }
    809   if(j)
    810     {
    811     for(i = 1; i <= ncols(I); i++)
    812       {
    813       for(k = 2; k <= size(I[i]); k++)
    814         {
    815         if((denominator(leadcoef(I[i][k])) mod p) == 0)
    816           {
    817           j = 0;
    818           break;
    819           }
    820         }
    821       if(!j)
    822         {
    823         break;
    824         }
    825       }
    826     }
    827   if(j)
    828     {
    829     if(!primeTest(I,p))
    830       {
    831       j = 0;
    832       }
    833     }
    834   }
    835 r[1] = p;
    836 def @R = ring(r);
    837 setring @R;
    838 ideal I = imap(R,I);
    839 ideal J = imap(R,J);
    840 attrib(J,"isSB",1);
    841 
    842 int t = timer;
    843 j = 1;
    844 if(isIncluded(I,J) == 0)
    845   {
    846   j = 0;
    847   }
    848 if(printlevel >= 11)
    849   {
    850   "isIncluded(I,J) takes "+string(timer - t)+" seconds";
    851   "j = "+string(j);
    852   }
    853 t = timer;
    854 if(j)
    855   {
    856   if(size(#) > 0)
    857     {
    858     ideal K = modpWalk(I,p,variant,#)[1];
    859     }
    860   else
    861     {
    862     ideal K = modpWalk(I,p,variant)[1];
    863     }
    864   t = timer;
    865   if(isIncluded(J,K) == 0)
    866     {
    867     j = 0;
    868     }
    869   if(printlevel >= 11)
    870     {
    871     "isIncluded(K,J) takes "+string(timer - t)+" seconds";
    872     "j = "+string(j);
    873     }
    874   }
    875 setring R;
    876 
    877 return(j);
    878 }
    879 example
    880 { "EXAMPLE:"; echo = 2;
    881    intvec L = 2,3,5;
    882    ring r = 0,(x,y,z),lp;
    883    ideal I = x+1,x+y+1;
    884    ideal J = x+1,y;
    885    primeTestSB(I,I,L,1);
    886    primeTestSB(I,J,L,1);
    887 }
    888 
    889 ////////////////////////////////////////////////////////////////////////////////
    890 static proc pTestSB(ideal I, ideal J, list L, int variant, list #)
    891 "USAGE:  pTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int
    892 RETURN:  1 (resp. 0) if for a randomly chosen prime p that is not in L
    893          J mod p is (resp. is not) a standard basis of I mod p
    894 EXAMPLE: example pTestSB; shows an example
    895 "
    896 {
    897    int i,j,k,p;
    898    def R = basering;
    899    list r = ringlist(R);
    900 
    901    while(!j)
    902    {
    903       j = 1;
    904       p = prime(random(1000000000,2134567879));
    905       for(i = 1; i <= size(L); i++)
    906       {
    907          if(p == L[i]) { j = 0; break; }
    908       }
    909       if(j)
    910       {
    911          for(i = 1; i <= ncols(I); i++)
    912          {
    913             for(k = 2; k <= size(I[i]); k++)
    914             {
    915                if((denominator(leadcoef(I[i][k])) mod p) == 0) { j = 0; break; }
    916             }
    917             if(!j){ break; }
    918          }
    919       }
    920       if(j)
    921       {
    922          if(!primeTest(I,p)) { j = 0; }
    923       }
    924    }
    925    r[1] = p;
    926    def @R = ring(r);
    927    setring @R;
    928    ideal I = imap(R,I);
    929    ideal J = imap(R,J);
    930    attrib(J,"isSB",1);
    931 
    932    int t = timer;
    933    j = 1;
    934    if(isIncluded(I,J) == 0) { j = 0; }
    935 
    936    if(printlevel >= 11)
    937    {
    938       "isIncluded(I,J) takes "+string(timer - t)+" seconds";
    939       "j = "+string(j);
    940    }
    941 
    942    t = timer;
    943    if(j)
    944    {
    945       if(size(#) > 0)
    946       {
    947          ideal K = modpStd(I,p,variant,#[1])[1];
    948       }
    949       else
    950       {
    951          ideal K = groebner(I);
    952       }
    953       t = timer;
    954       if(isIncluded(J,K) == 0) { j = 0; }
    955 
    956       if(printlevel >= 11)
    957       {
    958          "isIncluded(J,K) takes "+string(timer - t)+" seconds";
    959          "j = "+string(j);
    960       }
    961    }
    962    setring R;
    963    return(j);
    964 }
    965 example
    966 { "EXAMPLE:"; echo = 2;
    967    intvec L = 2,3,5;
    968    ring r = 0,(x,y,z),dp;
    969    ideal I = x+1,x+y+1;
    970    ideal J = x+1,y;
    971    pTestSB(I,I,L,2);
    972    pTestSB(I,J,L,2);
    973 }
    974 ////////////////////////////////////////////////////////////////////////////////
    975 static proc mixedTest()
    976 "USAGE:  mixedTest();
    977 RETURN:  1 if ordering of basering is mixed, 0 else
    978 EXAMPLE: example mixedTest(); shows an example
    979 "
    980 {
    981    int i,p,m;
    982    for(i = 1; i <= nvars(basering); i++)
    983    {
    984       if(var(i) > 1)
    985       {
    986          p++;
    987       }
    988       else
    989       {
    990          m++;
    991       }
    992    }
    993    if((p > 0) && (m > 0)) { return(1); }
    994    return(0);
    995 }
    996 example
    997 { "EXAMPLE:"; echo = 2;
    998    ring R1 = 0,(x,y,z),dp;
    999    mixedTest();
    1000    ring R2 = 31,(x(1..4),y(1..3)),(ds(4),lp(3));
    1001    mixedTest();
    1002    ring R3 = 181,x(1..9),(dp(5),lp(4));
    1003    mixedTest();
    1004 }
     139    int pertdeg = 3;
     140    ideal J = modrWalk(I,radius,pertdeg);
     141    J;
     142}
     143
     144proc modfWalk(ideal I, list #)
     145"USAGE:   modfWalk(I, [, v, w]); I ideal, v intvec, w intvec
     146RETURN:   a standard basis of I
     147NOTE:     The procedure computes a standard basis of I (over the rational
     148          numbers) by using modular methods.
     149SEE ALSO: modular
     150EXAMPLE:  example modfWalk; shows an example"
     151{
     152    /* read optional parameter */
     153    if (size(#) > 0) {
     154        if (size(#) == 1) {
     155            intvec w = #[1];
     156        }
     157        if (size(#) == 2) {
     158            intvec v = #[1];
     159            intvec w = #[2];
     160        }
     161        if (size(#) > 2 || typeof(#[1]) != "intvec") {
     162            ERROR("wrong optional parameter");
     163        }
     164    }
     165
     166    /* save options */
     167    intvec opt = option(get);
     168    option(redSB);
     169
     170    /* set additional parameters */
     171    int reduction = 1;
     172    int printout = 0;
     173
     174    /* call modular() */
     175    if (size(#) > 0) {
     176        I = modular("fwalk", list(I,reduction,printout,#));
     177    }
     178    else {
     179        I = modular("fwalk", list(I,reduction,printout));
     180    }
     181
     182    /* return the result */
     183    attrib(I, "isSB", 1);
     184    option(set, opt);
     185    return(I);
     186}
     187example
     188{
     189    "EXAMPLE:";
     190    echo = 2;
     191    ring R1 = 0, (x,y,z,t), dp;
     192    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
     193    I = std(I);
     194    ring R2 = 0, (x,y,z,t), lp;
     195    ideal I = fetch(R1, I);
     196    ideal J = modfWalk(I);
     197    J;
     198}
     199
     200proc modfrWalk(ideal I, int radius, list #)
     201"USAGE:   modfrWalk(I, radius [, v, w]); I ideal, radius int, v intvec, w intvec
     202RETURN:   a standard basis of I
     203NOTE:     The procedure computes a standard basis of I (over the rational
     204          numbers) by using modular methods.
     205SEE ALSO: modular
     206EXAMPLE:  example modfrWalk; shows an example"
     207{
     208    /* read optional parameter */
     209    if (size(#) > 0) {
     210        if (size(#) == 1) {
     211            intvec w = #[1];
     212        }
     213        if (size(#) == 2) {
     214            intvec v = #[1];
     215            intvec w = #[2];
     216        }
     217        if (size(#) > 2 || typeof(#[1]) != "intvec") {
     218            ERROR("wrong optional parameter");
     219        }
     220    }
     221
     222    /* save options */
     223    intvec opt = option(get);
     224    option(redSB);
     225
     226    /* set additional parameters */
     227    int reduction = 1;
     228    int printout = 0;
     229
     230    /* call modular() */
     231    if (size(#) > 0) {
     232        I = modular("frandwalk", list(I,radius,reduction,printout,#));
     233    }
     234    else {
     235        I = modular("frandwalk", list(I,radius,reduction,printout));
     236    }
     237
     238    /* return the result */
     239    attrib(I, "isSB", 1);
     240    option(set, opt);
     241    return(I);
     242}
     243example
     244{
     245    "EXAMPLE:";
     246    echo = 2;
     247    ring R1 = 0, (x,y,z,t), dp;
     248    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
     249    I = std(I);
     250    ring R2 = 0, (x,y,z,t), lp;
     251    ideal I = fetch(R1, I);
     252    int radius = 2;
     253    ideal J = modfrWalk(I,radius);
     254    J;
     255}
  • Singular/LIB/ncfactor.lib

    r38301d r47362dc  
    64986498    def Firstweyl = nc_algebra(1,1);
    64996499    setring Firstweyl;
    6500     ideal M = 0:nvars(@r);
     6500    ideal M = ideal(0:nvars(@r));
    65016501    M[varnumX]=var(2);
    65026502    M[varnumD]=var(1);
  • Singular/LIB/nctools.lib

    r38301d r47362dc  
    15971597  matrix U = UpOneMatrix(N);
    15981598  if (size(ideal(C-U)) <> 0) { return(notW); } // lt(xy)<>lt(yx)
    1599   ideal I = D;
     1599  ideal I = ideal(D);
    16001600  if (size(I) <> n) { return(notW); } // not n entries<>0
    16011601  I = simplify(I,4+2);
  • Singular/LIB/primdec.lib

    r38301d r47362dc  
    34133413
    34143414  ideal jkeep=@j;
    3415   if(ordstr(@P)[1]=="w")
     3415  if((ordstr(@P)[1]=="w")&&(size(ringlist(@P)[3])==2)) // weighted ordering
    34163416  {
    34173417    list gnir_l=ringlist(gnir);
  • Singular/LIB/rwalk.lib

    • Property mode changed from 100644 to 100755
    r4132ee r47362dc  
    1010rwalk(ideal,int,int[,intvec,intvec]);   standard basis of ideal via Random Walk algorithm
    1111rwalk(ideal,int[,intvec,intvec]);       standard basis of ideal via Random Perturbation Walk algorithm
    12 frwalk(ideal,int[,intvec,intvec]);      standard basis of ideal via Random Fractal Walk algorithm
     12frandwalk(ideal,int[,intvec,intvec]);      standard basis of ideal via Random Fractal Walk algorithm
    1313";
    1414
     
    141141 * Random Walk  *
    142142 ****************/
    143 proc rwalk(ideal Go, int radius, int pert_deg, list #)
     143proc rwalk(ideal Go, int radius, int pert_deg, int reduction, int printout, list #)
    144144"SYNTAX: rwalk(ideal i, int radius);
    145145         if size(#)>0 then rwalk(ideal i, int radius, intvec v, intvec w);
     146         intermediate Groebner bases are not reduced if reduction = 0
    146147TYPE:    ideal
    147148PURPOSE: compute the standard basis of the ideal, calculated via
     
    178179
    179180ideal G = fetch(xR, Go);
    180 G = system("Mrwalk", G, curr_weight, target_weight, radius, pert_deg, basering);
     181G = system("Mrwalk", G, curr_weight, target_weight, radius, pert_deg, reduction, printout);
    181182
    182183setring xR;
     
    196197  int radius = 1;
    197198  int perturb_deg = 2;
    198   rwalk(I,radius,perturb_deg);
     199  int reduction = 0;
     200  int printout = 1;
     201  rwalk(I,radius,perturb_deg,reduction,printout);
    199202}
    200203
     
    202205 * Perturbation Walk with random element *
    203206 *****************************************/
    204 proc prwalk(ideal Go, int radius, int o_pert_deg, int t_pert_deg, list #)
     207proc prwalk(ideal Go, int radius, int o_pert_deg, int t_pert_deg, int reduction, int printout, list #)
    205208"SYNTAX: rwalk(ideal i, int radius);
    206209         if size(#)>0 then rwalk(ideal i, int radius, intvec v, intvec w);
     
    227230  OSCTW= OrderStringalp_NP("al", #);
    228231  }
     232int nP = OSCTW[1];
    229233string ord_str = OSCTW[2];
    230234intvec curr_weight = OSCTW[3]; // original weight vector
     
    238242
    239243ideal G = fetch(xR, Go);
    240 G = system("Mprwalk", G, curr_weight, target_weight, radius, o_pert_deg, t_pert_deg, basering);
     244G = system("Mprwalk", G, curr_weight, target_weight, radius, o_pert_deg, t_pert_deg,
     245           nP, reduction, printout);
    241246
    242247setring xR;
     
    257262  int o_perturb_deg = 2;
    258263  int t_perturb_deg = 2;
    259   prwalk(I,radius,o_perturb_deg,t_perturb_deg);
     264  int reduction = 0;
     265  int printout = 2;
     266  prwalk(I,radius,o_perturb_deg,t_perturb_deg,reduction,printout);
    260267}
    261268
     
    263270 * Fractal Walk with random element *
    264271 ************************************/
    265 proc frandwalk(ideal Go, int radius, list #)
    266 "SYNTAX: frwalk(ideal i, int radius);
    267          frwalk(ideal i, int radius, intvec v, intvec w);
     272proc frandwalk(ideal Go, int radius, int reduction, int printout, list #)
     273"SYNTAX: frwalk(ideal i, int radius, int reduction, int printout);
     274         frwalk(ideal i, int radius, int reduction, int printout, intvec v, intvec w);
    268275TYPE:    ideal
    269276PURPOSE: compute the standard basis of the ideal, calculated via
     
    299306   ideal G = fetch(xR, Go);
    300307   int pert_deg = 2;
    301    G = system("Mfrwalk", G, curr_weight, target_weight, radius);
     308
     309   G = system("Mfrwalk", G, curr_weight, target_weight, radius, reduction, printout);
    302310
    303311   setring xR;
     
    314322    ring r = 0,(z,y,x), lp;
    315323    ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
    316     frandwalk(I,2);
    317 }
     324    int reduction = 0;
     325    frandwalk(I,2,0,1);
     326}
  • Singular/LIB/swalk.lib

    • Property mode changed from 100644 to 100755
  • Singular/blackbox.cc

    r38301d r47362dc  
    5454BOOLEAN blackbox_default_serialize(blackbox */*b*/, void */*d*/, si_link /*f*/)
    5555{
     56  WerrorS("blackbox_serialize is not implemented");
    5657  return TRUE;
    5758}
     
    5960BOOLEAN blackbox_default_deserialize(blackbox **/*b*/, void **/*d*/, si_link /*f*/)
    6061{
     62  WerrorS("blackbox_deserialize is not implemented");
    6163  return TRUE;
    6264}
  • Singular/cntrlc.cc

    r38301d r47362dc  
    55* ABSTRACT - interupt handling
    66*/
    7 #include <kernel/mod2.h>
    8 
    9 /* includes */
    10 #ifdef DecAlpha_OSF1
    11 #define _XOPEN_SOURCE_EXTENDED
    12 #endif /* MP3-Y2 0.022UF */
    13 
    14 #include <omalloc/omalloc.h>
    15 
    16 #include <reporter/si_signals.h>
    17 #include <Singular/fevoices.h>
    18 
    19 #include <Singular/tok.h>
    20 #include <Singular/ipshell.h>
    21 void sig_chld_hdl(int sig); /*#include <Singular/links/ssiLink.h>*/
    22 #include <Singular/cntrlc.h>
    23 #include <Singular/feOpt.h>
    24 #include <Singular/misc_ip.h>
    25 #include <Singular/links/silink.h>
    26 #include <Singular/links/ssiLink.h>
    27 
    287#include <stdio.h>
    298#include <stddef.h>
     
    3312#include <sys/types.h>
    3413#include <sys/wait.h>
     14
     15#include <kernel/mod2.h>
     16
     17#include <omalloc/omalloc.h>
     18
     19#include <reporter/si_signals.h>
     20#include <Singular/fevoices.h>
     21
     22#include <Singular/tok.h>
     23#include <Singular/ipshell.h>
     24#include <Singular/cntrlc.h>
     25#include <Singular/feOpt.h>
     26#include <Singular/misc_ip.h>
     27#include <Singular/links/silink.h>
     28#include <Singular/links/ssiLink.h>
    3529
    3630/* undef, if you don't want GDB to come up on error */
     
    352346      if (feOptValue(FE_OPT_EMACS) == NULL)
    353347      {
    354         fputs("abort after this command(a), abort immediately(r), print backtrace(b), continue(c) or quit Singular(q) ?",stderr);fflush(stderr);
     348        fputs("abort after this command(a), abort immediately(r), print backtrace(b), continue(c) or quit Singular(q) ?",stderr);
     349        fflush(stderr);fflush(stdin);
    355350        c = fgetc(stdin);
    356351      }
     
    371366                  fputs("** Warning: Singular should be restarted as soon as possible **\n",stderr);
    372367                  fflush(stderr);
     368                  extern void my_yy_flush();
     369                  my_yy_flush();
     370                  currentVoice=feInitStdin(NULL);
    373371                  longjmp(si_start_jmpbuf,1);
    374372                }
  • Singular/dyn_modules/gfanlib/bbcone.cc

    r38301d r47362dc  
    17391739}
    17401740
    1741 
     1741#if 0
     1742BOOLEAN bbcone_serialize(blackbox *b, void *d, si_link f)
     1743{
     1744  ssiInfo *dd = (ssiInfo *)f->data;
     1745  sleftv l;
     1746  memset(&l,0,sizeof(l));
     1747  l.rtyp=STRING_CMD;
     1748  l.data=(void*)"cone";
     1749  f->m->Write(f, &l);
     1750  gfan::ZCone *Z=((gfan::ZCone*) d;
     1751  /* AMBIENT_DIM */ fprintf(dd->f_write("%d ",Z->ambientDimension());
     1752  /* FACETS or INEQUALITIES */ fprintf(dd->f_write("%d ",Z->areFacetsKnown());
     1753  gfan::ZMatrix i=Z->getInequalities();
     1754....
     1755  /* LINEAR_SPAN or EQUATIONS */ fprintf(dd->f_write("%d ",Z->areImpliedEquationsKnown());
     1756  gfan::ZMatrix e=Z->getEquations();
     1757....
     1758  /* RAYS */
     1759  gfan::ZMatrix r=Z->extremeRays();
     1760....
     1761  /* LINEALITY_SPACE */
     1762  gfan::ZMatrix l=Z->generatorsOfLinealitySpace();
     1763....
     1764  return FALSE;
     1765}
     1766BOOLEAN bbcone_deserialize(blackbox **b, void **d, si_link f)
     1767{
     1768  ssiInfo *dd = (ssiInfo *)f->data;
     1769  gfan::ZCone *Z;
     1770  /* AMBIENT_DIM */ = s_readint(dd->f_read);
     1771  /* areFacetsKnown: */ = s_readint(dd->f_read);
     1772  if (areFacetsKnown)
     1773  ....FACETS
     1774  else
     1775  ....INEQUALITIES
     1776  /* areImpliedEquationsKnown*/ = s_readint(dd->f_read);
     1777  if(areImpliedEquationsKnown)
     1778  ....EQUATIONS
     1779  else
     1780  ...LINEAR_SPAN
     1781  ...RAYS
     1782  ...LINEALITY_SPACE
     1783  *d=Z;
     1784  return FALSE;
     1785}
     1786#endif
    17421787void bbcone_setup(SModulFunctions* p)
    17431788{
     
    17531798  b->blackbox_Assign=bbcone_Assign;
    17541799  b->blackbox_Op2=bbcone_Op2;
     1800  //b->blackbox_serialize=bbcone_serialize;
     1801  //b->blackbox_deserialize=bbcone_deserialize;
    17551802  p->iiAddCproc("","coneViaInequalities",FALSE,coneViaNormals);
    17561803  p->iiAddCproc("","coneViaPoints",FALSE,coneViaRays);
  • Singular/extra.cc

    r4132ee r47362dc  
    18641864    if (strcmp(sys_cmd, "Mwalk") == 0)
    18651865    {
    1866       const short t[]={4,IDEAL_CMD,INTVEC_CMD,INTVEC_CMD,RING_CMD};
     1866      const short t[]={6,IDEAL_CMD,INTVEC_CMD,INTVEC_CMD,RING_CMD,INT_CMD,INT_CMD};
    18671867      if (!iiCheckTypes(h,t,1)) return TRUE;
    18681868      if (((intvec*) h->next->Data())->length() != currRing->N &&
     
    18751875      ideal arg1 = (ideal) h->Data();
    18761876      intvec* arg2 = (intvec*) h->next->Data();
    1877       intvec* arg3   =  (intvec*) h->next->next->Data();
    1878       ring arg4   =  (ring) h->next->next->next->Data();
    1879       ideal result = (ideal) Mwalk(arg1, arg2, arg3,arg4);
     1877      intvec* arg3 = (intvec*) h->next->next->Data();
     1878      ring arg4 = (ring) h->next->next->next->Data();
     1879      int arg5 = (int) (long) h->next->next->next->next->Data();
     1880      int arg6 = (int) (long) h->next->next->next->next->next->Data();
     1881      ideal result = (ideal) Mwalk(arg1, arg2, arg3, arg4, arg5, arg6);
    18801882      res->rtyp = IDEAL_CMD;
    18811883      res->data =  result;
     
    19131915    if (strcmp(sys_cmd, "Mpwalk") == 0)
    19141916    {
    1915       const short t[]={6,IDEAL_CMD,INT_CMD,INT_CMD,INTVEC_CMD,INTVEC_CMD,INT_CMD};
     1917      const short t[]={8,IDEAL_CMD,INT_CMD,INT_CMD,INTVEC_CMD,INTVEC_CMD,INT_CMD,INT_CMD,INT_CMD};
    19161918      if (!iiCheckTypes(h,t,1)) return TRUE;
    19171919      if(((intvec*) h->next->next->next->Data())->length() != currRing->N &&
     
    19271929      intvec* arg5 = (intvec*) h->next->next->next->next->Data();
    19281930      int arg6 = (int) (long) h->next->next->next->next->next->Data();
    1929       ideal result = (ideal) Mpwalk(arg1, arg2, arg3, arg4, arg5,arg6);
     1931      int arg7 = (int) (long) h->next->next->next->next->next->next->Data();
     1932      int arg8 = (int) (long) h->next->next->next->next->next->next->next->Data();
     1933      ideal result = (ideal) Mpwalk(arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8);
    19301934      res->rtyp = IDEAL_CMD;
    19311935      res->data =  result;
     
    19391943    if (strcmp(sys_cmd, "Mrwalk") == 0)
    19401944    {
    1941       const short t[]={6,IDEAL_CMD,INTVEC_CMD,INTVEC_CMD,INT_CMD,INT_CMD,RING_CMD};
     1945      const short t[]={7,IDEAL_CMD,INTVEC_CMD,INTVEC_CMD,INT_CMD,INT_CMD,INT_CMD,INT_CMD};
    19421946      if (!iiCheckTypes(h,t,1)) return TRUE;
    1943       if((((intvec*) h->next->Data())->length() != currRing->N &&
    1944          ((intvec*) h->next->next->Data())->length() != currRing->N ) &&
    1945          (((intvec*) h->next->Data())->length() != (currRing->N)*(currRing->N) &&
    1946          ((intvec*) h->next->next->Data())->length() != (currRing->N)*(currRing->N) ))
     1947      if(((intvec*) h->next->Data())->length() != currRing->N &&
     1948         ((intvec*) h->next->Data())->length() != (currRing->N)*(currRing->N) &&
     1949         ((intvec*) h->next->next->Data())->length() != currRing->N &&
     1950         ((intvec*) h->next->next->Data())->length() != (currRing->N)*(currRing->N) )
    19471951      {
    19481952        Werror("system(\"Mrwalk\" ...) intvecs not of length %d or %d\n",
     
    19551959      int arg4 = (int)(long) h->next->next->next->Data();
    19561960      int arg5 = (int)(long) h->next->next->next->next->Data();
    1957       ring arg6 = (ring) h->next->next->next->next->next->Data();
    1958       ideal result = (ideal) Mrwalk(arg1, arg2, arg3, arg4, arg5, arg6);
     1961      int arg6 = (int)(long) h->next->next->next->next->next->Data();
     1962      int arg7 = (int)(long) h->next->next->next->next->next->next->Data();
     1963      ideal result = (ideal) Mrwalk(arg1, arg2, arg3, arg4, arg5, arg6, arg7);
    19591964      res->rtyp = IDEAL_CMD;
    19601965      res->data =  result;
     
    20182023    if (strcmp(sys_cmd, "Mfwalk") == 0)
    20192024    {
    2020       const short t[]={3,IDEAL_CMD,INTVEC_CMD,INTVEC_CMD};
     2025      const short t[]={5,IDEAL_CMD,INTVEC_CMD,INTVEC_CMD,INT_CMD,INT_CMD};
    20212026      if (!iiCheckTypes(h,t,1)) return TRUE;
    20222027      if (((intvec*) h->next->Data())->length() != currRing->N &&
     
    20252030        Werror("system(\"Mfwalk\" ...) intvecs not of length %d\n",
    20262031                 currRing->N);
    2027         return TRUE;
    2028       }
    2029       ideal arg1 = (ideal) h->Data();
    2030       intvec* arg2 = (intvec*) h->next->Data();
    2031       intvec* arg3   =  (intvec*) h->next->next->Data();
    2032       ideal result = (ideal) Mfwalk(arg1, arg2, arg3);
    2033       res->rtyp = IDEAL_CMD;
    2034       res->data =  result;
    2035       return FALSE;
    2036     }
    2037     else
    2038   #endif
    2039   /*==================== Mfrwalk =================*/
    2040   #ifdef HAVE_WALK
    2041     if (strcmp(sys_cmd, "Mfrwalk") == 0)
    2042     {
    2043       const short t[]={6,IDEAL_CMD,INTVEC_CMD,INTVEC_CMD,INT_CMD,INT_CMD,RING_CMD};
    2044       if (!iiCheckTypes(h,t,1)) return TRUE;
    2045       if (((intvec*) h->next->Data())->length() != currRing->N &&
    2046           ((intvec*) h->next->next->Data())->length() != currRing->N)
    2047       {
    2048         Werror("system(\"Mfrwalk\" ...) intvecs not of length %d\n",currRing->N);
    20492032        return TRUE;
    20502033      }
     
    20532036      intvec* arg3 = (intvec*) h->next->next->Data();
    20542037      int arg4 = (int)(long) h->next->next->next->Data();
    2055       ideal result = (ideal) Mfrwalk(arg1, arg2, arg3, arg4);
     2038      int arg5 = (int)(long) h->next->next->next->next->Data();
     2039      ideal result = (ideal) Mfwalk(arg1, arg2, arg3, arg4, arg5);
    20562040      res->rtyp = IDEAL_CMD;
    20572041      res->data =  result;
     
    20592043    }
    20602044    else
     2045  #endif
     2046  /*==================== Mfrwalk =================*/
     2047  #ifdef HAVE_WALK
     2048    if (strcmp(sys_cmd, "Mfrwalk") == 0)
     2049    {
     2050      const short t[]={6,IDEAL_CMD,INTVEC_CMD,INTVEC_CMD,INT_CMD,INT_CMD,INT_CMD};
     2051      if (!iiCheckTypes(h,t,1)) return TRUE;
     2052/*
     2053      if (((intvec*) h->next->Data())->length() != currRing->N &&
     2054          ((intvec*) h->next->next->Data())->length() != currRing->N)
     2055      {
     2056        Werror("system(\"Mfrwalk\" ...) intvecs not of length %d\n",currRing->N);
     2057        return TRUE;
     2058      }
     2059*/
     2060      if((((intvec*) h->next->Data())->length() != currRing->N &&
     2061         ((intvec*) h->next->next->Data())->length() != currRing->N ) &&
     2062         (((intvec*) h->next->Data())->length() != (currRing->N)*(currRing->N) &&
     2063         ((intvec*) h->next->next->Data())->length() != (currRing->N)*(currRing->N) ))
     2064      {
     2065        Werror("system(\"Mfrwalk\" ...) intvecs not of length %d or %d\n",
     2066               currRing->N,(currRing->N)*(currRing->N));
     2067        return TRUE;
     2068      }
     2069
     2070      ideal arg1 = (ideal) h->Data();
     2071      intvec* arg2 = (intvec*) h->next->Data();
     2072      intvec* arg3 = (intvec*) h->next->next->Data();
     2073      int arg4 = (int)(long) h->next->next->next->Data();
     2074      int arg5 = (int)(long) h->next->next->next->next->Data();
     2075      int arg6 = (int)(long) h->next->next->next->next->next->Data();
     2076      ideal result = (ideal) Mfrwalk(arg1, arg2, arg3, arg4, arg5, arg6);
     2077      res->rtyp = IDEAL_CMD;
     2078      res->data =  result;
     2079      return FALSE;
     2080    }
     2081    else
    20612082  /*==================== Mprwalk =================*/
    20622083    if (strcmp(sys_cmd, "Mprwalk") == 0)
    20632084    {
    2064       const short t[]={7,IDEAL_CMD,INTVEC_CMD,INTVEC_CMD,INT_CMD,INT_CMD,INT_CMD,RING_CMD};
     2085      const short t[]={9,IDEAL_CMD,INTVEC_CMD,INTVEC_CMD,INT_CMD,INT_CMD,INT_CMD,INT_CMD,INT_CMD,INT_CMD};
    20652086      if (!iiCheckTypes(h,t,1)) return TRUE;
    2066       if (((intvec*) h->next->Data())->length() != currRing->N &&
    2067           ((intvec*) h->next->next->Data())->length() != currRing->N )
    2068       {
    2069         Werror("system(\"Mrwalk\" ...) intvecs not of length %d\n",
    2070                currRing->N);
     2087      if((((intvec*) h->next->Data())->length() != currRing->N &&
     2088         ((intvec*) h->next->next->Data())->length() != currRing->N ) &&
     2089         (((intvec*) h->next->Data())->length() != (currRing->N)*(currRing->N) &&
     2090         ((intvec*) h->next->next->Data())->length() != (currRing->N)*(currRing->N) ))
     2091      {
     2092        Werror("system(\"Mrwalk\" ...) intvecs not of length %d or %d\n",
     2093               currRing->N,(currRing->N)*(currRing->N));
    20712094        return TRUE;
    20722095      }
     
    20772100      int arg5 = (int)(long) h->next->next->next->next->Data();
    20782101      int arg6 = (int)(long) h->next->next->next->next->next->Data();
    2079       ring arg7 = (ring) h->next->next->next->next->next->next->Data();
    2080       ideal result = (ideal) Mprwalk(arg1, arg2, arg3, arg4, arg5, arg6, arg7);
     2102      int arg7 = (int)(long) h->next->next->next->next->next->next->Data();
     2103      int arg8 = (int)(long) h->next->next->next->next->next->next->next->Data();
     2104      int arg9 = (int)(long) h->next->next->next->next->next->next->next->next->Data();
     2105      ideal result = (ideal) Mprwalk(arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9);
    20812106      res->rtyp = IDEAL_CMD;
    20822107      res->data =  result;
  • Singular/fevoices.cc

    r38301d r47362dc  
    1717#include <Singular/subexpr.h>
    1818#include <Singular/ipshell.h>
     19#include <Singular/sdb.h>
    1920
    2021#include <stdlib.h>
  • Singular/ipassign.cc

    r38301d r47362dc  
    844844  if (res->data!=NULL) idDelete((ideal*)&res->data);
    845845  matrix m=(matrix)a->CopyD(MATRIX_CMD);
     846  if (TEST_V_ALLWARN)
     847    if (MATROWS(m)>1)
     848      Warn("assign matrix with %d rows to an ideal in >>%s<<",MATROWS(m),my_yylinebuf);
    846849  IDELEMS((ideal)m)=MATROWS(m)*MATCOLS(m);
    847850  ((ideal)m)->rank=1;
  • Singular/ipshell.cc

    r38301d r47362dc  
    6565#include <Singular/subexpr.h>
    6666#include <Singular/fevoices.h>
     67#include <Singular/sdb.h>
    6768
    6869#include <math.h>
     
    10291030void iiDebug()
    10301031{
     1032#ifdef HAVE_SDB
     1033  sdb_flags=1;
     1034#endif
    10311035  Print("\n-- break point in %s --\n",VoiceName());
    10321036  if (iiDebugMarker) VoiceBackTrack();
  • Singular/scanner.cc

    r38301d r47362dc  
    23302330  //yy_flush_buffer((YY_BUFFER_STATE)oldb);
    23312331}
     2332
     2333void my_yy_flush() { YY_FLUSH_BUFFER;BEGIN(0); }
     2334
  • Singular/scanner.ll

    r38301d r47362dc  
    380380  //yy_flush_buffer((YY_BUFFER_STATE)oldb);
    381381}
     382
     383void my_yy_flush() { YY_FLUSH_BUFFER;BEGIN(0); }
  • Singular/sdb.cc

    r38301d r47362dc  
    244244          "p <var> - display type and value of the variable <var>\n"
    245245          "q <flags> - quit debugger, set debugger flags(0,1,2)\n"
     246          "   0: stop debug, 1:continue, 2: throw an error, return to toplevel\n"
    246247          "Q - quit Singular\n");
    247248          int i;
  • Singular/tesths.cc

    r38301d r47362dc  
    146146  {
    147147    if (feOptValue(FE_OPT_SORT)) On(SW_USE_NTL_SORT);
    148 #ifdef HAVE_SDB
    149     sdb_flags = 0;
    150 #endif
    151148    dup2(1,2);
    152149    /* alternative:
  • Singular/walk.cc

    • Property mode changed from 100644 to 100755
    r4132ee r47362dc  
    2727
    2828#define FIRST_STEP_FRACTAL // to define the first step of the fractal
    29 //#define MSTDCC_FRACTAL // apply Buchberger alg to compute a red GB, if
     29#define MSTDCC_FRACTAL // apply Buchberger alg to compute a red GB, if
    3030//                          tau doesn't stay in the correct cone
    3131
     
    4242#include <Singular/ipshell.h>
    4343#include <Singular/ipconv.h>
     44#include <coeffs/ffields.h>
    4445#include <coeffs/coeffs.h>
    4546#include <Singular/subexpr.h>
     47#include <polys/templates/p_Procs.h>
    4648
    4749#include <polys/monomials/maps.h>
     
    429431#endif
    430432
    431 #ifdef CHECK_IDEAL_MWALK
     433//#ifdef CHECK_IDEAL_MWALK
    432434static void idString(ideal L, const char* st)
    433435{
     
    441443  Print(" %s;", pString(L->m[nL-1]));
    442444}
    443 #endif
     445//#endif
    444446
    445447#if defined(CHECK_IDEAL_MWALK) || defined(ENDWALKS)
     
    511513
    512514//unused
    513 #if 0
     515//#if 0
    514516static void MivString(intvec* iva, intvec* ivb, intvec* ivc)
    515517{
     
    534536  Print("%d)", (*ivc)[nV]);
    535537}
    536 #endif
     538//#endif
    537539
    538540/********************************************************************
     
    558560  }
    559561  return p0;
     562}
     563
     564/*****************************************************************************
     565 * compute the gcd of the entries of the vectors curr_weight and diff_weight *
     566 *****************************************************************************/
     567static int simplify_gcd(intvec* curr_weight, intvec* diff_weight)
     568{
     569  int j;
     570  int nRing = currRing->N;
     571  int gcd_tmp = (*curr_weight)[0];
     572  for (j=1; j<nRing; j++)
     573  {
     574    gcd_tmp = gcd(gcd_tmp, (*curr_weight)[j]);
     575    if(gcd_tmp == 1)
     576    {
     577      break;
     578    }
     579  }
     580  if(gcd_tmp != 1)
     581  {
     582    for (j=0; j<nRing; j++)
     583    {
     584    gcd_tmp = gcd(gcd_tmp, (*diff_weight)[j]);
     585    if(gcd_tmp == 1)
     586      {
     587        break;
     588      }
     589    }
     590  }
     591  return gcd_tmp;
    560592}
    561593
     
    774806  for(i=nG-1; i>=0; i--)
    775807  {
    776     mi = MpolyInitialForm(G->m[i], iv);
    777     gi = G->m[i];
    778 
     808    mi = pHead(MpolyInitialForm(G->m[i], iv));
     809    //Print("\n **// test_w_in_ConeCC: lm(initial)= %s \n",pString(mi));
     810    gi = pHead(G->m[i]);
     811    //Print("\n **// test_w_in_ConeCC: lm(ideal)= %s \n",pString(gi));
    779812    if(mi == NULL)
    780813    {
     
    953986}
    954987
    955 /*****************************************************************************
    956 * create a weight matrix order as intvec of an extra weight vector (a(iv),lp)*
    957 ******************************************************************************/
     988/*********************************************************************************
     989* create a weight matrix order as intvec of an extra weight vector (a(iv),M(iw)) *
     990**********************************************************************************/
    958991intvec* MivMatrixOrderRefine(intvec* iv, intvec* iw)
    959992{
    960   assume(iv->length() == iw->length());
    961   int i, nR = iv->length();
    962 
     993  assume((iv->length())*(iv->length()) == iw->length());
     994  int i,j, nR = iv->length();
     995 
    963996  intvec* ivm = new intvec(nR*nR);
    964997
     
    966999  {
    9671000    (*ivm)[i] = (*iv)[i];
    968     (*ivm)[i+nR] = (*iw)[i];
    969   }
    970   for(i=2; i<nR; i++)
    971   {
    972     (*ivm)[i*nR+i-2] = 1;
     1001  }
     1002  for(i=1; i<nR; i++)
     1003  {
     1004    for(j=0; j<nR; j++)
     1005    {
     1006      (*ivm)[j+i*nR] = (*iw)[j+i*nR];
     1007    }
    9731008  }
    9741009  return ivm;
     
    16121647    (* result)[i] = mpz_get_si(pert_vector[i]);
    16131648  }
    1614 
     1649/*
    16151650  j = 0;
    1616   for(i=0; i<nV; i++)
     1651  for(i=0; i<niv; i++)
    16171652  {
    16181653    (* result1)[i] = mpz_get_si(pert_vector[i]);
     
    16241659    }
    16251660  }
    1626   if(j > nV - 1)
     1661  if(j > niv - 1)
    16271662  {
    16281663    // Print("\n//  MfPertwalk: geaenderter vector gleich Null! \n");
     
    16301665    goto CHECK_OVERFLOW;
    16311666  }
    1632 
     1667/*
    16331668// check that the perturbed weight vector lies in the Groebner cone
     1669  Print("\n========================================\n//** MfPertvector: test in Cone.\n");
    16341670  if(test_w_in_ConeCC(G,result1) != 0)
    16351671  {
     
    16471683    // Print("\n// Mfpertwalk: geaenderter vector liegt nicht in Groebnerkegel! \n");
    16481684  }
    1649 
     1685  Print("\n ========================================\n");*/
    16501686  CHECK_OVERFLOW:
    16511687
     
    18591895}
    18601896
     1897
     1898/**************************************************************
     1899 * Look for the position of the smallest absolut value in vec *
     1900 **************************************************************/
     1901static int MivAbsMaxArg(intvec* vec)
     1902{
     1903  int k = MivAbsMax(vec);
     1904  int i=0;
     1905  while(1)
     1906  {
     1907    if((*vec)[i] == k || (*vec)[i] == -k)
     1908    {
     1909      break;
     1910    }
     1911    i++;
     1912  }
     1913  return i;
     1914}
     1915
     1916
    18611917/**********************************************************************
    18621918 * Compute a next weight vector between curr_weight and target_weight *
     
    18731929
    18741930  int nRing = currRing->N;
    1875   int checkRed, j, kkk, nG = IDELEMS(G);
     1931  int checkRed, j, nG = IDELEMS(G);
    18761932  intvec* ivtemp;
    18771933
     
    19111967  mpz_init(dcw);
    19121968
    1913   //int tn0, tn1, tz1, ncmp, gcd_tmp, ntmp;
    19141969  int gcd_tmp;
    19151970  intvec* diff_weight = MivSub(target_weight, curr_weight);
     
    19171972  intvec* diff_weight1 = MivSub(target_weight, curr_weight);
    19181973  poly g;
    1919   //poly g, gw;
     1974
    19201975  for (j=0; j<nG; j++)
    19211976  {
     
    19341989        mpz_set_si(MwWd, MivDotProduct(ivtemp, curr_weight));
    19351990        mpz_sub(s_zaehler, deg_w0_p1, MwWd);
    1936 
    19371991        if(mpz_cmp(s_zaehler, t_null) != 0)
    19381992        {
    19391993          mpz_set_si(MwWd, MivDotProduct(ivtemp, diff_weight));
    19401994          mpz_sub(s_nenner, MwWd, deg_d0_p1);
    1941 
    19421995          // check for 0 < s <= 1
    19431996          if( (mpz_cmp(s_zaehler,t_null) > 0 &&
     
    19792032    }
    19802033  }
    1981 //Print("\n// Alloc Size = %d \n", nRing*sizeof(mpz_t));
     2034  //Print("\n// Alloc Size = %d \n", nRing*sizeof(mpz_t));
    19822035  mpz_t *vec=(mpz_t*)omAlloc(nRing*sizeof(mpz_t));
    19832036
    19842037
    1985   // there is no 0<t<1 and define the next weight vector that is equal to the current weight vector
     2038  // there is no 0<t<1 and define the next weight vector that is equal
     2039  // to the current weight vector
    19862040  if(mpz_cmp(t_nenner, t_null) == 0)
    19872041  {
     
    20542108#endif
    20552109
    2056   //  BOOLEAN isdwpos;
    2057 
    2058   // construct a new weight vector
     2110// construct a new weight vector and check whether vec[j] is overflow,
     2111// i.e. vec[j] > 2^31.
     2112// If vec[j] doesn't overflow, define a weight vector. Otherwise,
     2113// report that overflow appears. In the second case, test whether the
     2114// the correctness of the new vector plays an important role
     2115
    20592116  for (j=0; j<nRing; j++)
    20602117  {
     
    21002157    }
    21012158  }
    2102 
     2159  // reduce the vector with the gcd
     2160  if(mpz_cmp_si(ggt,1) != 0)
     2161  {
     2162    for (j=0; j<nRing; j++)
     2163    {
     2164      mpz_divexact(vec[j], vec[j], ggt);
     2165    }
     2166  }
    21032167#ifdef  NEXT_VECTORS_CC
    21042168  PrintS("\n// gcd of elements of the vector: ");
     
    21062170#endif
    21072171
    2108 /**********************************************************************
    2109 * construct a new weight vector and check whether vec[j] is overflow, *
    2110 * i.e. vec[j] > 2^31.                                                 *
    2111 * If vec[j] doesn't overflow, define a weight vector. Otherwise,      *
    2112 * report that overflow appears. In the second case, test whether the  *
    2113 * the correctness of the new vector plays an important role           *
    2114 **********************************************************************/
    2115   kkk=0;
    21162172  for(j=0; j<nRing; j++)
    21172173    {
     
    21292185
    21302186  REDUCTION:
     2187  checkRed = 1;
    21312188  for (j=0; j<nRing; j++)
    21322189  {
    2133     (*diff_weight)[j] = mpz_get_si(vec[j]);
    2134   }
    2135   while(MivAbsMax(diff_weight) >= 5)
    2136   {
    2137     for (j=0; j<nRing; j++)
    2138     {
    2139       if(mpz_cmp_si(ggt,1)==0)
    2140       {
    2141         (*diff_weight1)[j] = floor(0.1*(*diff_weight)[j] + 0.5);
    2142         // Print("\n//  vector[%d] = %d \n",j+1, (*diff_weight1)[j]);
    2143       }
    2144       else
    2145       {
    2146         mpz_divexact(vec[j], vec[j], ggt);
    2147         (*diff_weight1)[j] = floor(0.1*(*diff_weight)[j] + 0.5);
    2148         // Print("\n//  vector[%d] = %d \n",j+1, (*diff_weight1)[j]);
    2149       }
    2150 /*
    2151       if((*diff_weight1)[j] == 0)
    2152       {
    2153         kkk = kkk + 1;
    2154       }
    2155 */
    2156     }
    2157 
    2158 
    2159 /*
    2160     if(kkk > nRing - 1)
    2161     {
    2162       // diff_weight was reduced to zero
    2163       // Print("\n // MwalkNextWeightCC: geaenderter Vector gleich Null! \n");
    2164       goto TEST_OVERFLOW;
    2165     }
    2166 */
    2167 
    2168     if(test_w_in_ConeCC(G,diff_weight1) != 0)
    2169     {
    2170       Print("\n// MwalkNextWeightCC: geaenderter vector liegt in Groebnerkegel! \n");
    2171       for (j=0; j<nRing; j++)
    2172       {
    2173         (*diff_weight)[j] = (*diff_weight1)[j];
    2174       }
    2175       if(MivAbsMax(diff_weight) < 5)
    2176       {
    2177         checkRed = 1;
    2178         goto SIMPLIFY_GCD;
    2179       }
    2180     }
    2181     else
    2182     {
    2183       // Print("\n// MwalkNextWeightCC: geaenderter vector liegt nicht in Groebnerkegel! \n");
    2184       break;
    2185     }
     2190    (*diff_weight1)[j] = mpz_get_si(vec[j]);
     2191  }
     2192  while(test_w_in_ConeCC(G,diff_weight1))
     2193  {
     2194    for(j=0; j<nRing; j++)
     2195    {
     2196      (*diff_weight)[j] = (*diff_weight1)[j];
     2197      mpz_set_si(vec[j], (*diff_weight)[j]);     
     2198    }
     2199    for(j=0; j<nRing; j++)
     2200    {
     2201      (*diff_weight1)[j] = floor(0.1*(*diff_weight)[j] + 0.5);
     2202    }
     2203  }
     2204  if(MivAbsMax(diff_weight)>10000)
     2205  {
     2206    for(j=0; j<nRing; j++)
     2207    {
     2208      (*diff_weight1)[j] = (*diff_weight)[j];
     2209    }
     2210    j = 0;
     2211    while(test_w_in_ConeCC(G,diff_weight1))
     2212    {
     2213      (*diff_weight)[j] = (*diff_weight1)[j];
     2214      mpz_set_si(vec[j], (*diff_weight)[j]);
     2215      j = MivAbsMaxArg(diff_weight1);
     2216      (*diff_weight1)[j] = floor(0.1*(*diff_weight1)[j] + 0.5);
     2217    }
     2218    goto SIMPLIFY_GCD;
    21862219  }
    21872220
     
    22222255   mpz_clear(t_null);
    22232256
    2224 
    2225 
    22262257  if(Overflow_Error == FALSE)
    22272258  {
    22282259    Overflow_Error = nError;
    22292260  }
    2230  rComplete(currRing);
    2231   for(kkk=0; kkk<IDELEMS(G);kkk++)
    2232   {
    2233     poly p=G->m[kkk];
     2261  rComplete(currRing);
     2262  for(j=0; j<IDELEMS(G); j++)
     2263  {
     2264    poly p=G->m[j];
    22342265    while(p!=NULL)
    22352266    {
     
    22712302}
    22722303
    2273 /**************************************************************
     2304/********************************************************************
    22742305 * define and execute a new ring which order is (a(vb),a(va),lp,C)  *
    2275  * ************************************************************/
    2276 #if 0
    2277 // unused
     2306 * ******************************************************************/
    22782307static void VMrHomogeneous(intvec* va, intvec* vb)
    22792308{
     
    23532382  rChangeCurrRing(r);
    23542383}
    2355 #endif
     2384
    23562385
    23572386/**************************************************************
     
    24282457}
    24292458
    2430 static ring VMrDefault1(intvec* va)
    2431 {
    2432 
    2433   if ((currRing->ppNoether)!=NULL)
    2434   {
    2435     pDelete(&(currRing->ppNoether));
    2436   }
    2437   if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) ||
    2438       ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data))))
    2439   {
    2440     sLastPrinted.CleanUp();
    2441   }
    2442 
    2443   ring r = (ring) omAlloc0Bin(sip_sring_bin);
    2444   int i, nv = currRing->N;
    2445 
    2446   r->cf  = currRing->cf;
    2447   r->N   = currRing->N;
    2448 
    2449   int nb = 4;
    2450 
    2451   //names
    2452   char* Q; // In order to avoid the corrupted memory, do not change.
    2453   r->names = (char **) omAlloc0(nv * sizeof(char_ptr));
    2454   for(i=0; i<nv; i++)
    2455   {
    2456     Q = currRing->names[i];
    2457     r->names[i]  = omStrDup(Q);
    2458   }
    2459 
    2460   /*weights: entries for 3 blocks: NULL Made:???*/
    2461   r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
    2462   r->wvhdl[0] = (int*) omAlloc(nv*sizeof(int));
    2463   for(i=0; i<nv; i++)
    2464     r->wvhdl[0][i] = (*va)[i];
    2465 
    2466   /* order: a,lp,C,0 */
    2467   r->order = (int *) omAlloc(nb * sizeof(int *));
    2468   r->block0 = (int *)omAlloc0(nb * sizeof(int *));
    2469   r->block1 = (int *)omAlloc0(nb * sizeof(int *));
    2470 
    2471   // ringorder a for the first block: var 1..nv
    2472   r->order[0]  = ringorder_a;
    2473   r->block0[0] = 1;
    2474   r->block1[0] = nv;
    2475 
    2476   // ringorder lp for the second block: var 1..nv
    2477   r->order[1]  = ringorder_lp;
    2478   r->block0[1] = 1;
    2479   r->block1[1] = nv;
    2480 
    2481   // ringorder C for the third block
    2482   // it is very important within "idLift",
    2483   // especially, by ring syz_ring=rCurrRingAssure_SyzComp();
    2484   // therefore, nb  must be (nBlocks(currRing)  + 1)
    2485   r->order[2]  = ringorder_C;
    2486 
    2487   // the last block: everything is 0
    2488   r->order[3]  = 0;
    2489 
    2490   // polynomial ring
    2491   r->OrdSgn    = 1;
    2492 
    2493   // complete ring intializations
    2494 
    2495   rComplete(r);
    2496 
    2497   //rChangeCurrRing(r);
    2498   return r;
    2499 }
    2500 
    25012459/****************************************************************
    25022460 * define and execute a new ring with ordering (a(va),Wp(vb),C) *
    25032461 * **************************************************************/
    2504 
    25052462static ring VMrRefine(intvec* va, intvec* vb)
    25062463{
     
    25762533
    25772534  // complete ring intializations
    2578 
     2535 
    25792536  rComplete(r);
    25802537
     
    28062763}
    28072764
    2808 
    2809 /* define a ring with parameters und change to it */
    2810 /* DefRingPar and DefRingParlp corrupt still memory */
     2765/***************************************************
     2766* define a ring with parameters und change to it   *
     2767* DefRingPar and DefRingParlp corrupt still memory *
     2768****************************************************/
    28112769static void DefRingPar(intvec* va)
    28122770{
     
    29562914}
    29572915
    2958 //unused
    2959 /**************************************************************
    2960  * check wheather one or more components of a vector are zero *
    2961  **************************************************************/
    2962 #if 0
     2916/*************************************************************
     2917 * check whether one or more components of a vector are zero *
     2918 *************************************************************/
    29632919static int isNolVector(intvec* hilb)
    29642920{
     
    29732929  return 0;
    29742930}
    2975 #endif
     2931
     2932/*************************************************************
     2933 * check whether one or more components of a vector are <= 0 *
     2934 *************************************************************/
     2935static int isNegNolVector(intvec* hilb)
     2936{
     2937  int i;
     2938  for(i=hilb->length()-1; i>=0; i--)
     2939  {
     2940    if((* hilb)[i]<=0)
     2941    {
     2942      return 1;
     2943    }
     2944  }
     2945  return 0;
     2946}
     2947
     2948/**************************************************************************
     2949* Gomega is the initial ideal of G w. r. t. the current weight vector     *
     2950* curr_weight. Check whether curr_weight lies on a border of the Groebner *
     2951* cone, i. e. check whether a monomial is divisible by a leading monomial *
     2952***************************************************************************/
     2953static ideal middleOfCone(ideal G, ideal Gomega)
     2954{
     2955  BOOLEAN middle = FALSE;
     2956  int i,j,N = IDELEMS(Gomega);
     2957  poly p,lm,factor1,factor2;
     2958  //PrintS("\n//** idCopy\n");
     2959  ideal Go = idCopy(G);
     2960 
     2961  //PrintS("\n//** jetzt for-Loop!\n");
     2962
     2963  // check whether leading monomials of G and Gomega coincide
     2964  // and return NULL if not
     2965  for(i=0; i<N; i++)
     2966  {
     2967    p = pCopy(Gomega->m[i]);
     2968    lm = pCopy(pHead(G->m[i]));
     2969    if(!pIsConstant(pSub(p,lm)))
     2970    {
     2971      //pDelete(&p);
     2972      //pDelete(&lm);
     2973      idDelete(&Go);
     2974      return NULL;
     2975    }
     2976    //pDelete(&p);
     2977    //pDelete(&lm);
     2978  }
     2979  for(i=0; i<N; i++)
     2980  {
     2981    for(j=0; j<N; j++)
     2982    {
     2983      if(i!=j)
     2984      {
     2985        p = pCopy(Gomega->m[i]);
     2986        lm = pCopy(Gomega->m[j]);
     2987        pIter(p);
     2988        while(p!=NULL)
     2989        {
     2990          if(pDivisibleBy(lm,p))
     2991          {
     2992            if(middle == FALSE)
     2993            {
     2994              middle = TRUE;
     2995            }
     2996            factor1 = singclap_pdivide(pHead(p),lm,currRing);
     2997            factor2 = pMult(pCopy(factor1),pCopy(Go->m[j]));
     2998            pDelete(&factor1);
     2999            Go->m[i] = pAdd((Go->m[i]),pNeg(pCopy(factor2)));
     3000            pDelete(&factor2);
     3001          }
     3002          pIter(p);
     3003        }
     3004        pDelete(&lm);
     3005        pDelete(&p);
     3006      }
     3007    }
     3008  }
     3009 
     3010  //PrintS("\n//** jetzt Delete!\n");
     3011  //pDelete(&p);
     3012  //pDelete(&factor);
     3013  //pDelete(&lm);
     3014  if(middle == TRUE)
     3015  {
     3016    //PrintS("\n//** middle TRUE!\n");
     3017    return Go;
     3018  }
     3019  //PrintS("\n//** middle FALSE!\n");
     3020  idDelete(&Go);
     3021  return NULL;
     3022}
    29763023
    29773024/******************************  Februar 2002  ****************************
     
    31283175    else
    31293176    {
    3130       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung
     3177      rChangeCurrRing(VMrDefault(curr_weight));
    31313178    }
    31323179    newRing = currRing;
     
    32803327  }
    32813328  return 0;
     3329}
     3330
     3331/*****************************************
     3332 * return maximal polynomial length of G *
     3333 *****************************************/
     3334static int maxlengthpoly(ideal G)
     3335{
     3336  int i,k,length=0;
     3337  for(i=IDELEMS(G)-1; i>=0; i--)
     3338  {
     3339    k = pLength(G->m[i]);
     3340    if(k>length)
     3341    {
     3342      length = k;
     3343    }
     3344  }
     3345  return length;
    32823346}
    32833347
     
    38843948    else
    38853949    {
    3886       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung
     3950      rChangeCurrRing(VMrDefault(curr_weight));
    38873951    }
    38883952    newRing = currRing;
     
    41454209    else
    41464210    {
    4147       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung
     4211      rChangeCurrRing(VMrDefault(curr_weight));
    41484212    }
    41494213    newRing = currRing;
     
    42874351intvec* Xivlp;
    42884352
    4289 #if 0
     4353
    42904354/********************************
    42914355 * compute a next weight vector *
    42924356 ********************************/
    4293 static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight, intvec* target_weight, int weight_rad, int pert_deg)
     4357static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight,
     4358               intvec* target_weight, int weight_rad, int pert_deg)
    42944359{
    4295   int i, weight_norm;
    4296   int nV = currRing->N;
     4360  assume(currRing != NULL && curr_weight != NULL &&
     4361         target_weight != NULL && G->m[0] != NULL);
     4362
     4363  int i,weight_norm,nV = currRing->N;
    42974364  intvec* next_weight2;
    42984365  intvec* next_weight22 = new intvec(nV);
    4299   intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight, G);
    4300   if(MivComp(next_weight, target_weight) == 1)
    4301   {
    4302     return(next_weight);
    4303   }
    4304   else
    4305   {
    4306     //compute a perturbed next weight vector "next_weight1"
    4307     intvec* next_weight1 = MkInterRedNextWeight(MPertVectors(G, MivMatrixOrder(curr_weight), pert_deg), target_weight, G);
    4308     //Print("\n // size of next_weight1 = %d", sizeof((*next_weight1)));
    4309 
    4310     //compute a random next weight vector "next_weight2"
    4311     while(1)
    4312     {
    4313       weight_norm = 0;
    4314       while(weight_norm == 0)
     4366  intvec* result = new intvec(nV);
     4367
     4368  intvec* next_weight1 =MkInterRedNextWeight(curr_weight,target_weight,G);
     4369  //compute a random next weight vector "next_weight2"
     4370  while(1)
     4371  {
     4372    weight_norm = 0;
     4373    while(weight_norm == 0)
     4374    {
     4375
     4376      for(i=0; i<nV; i++)
     4377      {
     4378        (*next_weight22)[i] = rand() % 60000 - 30000;
     4379        weight_norm = weight_norm + (*next_weight22)[i]*(*next_weight22)[i];
     4380      }
     4381      weight_norm = 1 + floor(sqrt(weight_norm));
     4382    }
     4383
     4384    for(i=0; i<nV; i++)
     4385    {
     4386      if((*next_weight22)[i] < 0)
     4387      {
     4388        (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);
     4389      }
     4390      else
     4391      {
     4392        (*next_weight22)[i] = (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);
     4393      }
     4394    }
     4395
     4396    if(test_w_in_ConeCC(G, next_weight22) == 1)
     4397    {
     4398      next_weight2 = MkInterRedNextWeight(next_weight22,target_weight,G);
     4399      if(MivAbsMax(next_weight2)>1147483647)
    43154400      {
    43164401        for(i=0; i<nV; i++)
    43174402        {
    4318           //Print("\n//  next_weight[%d]  = %d", i, (*next_weight)[i]);
    4319           (*next_weight22)[i] = rand() % 60000 - 30000;
    4320           weight_norm = weight_norm + (*next_weight22)[i]*(*next_weight22)[i];
     4403          (*next_weight22)[i] = (*next_weight2)[i];
    43214404        }
    4322         weight_norm = 1 + floor(sqrt(weight_norm));
    4323       }
    4324 
    4325       for(i=nV-1; i>=0; i--)
    4326       {
    4327         if((*next_weight22)[i] < 0)
     4405        i = 0;
     4406        /* reduce the size of the maximal entry of the vector*/
     4407        while(test_w_in_ConeCC(G,next_weight22))
    43284408        {
    4329           (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);
     4409          (*next_weight2)[i] = (*next_weight22)[i];
     4410          i = MivAbsMaxArg(next_weight22);
     4411          (*next_weight22)[i] = floor(0.1*(*next_weight22)[i] + 0.5);
    43304412        }
    4331         else
    4332         {
    4333           (*next_weight22)[i] = (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);
    4334         }
    4335       //Print("\n//  next_weight22[%d]  = %d", i, (*next_weight22)[i]);
    4336       }
    4337 
    4338       if(test_w_in_ConeCC(G, next_weight22) == 1)
    4339       {
    4340         //Print("\n//MWalkRandomNextWeight: next_weight2 im Kegel\n");
    4341         next_weight2 = MkInterRedNextWeight(next_weight22, target_weight, G);
    4342         delete next_weight22;
    4343         break;
    4344       }
    4345     }
    4346     intvec* result = new intvec(nV);
    4347     ideal G_test = MwalkInitialForm(G, next_weight);
     4413      }
     4414      delete next_weight22;
     4415      break;
     4416    }
     4417  }
     4418 
     4419  // compute "usual" next weight vector
     4420  intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight, G);
     4421  ideal G_test = MwalkInitialForm(G, next_weight);
     4422  ideal G_test2 = MwalkInitialForm(G, next_weight2);
     4423
     4424  // compare next weights
     4425  if(Overflow_Error == FALSE)
     4426  {
    43484427    ideal G_test1 = MwalkInitialForm(G, next_weight1);
    4349     ideal G_test2 = MwalkInitialForm(G, next_weight2);
    4350 
    4351     // compare next_weights
    4352     if(IDELEMS(G_test1) < IDELEMS(G_test))
    4353     {
    4354       if(IDELEMS(G_test2) <= IDELEMS(G_test1)) // |G_test2| <= |G_test1| < |G_test|
     4428    if(G_test1->m[0] != NULL && maxlengthpoly(G_test1) < maxlengthpoly(G_test))//if(IDELEMS(G_test1) < IDELEMS(G_test))
     4429    {
     4430      if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test1))//if(IDELEMS(G_test2) < IDELEMS(G_test1))
    43554431      {
    43564432        for(i=0; i<nV; i++)
     
    43594435        }
    43604436      }
    4361       else // |G_test1| < |G_test|, |G_test1| < |G_test2|
     4437      else
    43624438      {
    43634439        for(i=0; i<nV; i++)
     
    43654441          (*result)[i] = (*next_weight1)[i];
    43664442        }
    4367       }
     4443      }   
    43684444    }
    43694445    else
    43704446    {
    4371       if(IDELEMS(G_test2) <= IDELEMS(G_test)) // |G_test2| <= |G_test| <= |G_test1|
     4447      if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test))//if(IDELEMS(G_test2) < IDELEMS(G_test)) // |G_test2| < |G_test| <= |G_test1|
    43724448      {
    43734449        for(i=0; i<nV; i++)
     
    43764452        }
    43774453      }
    4378       else // |G_test| <= |G_test1|, |G_test| < |G_test2|
     4454      else
    43794455      {
    43804456        for(i=0; i<nV; i++)
     
    43844460      }
    43854461    }
    4386     delete next_weight;
    4387     delete next_weight1;
    4388     idDelete(&G_test);
    43894462    idDelete(&G_test1);
    4390     idDelete(&G_test2);
    4391     if(test_w_in_ConeCC(G, result) == 1)
    4392     {
    4393       delete next_weight2;
    4394       return result;
     4463  }
     4464  else
     4465  {
     4466    Overflow_Error = FALSE;
     4467    if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test))//if(IDELEMS(G_test2) < IDELEMS(G_test))
     4468    {
     4469      for(i=1; i<nV; i++)
     4470      {
     4471        (*result)[i] = (*next_weight2)[i];
     4472      }
    43954473    }
    43964474    else
    43974475    {
    4398       delete result;
    4399       return next_weight2;
    4400     }
    4401   }
    4402 }
    4403 #endif
    4404 
    4405 /********************************
    4406  * compute a next weight vector *
    4407  ********************************/
    4408 static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight, intvec* target_weight, int weight_rad, int pert_deg)
    4409 {
    4410   int i, weight_norm;
    4411   //int randCount=0;
    4412   int nV = currRing->N;
    4413   intvec* next_weight2;
    4414   intvec* next_weight22 = new intvec(nV);
    4415   intvec* result = new intvec(nV);
    4416 
    4417   //compute a perturbed next weight vector "next_weight1"
    4418   //intvec* next_weight1 = MkInterRedNextWeight(MPertVectors(G,MivMatrixOrderRefine(curr_weight,target_weight),pert_deg),target_weight,G);
    4419   intvec* next_weight1 =MkInterRedNextWeight(curr_weight,target_weight,G);
    4420   //compute a random next weight vector "next_weight2"
    4421   while(1)
    4422   {
    4423     weight_norm = 0;
    4424     while(weight_norm == 0)
    4425     {
    44264476      for(i=0; i<nV; i++)
    44274477      {
    4428         (*next_weight22)[i] = rand() % 60000 - 30000;
    4429         weight_norm = weight_norm + (*next_weight22)[i]*(*next_weight22)[i];
    4430       }
    4431       weight_norm = 1 + floor(sqrt(weight_norm));
    4432     }
    4433     for(i=0; i<nV; i++)
    4434     {
    4435       if((*next_weight22)[i] < 0)
    4436       {
    4437         (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);
    4438       }
    4439       else
    4440       {
    4441         (*next_weight22)[i] = (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);
    4442       }
    4443     }
    4444     if(test_w_in_ConeCC(G, next_weight22) == 1)
    4445     {
    4446       next_weight2 = MkInterRedNextWeight(next_weight22,target_weight,G);
    4447       delete next_weight22;
    4448       break;
    4449     }
    4450   }
    4451   // compute "usual" next weight vector
    4452   intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight, G);
    4453   ideal G_test = MwalkInitialForm(G, next_weight);
    4454   ideal G_test2 = MwalkInitialForm(G, next_weight2);
    4455 
    4456   // compare next weights
    4457   if(Overflow_Error == FALSE)
    4458   {
    4459     ideal G_test1 = MwalkInitialForm(G, next_weight1);
    4460     if(IDELEMS(G_test1) < IDELEMS(G_test))
    4461     {
    4462       if(IDELEMS(G_test2) < IDELEMS(G_test1))
    4463       {
    4464         // |G_test2| < |G_test1| < |G_test|
    4465         for(i=0; i<nV; i++)
    4466         {
    4467           (*result)[i] = (*next_weight2)[i];
    4468         }
    4469       }
    4470       else
    4471       {
    4472         // |G_test1| < |G_test|, |G_test1| <= |G_test2|
    4473         for(i=0; i<nV; i++)
    4474         {
    4475           (*result)[i] = (*next_weight1)[i];
    4476         }
    4477       }
    4478     }
    4479     else
    4480     {
    4481       if(IDELEMS(G_test2) < IDELEMS(G_test)) // |G_test2| < |G_test| <= |G_test1|
    4482       {
    4483         for(i=0; i<nV; i++)
    4484         {
    4485           (*result)[i] = (*next_weight2)[i];
    4486         }
    4487       }
    4488       else
    4489       {
    4490         // |G_test| < |G_test1|, |G_test| <= |G_test2|
    4491         for(i=0; i<nV; i++)
    4492         {
    4493           (*result)[i] = (*next_weight)[i];
    4494         }
    4495       }
    4496     }
    4497     idDelete(&G_test1);
    4498   }
    4499   else
    4500   {
    4501     Overflow_Error = FALSE;
    4502     if(IDELEMS(G_test2) < IDELEMS(G_test))
    4503     {
    4504       for(i=1; i<nV; i++)
    4505       {
    4506         (*result)[i] = (*next_weight2)[i];
    4507       }
    4508     }
    4509     else
    4510     {
    4511       for(i=0; i<nV; i++)
    4512       {
    45134478        (*result)[i] = (*next_weight)[i];
    45144479      }
    45154480    }
    45164481  }
     4482  //PrintS("\n MWalkRandomNextWeight: Ende ok!\n");
    45174483  idDelete(&G_test);
    45184484  idDelete(&G_test2);
     
    45754541    else
    45764542    {
    4577       rChangeCurrRing(VMrDefault(orig_target_weight)); // Aenderung
     4543      rChangeCurrRing(VMrDefault(orig_target_weight));
    45784544    }
    45794545    TargetRing = currRing;
     
    46464612    else
    46474613    {
    4648       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung
     4614      rChangeCurrRing(VMrDefault(curr_weight));
    46494615    }
    46504616    newRing = currRing;
     
    47554721    else
    47564722    {
    4757       rChangeCurrRing(VMrDefault(orig_target_weight)); // Aenderung
     4723      rChangeCurrRing(VMrDefault(orig_target_weight));
    47584724    }
    47594725    F1 = idrMoveR(G, newRing,currRing);
     
    47864752      else
    47874753      {
    4788         rChangeCurrRing(VMrDefault(orig_target_weight)); // Aenderung
     4754        rChangeCurrRing(VMrDefault(orig_target_weight));
    47894755      }
    47904756    KSTD_Finish:
     
    48844850      tim = clock();
    48854851      /*
    4886         Print("\n// **** Grï¿œbnerwalk took %d steps and ", nwalk);
     4852        Print("\n// **** Groebnerwalk took %d steps and ", nwalk);
    48874853        PrintS("\n// **** call the rec. Pert. Walk to compute a red GB of:");
    48884854        idElements(Gomega, "G_omega");
     
    49144880      oldRing = currRing;
    49154881
    4916       /* create a new ring newRing */
     4882      // create a new ring newRing
    49174883       if (rParameter(currRing) != NULL)
    49184884       {
     
    49214887       else
    49224888       {
    4923          rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung
     4889         rChangeCurrRing(VMrDefault(curr_weight));
    49244890       }
    49254891      newRing = currRing;
     
    49474913      else
    49484914      {
    4949         rChangeCurrRing(VMrDefault(curr_weight));  //Aenderung
     4915        rChangeCurrRing(VMrDefault(curr_weight));
    49504916      }
    49514917      newRing = currRing;
     
    49594925      M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
    49604926      delete hilb_func;
    4961 #endif // BUCHBERGER_ALG
     4927#endif
    49624928      tstd = tstd + clock() - to;
    49634929
     
    49684934
    49694935      to = clock();
    4970       // compute a representation of the generators of submod (M) with respect to those of mod (Gomega). Gomega is a reduced Groebner basis w.r.t. the current ring.
     4936      // compute a representation of the generators of submod (M) with respect
     4937      // to those of mod (Gomega).
     4938      // Gomega is a reduced Groebner basis w.r.t. the current ring.
    49714939      F = MLifttwoIdeal(Gomega2, M1, G);
    49724940      tlift = tlift + clock() - to;
     
    50184986      else
    50194987      {
    5020         rChangeCurrRing(VMrDefault(target_weight)); // Aenderung
     4988        rChangeCurrRing(VMrDefault(target_weight));
    50214989      }
    50224990      F1 = idrMoveR(G, newRing,currRing);
     
    50655033 * THE GROEBNER WALK ALGORITHM *
    50665034 *******************************/
    5067 ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M, ring baseRing)
     5035ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M,
     5036            ring baseRing, int reduction, int printout)
    50685037{
    5069   BITSET save1 = si_opt_1; // save current options
    5070   //si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
    5071   //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
    5072   //si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB));
     5038  // save current options
     5039  BITSET save1 = si_opt_1;
     5040  if(reduction == 0)
     5041  {
     5042    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     5043    si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     5044  }
    50735045  Set_Error(FALSE);
    50745046  Overflow_Error = FALSE;
     5047  //BOOLEAN endwalks = FALSE;
    50755048#ifdef TIME_TEST
    50765049  clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
     
    50805053#endif
    50815054  nstep=0;
    5082   int i,nwalk,endwalks = 0;
     5055  int i,nwalk;
    50835056  int nV = baseRing->N;
    50845057
    5085   ideal Gomega, M, F, Gomega1, Gomega2, M1; //, F1;
     5058  ideal Gomega, M, F, FF, Gomega1, Gomega2, M1;
    50865059  ring newRing;
    50875060  ring XXRing = baseRing;
     5061  ring targetRing;
    50885062  intvec* ivNull = new intvec(nV);
    50895063  intvec* curr_weight = new intvec(nV);
    50905064  intvec* target_weight = new intvec(nV);
    50915065  intvec* exivlp = Mivlp(nV);
     5066/*
    50925067  intvec* tmp_weight = new intvec(nV);
    50935068  for(i=0; i<nV; i++)
    50945069  {
    5095     (*tmp_weight)[i] = (*target_M)[i];
    5096   }
     5070    (*tmp_weight)[i] = (*orig_M)[i];
     5071  }
     5072*/
    50975073  for(i=0; i<nV; i++)
    50985074  {
     
    51115087#endif
    51125088  rComplete(currRing);
    5113 #ifdef CHECK_IDEAL_MWALK
    5114     idString(Go,"Go");
    5115 #endif
     5089//#ifdef CHECK_IDEAL_MWALK
     5090  if(printout > 2)
     5091  {
     5092    idString(Go,"//** Mwalk: Go");
     5093  }
     5094//#endif
     5095
     5096  if(target_M->length() == nV)
     5097  {
     5098   // define the target ring
     5099    targetRing = VMrDefault(target_weight);
     5100  }
     5101  else
     5102  {
     5103    targetRing = VMatrDefault(target_M);
     5104  }
     5105  if(orig_M->length() == nV)
     5106  {
     5107    // define a new ring with ordering "(a(curr_weight),lp)
     5108    newRing = VMrDefault(curr_weight);
     5109  }
     5110  else
     5111  {
     5112    newRing = VMatrDefault(orig_M);
     5113  }
     5114  rChangeCurrRing(newRing);
     5115
    51165116#ifdef TIME_TEST
    51175117  to = clock();
    51185118#endif
    5119      if(orig_M->length() == nV)
    5120       {
    5121         newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
    5122       }
    5123       else
    5124       {
    5125         newRing = VMatrDefault(orig_M);
    5126       }
    5127   rChangeCurrRing(newRing);
     5119
    51285120  ideal G = MstdCC(idrMoveR(Go,baseRing,currRing));
    5129   baseRing = currRing;
     5121
    51305122#ifdef TIME_TEST
    51315123  tostd = clock()-to;
    51325124#endif
    51335125
     5126  baseRing = currRing;
    51345127  nwalk = 0;
     5128
    51355129  while(1)
    51365130  {
    51375131    nwalk ++;
    51385132    nstep ++;
     5133
    51395134#ifdef TIME_TEST
    51405135    to = clock();
    51415136#endif
    5142 #ifdef CHECK_IDEAL_MWALK
    5143     idString(G,"G");
    5144 #endif
    5145     Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
     5137    // compute an initial form ideal of <G> w.r.t. "curr_vector"
     5138    Gomega = MwalkInitialForm(G, curr_weight);
    51465139#ifdef TIME_TEST
    5147     tif = tif + clock()-to; //time for computing initial form ideal
    5148 #endif
    5149 #ifdef CHECK_IDEAL_MWALK
    5150     idString(Gomega,"Gomega");
    5151 #endif
     5140    tif = tif + clock()-to;
     5141#endif
     5142
     5143//#ifdef CHECK_IDEAL_MWALK
     5144    if(printout > 1)
     5145    {
     5146      idString(Gomega,"//** Mwalk: Gomega");
     5147    }
     5148//#endif
     5149
     5150    if(reduction == 0)
     5151    {
     5152      FF = middleOfCone(G,Gomega);
     5153      if( FF != NULL)
     5154      {
     5155        idDelete(&G);
     5156        G = idCopy(FF);
     5157        idDelete(&FF);
     5158        goto NEXT_VECTOR;
     5159      }
     5160    }
     5161
    51525162#ifndef  BUCHBERGER_ALG
    51535163    if(isNolVector(curr_weight) == 0)
    51545164    {
    51555165      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
    5156     }
     5166    }   
    51575167    else
    51585168    {
     
    51605170    }
    51615171#endif
     5172
    51625173    if(nwalk == 1)
    51635174    {
    51645175      if(orig_M->length() == nV)
    51655176      {
    5166         newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
     5177        // define a new ring with ordering "(a(curr_weight),lp)
     5178        newRing = VMrDefault(curr_weight);
    51675179      }
    51685180      else
     
    51755187     if(target_M->length() == nV)
    51765188     {
    5177        newRing = VMrRefine(curr_weight,target_weight); //define a new ring with ordering "(a(curr_weight),Wp(target_weight))"
     5189       //define a new ring with ordering "(a(curr_weight),lp)"
     5190       newRing = VMrDefault(curr_weight);
    51785191     }
    51795192     else
    51805193     {
     5194       //define a new ring with matrix ordering
    51815195       newRing = VMatrRefine(target_M,curr_weight);
    51825196     }
     
    51995213#endif
    52005214    idSkipZeroes(M);
    5201 #ifdef CHECK_IDEAL_MWALK
    5202     PrintS("\n//** Mwalk: computed M.\n");
    5203     idString(M, "M");
    5204 #endif
     5215//#ifdef CHECK_IDEAL_MWALK
     5216    if(printout > 2)
     5217    {
     5218      idString(M, "//** Mwalk: M");
     5219    }
     5220//#endif
    52055221    //change the ring to baseRing
    52065222    rChangeCurrRing(baseRing);
     
    52125228    to = clock();
    52135229#endif
    5214     // compute a representation of the generators of submod (M) with respect to those of mod (Gomega), where Gomega is a reduced Groebner basis w.r.t. the current ring
     5230    // compute a representation of the generators of submod (M) with respect to those of mod (Gomega),
     5231    // where Gomega is a reduced Groebner basis w.r.t. the current ring
    52155232    F = MLifttwoIdeal(Gomega2, M1, G);
     5233
    52165234#ifdef TIME_TEST
    52175235    tlift = tlift + clock() - to;
    52185236#endif
    5219 #ifdef CHECK_IDEAL_MWALK
    5220     idString(F, "F");
    5221 #endif
     5237//#ifdef CHECK_IDEAL_MWALK
     5238    if(printout > 2)
     5239    {
     5240      idString(F, "//** Mwalk: F");
     5241    }
     5242//#endif
    52225243    idDelete(&Gomega2);
    52235244    idDelete(&M1);
     5245
    52245246    rChangeCurrRing(newRing); // change the ring to newRing
    52255247    G = idrMoveR(F,baseRing,currRing);
    52265248    idDelete(&F);
     5249    idSkipZeroes(G);
     5250
     5251//#ifdef CHECK_IDEAL_MWALK
     5252    if(printout > 2)
     5253    {
     5254      idString(G, "//** Mwalk: G");
     5255    }
     5256//#endif
     5257
     5258    rChangeCurrRing(targetRing);
     5259    G = idrMoveR(G,newRing,currRing);
     5260    // test whether target cone is reached
     5261    if(reduction !=0 && test_w_in_ConeCC(G,curr_weight) == 1)
     5262    {
     5263      baseRing = currRing;
     5264      break;
     5265      //endwalks = TRUE;
     5266    }
     5267
     5268    rChangeCurrRing(newRing);
     5269    G = idrMoveR(G,targetRing,currRing);
    52275270    baseRing = currRing;
     5271
     5272/*
    52285273#ifdef TIME_TEST
    52295274    to = clock();
    52305275#endif
    5231     //G = kStd(F1,NULL,testHomog,NULL,NULL,0,0,NULL);
     5276
    52325277#ifdef TIME_TEST
    52335278    tstd = tstd + clock() - to;
    52345279#endif
    5235     idSkipZeroes(G);
    5236 #ifdef CHECK_IDEAL_MWALK
    5237     idString(G, "G");
    5238 #endif
     5280*/
     5281
     5282
    52395283#ifdef TIME_TEST
    52405284    to = clock();
    52415285#endif
     5286    NEXT_VECTOR:
    52425287    intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
    52435288#ifdef TIME_TEST
    52445289    tnw = tnw + clock() - to;
    52455290#endif
    5246 #ifdef PRINT_VECTORS
    5247     MivString(curr_weight, target_weight, next_weight);
    5248 #endif
    5249     if(MivComp(next_weight, ivNull) == 1 || MivComp(target_weight,curr_weight) == 1)// || test_w_in_ConeCC(G, target_weight) == 1 || MivComp(next_weight,curr_weight) == 1)
    5250     {
    5251 #ifdef CHECK_IDEAL_MWALK
    5252       PrintS("\n//** Mwalk: entering last cone.\n");
    5253 #endif
     5291//#ifdef PRINT_VECTORS
     5292    if(printout > 0)
     5293    {
     5294      MivString(curr_weight, target_weight, next_weight);
     5295    }
     5296//#endif
     5297    if(MivComp(target_weight,curr_weight) == 1)// || endwalks == TRUE)
     5298    {/*
     5299//#ifdef CHECK_IDEAL_MWALK
     5300      if(printout > 0)
     5301      {
     5302        PrintS("\n//** Mwalk: entering last cone.\n");
     5303      }
     5304//#endif
     5305
    52545306      Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
    52555307      if(target_M->length() == nV)
     
    52645316      Gomega1 = idrMoveR(Gomega, baseRing,currRing);
    52655317      idDelete(&Gomega);
    5266 #ifdef CHECK_IDEAL_MWALK
    5267       idString(Gomega1, "Gomega");
    5268 #endif
     5318//#ifdef CHECK_IDEAL_MWALK
     5319      if(printout > 1)
     5320      {
     5321        idString(Gomega1, "//** Mwalk: Gomega");
     5322      }
     5323      PrintS("\n //** Mwalk: kStd(Gomega)");
     5324//#endif
    52695325      M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL);
    5270 #ifdef CHECK_IDEAL_MWALK
    5271       idString(M,"M");
    5272 #endif
     5326//#ifdef CHECK_IDEAL_MWALK
     5327      if(printout > 1)
     5328      {
     5329        idString(M,"//** Mwalk: M");
     5330      }
     5331//#endif
    52735332      rChangeCurrRing(baseRing);
    52745333      M1 =  idrMoveR(M, newRing,currRing);
     
    52765335      Gomega2 = idrMoveR(Gomega1, newRing,currRing);
    52775336      idDelete(&Gomega1);
     5337      //PrintS("\n //** Mwalk: MLifttwoIdeal");
    52785338      F = MLifttwoIdeal(Gomega2, M1, G);
    5279 #ifdef CHECK_IDEAL_MWALK
    5280       idString(F,"F");
    5281 #endif
     5339//#ifdef CHECK_IDEAL_MWALK
     5340      if(printout > 2)
     5341      {
     5342        idString(F,"//** Mwalk: F");
     5343      }
     5344//#endif
    52825345      idDelete(&Gomega2);
    52835346      idDelete(&M1);
     
    52915354      to = clock();
    52925355#endif
    5293  //     if(si_opt_1 == (Sy_bit(OPT_REDSB)))
    5294   //    {
    5295         G = kInterRedCC(G,NULL); //reduce the Groebner basis <G> w.r.t. currRing, if option(redSB) is set
    5296   //    }
     5356      //PrintS("\n //**Mwalk: Interreduce");
     5357      //interreduce the Groebner basis <G> w.r.t. currRing
     5358      //G = kInterRedCC(G,NULL);
    52975359#ifdef TIME_TEST
    52985360      tred = tred + clock() - to;
    52995361#endif
    53005362      idSkipZeroes(G);
    5301       delete next_weight;
     5363      delete next_weight; */
    53025364      break;
    5303 #ifdef CHECK_IDEAL_MWALK
    5304       PrintS("\n//** Mwalk: last cone.\n");
    5305 #endif
    5306     }
    5307 #ifdef CHECK_IDEAL_MWALK
    5308     PrintS("\n//** Mwalk: update weight vectors.\n");
    5309 #endif
     5365    }
     5366
    53105367    for(i=nV-1; i>=0; i--)
    53115368    {
    5312       (*tmp_weight)[i] = (*curr_weight)[i];
     5369      //(*tmp_weight)[i] = (*curr_weight)[i];
    53135370      (*curr_weight)[i] = (*next_weight)[i];
    53145371    }
     
    53175374  rChangeCurrRing(XXRing);
    53185375  ideal result = idrMoveR(G,baseRing,currRing);
     5376  idDelete(&Go);
    53195377  idDelete(&G);
    5320 /*#ifdef CHECK_IDEAL_MWALK
    5321   pDelete(&p);
    5322 #endif*/
    5323   delete tmp_weight;
     5378  //delete tmp_weight;
    53245379  delete ivNull;
    53255380  delete exivlp;
     
    53285383#endif
    53295384#ifdef TIME_TEST
    5330   Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep);
    53315385  TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep);
    5332   Print("\n//** Mwalk: Ergebnis.\n");
    53335386  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
    53345387  //Print("\n// Overflow_Error? (%d)\n", Overflow_Error);
    53355388#endif
     5389  Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep);
    53365390  return(result);
    53375391}
    53385392
    5339 // 07.11.2012
    5340 // THE RANDOM WALK ALGORITHM  ideal Go, intvec* orig_M, intvec* target_M, ring baseRing
    5341 ideal Mrwalk(ideal Go, intvec* orig_M, intvec* target_M, int weight_rad, int pert_deg, ring baseRing)
     5393// THE RANDOM WALK ALGORITHM
     5394ideal Mrwalk(ideal Go, intvec* orig_M, intvec* target_M, int weight_rad, int pert_deg,
     5395             int reduction, int printout)
    53425396{
    53435397  BITSET save1 = si_opt_1; // save current options
    5344   //si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
    5345   //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
    5346   //si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB));
     5398  if(reduction == 0)
     5399  {
     5400    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     5401    si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     5402  }
    53475403  Set_Error(FALSE);
    53485404  Overflow_Error = FALSE;
     5405  //BOOLEAN endwalks = FALSE;
    53495406#ifdef TIME_TEST
    53505407  clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
     
    53545411#endif
    53555412  nstep=0;
    5356   int i,nwalk,endwalks = 0;
    5357   int nV = baseRing->N;
    5358 
    5359   ideal Gomega, M, F, Gomega1, Gomega2, M1; //, F1;
     5413  int i,polylength,nwalk;
     5414  int nV = currRing->N;
     5415
     5416  ideal Gomega, M, F,FF, Gomega1, Gomega2, M1;
    53605417  ring newRing;
    5361   ring XXRing = baseRing;
     5418  ring targetRing;
     5419  ring baseRing = currRing;
     5420  ring XXRing = currRing;
    53625421  intvec* ivNull = new intvec(nV);
    53635422  intvec* curr_weight = new intvec(nV);
    53645423  intvec* target_weight = new intvec(nV);
    53655424  intvec* exivlp = Mivlp(nV);
     5425  intvec* next_weight= new intvec(nV);
     5426/*
    53665427  intvec* tmp_weight = new intvec(nV);
    53675428  for(i=0; i<nV; i++)
     
    53695430    (*tmp_weight)[i] = (*target_M)[i];
    53705431  }
     5432*/
    53715433  for(i=0; i<nV; i++)
    53725434  {
     
    53855447#endif
    53865448  rComplete(currRing);
    5387 #ifdef CHECK_IDEAL_MWALK
    5388     idString(Go,"Go");
    5389 #endif
    53905449#ifdef TIME_TEST
    53915450  to = clock();
    53925451#endif
    5393      if(orig_M->length() == nV)
    5394       {
    5395         newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
    5396       }
    5397       else
    5398       {
    5399         newRing = VMatrDefault(orig_M);
    5400       }
     5452
     5453  if(target_M->length() == nV)
     5454  {
     5455   // define the target ring
     5456    targetRing = VMrDefault(target_weight);
     5457  }
     5458  else
     5459  {
     5460    targetRing = VMatrDefault(target_M);
     5461  }
     5462  if(orig_M->length() == nV)
     5463  {
     5464    newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
     5465  }
     5466  else
     5467  {
     5468    newRing = VMatrDefault(orig_M);
     5469  }
    54015470  rChangeCurrRing(newRing);
    54025471  ideal G = MstdCC(idrMoveR(Go,baseRing,currRing));
     
    54145483    to = clock();
    54155484#endif
    5416 #ifdef CHECK_IDEAL_MWALK
    5417     idString(G,"G");
    5418 #endif
     5485
    54195486    Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
     5487    //polylength = 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise
     5488    polylength = lengthpoly(Gomega);
    54205489#ifdef TIME_TEST
    54215490    tif = tif + clock()-to; //time for computing initial form ideal
    54225491#endif
    5423 #ifdef CHECK_IDEAL_MWALK
    5424     idString(Gomega,"Gomega");
    5425 #endif
     5492//#ifdef CHECK_IDEAL_MWALK
     5493    if(printout > 1)
     5494    {
     5495      idString(Gomega,"//** Mrwalk: Gomega");
     5496    }
     5497//#endif
     5498    // test whether target cone is reached
     5499/*    if(test_w_in_ConeCC(G,target_weight) == 1)
     5500    {
     5501      endwalks = TRUE;
     5502    }*/
     5503    if(reduction == 0)
     5504    {
     5505     
     5506      FF = middleOfCone(G,Gomega);
     5507     
     5508      if(FF != NULL)
     5509      {
     5510        idDelete(&G);
     5511        G = idCopy(FF);
     5512        idDelete(&FF);
     5513       
     5514        goto NEXT_VECTOR;
     5515      }
     5516    }
    54265517#ifndef  BUCHBERGER_ALG
    54275518    if(isNolVector(curr_weight) == 0)
    54285519    {
    54295520      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
    5430     }
     5521    }   
    54315522    else
    54325523    {
     
    54495540     if(target_M->length() == nV)
    54505541     {
    5451        newRing = VMrRefine(curr_weight,target_weight); //define a new ring with ordering "(a(curr_weight),Wp(target_weight))"
     5542       newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
     5543       //newRing = VMrRefine(curr_weight,target_weight); //define a new ring with ordering "(a(curr_weight),Wp(target_weight))"
    54525544     }
    54535545     else
     
    54735565#endif
    54745566    idSkipZeroes(M);
    5475 #ifdef CHECK_IDEAL_MWALK
    5476     PrintS("\n//** Mwalk: computed M.\n");
    5477     idString(M, "M");
    5478 #endif
     5567//#ifdef CHECK_IDEAL_MWALK
     5568    if(printout > 2)
     5569    {
     5570      idString(M, "//** Mrwalk: M");
     5571    }
     5572//#endif
    54795573    //change the ring to baseRing
    54805574    rChangeCurrRing(baseRing);
     
    54865580    to = clock();
    54875581#endif
    5488     // compute a representation of the generators of submod (M) with respect to those of mod (Gomega), where Gomega is a reduced Groebner basis w.r.t. the current ring
     5582    // compute a representation of the generators of submod (M) with respect to those of mod (Gomega),
     5583    // where Gomega is a reduced Groebner basis w.r.t. the current ring
    54895584    F = MLifttwoIdeal(Gomega2, M1, G);
    54905585#ifdef TIME_TEST
    54915586    tlift = tlift + clock() - to;
    54925587#endif
    5493 #ifdef CHECK_IDEAL_MWALK
    5494     idString(F, "F");
    5495 #endif
     5588//#ifdef CHECK_IDEAL_MWALK
     5589    if(printout > 2)
     5590    {
     5591      idString(F, "//** Mrwalk: F");
     5592    }
     5593//#endif
    54965594    idDelete(&Gomega2);
    54975595    idDelete(&M1);
     
    55025600#ifdef TIME_TEST
    55035601    to = clock();
    5504 #endif
    5505     //G = kStd(F1,NULL,testHomog,NULL,NULL,0,0,NULL);
    5506 #ifdef TIME_TEST
    55075602    tstd = tstd + clock() - to;
    55085603#endif
    55095604    idSkipZeroes(G);
    5510 #ifdef CHECK_IDEAL_MWALK
    5511     idString(G, "G");
    5512 #endif
     5605//#ifdef CHECK_IDEAL_MWALK
     5606    if(printout > 2)
     5607    {
     5608      idString(G, "//** Mrwalk: G");
     5609    }
     5610//#endif
     5611
     5612    rChangeCurrRing(targetRing);
     5613    G = idrMoveR(G,newRing,currRing);
     5614    // test whether target cone is reached
     5615    if(reduction !=0 && test_w_in_ConeCC(G,curr_weight) == 1)
     5616    {
     5617      baseRing = currRing;
     5618      break;
     5619    }
     5620
     5621    rChangeCurrRing(newRing);
     5622    G = idrMoveR(G,targetRing,currRing);
     5623    baseRing = currRing;
     5624
     5625
     5626    NEXT_VECTOR:
    55135627#ifdef TIME_TEST
    55145628    to = clock();
    55155629#endif
    5516     intvec* next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, pert_deg);//next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
     5630    next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
     5631    if(polylength > 0)
     5632    {
     5633      //there is a polynomial in Gomega with at least 3 monomials,
     5634      //low-dimensional facet of the cone
     5635      delete next_weight;
     5636      next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, pert_deg);
     5637    }
    55175638#ifdef TIME_TEST
    55185639    tnw = tnw + clock() - to;
    55195640#endif
    5520 #ifdef PRINT_VECTORS
    5521     MivString(curr_weight, target_weight, next_weight);
    5522 #endif
    5523     if(MivComp(next_weight, ivNull) == 1 || MivComp(target_weight,curr_weight) == 1)// || test_w_in_ConeCC(G, target_weight) == 1 || MivComp(next_weight,curr_weight) == 1)
    5524     {
    5525 #ifdef CHECK_IDEAL_MWALK
    5526       PrintS("\n//** Mwalk: entering last cone.\n");
    5527 #endif
     5641//#ifdef PRINT_VECTORS
     5642    if(printout > 0)
     5643    {
     5644      MivString(curr_weight, target_weight, next_weight);
     5645    }
     5646//#endif
     5647    if(MivComp(next_weight, ivNull) == 1 || MivComp(target_weight,curr_weight) == 1)// || endwalks == TRUE)
     5648    {/*
     5649//#ifdef CHECK_IDEAL_MWALK
     5650      if(printout > 0)
     5651      {
     5652        PrintS("\n//** Mrwalk: entering last cone.\n");
     5653      }
     5654//#endif
     5655
    55285656      Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
    55295657      if(target_M->length() == nV)
     
    55385666      Gomega1 = idrMoveR(Gomega, baseRing,currRing);
    55395667      idDelete(&Gomega);
    5540 #ifdef CHECK_IDEAL_MWALK
    5541       idString(Gomega1, "Gomega");
    5542 #endif
     5668//#ifdef CHECK_IDEAL_MWALK
     5669      if(printout > 1)
     5670      {
     5671        idString(Gomega1, "//** Mrwalk: Gomega");
     5672      }
     5673      PrintS("\n //** Mrwalk: kStd(Gomega)");
     5674//#endif
    55435675      M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL);
    5544 #ifdef CHECK_IDEAL_MWALK
    5545       idString(M,"M");
    5546 #endif
     5676//#ifdef CHECK_IDEAL_MWALK
     5677      if(printout > 1)
     5678      {
     5679        idString(M,"//** Mrwalk: M");
     5680      }
     5681//#endif
    55475682      rChangeCurrRing(baseRing);
    55485683      M1 =  idrMoveR(M, newRing,currRing);
     
    55505685      Gomega2 = idrMoveR(Gomega1, newRing,currRing);
    55515686      idDelete(&Gomega1);
     5687      //PrintS("\n //** Mrwalk: MLifttwoIdeal");
    55525688      F = MLifttwoIdeal(Gomega2, M1, G);
    5553 #ifdef CHECK_IDEAL_MWALK
    5554       idString(F,"F");
    5555 #endif
     5689//#ifdef CHECK_IDEAL_MWALK
     5690      if(printout > 2)
     5691      {
     5692        idString(F,"//** Mrwalk: F");
     5693      }
     5694//#endif
    55565695      idDelete(&Gomega2);
    55575696      idDelete(&M1);
     
    55655704      to = clock();
    55665705#endif
    5567  //     if(si_opt_1 == (Sy_bit(OPT_REDSB)))
    5568   //    {
    5569         //G = kInterRedCC(G,NULL); //reduce the Groebner basis <G> w.r.t. currRing, if option(redSB) is set
    5570   //    }
     5706      //PrintS("\n //**Mrwalk: Interreduce");
     5707      //interreduce the Groebner basis <G> w.r.t. currRing
     5708      //G = kInterRedCC(G,NULL);
    55715709#ifdef TIME_TEST
    55725710      tred = tred + clock() - to;
    55735711#endif
    55745712      idSkipZeroes(G);
    5575       delete next_weight;
     5713      delete next_weight;*/
    55765714      break;
    5577 #ifdef CHECK_IDEAL_MWALK
    5578       PrintS("\n//** Mwalk: last cone.\n");
    5579 #endif
    5580     }
    5581 #ifdef CHECK_IDEAL_MWALK
    5582     PrintS("\n//** Mwalk: update weight vectors.\n");
    5583 #endif
     5715    }
     5716
    55845717    for(i=nV-1; i>=0; i--)
    55855718    {
    5586       (*tmp_weight)[i] = (*curr_weight)[i];
     5719      //(*tmp_weight)[i] = (*curr_weight)[i];
    55875720      (*curr_weight)[i] = (*next_weight)[i];
    55885721    }
    55895722    delete next_weight;
    55905723  }
     5724  baseRing = currRing;
    55915725  rChangeCurrRing(XXRing);
    55925726  ideal result = idrMoveR(G,baseRing,currRing);
    55935727  idDelete(&G);
    5594 /*#ifdef CHECK_IDEAL_MWALK
    5595   pDelete(&p);
    5596 #endif*/
    5597   delete tmp_weight;
     5728  si_opt_1 = save1; //set original options, e. g. option(RedSB)
     5729  //delete tmp_weight;
    55985730  delete ivNull;
    55995731  delete exivlp;
     
    56015733  delete last_omega;
    56025734#endif
     5735  Print("\n//** Mrwalk: Groebner Walk took %d steps.\n", nstep);
    56035736#ifdef TIME_TEST
    5604   Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep);
    56055737  TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep);
    5606   Print("\n//** Mwalk: Ergebnis.\n");
    56075738  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
    56085739  //Print("\n// Overflow_Error? (%d)\n", Overflow_Error);
     
    57515882// use kStd, if nP = 0, else call LastGB
    57525883ideal Mpwalk(ideal Go, int op_deg, int tp_deg,intvec* curr_weight,
    5753              intvec* target_weight, int nP)
     5884             intvec* target_weight, int nP, int reduction, int printout)
    57545885{
     5886  BITSET save1 = si_opt_1; // save current options
     5887  if(reduction == 0)
     5888  {
     5889    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     5890    //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     5891    //si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB));
     5892  }
    57555893  Set_Error(FALSE  );
    57565894  Overflow_Error = FALSE;
     
    57665904  nstep = 0;
    57675905  int i, ntwC=1, ntestw=1,  nV = currRing->N;
    5768   int endwalks=0;
    5769 
    5770   ideal Gomega, M, F, G, Gomega1, Gomega2, M1,F1,Eresult,ssG;
     5906  BOOLEAN endwalks = FALSE;
     5907
     5908  ideal Gomega, M, F, FF, G, Gomega1, Gomega2, M1,F1,Eresult,ssG;
    57715909  ring newRing, oldRing, TargetRing;
    57725910  intvec* iv_M_dp;
     
    57905928  ring XXRing = currRing;
    57915929
    5792 
    57935930  to = clock();
    5794   /* perturbs the original vector */
     5931  // perturbs the original vector
    57955932  if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) := "dp"
    57965933  {
     
    58095946      DefRingPar(curr_weight);
    58105947    else
    5811       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 1
     5948      rChangeCurrRing(VMrDefault(curr_weight));
    58125949
    58135950    G = idrMoveR(Go, XXRing,currRing);
     
    58245961  ring HelpRing = currRing;
    58255962
    5826   /* perturbs the target weight vector */
     5963  // perturbs the target weight vector
    58275964  if(tp_deg > 1 && tp_deg <= nV)
    58285965  {
     
    58305967      DefRingPar(target_weight);
    58315968    else
    5832       rChangeCurrRing(VMrDefault(target_weight)); // Aenderung 2
     5969      rChangeCurrRing(VMrDefault(target_weight));
    58335970
    58345971    TargetRing = currRing;
     
    58375974    {
    58385975      iv_M_lp = MivMatrixOrderlp(nV);
    5839       //ivString(iv_M_lp, "iv_M_lp");
    5840       //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg);
    58415976      target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
    58425977    }
     
    58445979    {
    58455980      iv_M_lp = MivMatrixOrder(target_weight);
    5846       //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg);
    58475981      target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
    58485982    }
     
    58525986    G = idrMoveR(ssG, TargetRing,currRing);
    58535987  }
    5854   /*
    5855     Print("\n// Perturbationwalkalg. vom Gradpaar (%d,%d):",op_deg,tp_deg);
    5856     ivString(curr_weight, "new sigma");
    5857     ivString(target_weight, "new tau");
    5858   */
     5988    if(printout > 0)
     5989    {
     5990      Print("\n//** Mpwalk: Perturbation Walk of degree (%d,%d):",op_deg,tp_deg);
     5991      ivString(curr_weight, "//** Mpwalk: new current weight");
     5992      ivString(target_weight, "//** Mpwalk: new target weight");
     5993    }
    58595994  while(1)
    58605995  {
    58615996    nstep ++;
    58625997    to = clock();
    5863     /* compute an initial form ideal of <G> w.r.t. the weight vector
    5864        "curr_weight" */
     5998    // compute an initial form ideal of <G> w.r.t. the weight vector
     5999    // "curr_weight"
    58656000    Gomega = MwalkInitialForm(G, curr_weight);
    5866 
     6001//#ifdef CHECK_IDEAL_MWALK
     6002    if(printout > 1)
     6003    {
     6004      idString(Gomega,"//** Mpwalk: Gomega");
     6005    }
     6006//#endif
     6007    if(reduction == 0 && nstep > 1)
     6008    {
     6009      // check whether weight vector is in the interior of the cone
     6010      while(1)
     6011      {
     6012        FF = middleOfCone(G,Gomega);
     6013        if(FF != NULL)
     6014        {
     6015          idDelete(&G);
     6016          G = idCopy(FF);
     6017          idDelete(&FF);
     6018          next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
     6019          if(printout > 0)
     6020          {
     6021            MivString(curr_weight, target_weight, next_weight);
     6022          }
     6023        }
     6024        else
     6025        {
     6026          break;
     6027        }
     6028        for(i=nV-1; i>=0; i--)
     6029        {
     6030          (*curr_weight)[i] = (*next_weight)[i];
     6031        }
     6032        Gomega = MwalkInitialForm(G, curr_weight);
     6033        if(printout > 1)
     6034        {
     6035          idString(Gomega,"//** Mpwalk: Gomega");
     6036        }
     6037      }
     6038    }
    58676039
    58686040#ifdef ENDWALKS
    5869     if(endwalks == 1){
     6041    if(endwalks == TRUE)
     6042    {
    58706043      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
    58716044      idElements(G, "G");
    5872       // idElements(Gomega, "Gw");
    58736045      headidString(G, "G");
    5874       //headidString(Gomega, "Gw");
    58756046    }
    58766047#endif
     
    58916062      DefRingPar(curr_weight);
    58926063    else
    5893       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 3
     6064      rChangeCurrRing(VMrDefault(curr_weight));
    58946065
    58956066    newRing = currRing;
     
    58976068
    58986069#ifdef ENDWALKS
    5899     if(endwalks==1)
     6070    if(endwalks==TRUE)
    59006071    {
    59016072      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     
    59126083    tim = clock();
    59136084    to = clock();
    5914     /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
     6085    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
    59156086#ifdef  BUCHBERGER_ALG
    59166087    M = MstdhomCC(Gomega1);
     
    59186089    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
    59196090    delete hilb_func;
    5920 #endif // BUCHBERGER_ALG
    5921 
    5922     if(endwalks == 1){
     6091#endif
     6092//#ifdef CHECK_IDEAL_MWALK
     6093      if(printout > 2)
     6094      {
     6095        idString(M,"//** Mpwalk: M");
     6096      }
     6097//#endif
     6098
     6099    if(endwalks == TRUE){
    59236100      xtstd = xtstd+clock()-to;
    59246101#ifdef ENDWALKS
     
    59306107      tstd=tstd+clock()-to;
    59316108
    5932     /* change the ring to oldRing */
     6109    // change the ring to oldRing
    59336110    rChangeCurrRing(oldRing);
    59346111    M1 =  idrMoveR(M, newRing,currRing);
    59356112    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
    5936 
    5937     //if(endwalks==1)  PrintS("\n// Lifting is working:..");
    59386113
    59396114    to=clock();
     
    59426117       Gomega is a reduced Groebner basis w.r.t. the current ring */
    59436118    F = MLifttwoIdeal(Gomega2, M1, G);
    5944     if(endwalks != 1)
     6119    if(endwalks == FALSE)
    59456120      tlift = tlift+clock()-to;
    59466121    else
    59476122      xtlift=clock()-to;
    59486123
     6124//#ifdef CHECK_IDEAL_MWALK
     6125    if(printout > 2)
     6126    {
     6127      idString(F,"//** Mpwalk: F");
     6128    }
     6129//#endif
     6130
    59496131    idDelete(&M1);
    59506132    idDelete(&Gomega2);
    59516133    idDelete(&G);
    59526134
    5953     /* change the ring to newRing */
     6135    // change the ring to newRing
    59546136    rChangeCurrRing(newRing);
    5955     F1 = idrMoveR(F, oldRing,currRing);
    5956 
    5957     //if(endwalks==1)PrintS("\n// InterRed is working now:");
     6137    if(reduction == 0)
     6138    {
     6139      G = idrMoveR(F,oldRing,currRing);
     6140    }
     6141    else
     6142    {
     6143      F1 = idrMoveR(F, oldRing,currRing);
     6144      if(printout > 2)
     6145      {
     6146        PrintS("\n //** Mpwalk: reduce the Groebner basis.\n");
     6147      }
     6148      to=clock();
     6149      G = kInterRedCC(F1, NULL);
     6150      if(endwalks == FALSE)
     6151        tred = tred+clock()-to;
     6152      else
     6153        xtred=clock()-to;
     6154      idDelete(&F1);
     6155    }
     6156    if(endwalks == TRUE)
     6157      break;
    59586158
    59596159    to=clock();
    5960     /* reduce the Groebner basis <G> w.r.t. new ring */
    5961     G = kInterRedCC(F1, NULL);
    5962     if(endwalks != 1)
    5963       tred = tred+clock()-to;
    5964     else
    5965       xtred=clock()-to;
    5966 
    5967     idDelete(&F1);
    5968 
    5969     if(endwalks == 1)
    5970       break;
    5971 
    5972     to=clock();
    5973     /* compute a next weight vector */
     6160    // compute a next weight vector
    59746161    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
    59756162    tnw=tnw+clock()-to;
    5976 #ifdef PRINT_VECTORS
    5977     MivString(curr_weight, target_weight, next_weight);
    5978 #endif
     6163//#ifdef PRINT_VECTORS
     6164    if(printout > 0)
     6165    {
     6166      MivString(curr_weight, target_weight, next_weight);
     6167    }
     6168//#endif
    59796169
    59806170    if(Overflow_Error == TRUE)
     
    59946184    }
    59956185    if(MivComp(next_weight, target_weight) == 1)
    5996       endwalks = 1;
     6186      endwalks = TRUE;
    59976187
    59986188    for(i=nV-1; i>=0; i--)
     
    60006190
    60016191    delete next_weight;
    6002   }//while
     6192  }//end of while-loop
    60036193
    60046194  if(tp_deg != 1)
     
    60146204        DefRingPar(orig_target);
    60156205      else
    6016         rChangeCurrRing(VMrDefault(orig_target)); //Aenderung
     6206        rChangeCurrRing(VMrDefault(orig_target));
    60176207
    60186208    TargetRing=currRing;
     
    60306220    if( ntestw != 1 || ntwC == 0)
    60316221    {
    6032       /*
    6033       if(ntestw != 1){
     6222      if(ntestw != 1 && printout >2)
     6223      {
    60346224        ivString(pert_target_vector, "tau");
    60356225        PrintS("\n// ** perturbed target vector doesn't stay in cone!!");
    60366226        Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
    6037         idElements(F1, "G");
    6038       }
    6039       */
     6227        //idElements(F1, "G");
     6228      }
    60406229      // LastGB is "better" than the kStd subroutine
    60416230      to=clock();
     
    60686257    Eresult = idrMoveR(G, newRing,currRing);
    60696258  }
     6259  si_opt_1 = save1; //set original options, e. g. option(RedSB)
    60706260  delete ivNull;
    60716261  if(tp_deg != 1)
     
    60826272             tnw+xtnw);
    60836273
    6084   Print("\n// pSetm_Error = (%d)", ErrorCheck());
    6085   Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep,  Overflow_Error);
    6086 #endif
     6274  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
     6275  //Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep,  Overflow_Error);
     6276#endif
     6277  Print("\n//** Mpwalk: Perturbation Walk took %d steps.\n", nstep);
     6278  return(Eresult);
     6279}
     6280
     6281/*******************************************************
     6282 * THE PERTURBATION WALK ALGORITHM WITH RANDOM ELEMENT *
     6283 *******************************************************/
     6284ideal Mprwalk(ideal Go, intvec* orig_M, intvec* target_M, int weight_rad,
     6285              int op_deg, int tp_deg, int nP, int reduction, int printout)
     6286{
     6287  BITSET save1 = si_opt_1; // save current options
     6288  if(reduction == 0)
     6289  {
     6290    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     6291    //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     6292    //si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB));
     6293  }
     6294  Set_Error(FALSE);
     6295  Overflow_Error = FALSE;
     6296  //Print("// pSetm_Error = (%d)", ErrorCheck());
     6297
     6298  clock_t  tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
     6299  xtextra=0;
     6300  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0;
     6301  tinput = clock();
     6302
     6303  clock_t tim;
     6304
     6305  nstep = 0;
     6306  int i, ntwC=1, ntestw=1, polylength, nV = currRing->N;
     6307  BOOLEAN endwalks = FALSE;
     6308
     6309  ideal Gomega, M, F, FF, G, Gomega1, Gomega2, M1,F1,Eresult,ssG;
     6310  ring newRing, oldRing, TargetRing;
     6311  intvec* iv_M_dp;
     6312  intvec* iv_M_lp;
     6313  intvec* exivlp = Mivlp(nV);
     6314  intvec* curr_weight = new intvec(nV);
     6315  intvec* target_weight = new intvec(nV);
     6316  for(i=0; i<nV; i++)
     6317  {
     6318    (*curr_weight)[i] = (*orig_M)[i];
     6319    (*target_weight)[i] = (*target_M)[i];
     6320  }
     6321  intvec* orig_target = target_weight;
     6322  intvec* pert_target_vector = target_weight;
     6323  intvec* ivNull = new intvec(nV);
     6324  intvec* iv_dp = MivUnit(nV);// define (1,1,...,1)
     6325#ifndef BUCHBERGER_ALG
     6326  intvec* hilb_func;
     6327#endif
     6328  intvec* next_weight;
     6329
     6330  // to avoid (1,0,...,0) as the target vector
     6331  intvec* last_omega = new intvec(nV);
     6332  for(i=nV-1; i>0; i--)
     6333    (*last_omega)[i] = 1;
     6334  (*last_omega)[0] = 10000;
     6335
     6336  ring XXRing = currRing;
     6337
     6338  to = clock();
     6339  // perturbs the original vector
     6340  if(orig_M->length() == nV)
     6341  {
     6342    if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) := "dp"
     6343    {
     6344      G = MstdCC(Go);
     6345      tostd = clock()-to;
     6346      if(op_deg != 1)
     6347      {
     6348        iv_M_dp = MivMatrixOrderdp(nV);
     6349        curr_weight = MPertVectors(G, iv_M_dp, op_deg);
     6350      }
     6351    }
     6352    else
     6353    {
     6354      //define ring order := (a(curr_weight),lp);
     6355      if (rParameter(currRing) != NULL)
     6356        DefRingPar(curr_weight);
     6357      else
     6358        rChangeCurrRing(VMrDefault(curr_weight));
     6359
     6360      G = idrMoveR(Go, XXRing,currRing);
     6361      G = MstdCC(G);
     6362      tostd = clock()-to;
     6363      if(op_deg != 1)
     6364      {
     6365        iv_M_dp = MivMatrixOrder(curr_weight);
     6366        curr_weight = MPertVectors(G, iv_M_dp, op_deg);
     6367      }
     6368    }
     6369  }
     6370  else
     6371  {
     6372    rChangeCurrRing(VMatrDefault(orig_M));
     6373    G = idrMoveR(Go, XXRing,currRing);
     6374    G = MstdCC(G);
     6375    tostd = clock()-to;
     6376    if(op_deg != 1)
     6377    {
     6378      curr_weight = MPertVectors(G, orig_M, op_deg);
     6379    }
     6380  }
     6381
     6382  delete iv_dp;
     6383  if(op_deg != 1) delete iv_M_dp;
     6384
     6385  ring HelpRing = currRing;
     6386
     6387  // perturbs the target weight vector
     6388  if(target_M->length() == nV)
     6389  {
     6390    if(tp_deg > 1 && tp_deg <= nV)
     6391    {
     6392      if (rParameter(currRing) != NULL)
     6393        DefRingPar(target_weight);
     6394      else
     6395        rChangeCurrRing(VMrDefault(target_weight));
     6396
     6397      TargetRing = currRing;
     6398      ssG = idrMoveR(G,HelpRing,currRing);
     6399      if(MivSame(target_weight, exivlp) == 1)
     6400      {
     6401        iv_M_lp = MivMatrixOrderlp(nV);
     6402        target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
     6403      }
     6404      else
     6405      {
     6406        iv_M_lp = MivMatrixOrder(target_weight);
     6407        target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
     6408      }
     6409      delete iv_M_lp;
     6410      pert_target_vector = target_weight;
     6411      rChangeCurrRing(HelpRing);
     6412      G = idrMoveR(ssG, TargetRing,currRing);
     6413    }
     6414  }
     6415  else
     6416  {
     6417    if(tp_deg > 1 && tp_deg <= nV)
     6418    {
     6419      rChangeCurrRing(VMatrDefault(target_M));
     6420      TargetRing = currRing;
     6421      ssG = idrMoveR(G,HelpRing,currRing);
     6422      target_weight = MPertVectors(ssG, target_M, tp_deg);
     6423    }
     6424  }
     6425  if(printout > 0)
     6426  {
     6427    Print("\n//** Mprwalk: Random Perturbation Walk of degree (%d,%d):",op_deg,tp_deg);
     6428    ivString(curr_weight, "//** Mprwalk: new current weight");
     6429    ivString(target_weight, "//** Mprwalk: new target weight");
     6430  }
     6431  while(1)
     6432  {
     6433    nstep ++;
     6434    to = clock();
     6435    // compute an initial form ideal of <G> w.r.t. the weight vector
     6436    // "curr_weight"
     6437    Gomega = MwalkInitialForm(G, curr_weight);
     6438    polylength = lengthpoly(Gomega);
     6439//#ifdef CHECK_IDEAL_MWALK
     6440    if(printout > 1)
     6441    {
     6442      idString(Gomega,"//** Mprwalk: Gomega");
     6443    }
     6444//#endif
     6445
     6446    if(reduction == 0 && nstep > 1)
     6447    {
     6448      // check whether weight vector is in the interior of the cone
     6449      while(1)
     6450      {
     6451        FF = middleOfCone(G,Gomega);
     6452        if(FF != NULL)
     6453        {
     6454          idDelete(&G);
     6455          G = idCopy(FF);
     6456          idDelete(&FF);
     6457          next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
     6458          if(printout > 0)
     6459          {
     6460            MivString(curr_weight, target_weight, next_weight);
     6461          }
     6462        }
     6463        else
     6464        {
     6465          break;
     6466        }
     6467        for(i=nV-1; i>=0; i--)
     6468        {
     6469          (*curr_weight)[i] = (*next_weight)[i];
     6470        }
     6471        Gomega = MwalkInitialForm(G, curr_weight);
     6472        if(printout > 1)
     6473        {
     6474          idString(Gomega,"//** Mprwalk: Gomega");
     6475        }
     6476      }
     6477    }
     6478
     6479#ifdef ENDWALKS
     6480    if(endwalks == TRUE)
     6481    {
     6482      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6483      idElements(G, "G");
     6484      headidString(G, "G");
     6485    }
     6486#endif
     6487
     6488    tif = tif + clock()-to;
     6489
     6490#ifndef  BUCHBERGER_ALG
     6491    if(isNolVector(curr_weight) == 0)
     6492      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
     6493    else
     6494      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
     6495#endif // BUCHBERGER_ALG
     6496
     6497    oldRing = currRing;
     6498
     6499    if(target_M->length() == nV)
     6500    {
     6501      // define a new ring with ordering "(a(curr_weight),lp)
     6502      if (rParameter(currRing) != NULL)
     6503        DefRingPar(curr_weight);
     6504      else
     6505        rChangeCurrRing(VMrDefault(curr_weight));
     6506    }
     6507    else
     6508    {
     6509      rChangeCurrRing(VMatrRefine(target_M,curr_weight));
     6510    }
     6511    newRing = currRing;
     6512    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
     6513#ifdef ENDWALKS
     6514    if(endwalks == TRUE)
     6515    {
     6516      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6517      idElements(Gomega1, "Gw");
     6518      headidString(Gomega1, "headGw");
     6519      PrintS("\n// compute a rGB of Gw:\n");
     6520
     6521#ifndef  BUCHBERGER_ALG
     6522      ivString(hilb_func, "w");
     6523#endif
     6524    }
     6525#endif
     6526
     6527    tim = clock();
     6528    to = clock();
     6529    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
     6530#ifdef  BUCHBERGER_ALG
     6531    M = MstdhomCC(Gomega1);
     6532#else
     6533    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
     6534    delete hilb_func;
     6535#endif
     6536//#ifdef CHECK_IDEAL_MWALK
     6537    if(printout > 2)
     6538    {
     6539      idString(M,"//** Mprwalk: M");
     6540    }
     6541//#endif
     6542
     6543    if(endwalks == TRUE)
     6544    {
     6545      xtstd = xtstd+clock()-to;
     6546#ifdef ENDWALKS
     6547      Print("\n// time for the last std(Gw)  = %.2f sec\n",
     6548            ((double) clock())/1000000 -((double)tim) /1000000);
     6549#endif
     6550    }
     6551    else
     6552      tstd=tstd+clock()-to;
     6553
     6554    /* change the ring to oldRing */
     6555    rChangeCurrRing(oldRing);
     6556    M1 =  idrMoveR(M, newRing,currRing);
     6557    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
     6558
     6559    to=clock();
     6560    /* compute a representation of the generators of submod (M)
     6561       with respect to those of mod (Gomega).
     6562       Gomega is a reduced Groebner basis w.r.t. the current ring */
     6563    F = MLifttwoIdeal(Gomega2, M1, G);
     6564    if(endwalks == FALSE)
     6565      tlift = tlift+clock()-to;
     6566    else
     6567      xtlift=clock()-to;
     6568
     6569//#ifdef CHECK_IDEAL_MWALK
     6570    if(printout > 2)
     6571    {
     6572      idString(F,"//** Mprwalk: F");
     6573    }
     6574//#endif
     6575
     6576    idDelete(&M1);
     6577    idDelete(&Gomega2);
     6578    idDelete(&G);
     6579
     6580    // change the ring to newRing
     6581    rChangeCurrRing(newRing);
     6582    if(reduction == 0)
     6583    {
     6584      G = idrMoveR(F,oldRing,currRing);
     6585    }
     6586    else
     6587    {
     6588      F1 = idrMoveR(F, oldRing,currRing);
     6589      if(printout > 2)
     6590      {
     6591        PrintS("\n //** Mprwalk: reduce the Groebner basis.\n");
     6592      }
     6593      to=clock();
     6594      G = kInterRedCC(F1, NULL);
     6595      if(endwalks == FALSE)
     6596        tred = tred+clock()-to;
     6597      else
     6598        xtred=clock()-to;
     6599      idDelete(&F1);
     6600    }
     6601
     6602    if(endwalks == TRUE)
     6603      break;
     6604
     6605    to=clock();
     6606
     6607    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
     6608    if(polylength > 0)
     6609    {
     6610      if(printout > 2)
     6611      {
     6612        Print("\n //**Mprwalk: there is a polynomial in Gomega with at least 3 monomials, low-dimensional facet of the cone.\n");
     6613      }
     6614      delete next_weight;
     6615      next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, op_deg);
     6616    }
     6617    tnw=tnw+clock()-to;
     6618//#ifdef PRINT_VECTORS
     6619    if(printout > 0)
     6620    {
     6621      MivString(curr_weight, target_weight, next_weight);
     6622    }
     6623//#endif
     6624
     6625    if(Overflow_Error == TRUE)
     6626    {
     6627      ntwC = 0;
     6628      //ntestomega = 1;
     6629      //Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6630      //idElements(G, "G");
     6631      delete next_weight;
     6632      goto FINISH_160302;
     6633    }
     6634    if(MivComp(next_weight, ivNull) == 1){
     6635      newRing = currRing;
     6636      delete next_weight;
     6637      //Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6638      break;
     6639    }
     6640    if(MivComp(next_weight, target_weight) == 1)
     6641      endwalks = TRUE;
     6642
     6643    for(i=nV-1; i>=0; i--)
     6644      (*curr_weight)[i] = (*next_weight)[i];
     6645
     6646    delete next_weight;
     6647  }//while
     6648
     6649  if(tp_deg != 1)
     6650  {
     6651    FINISH_160302:
     6652    if(target_M->length() == nV)
     6653    {
     6654      if(MivSame(orig_target, exivlp) == 1)
     6655        if (rParameter(currRing) != NULL)
     6656          DefRingParlp();
     6657        else
     6658          VMrDefaultlp();
     6659      else
     6660        if (rParameter(currRing) != NULL)
     6661          DefRingPar(orig_target);
     6662        else
     6663          rChangeCurrRing(VMrDefault(orig_target));
     6664    }
     6665    else
     6666    {
     6667      rChangeCurrRing(VMatrDefault(target_M));
     6668    }
     6669    TargetRing=currRing;
     6670    F1 = idrMoveR(G, newRing,currRing);
     6671#ifdef CHECK_IDEAL
     6672      headidString(G, "G");
     6673#endif
     6674
     6675    // check whether the pertubed target vector stays in the correct cone
     6676    if(ntwC != 0)
     6677    {
     6678      ntestw = test_w_in_ConeCC(F1, pert_target_vector);
     6679    }
     6680    if(ntestw != 1 || ntwC == 0)
     6681    {
     6682      if(ntestw != 1 && printout > 2)
     6683      {
     6684        ivString(pert_target_vector, "tau");
     6685        PrintS("\n// **Mprwalk: perturbed target vector doesn't stay in cone.");
     6686        Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6687        //idElements(F1, "G");
     6688      }
     6689      // LastGB is "better" than the kStd subroutine
     6690      to=clock();
     6691      ideal eF1;
     6692      if(nP == 0 || tp_deg == 1 || MivSame(orig_target, exivlp) != 1 || target_M->length() != nV)
     6693      {
     6694        if(printout > 2)
     6695        {
     6696          PrintS("\n// ** Mprwalk: Call \"std\" to compute a Groebner basis.\n");
     6697        }
     6698        eF1 = MstdCC(F1);
     6699        idDelete(&F1);
     6700      }
     6701      else
     6702      {
     6703        if(printout > 2)
     6704        {
     6705          PrintS("\n// **Mprwalk: Call \"LastGB\" to compute a Groebner basis.\n");
     6706        }
     6707        rChangeCurrRing(newRing);
     6708        ideal F2 = idrMoveR(F1, TargetRing,currRing);
     6709        eF1 = LastGB(F2, curr_weight, tp_deg-1);
     6710        F2=NULL;
     6711      }
     6712      xtextra=clock()-to;
     6713      ring exTargetRing = currRing;
     6714
     6715      rChangeCurrRing(XXRing);
     6716      Eresult = idrMoveR(eF1, exTargetRing,currRing);
     6717    }
     6718    else{
     6719      rChangeCurrRing(XXRing);
     6720      Eresult = idrMoveR(F1, TargetRing,currRing);
     6721    }
     6722  }
     6723  else {
     6724    rChangeCurrRing(XXRing);
     6725    Eresult = idrMoveR(G, newRing,currRing);
     6726  }
     6727  si_opt_1 = save1; //set original options, e. g. option(RedSB)
     6728  delete ivNull;
     6729  if(tp_deg != 1)
     6730    delete target_weight;
     6731
     6732  if(op_deg != 1 )
     6733    delete curr_weight;
     6734
     6735  delete exivlp;
     6736  delete last_omega;
     6737
     6738#ifdef TIME_TEST
     6739  TimeStringFractal(tinput, tostd, tif+xtif, tstd+xtstd,0, tlift+xtlift, tred+xtred,
     6740             tnw+xtnw);
     6741
     6742  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
     6743  //Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep,  Overflow_Error);
     6744#endif
     6745  Print("\n//** Mprwalk: Perturbation Walk took %d steps.\n", nstep);
    60876746  return(Eresult);
    60886747}
     
    61106769 * Perturb the start weight vector at the top level, i.e. nlev = 1     *
    61116770 ***********************************************************************/
    6112 static ideal rec_fractal_call(ideal G, int nlev, intvec* omtmp)
     6771static ideal rec_fractal_call(ideal G, int nlev, intvec* ivtarget,
     6772             int reduction, int printout)
    61136773{
    61146774  Overflow_Error =  FALSE;
    6115   //Print("\n\n// Entering the %d-th recursion:", nlev);
    6116 
     6775  if(printout >0)
     6776  {
     6777    Print("\n\n// Entering the %d-th recursion:", nlev);
     6778  }
    61176779  int i, nV = currRing->N;
    61186780  ring new_ring, testring;
    61196781  //ring extoRing;
    6120   ideal Gomega, Gomega1, Gomega2, F, F1, Gresult, Gresult1, G1, Gt;
     6782  ideal Gomega, Gomega1, Gomega2, FF, F, F1, Gresult, Gresult1, G1, Gt;
    61216783  int nwalks = 0;
    61226784  intvec* Mwlp;
     
    61276789  intvec* next_vect;
    61286790  intvec* omega2 = new intvec(nV);
    6129   intvec* altomega = new intvec(nV);
    6130 
     6791  intvec* omtmp = new intvec(nV);
     6792  //intvec* altomega = new intvec(nV);
     6793
     6794  for(i = nV -1; i>0; i--)
     6795  {
     6796    (*omtmp)[i] = (*ivtarget)[i];
     6797  }
    61316798  //BOOLEAN isnewtarget = FALSE;
    61326799
     
    61696836    NEXT_VECTOR_FRACTAL:
    61706837    to=clock();
    6171     /* determine the next border */
     6838    // determine the next border
    61726839    next_vect = MkInterRedNextWeight(omega,omega2,G);
    61736840    xtnw=xtnw+clock()-to;
    6174 #ifdef PRINT_VECTORS
    6175     MivString(omega, omega2, next_vect);
    6176 #endif
     6841
    61776842    oRing = currRing;
    61786843
    6179     /* We only perturb the current target vector at the recursion  level 1 */
     6844    // We only perturb the current target vector at the recursion level 1
    61806845    if(Xngleich == 0 && nlev == 1) //(ngleich == 0) important, e.g. ex2, ex3
    6181       if (MivComp(next_vect, omega2) != 1)
    6182       {
    6183         /* to dispense with taking initial (and lifting/interreducing
    6184            after the call of recursion */
    6185         //Print("\n\n// ** Perturb the both vectors with degree %d with",nlev);
    6186         //idElements(G, "G");
     6846      if (MivComp(next_vect, omega2) == 1)
     6847      {
     6848        // to dispense with taking initial (and lifting/interreducing
     6849        // after the call of recursion
     6850        if(printout > 0)
     6851        {
     6852          Print("\n//** rec_fractal_call: Perturb the both vectors with degree %d.",nlev);
     6853          //idElements(G, "G");
     6854        }
    61876855
    61886856        Xngleich = 1;