Changeset 4052dc in git

Timestamp:

Feb 6, 2015, 10:25:18 AM (9 years ago)

Author:

Stephan Oberfranz <oberfran@…>

Branches:

(u'spielwiese', '17f1d200f27c5bd38f5dfc6e8a0879242279d1d8')

Children:

99f4c6811eadeb02db7dd8cfde8940b9519df3c1e05649bb049ee8f1b6dfccd3446c63020e099649

Parents:

dbf60932968d7dd0d5f400c4209b39e40b442f58

Message:

---

Location:

Singular

Files:

• LIB/grwalk.lib (modified) (3 diffs)
• walk.cc (modified) (65 diffs)
• walk.h (modified) (1 diff)

Singular/LIB/grwalk.lib

@@ -274 +274 @@
    string ord_str =   OSCTW[2];
    intvec curr_weight   =   OSCTW[3]; /* original weight vector */
-   intvec target_weight =   OSCTW[4]; /* terget weight vector */
+   intvec target_weight =   OSCTW[4]; /* target weight vector */
    kill OSCTW;
    option(redSB);
@@ -299 +299 @@
   //** compute a Groebner basis of I w.r.t. lp.
   ring r = 32003,(z,y,x), lp;
-  ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
+  ideal I = zy2+yx2+yx+3,
+            z3x+y3+zyx-yx2-yx-3,
+            z2yx3-y5+z2yx2+y3x2+y2x3+y3x+y2x2+3z2x+3y2+3yx,
+            zyx5+y6-y4x2-y3x3+2zyx4-y4x-y3x2+zyx3-3z2yx+3zx3-3y3-3y2x+3zx2,
+            yx7-y7+y5x2+y4x3+3yx6+y5x+y4x2+3yx5-6zyx3+yx4+3x5+3y4+3y3x-6zyx2+6x4+3x3-9zx;
   gwalk(I,0,1);
 }
@@ -479 +483 @@
   ideal G = fetch(xR, Go);
 
-<<<<<<< HEAD
   G = system("Mpwalk",G,n1,n2,curr_weight,target_weight,nP,reduction,printout);
  
-=======
-  G = system("Mpwalk", G, n1, n2, curr_weight, target_weight,nP);
-
->>>>>>> f533f6f7667328bccb271b19b2f603aaebe41596
   setring xR;
-  //kill Go;
+  //kill Go; //unused
 
   keepring basering;

Singular/walk.cc

@@ -27 +27 @@
 
 #define FIRST_STEP_FRACTAL // to define the first step of the fractal
-//#define MSTDCC_FRACTAL // apply Buchberger alg to compute a red GB, if
+#define MSTDCC_FRACTAL // apply Buchberger alg to compute a red GB, if
 //                          tau doesn't stay in the correct cone
 
@@ -42 +42 @@
 #include <Singular/ipshell.h>
 #include <Singular/ipconv.h>
+#include <coeffs/ffields.h>
 #include <coeffs/coeffs.h>
 #include <Singular/subexpr.h>
+#include <polys/templates/p_Procs.h>
 
 #include <polys/monomials/maps.h>
@@ -988 +990 @@
 intvec* MivMatrixOrderRefine(intvec* iv, intvec* iw)
 {
-<<<<<<< HEAD
   assume((iv->length())*(iv->length()) == iw->length());
   int i,j, nR = iv->length();
   
-=======
-  assume(iv->length() == iw->length());
-  int i, nR = iv->length();
-
->>>>>>> f533f6f7667328bccb271b19b2f603aaebe41596
   intvec* ivm = new intvec(nR*nR);
 
@@ -2460 +2456 @@
   //rChangeCurrRing(r);
 }
