Changeset a14482 in git

Timestamp:

Feb 26, 2015, 6:55:30 PM (9 years ago)

Author:

Stephan Oberfranz <oberfran@…>

Branches:

(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '38dfc5131670d387a89455159ed1e071997eec94')

Children:

e57a75c6d0050e3db7e1a87294a311d4784df43c

Parents:

eca8d9cc9b3584a2f5741c75f69099751c4cbe8f

Message:

update

Location:

Singular

Files:

• LIB/modwalk.lib (modified) (1 diff)
• walk.cc (modified) (24 diffs)

Singular/LIB/modwalk.lib

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

Singular/walk.cc

@@ -1106 +1106 @@
   if(pdeg > nV || pdeg <= 0)
   {
-    WerrorS("//** The perturbed degree is wrong!!");
+    WerrorS("//**MPertVectors: The perturbed degree is wrong!!");
     return v_null;
   }
@@ -1116 +1116 @@
   }
   mpz_t *pert_vector = (mpz_t*)omAlloc(nV*sizeof(mpz_t));
-  //mpz_t *pert_vector1 = (mpz_t*)omAlloc(nV*sizeof(mpz_t));
 
   for(i=0; i<nV; i++)
@@ -1123 +1122 @@
    // mpz_init_set_si(pert_vector1[i], (*ivtarget)[i]);
   }
-  // Calculate max1 = Max(A2)+Max(A3)+...+Max(Apdeg),
+  // Compute max1 = Max(A2)+Max(A3)+...+Max(Apdeg),
   // where the Ai are the i-te rows of the matrix target_ord.
   int ntemp, maxAi, maxA=0;
@@ -1148 +1147 @@
   }
 
-  // Calculate inveps = 1/eps, where 1/eps > totaldeg(p)*max1 for all p in G.
+  // Compute inveps = 1/eps, where 1/eps > totaldeg(p)*max1 for all p in G.
 
   intvec* ivUnit = Mivdp(nV);
@@ -1180 +1179 @@
   }
 #else
-  // PrintS("\n// the \"big\" inverse epsilon: ");
+  PrintS("\n //**MPertVectors: the \"big\" inverse epsilon:");
   mpz_out_str(stdout, 10, inveps);
 #endif
@@ -1234 +1233 @@
   if(j > nV - 1)
   {
-    // Print("\n//  MPertVectors: geaenderter vector gleich Null! \n");
+    //perturbed vector equals zero
     delete pert_vector1;
     goto CHECK_OVERFLOW;
   }
 
-// check that the perturbed weight vector lies in the Groebner cone
+  // check that the perturbed weight vector lies in the Groebner cone
   if(test_w_in_ConeCC(G,pert_vector1) != 0)
   {
-    // Print("\n//  MPertVectors: geaenderter vector liegt in Groebnerkegel! \n");
     for(i=0; i<nV; i++)
     {
@@ -1248 +1246 @@
     }
   }
-  else
-  {
-    //Print("\n// MpertVectors: geaenderter vector liegt nicht in Groebnerkegel! \n");
-  }
+
   delete pert_vector1;
 
@@ -1257 +1252 @@
   intvec* result = new intvec(nV);
 
-  /* 2147483647 is max. integer representation in SINGULAR */
+  // 2147483647 is max. integer representation in SINGULAR
   mpz_t sing_int;
   mpz_init_set_ui(sing_int,  2147483647);
@@ -1644 +1639 @@
   {
     mpz_divexact(pert_vector[i], pert_vector[i], ztmp);
-    (* result)[i] = mpz_get_si(pert_vector[i]);
-  }
-
+    (* result1)[i] = mpz_get_si(pert_vector[i]);
+  }
   j = 0;