-//unused
-#if 0
-static ring VMrDefault1(intvec* va)
-{
-
-  if ((currRing->ppNoether)!=NULL)
-  {
-    pDelete(&(currRing->ppNoether));
-  }
-  if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) ||
-      ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data))))
-  {
-    sLastPrinted.CleanUp();
-  }
-
-  ring r = (ring) omAlloc0Bin(sip_sring_bin);
-  int i, nv = currRing->N;
-
-  r->cf  = currRing->cf;
-  r->N   = currRing->N;
-
-  int nb = 4;
-
-  //names
-  char* Q; // In order to avoid the corrupted memory, do not change.
-  r->names = (char **) omAlloc0(nv * sizeof(char_ptr));
-  for(i=0; i<nv; i++)
-  {
-    Q = currRing->names[i];
-    r->names[i]  = omStrDup(Q);
-  }
-
-  /*weights: entries for 3 blocks: NULL Made:???*/
-  r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
-  r->wvhdl[0] = (int*) omAlloc(nv*sizeof(int));
-  for(i=0; i<nv; i++)
-    r->wvhdl[0][i] = (*va)[i];
-
-  /* order: a,lp,C,0 */
-  r->order = (int *) omAlloc(nb * sizeof(int *));
-  r->block0 = (int *)omAlloc0(nb * sizeof(int *));
-  r->block1 = (int *)omAlloc0(nb * sizeof(int *));
-
-  // ringorder a for the first block: var 1..nv
-  r->order[0]  = ringorder_a;
-  r->block0[0] = 1;
-  r->block1[0] = nv;
-
-  // ringorder lp for the second block: var 1..nv
-  r->order[1]  = ringorder_lp;
-  r->block0[1] = 1;
-  r->block1[1] = nv;
-
-  // ringorder C for the third block
-  // it is very important within "idLift",
-  // especially, by ring syz_ring=rCurrRingAssure_SyzComp();
-  // therefore, nb  must be (nBlocks(currRing)  + 1)
-  r->order[2]  = ringorder_C;
-
-  // the last block: everything is 0
-  r->order[3]  = 0;
-
-  // polynomial ring
-  r->OrdSgn    = 1;
-
-  // complete ring intializations
-
-  rComplete(r);
-
-  //rChangeCurrRing(r);
-  return r;
-}
-#endif
+
 /****************************************************************
  * define and execute a new ring with ordering (a(va),Wp(vb),C) *
  * **************************************************************/
-
 static ring VMrRefine(intvec* va, intvec* vb)
 {
@@ -2610 +2533 @@
 
   // complete ring intializations
-
+  
   rComplete(r);
 
@@ -2840 +2763 @@
 }
 
-
-/* define a ring with parameters und change to it */
-/* DefRingPar and DefRingParlp corrupt still memory */
+/***************************************************
+* define a ring with parameters und change to it   *
+* DefRingPar and DefRingParlp corrupt still memory *
+****************************************************/
 static void DefRingPar(intvec* va)
 {
@@ -2990 +2914 @@
 }
 
-//unused
-/**************************************************************
- * check wheather one or more components of a vector are zero *
- **************************************************************/
-//#if 0
+/*************************************************************
+ * check whether one or more components of a vector are zero *
+ *************************************************************/
 static int isNolVector(intvec* hilb)
 {
@@ -3007 +2929 @@
   return 0;
 }
-//#endif
+
+/*************************************************************
+ * check whether one or more components of a vector are <= 0 *
+ *************************************************************/
+static int isNegNolVector(intvec* hilb)
+{
+  int i;
+  for(i=hilb->length()-1; i>=0; i--)
+  {
+    if((* hilb)[i]<=0)
+    {
+      return 1;
+    }
+  }
+  return 0;
+}
+
+/**************************************************************************
+* Gomega is the initial ideal of G w. r. t. the current weight vector     *
+* curr_weight. Check whether curr_weight lies on a border of the Groebner *
+* cone, i. e. check whether a monomial is divisible by a leading monomial *
+***************************************************************************/
+static ideal middleOfCone(ideal G, ideal Gomega)
+{
+  BOOLEAN middle = FALSE;
+  int i,j,N = IDELEMS(Gomega);
+  poly p,lm,factor1,factor2;
+  PrintS("\n//** idCopy\n");
+  ideal Go = idCopy(G);
+  
+  PrintS("\n//** jetzt for-Loop!\n");
+  for(i=0; i<N; i++)
+  {
+    for(j=0; j<N; j++)
+    {
+      if(i!=j)
+      {
+        p = pCopy(Gomega->m[i]);
+        lm = pCopy(Gomega->m[j]);
+        pIter(p);
+        while(p!=NULL)
+        {
+          if(pDivisibleBy(lm,p))
+          {
+            if(middle == FALSE)
+            {
+              middle = TRUE;
+            }
+            factor1 = singclap_pdivide(pHead(p),lm,currRing);
+            factor2 = pMult(pCopy(factor1),pCopy(Go->m[j]));
+            pDelete(&factor1);
+            Go->m[i] = pAdd((Go->m[i]),pNeg(pCopy(factor2)));
+            pDelete(&factor2);
+          }
+          pIter(p);
+        }
+        pDelete(&lm);
+        pDelete(&p);
+      }
+    }
+  }
+  
+  //PrintS("\n//** jetzt Delete!\n");
+  //pDelete(&p);
+  //pDelete(&factor);
+  //pDelete(&lm);
+  if(middle == TRUE)
+  {
+    PrintS("\n//** middle TRUE!\n");
+    return Go;
+  }
+  PrintS("\n//** middle FALSE!\n");
+  idDelete(&Go);
+  return NULL; 
+}
 
 /******************************  Februar 2002  ****************************
@@ -3314 +3310 @@
   }
   return 0;
+}
+
+/*****************************************
+ * return maximal polynomial length of G *
+ *****************************************/
+static int maxlengthpoly(ideal G)
+{
+  int i,k,length=0;
+  for(i=IDELEMS(G)-1; i>=0; i--)
+  {
+    k = pLength(G->m[i]);
+    if(k>length)
+    {
+      length = k;
+    }
+  }
+  return length;
 }
 
@@ -4325 +4338 @@
  * compute a next weight vector *
  ********************************/
-static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight, intvec* target_weight, int weight_rad, int pert_deg)
-{
-  int i, weight_norm;
-  //int randCount=0;
-  int nV = currRing->N;
+static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight,
+               intvec* target_weight, int weight_rad, int pert_deg)
+{
+  assume(currRing != NULL && curr_weight != NULL &&
+         target_weight != NULL && G->m[0] != NULL);
+
+  int i,weight_norm,nV = currRing->N;
   intvec* next_weight2;
   intvec* next_weight22 = new intvec(nV);
   intvec* result = new intvec(nV);
 
-  //compute a perturbed next weight vector "next_weight1"
   intvec* next_weight1 =MkInterRedNextWeight(curr_weight,target_weight,G);
   //compute a random next weight vector "next_weight2"
@@ -4342 +4356 @@
     while(weight_norm == 0)
     {
+
       for(i=0; i<nV; i++)
       {
@@ -4349 +4364 @@
       weight_norm = 1 + floor(sqrt(weight_norm));
     }
+
     for(i=0; i<nV; i++)
     {
@@ -4360 +4376 @@
       }
     }
+
     if(test_w_in_ConeCC(G, next_weight22) == 1)
     {
@@ -4370 +4387 @@
         }
         i = 0;
+        /* reduce the size of the maximal entry of the vector*/
         while(test_w_in_ConeCC(G,next_weight22))
         {
@@ -4381 +4399 @@
     }
   }
+  
   // compute "usual" next weight vector
   intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight, G);
@@ -4390 +4409 @@
   {
     ideal G_test1 = MwalkInitialForm(G, next_weight1);
-    if(IDELEMS(G_test1) < IDELEMS(G_test))
-    {
-      if(IDELEMS(G_test2) < IDELEMS(G_test1))
+    if(G_test1->m[0] != NULL && maxlengthpoly(G_test1) < maxlengthpoly(G_test))//if(IDELEMS(G_test1) < IDELEMS(G_test))
+    {
+      if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test1)) //if(IDELEMS(G_test2) < IDELEMS(G_test1))
       {
         // |G_test2| < |G_test1| < |G_test|
@@ -4407 +4426 @@
           (*result)[i] = (*next_weight1)[i];
         }
-      }
+      }    
     }
     else
     {
-      if(IDELEMS(G_test2) < IDELEMS(G_test)) // |G_test2| < |G_test| <= |G_test1|
+      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|
       {
         for(i=0; i<nV; i++)
@@ -4432 +4451 @@
   {
     Overflow_Error = FALSE;
-    if(IDELEMS(G_test2) < IDELEMS(G_test))
+    if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test))//if(IDELEMS(G_test2) < IDELEMS(G_test))
     {
       for(i=1; i<nV; i++)
@@ -4447 +4466 @@
     }
   }
+  PrintS("\n MWalkRandomNextWeight: Ende ok!\n");
   idDelete(&G_test);
   idDelete(&G_test2);
@@ -5023 +5043 @@
   int nV = baseRing->N;
 
-  ideal Gomega, M, F, Gomega1, Gomega2, M1; //, F1;
+  ideal Gomega, M, F, FF, Gomega1, Gomega2, M1; //, F1;
   ring newRing;
   ring XXRing = baseRing;
@@ -5070 +5090 @@
   }
   rChangeCurrRing(newRing);
-  ideal G = MstdCC(idrMoveR(Go,baseRing,currRing));
+  //ideal G = MstdCC(idrMoveR(Go,baseRing,currRing));
+  ideal G = idrMoveR(Go,baseRing,currRing);
   baseRing = currRing;
 #ifdef TIME_TEST
@@ -5102 +5123 @@
     }
 //#endif