-  for(i=0; i<nV; i++)
-  {
-    (* result1)[i] = mpz_get_si(pert_vector[i]);
-    (* result1)[i] = 0.1*(* result1)[i];
-    (* result1)[i] = floor((* result1)[i] + 0.5);
-    if((* result1)[i] == 0)
-    {
-      j++;
-    }
-  }
-  if(j > nV - 1)
-  {
-    // Print("\n//  MfPertwalk: geaenderter vector gleich Null! \n");
-    delete result1;
-    goto CHECK_OVERFLOW;
-  }
-
-// check that the perturbed weight vector lies in the Groebner cone
-  if(test_w_in_ConeCC(G,result1) != 0)
-  {
-    // Print("\n//  MfPertwalk: geaenderter vector liegt in Groebnerkegel! \n");
-    delete result;
-    result = result1;
+  while(test_w_in_ConeCC(G,result1) && j<nV)
+  {
+    j = 0;
     for(i=0; i<nV; i++)
     {
+      (* result)[i] = (* result1)[i];
       mpz_set_si(pert_vector[i], (*result1)[i]);
-    }
-  }
-  else
-  {
-    delete result1;
-    // Print("\n// Mfpertwalk: geaenderter vector liegt nicht in Groebnerkegel! \n");
-  }
+      (* result1)[i] = floor(0.1*(mpz_get_si(pert_vector[i])) + 0.5);
+      if((* result1)[i] == 0)
+      {
+        j++;
+      }
+    }
+  }
+  delete result1;
 
   CHECK_OVERFLOW:
@@ -4343 +4319 @@
   assume(currRing != NULL && curr_weight != NULL &&
          target_weight != NULL && G->m[0] != NULL);
-
-  int i,weight_norm,nV = currRing->N;
+  //PrintS("\n //**MWalkRandomNextWeight: Anfang ok!\n");
+  int i,k,nV = currRing->N;
+  long weight_norm;
   intvec* next_weight2;
-  intvec* next_weight22 = new intvec(nV);
+  //intvec* next_weight22 = new intvec(nV);
   intvec* result = new intvec(nV);
+  ideal G_test;
+  ideal G_test2;
+  BOOLEAN random = FALSE;
 
   intvec* next_weight1 =MkInterRedNextWeight(curr_weight,target_weight,G);
   //compute a random next weight vector "next_weight2"
-  while(1)
-  {
+  k = 1;
+  while(k<4)
+  {
+    k++;
     weight_norm = 0;
+    intvec* next_weight22 = new intvec(nV);
     while(weight_norm == 0)
     {
-
+      //PrintS("\nWhile WeightNorm\n");
       for(i=0; i<nV; i++)
       {
-        (*next_weight22)[i] = rand() % 60000 - 30000;
+        (*next_weight22)[i] = rand() % 60000 ;
         weight_norm = weight_norm + (*next_weight22)[i]*(*next_weight22)[i];
       }
+      //Print("\n// weight_norm vor Wurzel = %d\n", weight_norm);
       weight_norm = 1 + floor(sqrt(weight_norm));
-    }
-
+      
+    }
+//Print("\n// weight_norm nach Wurzel = %d\n", weight_norm);
     for(i=0; i<nV; i++)
     {
       if((*next_weight22)[i] < 0)
       {
-        (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);
+        (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor((weight_rad)*(*next_weight22)[i]/weight_norm);
       }
       else
       {
-        (*next_weight22)[i] = (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);
+        (*next_weight22)[i] = (*curr_weight)[i] + floor((weight_rad)*(*next_weight22)[i]/weight_norm);
       }
     }
@@ -4379 +4364 @@
     if(test_w_in_ConeCC(G, next_weight22) == 1)
     {
+      //PrintS("\n //** MWalkRandomNextWeight: test ob in Kegel.\n");
       next_weight2 = MkInterRedNextWeight(next_weight22,target_weight,G);
       if(MivAbsMax(next_weight2)>1147483647)
@@ -4395 +4381 @@
         }
       }
-      delete next_weight22;
+      random = TRUE;
       break;
     }
-  }
-  
+    delete next_weight22;
+  }
+  //PrintS("\n //**MWalkRandomNextWeight: Ende while-Schleife!\n");
+
   // compute "usual" next weight vector
   intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight, G);