+
+
+    if(reduction == 0)
+    {
+      PrintS("\n//** Mwalk: test middle of cone!\n");
+      FF = middleOfCone(G,Gomega);
+      PrintS("\n//** Mwalk: Test F!\n");
+      if( FF != NULL)
+      {
+        idDelete(&G);
+        G = idCopy(FF);
+        idDelete(&FF);
+        PrintS("\n//** Mwalk: FF nicht NULL! Compue next vector.\n");
+        goto NEXT_VECTOR;
+      } 
+    }
+
 #ifndef  BUCHBERGER_ALG
     if(isNolVector(curr_weight) == 0)
     {
       hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
-    }
+    }   
     else
     {
@@ -5204 +5242 @@
     to = clock();
 #endif
+    NEXT_VECTOR:
+    PrintS("\n//** Mwalk: NEXT_VECTOR!\n");
     intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
 #ifdef TIME_TEST
@@ -5239 +5279 @@
         idString(Gomega1, "//** Mwalk: Gomega");
       }
+      PrintS("\n //** Mwalk: kStd(Gomega)");
 //#endif
       M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL);
@@ -5252 +5293 @@
       Gomega2 = idrMoveR(Gomega1, newRing,currRing);
       idDelete(&Gomega1);
+      PrintS("\n //** Mwalk: MLifttwoIdeal");
       F = MLifttwoIdeal(Gomega2, M1, G);
 //#ifdef CHECK_IDEAL_MWALK
@@ -5270 +5312 @@
       to = clock();
 #endif
+      //PrintS("\n //**Mwalk: Interreduce");
       //interreduce the Groebner basis <G> w.r.t. currRing
-      G = kInterRedCC(G,NULL);
+      //G = kInterRedCC(G,NULL);
 #ifdef TIME_TEST
       tred = tred + clock() - to;
@@ -5369 +5412 @@
   }
   rChangeCurrRing(newRing);
-  ideal G = MstdCC(idrMoveR(Go,baseRing,currRing));
+  //ideal G = MstdCC(idrMoveR(Go,baseRing,currRing));
+  ideal G = idrMoveR(Go,baseRing,currRing);
   baseRing = currRing;
 #ifdef TIME_TEST
@@ -5405 +5449 @@
     {
       hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
-    }
+    }   
     else
     {
@@ -5514 +5558 @@
     }
 //#endif
-    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)
+    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)
     {
 #ifdef CHECK_IDEAL_MWALK
@@ -5565 +5609 @@
         idString(G,"//** Mrwalk: G");
       }
-/*#endif
+///#endif
       if(si_opt_1 == (Sy_bit(OPT_REDSB)))
-      {*/
-        G = kInterRedCC(G,NULL); //interreduce the Groebner basis <G> w.r.t. currRing
-//      }
+      {
+        //G = kInterRedCC(G,NULL); //interreduce the Groebner basis <G> w.r.t. currRing
+      }
 #ifdef TIME_TEST
       tred = tred + clock() - to;
@@ -5861 +5905 @@
     nstep ++;
     to = clock();
-    /* compute an initial form ideal of <G> w.r.t. the weight vector
-       "curr_weight" */
+    // compute an initial form ideal of <G> w.r.t. the weight vector
+    // "curr_weight"
     Gomega = MwalkInitialForm(G, curr_weight);
 //#ifdef CHECK_IDEAL_MWALK
-      if(printout > 1)
-      {
-        idString(Gomega,"//** Mpwalk: Gomega");
-      }
+    if(printout > 1)
+    {
+      idString(Gomega,"//** Mpwalk: Gomega");
+    }
 //#endif
 
@@ -5916 +5960 @@
     tim = clock();
     to = clock();
-    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
+    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
 #ifdef  BUCHBERGER_ALG
     M = MstdhomCC(Gomega1);
@@ -5940 +5984 @@
       tstd=tstd+clock()-to;
 
-    /* change the ring to oldRing */
+    // change the ring to oldRing
     rChangeCurrRing(oldRing);
     M1 =  idrMoveR(M, newRing,currRing);
@@ -5966 +6010 @@
     idDelete(&G);
 
-    /* change the ring to newRing */
+    // change the ring to newRing
     rChangeCurrRing(newRing);
     F1 = idrMoveR(F, oldRing,currRing);
-
+    //G = idrMoveR(F,oldRing,currRing);
     to=clock();
-    /* reduce the Groebner basis <G> w.r.t. new ring */
+    PrintS("\n //** Mpwalk: reduce the Groebner basis.\n");
     G = kInterRedCC(F1, NULL);
     if(endwalks != 1)
@@ -5988 +6032 @@
     tnw=tnw+clock()-to;
 //#ifdef PRINT_VECTORS