-  ideal G_test = MwalkInitialForm(G, next_weight);
-  ideal G_test2 = MwalkInitialForm(G, next_weight2);
-
+  G_test = MwalkInitialForm(G, next_weight);
+  if(random == TRUE)
+  {
+    G_test2 = MwalkInitialForm(G, next_weight2);
+  }
   // compare next weights
   if(Overflow_Error == FALSE)
   {
     ideal G_test1 = MwalkInitialForm(G, next_weight1);
-    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|
+    if(G_test1->m[0] != NULL && maxlengthpoly(G_test1) < maxlengthpoly(G_test))
+    {
+      if(random == TRUE &&
+         G_test2->m[0] != NULL &&
+         maxlengthpoly(G_test2) < maxlengthpoly(G_test1))
+      {
         for(i=0; i<nV; i++)
         {
@@ -4421 +4412 @@
       else
       {
-        // |G_test1| < |G_test|, |G_test1| <= |G_test2|
         for(i=0; i<nV; i++)
         {
@@ -4430 +4420 @@
     else
     {
-      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|
+      if(random == TRUE &&
+         G_test2->m[0] != NULL &&
+         maxlengthpoly(G_test2) < maxlengthpoly(G_test))
       {
         for(i=0; i<nV; i++)
@@ -4439 +4431 @@
       else
       {
-        // |G_test| < |G_test1|, |G_test| <= |G_test2|
         for(i=0; i<nV; i++)
         {
@@ -4451 +4442 @@
   {
     Overflow_Error = FALSE;
-    if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test))//if(IDELEMS(G_test2) < IDELEMS(G_test))
+    if(random == TRUE && 
+       G_test2->m[0] != NULL &&
+       maxlengthpoly(G_test2) < maxlengthpoly(G_test))
     {
       for(i=1; i<nV; i++)
@@ -4466 +4459 @@
     }
   }
-  PrintS("\n MWalkRandomNextWeight: Ende ok!\n");
+  //PrintS("\n MWalkRandomNextWeight: Ende ok!\n");
   idDelete(&G_test);
-  idDelete(&G_test2);
+  if(random == TRUE)
+  {
+    idDelete(&G_test2);
+  }
   if(test_w_in_ConeCC(G, result) == 1)
   {
-    delete next_weight2;
+    if(random == TRUE)
+    {
+      delete next_weight2;
+    }
     delete next_weight;
     delete next_weight1;
@@ -4479 +4478 @@
   {
     delete result;
-    delete next_weight2;
+    if(random == TRUE)
+    {
+      delete next_weight2;
+    }
     delete next_weight1;
     return next_weight;
@@ -6029 +6031 @@
     to=clock();
     // compute a next weight vector
-    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
+    next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);//MkInterRedNextWeight(curr_weight,target_weight, G);
     tnw=tnw+clock()-to;
 //#ifdef PRINT_VECTORS
@@ -7095 +7097 @@
       delete next_vect;
       next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,nlev);
-      if(isNegNolVector(next_vect)==1)
+      if(isNegNolVector(next_vect)==1 || MivSame(omega,next_vect) == 1)
       {
         delete next_vect;
@@ -7180 +7182 @@
           delete next_vect;
           next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,nlev);
-          if(isNegNolVector(next_vect)==1)
+          if(isNegNolVector(next_vect)==1 || MivSame(omega,next_vect) == 1)
           {
             delete next_vect;
@@ -7239 +7241 @@
         Print("\n//** rec_r_fractal_call: Leaving the %d-th recursion with %d steps.\n",
               nlev, nwalks);
-        //Print(" ** Overflow_Error? (%d)", Overflow_Error);
+      /*
+        Print(" ** Overflow_Error? (%d)", Overflow_Error);
+        idString(G1,"//** rec_r_fractal_call: G1");
+      */
       }
       nnflow ++;
@@ -7279 +7284 @@
           Print("\n//** rec_r_fractal_call: Leaving the %d-th recursion with %d steps.\n",
                 nlev, nwalks);
-          //Print(" ** Overflow_Error? (%d)", Overflow_Error);
+        /*
+          Print(" ** Overflow_Error? (%d)", Overflow_Error);
+          idString(Gt,"//** rec_r_fractal_call: Gt");
+        */
         }
         return (Gt);
@@ -7382 +7390 @@
           Print("\n//** rec_r_fractal_call: Leaving the %d-th recursion with %d steps.\n",
                 nlev,nwalks);
-          //Print(" ** Overflow_Error? (%d)", Overflow_Error);
+         /*
+          Print(" ** Overflow_Error? (%d)", Overflow_Error);
+          idString(G,"//** rec_r_fractal_call: G");
+         */
         }
         if(Overflow_Error == TRUE)