-    if(printout > 2)
+    if(printout > 0)
     {
       MivString(curr_weight, target_weight, next_weight);
@@ -6016 +6060 @@
 
     delete next_weight;
-  }//while
+  }//end of while-loop
 
   if(tp_deg != 1)
@@ -6264 +6308 @@
     nstep ++;
     to = clock();
-    /* compute an initial form ideal of <G> w.r.t. the weight vector
-       "curr_weight" */
+    // compute an initial form ideal of <G> w.r.t. the weight vector
+    // "curr_weight"
     Gomega = MwalkInitialForm(G, curr_weight);
 //#ifdef CHECK_IDEAL_MWALK
-    if(printout > 1)
+    if(printout > 0)
     {
       idString(Gomega,"//** Mprwalk: Gomega");
@@ -6308 +6352 @@
     newRing = currRing;
     Gomega1 = idrMoveR(Gomega, oldRing,currRing);
-
+    PrintS("\n //** Mprwalk: Changed ring.\n");
 #ifdef ENDWALKS
     if(endwalks==1)
@@ -6325 +6369 @@
     tim = clock();
     to = clock();
-    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
+    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
 #ifdef  BUCHBERGER_ALG
     M = MstdhomCC(Gomega1);
@@ -6333 +6377 @@
 #endif
 //#ifdef CHECK_IDEAL_MWALK
+    Print("\n //** Mprwalk: M");
     if(printout > 2)
     {
@@ -6364 +6409 @@
     else
       xtlift=clock()-to;
-
+    Print("\n //** Mprwalk: F.\n");
 //#ifdef CHECK_IDEAL_MWALK
     if(printout > 2)
@@ -6376 +6421 @@
     idDelete(&G);
 
-    /* change the ring to newRing */
+    // change the ring to newRing
     rChangeCurrRing(newRing);
-    F1 = idrMoveR(F, oldRing,currRing);
-
+    //F1 = idrMoveR(F, oldRing,currRing);
+    F1 = idrMoveR(F,oldRing,currRing);
     to=clock();
-    /* reduce the Groebner basis <G> w.r.t. new ring */
+    Print("\n //** Mprwalk: reducing the Groebner basis.\n");
     G = kInterRedCC(F1, NULL);
     if(endwalks != 1)
@@ -6394 +6439 @@
 
     to=clock();
-    // compute a next weight vector
+    Print("\n //** Mprwalk: compute a next weight vector.\n");
     next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
     if(polylength > 0)
     {
-      //there is a polynomial in Gomega with at least 3 monomials,
-      //low-dimensional facet of the cone
+      Print("\n //**Mprwalk: there is a polynomial in Gomega with at least 3 monomials, low-dimensional facet of the cone.\n");
       delete next_weight;
       next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, op_deg);
@@ -6405 +6449 @@
     tnw=tnw+clock()-to;
 //#ifdef PRINT_VECTORS
-    if(printout > 2)
+    if(printout > 0)
     {
       MivString(curr_weight, target_weight, next_weight);
@@ -6763 +6807 @@
       Gt = idrMoveR(G, oRing,currRing);
 
-      if(test_w_in_ConeCC(Gt, omega2) == 1) {
+      if(test_w_in_ConeCC(Gt, omega2) == 1)
+      {
         delete omega2;
         delete next_vect;
@@ -6771 +6816 @@
           Print("\n//** rec_fractal_call: Leaving the %d-th recursion with %d steps.\n",
               nlev, nwalks);
-          //Print(" ** Overflow_Error? (%d)", Overflow_Error);
         }
         return (Gt);
@@ -6777 +6821 @@
       else
       {
-        //ivString(omega2, "tau'");
-        //Print("\n//  tau' doesn't stay in the correct cone!!");
-
+        if(printout > 0)
+        {
+          Print("\n//** rec_fractal_call:  tau' doesn't stay in the correct cone!!");
+        }
 #ifndef  MSTDCC_FRACTAL
-        //07.08.03
-        //ivString(Xtau, "old Xtau");
         intvec* Xtautmp;
         if(ivtarget->length() == nV)
@@ -6799 +6842 @@
         if(MivSame(Xtau, Xtautmp) == 1)
         {
-          //PrintS("\n// Update vectors are equal to the old vectors!!");
+          if(printout > 0)
+          {
+            Print("\n//** rec_fractal_call: Update vectors are equal to the old vectors!!");
+          }
           delete Xtautmp;
           goto FRACTAL_MSTDCC;
@@ -6806 +6852 @@
         Xtau = Xtautmp;
         Xtautmp = NULL;
-        //ivString(Xtau, "new  Xtau");
 
         for(i=nV-1; i>=0; i--)
           (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i];
 
-        //Print("\n//  ring tau = %s;", rString(currRing));
         rChangeCurrRing(oRing);
         G = idrMoveR(Gt, testring,currRing);
@@ -6851 +6895 @@
         Gt = idrMoveR(G, oRing,currRing);
 
-        delete Xtau;
-        Xtau = NewVectorlp(Gt);
+        // perturb the original target vector w.r.t. the current GB
+        if(ivtarget->length() == nV)
+        {
+          delete Xtau;
+          Xtau = NewVectorlp(Gt);
+        }
+        else
+        {
+          delete Xtau;
+          Xtau = Mfpertvector(Gt,ivtarget);
+        }
 
         rChangeCurrRing(oRing);
@@ -6893 +6946 @@
     else
       hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
-#endif // BUCHBERGER_ALG
+#endif
+
     if(ivtarget->length() == nV)
     {
@@ -6910 +6964 @@
     // Fractal walk with the alternative recursion
     // alternative recursion
-    // if(nlev == nV || lengthpoly(Gomega1) == 0)
     if(nlev == Xnlev || lengthpoly(Gomega1) == 0)
-      //if(nlev == nV) // blind recursion
     {
       to=clock();
@@ -6920 +6972 @@
       Gresult =kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,omega);
       delete hilb_func;
-#endif // BUCHBERGER_ALG
+#endif
       xtstd=xtstd+clock()-to;
     }
@@ -6953 +7005 @@
 
     rChangeCurrRing(new_ring);
-    F1 = idrMoveR(F, oRing,currRing);
-
+    //F1 = idrMoveR(F, oRing,currRing);
+    G = idrMoveR(F,oRing,currRing);
     to=clock();
     /* Interreduce G */
-    G = kInterRedCC(F1, NULL);
+   // G = kInterRedCC(F1, NULL);
     xtred=xtred+clock()-to;
-    idDelete(&F1);
+    //idDelete(&F1);
   }
 }
@@ -7017 +7069 @@
   {
 #ifdef FIRST_STEP_FRACTAL
-    // perturb the current weight vector only on the top level or
-    // after perturbation of the both vectors, nlev = 2 as the top level
+    /*
+    perturb the current weight vector only on the top level or
+    after perturbation of the both vectors, nlev = 2 as the top level
+    */
     if((nlev == 1 && Xcall == 0) || (nlev == 2 && Xngleich == 1))
       if(islengthpoly2(G) == 1)
@@ -7033 +7087 @@
     /* determine the next border */
     next_vect = MkInterRedNextWeight(omega,omega2,G);
-    if(polylength > 0)
-    {
-      //there is a polynomial in Gomega with at least 3 monomials,
-      //low-dimensional facet of the cone
+    if(polylength > 0 && G->m[0] != NULL)
+    {
+      /*
+      there is a polynomial in Gomega with at least 3 monomials,
+      low-dimensional facet of the cone
+      */
       delete next_vect;
-      next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,1+nlev);
-      if(isNolVector(next_vect))
+      next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,nlev);
+      if(isNegNolVector(next_vect)==1)
       {
         delete next_vect;
@@ -7057 +7113 @@
         if(printout > 0)
         {
-          Print("\n//** rec_r_fractal_call: Perturb the both vectors with degree %d.",nlev);
+          Print("\n//** rec_r_fractal_call: Perturb both vectors with degree %d.",nlev);
           //idElements(G, "G");
         }
@@ -7109 +7165 @@
         to=clock();
 
-        /* to avoid the value of Overflow_Error that occur in Mfpertvector*/
+        /*
+        to avoid the value of Overflow_Error that occur in
+        Mfpertvector
+        */
         Overflow_Error = FALSE;
 
         next_vect = MkInterRedNextWeight(omega,omega2,G);
-        if(polylength > 0)
+        if(G->m[0] != NULL && polylength > 0)
         {
-          //there is a polynomial in Gomega with at least 3 monomials,
-          //low-dimensional facet of the cone
+          /*
+          there is a polynomial in Gomega with at least 3 monomials,
+          low-dimensional facet of the cone
+          */
           delete next_vect;
-          next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,1+nlev);
-          if(isNolVector(next_vect))
+          next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,nlev);
+          if(isNegNolVector(next_vect)==1)
           {
             delete next_vect;
@@ -7228 +7289 @@
           Print("\n//** rec_r_fractal_call: target weight doesn't stay in the correct cone.\n");
         }
+
 #ifndef  MSTDCC_FRACTAL
         //ivString(Xtau, "old Xtau");
@@ -7268 +7330 @@
         if(printout > 0)
         {
-          Print("\n//** rec_r_fractal_call: apply Buchberge's algorithm in ring = %s.\n",
+          Print("\n//** rec_r_fractal_call: apply Buchberger's algorithm in ring = %s.\n",
                 rString(currRing));
         }
@@ -7298 +7360 @@
         Gt = idrMoveR(G, oRing,currRing);
 
-        delete Xtau;
-        Xtau = NewVectorlp(Gt);
+        // perturb the original target vector w.r.t. the current GB
+        if(ivtarget->length() == nV)
+        {
+          delete Xtau;
+          Xtau = NewVectorlp(Gt);
+        }
+        else
+        {
+          delete Xtau;
+          Xtau = Mfpertvector(Gt,ivtarget);
+        }
 
         rChangeCurrRing(oRing);
@@ -7319 +7390 @@
         return(G);
       }
-    }
-
-    for(i=nV-1; i>=0; i--) {
+    } //end of if(MivComp(next_vect, XivNull) == 1)
+
+    for(i=nV-1; i>=0; i--)
+    {
       (*altomega)[i] = (*omega)[i];
       (*omega)[i] = (*next_vect)[i];
@@ -7374 +7446 @@
       rChangeCurrRing(oRing);
       Gomega1 = idrMoveR(Gomega1, oRing,currRing);
-      Gresult = rec_fractal_call(idCopy(Gomega1),nlev+1,omega,printout);
+      Gresult = rec_r_fractal_call(idCopy(Gomega1),nlev+1,omega,weight_rad,printout);
     }
     if(printout > 2)
@@ -7402 +7474 @@
 
     rChangeCurrRing(new_ring);
-    F1 = idrMoveR(F, oRing,currRing);
-
+    //F1 = idrMoveR(F, oRing,currRing);
+    G = idrMoveR(F,oRing,currRing);
     to=clock();
     // Interreduce G
-    G = kInterRedCC(F1, NULL);
+    //G = kInterRedCC(F1, NULL);
     xtred=xtred+clock()-to;
-    idDelete(&F1);
+    //idDelete(&F1);
   }
 }
@@ -8916 +8988 @@
 
   //Print("\n// \"Mpwalk\" (1,%d) took %d steps and %.2f sec. Overflow_Error (%d)", tp_deg, nwalk, ((double) clock()-tinput)/1000000, nOverflow_Error);
-<<<<<<< HEAD
-=======
-
-  return(result);
-}
-
-/*******************************************************
- * THE PERTURBATION WALK ALGORITHM WITH RANDOM ELEMENT *
- *******************************************************/
-ideal Mprwalk(ideal Go, intvec* curr_weight, intvec* target_weight, int weight_rad, int op_deg, int tp_deg, ring baseRing)
-{
-  BITSET save1 = si_opt_1; // save current options
-  si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
-  Set_Error(FALSE);
-  Overflow_Error = FALSE;
-#ifdef TIME_TEST
-  clock_t tinput=0, tostd=0, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
-  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0;
-  tinput = clock();
-  clock_t tim;
-#endif
-  int i,nwalk,nV = baseRing->N;
-
-  ideal G, Gomega, M, F, Gomega1, Gomega2, M1;
-  ring newRing;
-  ring XXRing = baseRing;
-  intvec* exivlp = Mivlp(nV);
-  intvec* orig_target = target_weight;
-  intvec* pert_target_vector = target_weight;
-  intvec* ivNull = new intvec(nV);
-  intvec* tmp_weight = new intvec(nV);
-#ifdef CHECK_IDEAL_MWALK
-  poly p;
-#endif
-  for(i=0; i<nV; i++)
-  {
-    (*tmp_weight)[i] = (*curr_weight)[i];
-  }
-#ifndef BUCHBERGER_ALG
-  intvec* hilb_func;
-   // to avoid (1,0,...,0) as the target vector
-  intvec* last_omega = new intvec(nV);
-  for(i=0 i<nV; i++)
-  {
-    (*last_omega)[i] = 1;
-  }
-  (*last_omega)[0] = 10000;
-#endif
-  baseRing = currRing;
-  newRing = VMrDefault(curr_weight);
-  rChangeCurrRing(newRing);
-  G = idrMoveR(Go,baseRing,currRing);
-#ifdef TIME_TEST
-  to = clock();
-#endif
-  G = kStd(G,NULL,testHomog,NULL,NULL,0,0,NULL);
-  idSkipZeroes(G);
-#ifdef TIME_TEST
-  tostd = tostd + to - clock();
-#endif
-#ifdef CHECK_IDEAL_MWALK
-  idString(G,"G");
-#endif
-  if(op_deg >1)
-  {
-    if(MivComp(curr_weight,MivUnit(nV)) == 1) //ring order is "dp"
-    {
-      curr_weight = MPertVectors(G, MivMatrixOrderdp(nV), op_deg);
-    }
-    else //ring order is not "dp"
-    {
-      curr_weight = MPertVectors(G, MivMatrixOrder(curr_weight), op_deg);
-    }
-  }
-  baseRing = currRing;
-  if(tp_deg > 1 && tp_deg <= nV)
-  {
-    pert_target_vector = target_weight;
-  }
-#ifdef CHECK_IDEAL_MWALK
-  ivString(curr_weight, "new curr_weight");
-  ivString(target_weight, "new target_weight");
-#endif
-  nwalk = 0;
-  while(1)
-  {
-    nwalk ++;
-#ifdef TIME_TEST
-    to = clock();
-#endif
-    Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
-#ifdef TIME_TEST
-    tif = tif + clock()-to; //time for computing initial form ideal
-#endif
-#ifdef CHECK_IDEAL_MWALK
-    idString(Gomega,"Gomega");
-#endif
-#ifndef  BUCHBERGER_ALG
-    if(isNolVector(curr_weight) == 0)
-    {
-      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
-    }
-    else
-    {
-      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
-    }
-#endif
-    if(nwalk == 1)
-    {
-      newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
-    }
-    else
-    {
-      newRing = VMrRefine(curr_weight,target_weight); //define a new ring with ordering "(a(curr_weight),Wp(target_weight))"
-    }
-    rChangeCurrRing(newRing);
-    Gomega1 = idrMoveR(Gomega, baseRing,currRing);
-    idDelete(&Gomega);
-    // compute a Groebner basis of <Gomega> w.r.t. "newRing"
-#ifdef TIME_TEST
-    to = clock();
-#endif
-#ifndef  BUCHBERGER_ALG
-    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
-    delete hilb_func;
-#else
-    M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL);
-#endif
-    idSkipZeroes(M);
-#ifdef TIME_TEST
-    tstd = tstd + clock() - to;
-#endif
-#ifdef CHECK_IDEAL_MWALK
-    idString(M, "M");
-#endif
-    //change the ring to baseRing
-    rChangeCurrRing(baseRing);
-    M1 =  idrMoveR(M, newRing,currRing);
-    idDelete(&M);
-    Gomega2 = idrMoveR(Gomega1, newRing,currRing);
-    idDelete(&Gomega1);
-    to = clock();
-    // 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
-    F = MLifttwoIdeal(Gomega2, M1, G);
-    idSkipZeroes(F);
-#ifdef TIME_TEST
-    tlift = tlift + clock() - to;
-#endif
-#ifdef CHECK_IDEAL_MWALK
-    idString(F,"F");
-#endif
-    rChangeCurrRing(newRing); // change the ring to newRing
-    G = idrMoveR(F,baseRing,currRing);
-    idDelete(&F);
-    baseRing = currRing; // set baseRing equal to newRing
-#ifdef CHECK_IDEAL_MWALK
-    idString(G,"G");
-#endif
-#ifdef TIME_TEST
-    to = clock();
-#endif
-    intvec* next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, op_deg);
-#ifdef TIME_TEST
-    tnw = tnw + clock() - to;
-#endif
-#ifdef PRINT_VECTORS
-    MivString(curr_weight, target_weight, next_weight);
-#endif
-    if(Overflow_Error == TRUE)
-    {
-      PrintS("\n//**Mprwalk: OVERFLOW: The computed vector does not stay in cone, the result may be wrong.\n");
-      delete next_weight;
-      break;
-    }
-
-    if(test_w_in_ConeCC(G,target_weight) == 1 || MivComp(next_weight, ivNull) == 1)
-    {
-      delete next_weight;
-      break;
-    }
-    //update tmp_weight and curr_weight
-    for(i=nV-1; i>=0; i--)
-    {
-      (*tmp_weight)[i] = (*curr_weight)[i];
-      (*curr_weight)[i] = (*next_weight)[i];
-    }
-    delete next_weight;
-  } //end of while-loop
->>>>>>> f533f6f7667328bccb271b19b2f603aaebe41596
   Print("\n// Mprwalk took %d steps. Ring= %s;\n", nwalk, rString(currRing));
   return(result);

Singular/walk.h

@@ -52 +52 @@
 
 /* Okt -- Nov'01 */
-// compute a Groebner basis of an ideal G w.r.t. lexicographic order
+// compute a Groebner basis of an ideal G w.r.t. lexicographic order 
 //ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M);
 ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M, ring baseRing, int reduction, int printout);