Changeset 2c2d2e8 in git


Ignore:
Timestamp:
Jun 3, 2015, 4:34:15 PM (8 years ago)
Author:
Stephan Oberfranz <oberfran@…>
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
20c1a357575a4f15b08754d8cab675b51b6d79b0
Parents:
5638b97590e59707679f6e18a166d55e968404c922c3f2c3748d3c6aa16cfcf2ab75e3fd23da838d
Message:
update

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

Legend:

Unmodified
Added
Removed
  • Singular/LIB/gradedModules.lib

    r5638b9 r2c2d2e8  
    2323    grobj(M,w[,d])  construct a graded object (map) given by matrix M
    2424    grtest(A)       check whether A is a valid graded object
    25     grisequal(A,B)  check whether A is exactly eqal to B? TODO: isomorphic!
    2625    grdeg(M)        compute graded degrees of columns of the map M
    2726    grview(M)       view the graded structure of map M
     
    5756    matrixpres(a)      matrix presentation of direct sum of Omega^{a[i]}(i)
    5857";
     58   
     59//    grisequal(A,B)  check whether A is exactly eqal to B? TODO: isomorphic!
    5960
    6061LIB "matrix.lib"; // ?
     
    7374"direct sum_i=1^size R(-v[i]), for source and targets of graded objects"
    7475{
     76  int n = size(v);
     77  if (n == 0) { return ("0"); }
     78  ASSUME(0, n > 0 );
     79
    7580  if (R == "")
    7681  {
     
    8085  ASSUME(0, defined(R) && (R != "") );
    8186
    82   int n = size(v);
    83 
    84   if (n == 0) { return ("0"); }
    85 
    86   ASSUME(0, n > 0 );
    87  
    8887  v = -v; // NOTE: due to Mathematical meanings of Singular data
    8988 
     
    178177  string R = nameof(basering);
    179178
    180   if( typeof( N ) == "list" )
    181   {
    182     msg = msg + " resolution";  // TODO: what about chain maps?
     179  if( typeof( N ) == "list" ) // TODO: find a better DS for graded resolutions / chain map !?
     180  {
     181    int n = size(N); ASSUME(0, n > 0);
     182   
     183    string msg1 = "";
    183184    if( size(R) >= 2 )
    184185    {
    185       msg = msg + "(let R:="+R+")";
    186       R = "R";
     186      msg1 = msg1 + "(let R:="+R+")";
     187      R = "R"; // !!!
     188    }   
     189    msg1 = msg1 + ": " ;
     190
     191   
     192
     193    list arr; arr[n] = list();
     194    int exact = (1==1);
     195   
     196    int i = 1;
     197   
     198    ASSUME(1, grtest(N[i]));
     199   
     200    string dst = grsumstr(R, grrange(N[i]));
     201    string src = grsumstr(R, grdeg(N[i]));
     202   
     203    arr[i] = list(dst,  src);
     204
     205    i = i + 1;
     206   
     207    while( i <= n )
     208    {
     209      ASSUME(1, grtest(N[i]));
     210     
     211      dst = grsumstr(R, grrange(N[i]));
     212     
     213      if( exact && (src != dst) )
     214      {
     215//        "src: [" + src+ "] != [" + dst + "] :(!!";
     216        exact = (1==0);
     217      }
     218
     219      src = grsumstr(R, grdeg(N[i]));
     220     
     221      arr[i] = list(dst,  src);
     222
     223      i = i + 1;
     224    };
     225
     226    string o = "";
     227
     228    if( exact )
     229    { // complex?
     230      msg = msg + " resolution" + msg1;
     231
     232      o = "d";
     233
     234      for( i = 1; i <= n; i++ )
     235      {
     236        msg = msg + newline + arr[i][1] + " <-- "+o+"_" + string(i) + " --";
     237      };
     238     
     239      msg =  msg + newline + arr[n][2];
     240      msg = msg + ", given by maps: ";
     241    } else
     242    {
     243//      print(arr);
     244
     245      msg = msg + "-object collection";     
     246      o = "o";
     247     
     248//      for( i = 1; i <= n; i++ )
     249//      {
     250//        msg = msg + newline + arr[i][1] + " <-- "+o+"_" + string(i) + " -- " + arr[i][2];
     251//      };     
     252      msg = msg + ", given by the following maps (named here as "+o+"_[1 .. "+string(n)+"]): ";
    187253    }
    188254
    189 
    190     int i = 1;  int n = size(N);
    191     string dst;
    192    
    193     msg = msg + ": " + newline + grsumstr(R, grrange(N[i])) + " <-- d_" + string(i) + " -- " ;
    194     for( ; i < n; i++ )
    195     {
    196       dst = grsumstr(R, grdeg(N[i])) +  " <-- d_" + string(i+1) + " --";
    197       msg = msg + newline + dst;
    198     };
    199 
    200     msg = msg + newline + grsumstr(R, grdeg(N[i])) + ", given by maps: ";
    201255
    202256    print(msg);
     
    204258    for( i = 1; i <= size(N); i++ )
    205259    {
    206       "d_" + string(i) + " :"; grview(N[i]);
     260      print( o+"_" + string(i) + " :" );
     261      grview( N[i] );
    207262    };
    208263
     
    393448  }
    394449
    395 /*
    396   if( flag ) // output details in case of non-identical permutation?
    397   {
    398     // grades & ordering permutation : gr[pivot] should be sorted:
    399     "s: ", s;
    400     "gr: ", gr;
    401     "pivot: ",    pivot;
    402   }
    403 */
    404450  ASSUME(1, issorted(intvec(gr[pivot]), 1));
    405451
     
    407453}
    408454
    409 // Q@Wolfram: what should be done with zero gens!?
    410455proc grdeg(M)
    411456"USAGE:  grdeg(M), graded object M
    412457RETURN:  intvec of degrees
    413 PURPOSE: graded degrees of columns (generators) of M
     458PURPOSE: graded degrees of columns (generators) of M, describing the source of M
    414459ASSUME:  M must be a graded object (matrix/module/ideal/mapping)
    415460NOTE:    if M has zero cols it shoud have attrib(M,'degHomog') set.
     
    423468    def t = attrib(M, "degHomog"); // graded degrees
    424469
    425     if( ncols(M) == 0 )    {     return (t);    } // for now
     470    if( size(t) == 0 ){ return (t); } // ZERO!
    426471   
    427     ASSUME(0, ncols(M) == size(t) );
     472    ASSUME(2, ncols(M) == size(t) );
    428473    return (t);
    429474  }
    430  
     475
     476  if( ncols(M) == 0 ) { return (0:0); } // FIXME: Why???
     477
    431478  ASSUME(0, ncols(M) > 0);
    432 
    433479  ASSUME(0, ncols(M) == size(M) );
    434480
     
    747793/////////////////////////////////////////////////////////
    748794
    749 // Q@Woflram?
    750795proc grzero()
    751796"USAGE:  grzero()
     
    755800"
    756801{
    757  return ( grobj(module([0]), intvec(0), intvec(0)) );
     802 return ( grobj(freemodule(0), intvec(0:0), intvec(0:0)) );
    758803}
    759804example
     
    762807  ring r=32003,(x,y,z),dp;
    763808
    764   def M = grpower( grshift( grzero(), 10), 5 );
    765 
    766   print(M);
    767 
    768   grview(M);
     809
     810  grview( grobj(freemodule(0), intvec(0:0), intvec(0:0)) );
     811  grview( grobj(freemodule(0), intvec(0:0)) );
     812
     813  grview( grzero() );
     814
     815//  def M = grpower( grshift( grzero(), 3), 2 ); grview(M);
     816}
     817
     818proc grtwists(intvec v)
     819"USAGE:  grtwists(v), intvec v
     820RETURN:  graded object representing S(v[1]) + ... + S(v[size(v)])
     821PURPOSE: compute presentation of S(v[1]) + ... + S(v[size(v)])
     822EXAMPLE: example grtwists; shows an example
     823"
     824{
     825  matrix m[size(v)][0];
     826  return( grobj(m, -v) ); // will set the rank as well
     827}
     828example
     829{ "EXAMPLE:"; echo = 2;
     830
     831  ring r=32003,(x,y,z),dp;
     832 
     833  grview( grtwists ( intvec(-4, 1, 6 )) );
     834
     835  grview( grtwists ( intvec(0:0) ) );
     836}
     837
     838proc grtwist(int a, int d)
     839"USAGE:  grtwist(a,d), int a, d
     840RETURN:  graded object representing S(d)^a
     841PURPOSE: compute presentation of S(d)^a
     842EXAMPLE: example grtwist; shows an example
     843"
     844{
     845  ASSUME(0, a > 0); 
     846  def Z = grtwists( intvec(d:a) ); // will set the rank as well
     847//  ASSUME(2, grisequal(Z, grpower( grshift(grzero(), d), a ) )); // optional check
     848  return(Z);
     849}
     850example
     851{ "EXAMPLE:"; echo = 2;
     852
     853  ring r=32003,(x,y,z),dp;
     854
     855//  grview(grpower( grshift(grzero(), 10), 5 ) );
     856
     857  grview( grtwist (5, 10) );
    769858}
    770859
     
    827916  intvec b = grrange(B);
    828917  intvec c = a,b;
    829   int r = nrows(A);
    830 
    831   module T = align(module(B), r); //  T;  print(T);  nrows(T); // BUG!!!!
    832   module S = module(A), T;
     918
     919  if( (ncols(B)>0) && (size(B)>0) )
     920  {
     921    int r = nrows(A);
     922    module T = align(module(B), r); //  T;  print(T);  nrows(T); // BUG!!!!
     923    module S = module(A), T;
     924  }
     925  else { def S = A; }
    833926
    834927  intvec da = grdeg(A);
     
    836929  intvec dc = da, db;
    837930
    838   return(grobj(S, c, dc));
     931
     932  def SS = grobj(S, c, dc);
     933 
     934  ASSUME(0, size( grrange(SS) ) == (size(a) + size(b)) );
     935  ASSUME(0, size( grdeg(SS) ) == (size(da) + size(db)) );
     936  ASSUME(0, ncols( SS ) == size( grdeg(SS) ) );
     937  ASSUME(0, nrows( SS ) == size( grrange(SS) ) );
     938 
     939  return(SS);
    839940}
    840941example
     
    883984{
    884985  ASSUME(1, grtest(M) );
     986           
     987           
    885988
    886989  intvec a = grrange(M);
     990  intvec t = grdeg(M);
     991
     992  if( size(a) == 0 && size(t) == 0 )
     993  {
     994    "!! Warning: shifting '0 <- 0' leaves it as it unchanged!";
     995    return (M);   
     996  }
     997 
    887998  a = a - intvec(d:size(a));
    888 
    889   intvec t = attrib(M, "degHomog");
    890999  t = t - intvec(d:size(t));
    8911000
     
    9041013
    9051014  grview(S);
     1015
     1016  grview( grshift( grzero(), 100 ) ); // does nothing...
    9061017}
    9071018
     
    9341045example
    9351046{ "EXAMPLE:"; echo = 2;
    936   TODO
    937 }
    938 
    939 proc grtwist(int a, int d)
    940 "USAGE:  grtwist(a,d), int a, d
    941 RETURN:  graded object representing S(d)^a
    942 PURPOSE: compute presentation of S(d)^a
    943 EXAMPLE: example grtwist; shows an example
    944 "
    945 {
    946   ASSUME(0, a > 0);
    947 
    948   module Z; Z[a] = [0];
    949   intvec w = -intvec(d:a);
    950   Z = grobj(Z, w, w); // will set the rank as well
    951   ASSUME(2, grisequal(Z, grpower( grshift(grzero(), d), a ) )); // optional check
    952   return(Z);
    953 }
    954 example
    955 { "EXAMPLE:"; echo = 2;
    956 
    957   ring r=32003,(x,y,z),dp;
    958 
    959   grview(grpower( grshift(grzero(), 10), 5 ) );
    960 
    961   grview( grtwist (5, 10) );
     1047//  TODO
    9621048}
    9631049
     
    9691055"
    9701056{
    971   if( ncols(A) == 0 )
    972   {
    973     ASSUME(0, size(w) == nrows(A) );
    974     attrib(A, "isHomog", w);
    975     attrib(A, "degHomog", intvec(0)); // FIXME: just for now: no intvec of 0-size :(
    976     ASSUME(0, grtest(A) );
    977     return (A);
    978   }
     1057  ASSUME(0, size(w) >= nrows(A) );
    9791058 
    9801059  module M = module(A);
    981   ASSUME(0, size(w) >= nrows(M) );
    9821060
    9831061  attrib( M, "rank", size(w) );
    9841062  attrib( M, "isHomog", w );
    9851063
     1064  intvec @ww = 0:0;
     1065 
    9861066  if( size(#) > 0 )
    9871067  {
    9881068    ASSUME(0, typeof(#[1]) == "intvec" );
    989     ASSUME(0, size(#[1]) == ncols(M) );
    990     attrib(M, "degHomog", #[1]);
     1069   
     1070    @ww = intvec( #[1] );
     1071
     1072    if( size(@ww) != ncols(M) )
     1073    {
     1074      if( (size(M) == 0) && (ncols(M) <= 1) && (size(w) == 0) && (size(@ww) > 0) )
     1075      {
     1076        matrix m[size(w)][size(@ww)]; M = module(m); attrib( M, "isHomog", w );
     1077      }
     1078    }
     1079   
     1080    ASSUME(0, size(@ww) == ncols(M) );
    9911081  }
    9921082  else
    9931083  {
    994     ASSUME(0, /* no zero cols please! */ size(M) == ncols(M) );
    995     attrib(M, "degHomog", grdeg(M));
    996   }
     1084    if( size(M) == ncols(M) ) /* no zero cols? */
     1085    {
     1086      @ww = grdeg(M); // let us compute them all :)
     1087    } else
     1088    {
     1089      if( (ncols(M) == 1) && (size(M) == 0) )
     1090      {
     1091        M = freemodule(0);
     1092      }
     1093     
     1094      attrib( M, "rank", size(w) );
     1095      attrib( M, "isHomog", w );
     1096     
     1097//      ASSUME(0, /* PROBLEM WITH ZERO COLUMNS / THEIR DEGREES! */ (ncols(M) == 0) );
     1098    }
     1099  }
     1100 
     1101//  type(@ww);  type(M);
     1102 
     1103  ASSUME(0, size(@ww) == ncols(M) ); // ?!
     1104 
     1105  attrib(M, "degHomog", @ww);
    9971106
    9981107  ASSUME(0, grtest(M) );
     
    10201129
    10211130  // Zero object:
    1022   matrix z[3][0];
    1023   def zz = grobj( z, intvec(1,2,3) ); // intvec() ??
    1024   grview(zz); 
     1131  matrix z[3][0];  grview( grobj( z, intvec(1,2,3) ) );
     1132  grview( grobj( freemodule(0), intvec(1,2,3) ) );
     1133
     1134  matrix z1[0][3]; grview( grobj( z1, 0:0, intvec(1,2,3) ) );
     1135  grview( grobj( freemodule(0), 0:0, intvec(1,2,3) ) );
     1136
     1137  matrix z0[0][0]; grview( grobj( z0, 0:0 ) );
     1138  grview( grobj( freemodule(0), 0:0 ) );
     1139 
     1140 
    10251141
    10261142}
     
    10671183  {
    10681184    intvec T = attrib(N, "degHomog"); // graded degrees
     1185
     1186    if( (ncols(N) == 1) && (size(T) == 0) && (size(N) == 0) )
     1187    {
     1188      //  if(b) { "Input seems to be a valid graded ZERO-arrow!"; };
     1189      return (1);
     1190    }
    10691191   
    10701192    if ( ncols(N) != size(T) )
     
    11201242  // the following calls will fail due to tests in grtest:
    11211243
     1244 grobj( module([x+y, x, 0, 0], [0, x+y, y, 0]), intvec(0,0,0,0) ); // enough row weights
    11221245// grobj( module([x+y, x, 0, 0], [0, x+y, y, 0]), intvec(0,0) ); // not enough row weights
    11231246// grobj( module([x,0], [0,0,0], [0, y]), intvec(1,2,3) ); // zero column needs (otherwise optional) total degrees
    1124 // grobj( module([x,0], [0,0,0], [0, y]), intvec(1,2,3), intvec(1, 1, 1) ); // incompatible total degrees (on non-zero columns)
     1247 grobj( module([x,0], [0,0,0], [0, y]), intvec(1,2,3), intvec(2, 10, 3) ); // compatible total degrees (on non-zero columns)
     1248// grobj( module([x,0], [0,0,0], [0, y]), intvec(1,2,3), intvec(2-1, 10, 3+1) ); // incompatible total degrees (on both non-zero columns)
    11251249
    11261250}
     
    11431267    return (B);   
    11441268  }
    1145 
    11461269 
    11471270  module T; T[d] = 0;
     
    11511274example
    11521275{ "EXAMPLE:"; echo = 2;
     1276
    11531277  ring r;
    1154   matrix m[1][0];
    1155   align(m, 3);
     1278  matrix m[2][0];
     1279  type( align(m, 3) );
     1280
     1281  matrix m[0][2];
     1282  type( align(m, 3) );
    11561283
    11571284
     
    11801307}
    11811308
    1182 
    1183 proc grtwists(intvec v)
    1184 "USAGE:  grtwists(v), intvec v
    1185 RETURN:  graded object representing S(v[1]) + ... + S(v[size(v)])
    1186 PURPOSE: compute presentation of S(v[1]) + ... + S(v[size(v)])
    1187 EXAMPLE: example grtwists; shows an example
    1188 "
    1189 {
    1190   int l = size(v);
    1191   module Z; Z[l] = [0];
    1192   Z = grobj(Z, -v, -v); // will set the rank as well
    1193   return(Z);
    1194 }
    1195 example
    1196 { "EXAMPLE:"; echo = 2;
    1197 
    1198   ring r=32003,(x,y,z),dp;
    1199  
    1200   grview( grtwists ( intvec(-4, 1, 6 )) );
    1201 }
    1202 
    1203 
    12041309proc grsyz(A)
    12051310"USAGE:  grsyz(M), graded object M
     
    12101315{
    12111316  ASSUME(1, grtest(A));
    1212   intvec v = grdeg(A);
    1213 
    1214 //  grrange(A);
    1215  
    1216   module M = syz(A);
    1217 //  print(M);
    1218  
    1219   if( size(M) > 0 ) { return( grobj( M, v ) ); }
    1220 
    1221   // zero syzygy?
    1222   return( grtwists(-v) ); // ???
     1317  return( grobj( syz(A), grdeg(A) ) );
     1318
     1319//   if( size(syz(A)) == 0 ) : zero syzygy? //  return( grtwists( -grdeg(A) ) );
    12231320}
    12241321example
     
    12271324  ring r=32003,(x,y,z),dp;
    12281325
    1229   module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) );
     1326  module A = grgroebner( grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ) );
    12301327  grview(A);
    12311328 
    1232   module B = grgroebner(A);
    1233   grview(B);
    1234 
    1235   module C = grsyz(B); 
    1236   grview(C);
     1329  grview(grsyz(A));
     1330 
     1331  module X = grgroebner( grobj( module([x]), intvec(2) ) );
     1332  grview(X);
     1333
     1334  // syzygy module should be zero!
     1335  grview(grsyz(X));
     1336 
     1337 
    12371338}
    12381339
     
    12481349  ASSUME(1, grtest(B));
    12491350
    1250   // TODO: ASSUME(0, grdeg() == grrange());!!!
     1351  intvec a = grdeg(A);
     1352  intvec b = grrange(B);
     1353 
     1354  ASSUME(0, (size(a) == size(b)) && (a == b));  // == for intvec :(
    12511355
    12521356  return ( grobj( A*B, grrange(A), grdeg(B) ) );
     
    12761380
    12771381
     1382// TODO: think about a proper data structure for a graded resolution!?
    12781383proc grres(def A, int l, list #)
    12791384"USAGE:  grres(M, l[, b]), graded object M, int l, int b
     
    12911396  if(b) { list r = res(A, l, #[1]); } else { list r = res(A, l); }
    12921397
    1293 //  r;  v;
    1294 
    12951398  l = size(r);
    12961399 
    1297   int i; module m;
     1400  int i;
    12981401 
    12991402  for ( i = 1; i <= l; i++ )
    13001403  {
    1301     if( size(r[i]) == 0 ){ r[i] = grtwists(-v); i++; break;   }
     1404    if( size(r[i]) == 0 )
     1405    {
     1406      r[i] = grobj(freemodule(0), v); // grtwists(-v);
     1407      i++;
     1408      break;
     1409    }
    13021410
    13031411    r[i] = grobj(r[i], v); v = grdeg(r[i]);
    13041412  }
     1413 
    13051414  i = i-1;
    13061415
     
    13261435
    13271436proc grtranspose(def M)
    1328 "
    1329 USAGE:   grtranspose(M), graded object M
     1437"USAGE:   grtranspose(M), graded object M
    13301438RETURN:  graded object
    13311439PURPOSE: graded transpose of M
     
    13511459 
    13521460
    1353   module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); A = grgroebner(A); grview(A);
    1354 
    1355   "graded transpose: "; module B = grtranspose(A); grview( B ); print(B);
    1356 
    1357   "... syzygy: "; module C = grsyz(B); grview(C);
    1358 
    1359   "... transposed: "; module D = grtranspose(C); grview( D ); print (D);
    1360 
    1361   "... and back to presentation: "; module E = grsyz( D ); grview(E); print(E);
    1362 
    1363   module F = grgens( E ); grview(F); print(F);
    1364  
    1365   module G = grpres( F ); grview(G); print(G);
     1461  // Corner cases: 0 <- 0!
     1462  module Z = grzero(); grview(Z);
     1463  grview( grtranspose( Z ) );
     1464
     1465 
     1466  // Corner cases: * <- 0
     1467  matrix M1[3][0];
     1468 
     1469  module Z1 = grobj( M1, intvec(-1, 0, 1) ); grview(Z1);
     1470  grview( grtranspose( Z1 ) );
     1471 
     1472 
     1473  // Corner cases: 0 <- *
     1474  matrix M2[0][3];
     1475
     1476  module Z2 = grobj( M2, 0:0, intvec(-1, 0, 1) ); grview(Z2);
     1477  grview( grtranspose( Z2 ) );
     1478 
    13661479}
    13671480
    13681481
    13691482proc grgens(def M)
    1370 "
    1371 USAGE:   grgens(M), graded object M (map)
     1483"USAGE:   grgens(M), graded object M (map)
    13721484RETURN:  graded object 
    13731485PURPOSE: try compute graded generators of coker(M) and return them as columns
     
    14041516
    14051517  grview( B ); print(B); // Ups :( != A
     1518
     1519  grview( grgens( grzero() ) );
    14061520 
    14071521}
     
    14091523
    14101524proc grpres(def M)
    1411 "
    1412 USAGE:   grpres(M), graded object M (submodule gens)
     1525"USAGE:   grpres(M), graded object M (submodule gens)
    14131526RETURN:  graded module (via coker)
    14141527PURPOSE: compute graded presentation matrix of submodule generated by columns of M
     
    14181531  ASSUME(1, grtest(M) );
    14191532
    1420   module N = grsyz(M);
     1533  def N = grsyz(M);
    14211534
    14221535//  ASSUME(3, grisequal( M, grgens( N ) ) );
     
    14291542  ring r=32003,(x,y,z),dp;
    14301543
    1431   module M = grtwists( intvec(-2, 0, 4, 4) ); grview(M);
    1432 
    1433   module N = grgens(M); grview( N ); print(N);
    1434 
    1435   module L = grpres( N ); grview( L ); print(L);
    1436 
    1437 
    1438   module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) );
    1439 
    1440   A = grgroebner(A); grview(A);
    1441 
    1442   module B = grgens(A); grview( B ); print(B);
    1443 
    1444   module C = grpres( B ); grview( C ); print(C);
     1544  def A = grgroebner( grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ) );
     1545  grview(A);
     1546
     1547  "graded transpose: "; def B = grtranspose(A); grview( B ); print(B);
     1548
     1549  "... syzygy: "; def C = grsyz(B); grview(C);
     1550
     1551  "... transposed: "; def D = grtranspose(C); grview( D ); print (D);
     1552
     1553  "... and back to presentation: "; def E = grsyz( D ); grview(E); print(E);
     1554
     1555  def F = grgens( E ); grview(F); print(F);
     1556  def G = grpres( F ); grview(G); print(G);
     1557
     1558
     1559  def M = grtwists( intvec(-2, 0, 4, 4) ); grview(M);
     1560 
     1561  def N = grgens(M); grview( N ); print(N);
     1562 
     1563  def L = grpres( N ); grview( L ); print(L); 
    14451564}
    14461565
     
    15041623}
    15051624
    1506 //
    1507 
    1508 
     1625// TODO: remove the following?
    15091626proc KeneshlouMatrixPresentation(intvec a)
    1510 "USAGE: intvec a.
     1627"USAGE: KeneshlouMatrixPresentation(intvec a), intvec a.
    15111628RETURN: graded object
    15121629PURPOSE: matrix presentation for direct sum of omega^a[i](i) in form of a graded object
     
    16751792
    16761793proc grrange(M)
     1794"USAGE: grrange(M), graded object M
     1795RETURN: intvec
     1796PURPOSE: get weights of module units, thus describing the target of M
     1797EXAMPLE: example grrange; shows an example
     1798"
    16771799{
    16781800//  ASSUME(1, grtest(M)); // Leads to recursive call due to grtest...
    16791801  return( attrib(M, "isHomog") );
    16801802}
     1803example
     1804{ "EXAMPLE:"; echo = 2;
     1805
     1806  ring r=32003,(x,y,z),dp;
     1807
     1808  module Z = grobj(freemodule(0),intvec(0:0),intvec(0:0));
     1809 
     1810  grrange(Z);
     1811  grdeg(Z);
     1812 
     1813  grview(Z);
     1814
     1815  module P=grobj(module([xy,0,xz]),intvec(0,1,0));
     1816
     1817  grrange(P);
     1818  grdeg(P);
     1819
     1820  grview(P);
     1821}
    16811822
    16821823proc grlift0(M, N, alpha1)
    1683 "PURPOSE: generic random alpha0 : coker(M) -> coker(N) from random alpha1
     1824"USAGE: grlift0(M, N, alpha1) TODO!
     1825PURPOSE: generic random alpha0 : coker(M) -> coker(N) from random alpha1
    16841826NOTE: this proc can work only if some assumptions are fulfilled (due
    16851827to Wolfram)! e.g. at the end of a resolution for the source module...
     
    17331875
    17341876proc grlifting(M,N)
    1735 "USAGE: graded objects M and N
     1877"USAGE: grlifting(M,N), graded objects M and N
    17361878RETURN: map of chain complexes (as a list)
    17371879PURPOSE: construct a map of chain complexes between free resolutions of Img(M) and Img(N).
     
    18231965
    18241966proc mappingcone(M,N)
    1825 "USAGE:M,N graded objects
     1967"USAGE: mappingcone(M,N), M,N graded objects
    18261968RETURN: chain complex (as a list)
    18271969PURPOSE: construct a free resolution of the cokernel of a random map between Img(M), and Img(N).
     
    18852027// correct
    18862028proc grrndmap(def S, def D, list #)
    1887 "USAGE: (S,D), graded objects S and D
     2029"USAGE: grrndmap(S,D), graded objects S and D
    18882030RETURN: graded object
    18892031PURPOSE: construct a random 0-deg graded homomorphism src(S) -> src(D)
     
    19192061
    19202062proc grrndmap2(def D, def S, list #)
    1921 "USAGE: (D,S), graded objects S and D
     2063"USAGE: grrndmap2(D,S), graded objects S and D
    19222064RETURN: graded object
    19232065PURPOSE: construct a random 0-deg graded homomorphism between target of D and S.
     
    19552097//
    19562098proc grlifting2(A,B)
    1957 "USAGE: (A,B), graded objects A and B (matrices defining maps)
     2099"USAGE: grlifting2(A,B), graded objects A and B (matrices defining maps)
    19582100RETURN: map of chain complexes (as a list)
    19592101PURPOSE: construct a map of chain complexes between free resolution of
     
    20282170// G0<-----F0+G1<------F1+G2<-------F2+G3<-----
    20292171proc mappingcone2(A,B)
    2030 "USAGE: (A,B), graded objects A and B (matrices defining maps)
     2172"USAGE: mappingcone2(A,B), graded objects A and B (matrices defining maps)
    20312173RETURN: chain complex (as a list)
    20322174PURPOSE: construct the free resolution of a cokernel of a random map between M=coker(A), and N=coker(B)
     
    21982340//G0<-----F0+G1<------F1+G2<-------F2+G3<-----
    21992341proc mappingcone3(A,B)
    2200 "USAGE: (A,B), graded objects A and B (matrices defining maps)
     2342"USAGE: mappingcone3(A,B), graded objects A and B (matrices defining maps)
    22012343RETURN: chain complex (as a list)
    22022344PURPOSE: construct a free resolution of the cokernel of a random map between M=coker(A), and N=coker(B)
     
    22902432// TODO: Please decide between KeneshlouMatrixPresentation and matrixpres, and replace one with the other!
    22912433proc matrixpres(intvec a)
    2292 "USAGE:  intvec a
     2434"USAGE:  matrixpres(a), intvec a
    22932435RETURN:  graded object
    22942436PURPOSE: matrix presentation for direct sum of omega^a[i](i) in form of a graded object
  • Singular/LIB/grwalk.lib

    r22c3f2 r2c2d2e8  
    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/modwalk.lib

    • Property mode changed from 100644 to 100755
    r22c3f2 r2c2d2e8  
    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/rwalk.lib

    • Property mode changed from 100644 to 100755
    r22c3f2 r2c2d2e8  
    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/dyn_modules/syzextra/mod_main.cc

    r5638b9 r2c2d2e8  
    6767BEGIN_NAMESPACE_NONAME
    6868
     69// returns TRUE, if idRankFreeModule(m) > 0 ???
     70/// test whether this input has vectors among entries or no enties
     71/// result must be FALSE for only 0-entries
     72static BOOLEAN id_IsModule(ideal id, ring r)
     73{
     74  id_Test(id, r);
     75
     76  if( id->rank != 1 ) return TRUE;
     77
     78  if (rRing_has_Comp(r))
     79  {
     80    const int l = IDELEMS(id);
     81
     82    for (int j=0; j<l; j++)
     83      if (id->m[j] != NULL && p_GetComp(id->m[j], r) > 0)
     84        return TRUE;   
     85
     86    return FALSE; // rank: 1, only zero or no entries? can be an ideal OR module... BUT in the use-case should better be an ideal!
     87  }
     88
     89  return FALSE;
     90}
     91
     92
     93   
    6994
    7095static inline void NoReturn(leftv& res)
     
    15251550  const int iLimit = r->typ[pos].data.is.limit;
    15261551  const ideal F = r->typ[pos].data.is.F;
     1552 
    15271553  ideal FF = id_Copy(F, r);
    1528 
    1529 
    15301554
    15311555  lists l=(lists)omAllocBin(slists_bin);
     
    15381562  //        l->m[1].rtyp = MODUL_CMD;
    15391563
    1540   if( idIsModule(FF, r) )
     1564  if( id_IsModule(FF, r) ) // ???
    15411565  {
    15421566    l->m[1].rtyp = MODUL_CMD;
     
    18241848  {
    18251849    iComp = (int)((long)(h->Data()));
    1826   } else
     1850  }
     1851  else
    18271852  {
    18281853      if( (!isSyz) && (-1 == posIS) )
  • Singular/extra.cc

    r5638b9 r2c2d2e8  
    37813781        bigintmat *b=(bigintmat *)h->Data();
    37823782        // just for tests: simply transpose
    3783         bigintmat *bb=b->transpose();
    3784         // return the result:
     3783        bigintmat *bb=b->transpose();
     3784        // return the result:
    37853785        res->rtyp=CMATRIX_CMD;
    3786         res->data=(char*)bb;
     3786        res->data=(char*)bb;
    37873787        return FALSE;
    37883788      }
     
    37903790      {
    37913791        WerrorS("system(\"LU\",<cmatrix>) expected");
    3792         return TRUE;
     3792        return TRUE;
    37933793      }
    37943794    }
  • Singular/iparith.cc

    r5638b9 r2c2d2e8  
    10941094  if (res->data==NULL)
    10951095  {
    1096      Werror("matrix size not compatible(%dx%d, %dx%d)",
     1096     Werror("matrix size not compatible(%dx%d, %dx%d) in *",
    10971097             MATROWS(A),MATCOLS(A),MATROWS(B),MATCOLS(B));
    10981098     return TRUE;
     
    14821482  poly r=p; // pointer to the beginning of component i
    14831483  poly o=NULL;
    1484   unsigned i=(unsigned)(long)v->Data();
     1484  int i=(int)(long)v->Data();
    14851485  while (p!=NULL)
    14861486  {
     
    53305330  return TRUE;
    53315331}
     5332static int WerrorS_dummy_cnt=0;
     5333static void WerrorS_dummy(const char *)
     5334{
     5335  WerrorS_dummy_cnt++;
     5336}
    53325337BOOLEAN jjLOAD_TRY(const char *s)
    53335338{
    5334   WerrorS("not yet");
    5335   return TRUE;
     5339  void (*WerrorS_save)(const char *s) = WerrorS_callback;
     5340  WerrorS_callback=WerrorS_dummy;
     5341  WerrorS_dummy_cnt=0;
     5342  BOOLEAN bo=jjLOAD(s,TRUE);
     5343  if (bo || (WerrorS_dummy_cnt>0)) Print("loading of >%s< failed\n",s);
     5344  WerrorS_callback=WerrorS_save;
     5345  errorreported=0;
     5346  return FALSE;
    53365347}
    53375348
     
    81148125{
    81158126  memset(res,0,sizeof(sleftv));
    8116   BOOLEAN call_failed=FALSE;
    81178127
    81188128  if (!errorreported)
     
    83028312{
    83038313  memset(res,0,sizeof(sleftv));
    8304   BOOLEAN call_failed=FALSE;
    83058314
    83068315  if (!errorreported)
     
    83468355    }
    83478356
    8348     BOOLEAN failed=FALSE;
    83498357    iiOp=op;
    83508358    int i=iiTabIndex(dArithTab1,JJTAB1LEN,op);
  • Singular/ipshell.cc

    r5638b9 r2c2d2e8  
    781781        {
    782782          v->rtyp=IDEAL_CMD;
     783          char *tmp = theMap->preimage;
     784          theMap->preimage=(char*)1L;
     785          // map gets 1 as its rank (as an ideal)
    783786          v->data=fast_map(IDIDEAL(w), src_ring, (ideal)theMap, currRing);
     787          theMap->preimage=tmp; // map gets its preimage back
    784788        }
    785789        else
  • Singular/maps_ip.cc

    r5638b9 r2c2d2e8  
    4343#include <polys/monomials/maps.h>
    4444#endif
    45 
    4645
    4746/*2
     
    391390    ideal src_id=idInit(1,1);
    392391    src_id->m[0]=p;
     392
     393    char *tmp = theMap->preimage;
     394    theMap->preimagei=(char*)1L; // map gets 1 as its rank (as an ideal)
    393395    ideal res_id=fast_map(src_id,currRing,(ideal)theMap,currRing);
     396    theMap->preimage=tmp; // map gets its preimage back
     397
    394398    res=res_id->m[0];
    395399    res_id->m[0]=NULL; idDelete(&res_id);
  • Singular/table.h

    r5638b9 r2c2d2e8  
    6161,{D(jjCOLS_IV),    COLS_CMD,        INT_CMD,        INTMAT_CMD    , ALLOW_PLURAL |ALLOW_RING}
    6262,{D(jjCOLS_BIM),   COLS_CMD,        INT_CMD,        BIGINTMAT_CMD , ALLOW_PLURAL |ALLOW_RING}
     63#ifdef SINGULAR_4_1
     64,{D(jjCOLS_BIM),   COLS_CMD,        INT_CMD,        CMATRIX_CMD , ALLOW_PLURAL |ALLOW_RING}
     65#endif
    6366,{  jjWRONG ,      COLS_CMD,        0,              INTVEC_CMD    , ALLOW_PLURAL |ALLOW_RING}
    6467,{D(jjCONTENT),    CONTENT_CMD,     POLY_CMD,       POLY_CMD      , ALLOW_PLURAL |ALLOW_RING}
     
    219222,{D(jjROWS_IV),    ROWS_CMD,        INT_CMD,        INTMAT_CMD    , ALLOW_PLURAL |ALLOW_RING}
    220223,{D(jjROWS_BIM),   ROWS_CMD,        INT_CMD,        BIGINTMAT_CMD , ALLOW_PLURAL |ALLOW_RING}
     224#ifdef SINGULAR_4_1
     225,{D(jjROWS_BIM),   ROWS_CMD,        INT_CMD,        CMATRIX_CMD , ALLOW_PLURAL |ALLOW_RING}
     226#endif
    221227,{D(jjCOUNT_IV),   ROWS_CMD,        INT_CMD,        INTVEC_CMD    , ALLOW_PLURAL |ALLOW_RING}
    222228,{D(jjSBA),        SBA_CMD,         IDEAL_CMD,      IDEAL_CMD     , ALLOW_PLURAL |ALLOW_RING}
  • Singular/walk.cc

    • Property mode changed from 100644 to 100755
    r22c3f2 r2c2d2e8  
    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
     
    953985}
    954986
    955 /*****************************************************************************
    956 * create a weight matrix order as intvec of an extra weight vector (a(iv),lp)*
    957 ******************************************************************************/
     987/*********************************************************************************
     988* create a weight matrix order as intvec of an extra weight vector (a(iv),M(iw)) *
     989**********************************************************************************/
    958990intvec* MivMatrixOrderRefine(intvec* iv, intvec* iw)
    959991{
    960   assume(iv->length() == iw->length());
    961   int i, nR = iv->length();
    962 
     992  assume((iv->length())*(iv->length()) == iw->length());
     993  int i,j, nR = iv->length();
     994 
    963995  intvec* ivm = new intvec(nR*nR);
    964996
     
    966998  {
    967999    (*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;
     1000  }
     1001  for(i=1; i<nR; i++)
     1002  {
     1003    for(j=0; j<nR; j++)
     1004    {
     1005      (*ivm)[j+i*nR] = (*iw)[j+i*nR];
     1006    }
    9731007  }
    9741008  return ivm;
     
    18591893}
    18601894
     1895
     1896/**************************************************************
     1897 * Look for the position of the smallest absolut value in vec *
     1898 **************************************************************/
     1899static int MivAbsMaxArg(intvec* vec)
     1900{
     1901  int k = MivAbsMax(vec);
     1902  int i=0;
     1903  while(1)
     1904  {
     1905    if((*vec)[i] == k || (*vec)[i] == -k)
     1906    {
     1907      break;
     1908    }
     1909    i++;
     1910  }
     1911  return i;
     1912}
     1913
     1914
    18611915/**********************************************************************
    18621916 * Compute a next weight vector between curr_weight and target_weight *
     
    18731927
    18741928  int nRing = currRing->N;
    1875   int checkRed, j, kkk, nG = IDELEMS(G);
     1929  int checkRed, j, nG = IDELEMS(G);
    18761930  intvec* ivtemp;
    18771931
     
    19111965  mpz_init(dcw);
    19121966
    1913   //int tn0, tn1, tz1, ncmp, gcd_tmp, ntmp;
    19141967  int gcd_tmp;
    19151968  intvec* diff_weight = MivSub(target_weight, curr_weight);
     
    19171970  intvec* diff_weight1 = MivSub(target_weight, curr_weight);
    19181971  poly g;
    1919   //poly g, gw;
     1972
    19201973  for (j=0; j<nG; j++)
    19211974  {
     
    19341987        mpz_set_si(MwWd, MivDotProduct(ivtemp, curr_weight));
    19351988        mpz_sub(s_zaehler, deg_w0_p1, MwWd);
    1936 
    19371989        if(mpz_cmp(s_zaehler, t_null) != 0)
    19381990        {
    19391991          mpz_set_si(MwWd, MivDotProduct(ivtemp, diff_weight));
    19401992          mpz_sub(s_nenner, MwWd, deg_d0_p1);
    1941 
    19421993          // check for 0 < s <= 1
    19431994          if( (mpz_cmp(s_zaehler,t_null) > 0 &&
     
    19792030    }
    19802031  }
    1981 //Print("\n// Alloc Size = %d \n", nRing*sizeof(mpz_t));
     2032  //Print("\n// Alloc Size = %d \n", nRing*sizeof(mpz_t));
    19822033  mpz_t *vec=(mpz_t*)omAlloc(nRing*sizeof(mpz_t));
    19832034
    19842035
    1985   // there is no 0<t<1 and define the next weight vector that is equal to the current weight vector
     2036  // there is no 0<t<1 and define the next weight vector that is equal
     2037  // to the current weight vector
    19862038  if(mpz_cmp(t_nenner, t_null) == 0)
    19872039  {
     
    20542106#endif
    20552107
    2056   //  BOOLEAN isdwpos;
    2057 
    2058   // construct a new weight vector
     2108// construct a new weight vector and check whether vec[j] is overflow,
     2109// i.e. vec[j] > 2^31.
     2110// If vec[j] doesn't overflow, define a weight vector. Otherwise,
     2111// report that overflow appears. In the second case, test whether the
     2112// the correctness of the new vector plays an important role
     2113
    20592114  for (j=0; j<nRing; j++)
    20602115  {
     
    21002155    }
    21012156  }
    2102 
     2157  // reduce the vector with the gcd
     2158  if(mpz_cmp_si(ggt,1) != 0)
     2159  {
     2160    for (j=0; j<nRing; j++)
     2161    {
     2162      mpz_divexact(vec[j], vec[j], ggt);
     2163    }
     2164  }
    21032165#ifdef  NEXT_VECTORS_CC
    21042166  PrintS("\n// gcd of elements of the vector: ");
     
    21062168#endif
    21072169
    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;
    21162170  for(j=0; j<nRing; j++)
    21172171    {
     
    21292183
    21302184  REDUCTION:
     2185  checkRed = 1;
    21312186  for (j=0; j<nRing; j++)
    21322187  {
    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     }
     2188    (*diff_weight1)[j] = mpz_get_si(vec[j]);
     2189  }
     2190  while(test_w_in_ConeCC(G,diff_weight1))
     2191  {
     2192    for(j=0; j<nRing; j++)
     2193    {
     2194      (*diff_weight)[j] = (*diff_weight1)[j];
     2195      mpz_set_si(vec[j], (*diff_weight)[j]);     
     2196    }
     2197    for(j=0; j<nRing; j++)
     2198    {
     2199      (*diff_weight1)[j] = floor(0.1*(*diff_weight)[j] + 0.5);
     2200    }
     2201  }
     2202  if(MivAbsMax(diff_weight)>10000)
     2203  {
     2204    for(j=0; j<nRing; j++)
     2205    {
     2206      (*diff_weight1)[j] = (*diff_weight)[j];
     2207    }
     2208    j = 0;
     2209    while(test_w_in_ConeCC(G,diff_weight1))
     2210    {
     2211      (*diff_weight)[j] = (*diff_weight1)[j];
     2212      mpz_set_si(vec[j], (*diff_weight)[j]);
     2213      j = MivAbsMaxArg(diff_weight1);
     2214      (*diff_weight1)[j] = floor(0.1*(*diff_weight1)[j] + 0.5);
     2215    }
     2216    goto SIMPLIFY_GCD;
    21862217  }
    21872218
     
    22222253   mpz_clear(t_null);
    22232254
    2224 
    2225 
    22262255  if(Overflow_Error == FALSE)
    22272256  {
    22282257    Overflow_Error = nError;
    22292258  }
    2230  rComplete(currRing);
    2231   for(kkk=0; kkk<IDELEMS(G);kkk++)
    2232   {
    2233     poly p=G->m[kkk];
     2259  rComplete(currRing);
     2260  for(j=0; j<IDELEMS(G); j++)
     2261  {
     2262    poly p=G->m[j];
    22342263    while(p!=NULL)
    22352264    {
     
    22712300}
    22722301
    2273 /**************************************************************
     2302/********************************************************************
    22742303 * define and execute a new ring which order is (a(vb),a(va),lp,C)  *
    2275  * ************************************************************/
    2276 #if 0
    2277 // unused
     2304 * ******************************************************************/
    22782305static void VMrHomogeneous(intvec* va, intvec* vb)
    22792306{
     
    23532380  rChangeCurrRing(r);
    23542381}
    2355 #endif
     2382
    23562383
    23572384/**************************************************************
     
    24282455}
    24292456
    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 
    25012457/****************************************************************
    25022458 * define and execute a new ring with ordering (a(va),Wp(vb),C) *
    25032459 * **************************************************************/
    2504 
    25052460static ring VMrRefine(intvec* va, intvec* vb)
    25062461{
     
    25762531
    25772532  // complete ring intializations
    2578 
     2533 
    25792534  rComplete(r);
    25802535
     
    28062761}
    28072762
    2808 
    2809 /* define a ring with parameters und change to it */
    2810 /* DefRingPar and DefRingParlp corrupt still memory */
     2763/***************************************************
     2764* define a ring with parameters und change to it   *
     2765* DefRingPar and DefRingParlp corrupt still memory *
     2766****************************************************/
    28112767static void DefRingPar(intvec* va)
    28122768{
     
    29562912}
    29572913
    2958 //unused
    2959 /**************************************************************
    2960  * check wheather one or more components of a vector are zero *
    2961  **************************************************************/
    2962 #if 0
     2914/*************************************************************
     2915 * check whether one or more components of a vector are zero *
     2916 *************************************************************/
    29632917static int isNolVector(intvec* hilb)
    29642918{
     
    29732927  return 0;
    29742928}
    2975 #endif
     2929
     2930/*************************************************************
     2931 * check whether one or more components of a vector are <= 0 *
     2932 *************************************************************/
     2933static int isNegNolVector(intvec* hilb)
     2934{
     2935  int i;
     2936  for(i=hilb->length()-1; i>=0; i--)
     2937  {
     2938    if((* hilb)[i]<=0)
     2939    {
     2940      return 1;
     2941    }
     2942  }
     2943  return 0;
     2944}
     2945
     2946/**************************************************************************
     2947* Gomega is the initial ideal of G w. r. t. the current weight vector     *
     2948* curr_weight. Check whether curr_weight lies on a border of the Groebner *
     2949* cone, i. e. check whether a monomial is divisible by a leading monomial *
     2950***************************************************************************/
     2951static ideal middleOfCone(ideal G, ideal Gomega)
     2952{
     2953  BOOLEAN middle = FALSE;
     2954  int i,j,N = IDELEMS(Gomega);
     2955  poly p,lm,factor1,factor2;
     2956  //PrintS("\n//** idCopy\n");
     2957  ideal Go = idCopy(G);
     2958 
     2959  //PrintS("\n//** jetzt for-Loop!\n");
     2960
     2961  // check whether leading monomials of G and Gomega coincide
     2962  // and return NULL if not
     2963  for(i=0; i<N; i++)
     2964  {
     2965    p = pCopy(Gomega->m[i]);
     2966    lm = pCopy(pHead(G->m[i]));
     2967    if(!pIsConstant(pSub(p,lm)))
     2968    {
     2969      //pDelete(&p);
     2970      //pDelete(&lm);
     2971      idDelete(&Go);
     2972      return NULL;
     2973    }
     2974    //pDelete(&p);
     2975    //pDelete(&lm);
     2976  }
     2977  for(i=0; i<N; i++)
     2978  {
     2979    for(j=0; j<N; j++)
     2980    {
     2981      if(i!=j)
     2982      {
     2983        p = pCopy(Gomega->m[i]);
     2984        lm = pCopy(Gomega->m[j]);
     2985        pIter(p);
     2986        while(p!=NULL)
     2987        {
     2988          if(pDivisibleBy(lm,p))
     2989          {
     2990            if(middle == FALSE)
     2991            {
     2992              middle = TRUE;
     2993            }
     2994            factor1 = singclap_pdivide(pHead(p),lm,currRing);
     2995            factor2 = pMult(pCopy(factor1),pCopy(Go->m[j]));
     2996            pDelete(&factor1);
     2997            Go->m[i] = pAdd((Go->m[i]),pNeg(pCopy(factor2)));
     2998            pDelete(&factor2);
     2999          }
     3000          pIter(p);
     3001        }
     3002        pDelete(&lm);
     3003        pDelete(&p);
     3004      }
     3005    }
     3006  }
     3007 
     3008  //PrintS("\n//** jetzt Delete!\n");
     3009  //pDelete(&p);
     3010  //pDelete(&factor);
     3011  //pDelete(&lm);
     3012  if(middle == TRUE)
     3013  {
     3014    //PrintS("\n//** middle TRUE!\n");
     3015    return Go;
     3016  }
     3017  //PrintS("\n//** middle FALSE!\n");
     3018  idDelete(&Go);
     3019  return NULL;
     3020}
    29763021
    29773022/******************************  Februar 2002  ****************************
     
    31283173    else
    31293174    {
    3130       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung
     3175      rChangeCurrRing(VMrDefault(curr_weight));
    31313176    }
    31323177    newRing = currRing;
     
    32803325  }
    32813326  return 0;
     3327}
     3328
     3329/*****************************************
     3330 * return maximal polynomial length of G *
     3331 *****************************************/
     3332static int maxlengthpoly(ideal G)
     3333{
     3334  int i,k,length=0;
     3335  for(i=IDELEMS(G)-1; i>=0; i--)
     3336  {
     3337    k = pLength(G->m[i]);
     3338    if(k>length)
     3339    {
     3340      length = k;
     3341    }
     3342  }
     3343  return length;
    32823344}
    32833345
     
    38843946    else
    38853947    {
    3886       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung
     3948      rChangeCurrRing(VMrDefault(curr_weight));
    38873949    }
    38883950    newRing = currRing;
     
    41454207    else
    41464208    {
    4147       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung
     4209      rChangeCurrRing(VMrDefault(curr_weight));
    41484210    }
    41494211    newRing = currRing;
     
    42874349intvec* Xivlp;
    42884350
    4289 #if 0
     4351
    42904352/********************************
    42914353 * compute a next weight vector *
    42924354 ********************************/
    4293 static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight, intvec* target_weight, int weight_rad, int pert_deg)
     4355static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight,
     4356               intvec* target_weight, int weight_rad, int pert_deg)
    42944357{
    4295   int i, weight_norm;
    4296   int nV = currRing->N;
     4358  assume(currRing != NULL && curr_weight != NULL &&
     4359         target_weight != NULL && G->m[0] != NULL);
     4360
     4361  int i,weight_norm,nV = currRing->N;
    42974362  intvec* next_weight2;
    42984363  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)
     4364  intvec* result = new intvec(nV);
     4365
     4366  intvec* next_weight1 =MkInterRedNextWeight(curr_weight,target_weight,G);
     4367  //compute a random next weight vector "next_weight2"
     4368  while(1)
     4369  {
     4370    weight_norm = 0;
     4371    while(weight_norm == 0)
     4372    {
     4373
     4374      for(i=0; i<nV; i++)
     4375      {
     4376        (*next_weight22)[i] = rand() % 60000 - 30000;
     4377        weight_norm = weight_norm + (*next_weight22)[i]*(*next_weight22)[i];
     4378      }
     4379      weight_norm = 1 + floor(sqrt(weight_norm));
     4380    }
     4381
     4382    for(i=0; i<nV; i++)
     4383    {
     4384      if((*next_weight22)[i] < 0)
     4385      {
     4386        (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);
     4387      }
     4388      else
     4389      {
     4390        (*next_weight22)[i] = (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);
     4391      }
     4392    }
     4393
     4394    if(test_w_in_ConeCC(G, next_weight22) == 1)
     4395    {
     4396      next_weight2 = MkInterRedNextWeight(next_weight22,target_weight,G);
     4397      if(MivAbsMax(next_weight2)>1147483647)
    43154398      {
    43164399        for(i=0; i<nV; i++)
    43174400        {
    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];
     4401          (*next_weight22)[i] = (*next_weight2)[i];
    43214402        }
    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)
     4403        i = 0;
     4404        /* reduce the size of the maximal entry of the vector*/
     4405        while(test_w_in_ConeCC(G,next_weight22))
    43284406        {
    4329           (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);
     4407          (*next_weight2)[i] = (*next_weight22)[i];
     4408          i = MivAbsMaxArg(next_weight22);
     4409          (*next_weight22)[i] = floor(0.1*(*next_weight22)[i] + 0.5);
    43304410        }
    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);
     4411      }
     4412      delete next_weight22;
     4413      break;
     4414    }
     4415  }
     4416 
     4417  // compute "usual" next weight vector
     4418  intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight, G);
     4419  ideal G_test = MwalkInitialForm(G, next_weight);
     4420  ideal G_test2 = MwalkInitialForm(G, next_weight2);
     4421
     4422  // compare next weights
     4423  if(Overflow_Error == FALSE)
     4424  {
    43484425    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|
     4426    if(G_test1->m[0] != NULL && maxlengthpoly(G_test1) < maxlengthpoly(G_test))//if(IDELEMS(G_test1) < IDELEMS(G_test))
     4427    {
     4428      if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test1))//if(IDELEMS(G_test2) < IDELEMS(G_test1))
    43554429      {
    43564430        for(i=0; i<nV; i++)
     
    43594433        }
    43604434      }
    4361       else // |G_test1| < |G_test|, |G_test1| < |G_test2|
     4435      else
    43624436      {
    43634437        for(i=0; i<nV; i++)
     
    43654439          (*result)[i] = (*next_weight1)[i];
    43664440        }
    4367       }
     4441      }   
    43684442    }
    43694443    else
    43704444    {
    4371       if(IDELEMS(G_test2) <= IDELEMS(G_test)) // |G_test2| <= |G_test| <= |G_test1|
     4445      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|
    43724446      {
    43734447        for(i=0; i<nV; i++)
     
    43764450        }
    43774451      }
    4378       else // |G_test| <= |G_test1|, |G_test| < |G_test2|
     4452      else
    43794453      {
    43804454        for(i=0; i<nV; i++)
     
    43844458      }
    43854459    }
    4386     delete next_weight;
    4387     delete next_weight1;
    4388     idDelete(&G_test);
    43894460    idDelete(&G_test1);
    4390     idDelete(&G_test2);
    4391     if(test_w_in_ConeCC(G, result) == 1)
    4392     {
    4393       delete next_weight2;
    4394       return result;
     4461  }
     4462  else
     4463  {
     4464    Overflow_Error = FALSE;
     4465    if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test))//if(IDELEMS(G_test2) < IDELEMS(G_test))
     4466    {
     4467      for(i=1; i<nV; i++)
     4468      {
     4469        (*result)[i] = (*next_weight2)[i];
     4470      }
    43954471    }
    43964472    else
    43974473    {
    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     {
    44264474      for(i=0; i<nV; i++)
    44274475      {
    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       {
    45134476        (*result)[i] = (*next_weight)[i];
    45144477      }
    45154478    }
    45164479  }
     4480  //PrintS("\n MWalkRandomNextWeight: Ende ok!\n");
    45174481  idDelete(&G_test);
    45184482  idDelete(&G_test2);
     
    45754539    else
    45764540    {
    4577       rChangeCurrRing(VMrDefault(orig_target_weight)); // Aenderung
     4541      rChangeCurrRing(VMrDefault(orig_target_weight));
    45784542    }
    45794543    TargetRing = currRing;
     
    46464610    else
    46474611    {
    4648       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung
     4612      rChangeCurrRing(VMrDefault(curr_weight));
    46494613    }
    46504614    newRing = currRing;
     
    47554719    else
    47564720    {
    4757       rChangeCurrRing(VMrDefault(orig_target_weight)); // Aenderung
     4721      rChangeCurrRing(VMrDefault(orig_target_weight));
    47584722    }
    47594723    F1 = idrMoveR(G, newRing,currRing);
     
    47864750      else
    47874751      {
    4788         rChangeCurrRing(VMrDefault(orig_target_weight)); // Aenderung
     4752        rChangeCurrRing(VMrDefault(orig_target_weight));
    47894753      }
    47904754    KSTD_Finish:
     
    48844848      tim = clock();
    48854849      /*
    4886         Print("\n// **** Grï¿œbnerwalk took %d steps and ", nwalk);
     4850        Print("\n// **** Groebnerwalk took %d steps and ", nwalk);
    48874851        PrintS("\n// **** call the rec. Pert. Walk to compute a red GB of:");
    48884852        idElements(Gomega, "G_omega");
     
    49144878      oldRing = currRing;
    49154879
    4916       /* create a new ring newRing */
     4880      // create a new ring newRing
    49174881       if (rParameter(currRing) != NULL)
    49184882       {
     
    49214885       else
    49224886       {
    4923          rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung
     4887         rChangeCurrRing(VMrDefault(curr_weight));
    49244888       }
    49254889      newRing = currRing;
     
    49474911      else
    49484912      {
    4949         rChangeCurrRing(VMrDefault(curr_weight));  //Aenderung
     4913        rChangeCurrRing(VMrDefault(curr_weight));
    49504914      }
    49514915      newRing = currRing;
     
    49594923      M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
    49604924      delete hilb_func;
    4961 #endif // BUCHBERGER_ALG
     4925#endif
    49624926      tstd = tstd + clock() - to;
    49634927
     
    49684932
    49694933      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.
     4934      // compute a representation of the generators of submod (M) with respect
     4935      // to those of mod (Gomega).
     4936      // Gomega is a reduced Groebner basis w.r.t. the current ring.
    49714937      F = MLifttwoIdeal(Gomega2, M1, G);
    49724938      tlift = tlift + clock() - to;
     
    50184984      else
    50194985      {
    5020         rChangeCurrRing(VMrDefault(target_weight)); // Aenderung
     4986        rChangeCurrRing(VMrDefault(target_weight));
    50214987      }
    50224988      F1 = idrMoveR(G, newRing,currRing);
     
    50655031 * THE GROEBNER WALK ALGORITHM *
    50665032 *******************************/
    5067 ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M, ring baseRing)
     5033ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M,
     5034            ring baseRing, int reduction, int printout)
    50685035{
    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));
     5036  // save current options
     5037  BITSET save1 = si_opt_1;
     5038  if(reduction == 0)
     5039  {
     5040    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     5041    si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     5042  }
    50735043  Set_Error(FALSE);
    50745044  Overflow_Error = FALSE;
     5045  BOOLEAN endwalks = FALSE;
    50755046#ifdef TIME_TEST
    50765047  clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
     
    50805051#endif
    50815052  nstep=0;
    5082   int i,nwalk,endwalks = 0;
     5053  int i,nwalk;
    50835054  int nV = baseRing->N;
    50845055
    5085   ideal Gomega, M, F, Gomega1, Gomega2, M1; //, F1;
     5056  ideal Gomega, M, F, FF, Gomega1, Gomega2, M1;
    50865057  ring newRing;
    50875058  ring XXRing = baseRing;
     
    50905061  intvec* target_weight = new intvec(nV);
    50915062  intvec* exivlp = Mivlp(nV);
     5063/*
    50925064  intvec* tmp_weight = new intvec(nV);
    50935065  for(i=0; i<nV; i++)
     
    50955067    (*tmp_weight)[i] = (*target_M)[i];
    50965068  }
     5069*/
    50975070  for(i=0; i<nV; i++)
    50985071  {
     
    51115084#endif
    51125085  rComplete(currRing);
    5113 #ifdef CHECK_IDEAL_MWALK
    5114     idString(Go,"Go");
    5115 #endif
     5086//#ifdef CHECK_IDEAL_MWALK
     5087  if(printout > 2)
     5088  {
     5089    idString(Go,"//** Mwalk: Go");
     5090  }
     5091//#endif
    51165092#ifdef TIME_TEST
    51175093  to = clock();
    51185094#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       }
     5095  if(orig_M->length() == nV)
     5096  {
     5097    // define a new ring with ordering "(a(curr_weight),lp)
     5098    newRing = VMrDefault(curr_weight);
     5099  }
     5100  else
     5101  {
     5102    newRing = VMatrDefault(orig_M);
     5103  }
    51275104  rChangeCurrRing(newRing);
    51285105  ideal G = MstdCC(idrMoveR(Go,baseRing,currRing));
     5106  //ideal G = idrMoveR(Go,baseRing,currRing);
    51295107  baseRing = currRing;
    51305108#ifdef TIME_TEST
     
    51405118    to = clock();
    51415119#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"
     5120//#ifdef CHECK_IDEAL_MWALK
     5121    if(printout > 2)
     5122    {
     5123      idString(G,"//** Mwalk: G");
     5124    }
     5125//#endif
     5126    // test whether target cone is reached
     5127    if(reduction !=0 && test_w_in_ConeCC(G,target_weight) == 1)
     5128    {
     5129      endwalks = TRUE;
     5130    }
     5131    // compute an initial form ideal of <G> w.r.t. "curr_vector"
     5132    Gomega = MwalkInitialForm(G, curr_weight);
    51465133#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
     5134    //time for computing initial form ideal
     5135    tif = tif + clock()-to;
     5136#endif
     5137//#ifdef CHECK_IDEAL_MWALK
     5138    if(printout > 1)
     5139    {
     5140      idString(Gomega,"//** Mwalk: Gomega");
     5141    }
     5142//#endif
     5143    if(reduction == 0)
     5144    {
     5145      //PrintS("\n//** Mwalk: test middle of cone!\n");
     5146      FF = middleOfCone(G,Gomega);
     5147      //PrintS("\n//** Mwalk: Test F!\n");
     5148      if( FF != NULL)
     5149      {
     5150        idDelete(&G);
     5151        G = idCopy(FF);
     5152        idDelete(&FF);
     5153        //PrintS("\n//** Mwalk: FF nicht NULL! Compue next vector.\n");
     5154        goto NEXT_VECTOR;
     5155      }
     5156    }
    51525157#ifndef  BUCHBERGER_ALG
    51535158    if(isNolVector(curr_weight) == 0)
    51545159    {
    51555160      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
    5156     }
     5161    }   
    51575162    else
    51585163    {
     
    51645169      if(orig_M->length() == nV)
    51655170      {
    5166         newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
     5171        // define a new ring with ordering "(a(curr_weight),lp)
     5172        newRing = VMrDefault(curr_weight);
    51675173      }
    51685174      else
     
    51755181     if(target_M->length() == nV)
    51765182     {
    5177        newRing = VMrRefine(curr_weight,target_weight); //define a new ring with ordering "(a(curr_weight),Wp(target_weight))"
     5183       //define a new ring with ordering "(a(curr_weight),Wp(target_weight))"
     5184       newRing = VMrRefine(curr_weight,target_weight);
    51785185     }
    51795186     else
    51805187     {
     5188       //define a new ring with matrix ordering
    51815189       newRing = VMatrRefine(target_M,curr_weight);
    51825190     }
     
    51995207#endif
    52005208    idSkipZeroes(M);
    5201 #ifdef CHECK_IDEAL_MWALK
    5202     PrintS("\n//** Mwalk: computed M.\n");
    5203     idString(M, "M");
    5204 #endif
     5209//#ifdef CHECK_IDEAL_MWALK
     5210    if(printout > 2)
     5211    {
     5212      idString(M, "//** Mwalk: M");
     5213    }
     5214//#endif
    52055215    //change the ring to baseRing
    52065216    rChangeCurrRing(baseRing);
     
    52125222    to = clock();
    52135223#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
     5224    // compute a representation of the generators of submod (M) with respect to those of mod (Gomega),
     5225    // where Gomega is a reduced Groebner basis w.r.t. the current ring
    52155226    F = MLifttwoIdeal(Gomega2, M1, G);
    52165227#ifdef TIME_TEST
    52175228    tlift = tlift + clock() - to;
    52185229#endif
    5219 #ifdef CHECK_IDEAL_MWALK
    5220     idString(F, "F");
    5221 #endif
     5230//#ifdef CHECK_IDEAL_MWALK
     5231    if(printout > 2)
     5232    {
     5233      idString(F, "//** Mwalk: F");
     5234    }
     5235//#endif
    52225236    idDelete(&Gomega2);
    52235237    idDelete(&M1);
     
    52295243    to = clock();
    52305244#endif
    5231     //G = kStd(F1,NULL,testHomog,NULL,NULL,0,0,NULL);
     5245
    52325246#ifdef TIME_TEST
    52335247    tstd = tstd + clock() - to;
    52345248#endif
    52355249    idSkipZeroes(G);
    5236 #ifdef CHECK_IDEAL_MWALK
    5237     idString(G, "G");
    5238 #endif
     5250//#ifdef CHECK_IDEAL_MWALK
     5251    if(printout > 2)
     5252    {
     5253      idString(G, "//** Mwalk: G");
     5254    }
     5255//#endif
    52395256#ifdef TIME_TEST
    52405257    to = clock();
    52415258#endif
     5259    NEXT_VECTOR:
    52425260    intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
    52435261#ifdef TIME_TEST
    52445262    tnw = tnw + clock() - to;
    52455263#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
     5264//#ifdef PRINT_VECTORS
     5265    if(printout > 0)
     5266    {
     5267      MivString(curr_weight, target_weight, next_weight);
     5268    }
     5269//#endif
     5270    if(MivComp(target_weight,curr_weight) == 1 || endwalks == TRUE)
     5271    {
     5272
     5273//#ifdef CHECK_IDEAL_MWALK
     5274      if(printout > 0)
     5275      {
     5276        PrintS("\n//** Mwalk: entering last cone.\n");
     5277      }
     5278//#endif
     5279
    52545280      Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
    52555281      if(target_M->length() == nV)
     
    52645290      Gomega1 = idrMoveR(Gomega, baseRing,currRing);
    52655291      idDelete(&Gomega);
    5266 #ifdef CHECK_IDEAL_MWALK
    5267       idString(Gomega1, "Gomega");
    5268 #endif
     5292//#ifdef CHECK_IDEAL_MWALK
     5293      if(printout > 1)
     5294      {
     5295        idString(Gomega1, "//** Mwalk: Gomega");
     5296      }
     5297      PrintS("\n //** Mwalk: kStd(Gomega)");
     5298//#endif
    52695299      M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL);
    5270 #ifdef CHECK_IDEAL_MWALK
    5271       idString(M,"M");
    5272 #endif
     5300//#ifdef CHECK_IDEAL_MWALK
     5301      if(printout > 1)
     5302      {
     5303        idString(M,"//** Mwalk: M");
     5304      }
     5305//#endif
    52735306      rChangeCurrRing(baseRing);
    52745307      M1 =  idrMoveR(M, newRing,currRing);
     
    52765309      Gomega2 = idrMoveR(Gomega1, newRing,currRing);
    52775310      idDelete(&Gomega1);
     5311      //PrintS("\n //** Mwalk: MLifttwoIdeal");
    52785312      F = MLifttwoIdeal(Gomega2, M1, G);
    5279 #ifdef CHECK_IDEAL_MWALK
    5280       idString(F,"F");
    5281 #endif
     5313//#ifdef CHECK_IDEAL_MWALK
     5314      if(printout > 2)
     5315      {
     5316        idString(F,"//** Mwalk: F");
     5317      }
     5318//#endif
    52825319      idDelete(&Gomega2);
    52835320      idDelete(&M1);
     
    52915328      to = clock();
    52925329#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   //    }
     5330      /*PrintS("\n //**Mwalk: Interreduce");
     5331      //interreduce the Groebner basis <G> w.r.t. currRing
     5332      G = kInterRedCC(G,NULL);*/
    52975333#ifdef TIME_TEST
    52985334      tred = tred + clock() - to;
     
    53015337      delete next_weight;
    53025338      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
     5339    }
     5340
    53105341    for(i=nV-1; i>=0; i--)
    53115342    {
    5312       (*tmp_weight)[i] = (*curr_weight)[i];
     5343      //(*tmp_weight)[i] = (*curr_weight)[i];
    53135344      (*curr_weight)[i] = (*next_weight)[i];
    53145345    }
     
    53185349  ideal result = idrMoveR(G,baseRing,currRing);
    53195350  idDelete(&G);
    5320 /*#ifdef CHECK_IDEAL_MWALK
    5321   pDelete(&p);
    5322 #endif*/
    5323   delete tmp_weight;
     5351  //delete tmp_weight;
    53245352  delete ivNull;
    53255353  delete exivlp;
     
    53285356#endif
    53295357#ifdef TIME_TEST
    5330   Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep);
    53315358  TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep);
    5332   Print("\n//** Mwalk: Ergebnis.\n");
    53335359  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
    53345360  //Print("\n// Overflow_Error? (%d)\n", Overflow_Error);
    53355361#endif
     5362  Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep);
    53365363  return(result);
    53375364}
    53385365
    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)
     5366// THE RANDOM WALK ALGORITHM
     5367ideal Mrwalk(ideal Go, intvec* orig_M, intvec* target_M, int weight_rad, int pert_deg,
     5368             int reduction, int printout)
    53425369{
    53435370  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));
     5371  if(reduction == 0)
     5372  {
     5373    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     5374    si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     5375  }
    53475376  Set_Error(FALSE);
    53485377  Overflow_Error = FALSE;
     5378  BOOLEAN endwalks = FALSE;
    53495379#ifdef TIME_TEST
    53505380  clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
     
    53545384#endif
    53555385  nstep=0;
    5356   int i,nwalk,endwalks = 0;
    5357   int nV = baseRing->N;
    5358 
    5359   ideal Gomega, M, F, Gomega1, Gomega2, M1; //, F1;
     5386  int i,polylength,nwalk;
     5387  int nV = currRing->N;
     5388
     5389  ideal Gomega, M, F,FF, Gomega1, Gomega2, M1;
    53605390  ring newRing;
    5361   ring XXRing = baseRing;
     5391  ring baseRing = currRing;
     5392  ring XXRing = currRing;
    53625393  intvec* ivNull = new intvec(nV);
    53635394  intvec* curr_weight = new intvec(nV);
    53645395  intvec* target_weight = new intvec(nV);
    53655396  intvec* exivlp = Mivlp(nV);
     5397  intvec* next_weight= new intvec(nV);
     5398/*
    53665399  intvec* tmp_weight = new intvec(nV);
    53675400  for(i=0; i<nV; i++)
     
    53695402    (*tmp_weight)[i] = (*target_M)[i];
    53705403  }
     5404*/
    53715405  for(i=0; i<nV; i++)
    53725406  {
     
    53855419#endif
    53865420  rComplete(currRing);
    5387 #ifdef CHECK_IDEAL_MWALK
    5388     idString(Go,"Go");
    5389 #endif
    53905421#ifdef TIME_TEST
    53915422  to = clock();
    53925423#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       }
     5424  if(orig_M->length() == nV)
     5425  {
     5426    newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
     5427  }
     5428  else
     5429  {
     5430    newRing = VMatrDefault(orig_M);
     5431  }
    54015432  rChangeCurrRing(newRing);
    54025433  ideal G = MstdCC(idrMoveR(Go,baseRing,currRing));
     
    54145445    to = clock();
    54155446#endif
    5416 #ifdef CHECK_IDEAL_MWALK
    5417     idString(G,"G");
    5418 #endif
     5447//#ifdef CHECK_IDEAL_MWALK
     5448    if(printout > 2)
     5449    {
     5450      idString(G,"//** Mrwalk: G");
     5451    }
     5452//#endif
    54195453    Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
     5454    //polylength = 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise
     5455    polylength = lengthpoly(Gomega);
    54205456#ifdef TIME_TEST
    54215457    tif = tif + clock()-to; //time for computing initial form ideal
    54225458#endif
    5423 #ifdef CHECK_IDEAL_MWALK
    5424     idString(Gomega,"Gomega");
    5425 #endif
     5459//#ifdef CHECK_IDEAL_MWALK
     5460    if(printout > 1)
     5461    {
     5462      idString(Gomega,"//** Mrwalk: Gomega");
     5463    }
     5464//#endif
     5465    // test whether target cone is reached
     5466    if(test_w_in_ConeCC(G,target_weight) == 1)
     5467    {
     5468      endwalks = TRUE;
     5469    }
     5470    if(reduction == 0)
     5471    {
     5472      //PrintS("\n//** Mrwalk: test middle of cone!\n");
     5473      FF = middleOfCone(G,Gomega);
     5474      //PrintS("\n//** Mrwalk: Test F!\n");
     5475      if(FF != NULL)
     5476      {
     5477        idDelete(&G);
     5478        G = idCopy(FF);
     5479        idDelete(&FF);
     5480        //PrintS("\n//** Mrwalk: FF nicht NULL! Compue next vector.\n");
     5481        goto NEXT_VECTOR;
     5482      }
     5483    }
    54265484#ifndef  BUCHBERGER_ALG
    54275485    if(isNolVector(curr_weight) == 0)
    54285486    {
    54295487      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
    5430     }
     5488    }   
    54315489    else
    54325490    {
     
    54735531#endif
    54745532    idSkipZeroes(M);
    5475 #ifdef CHECK_IDEAL_MWALK
    5476     PrintS("\n//** Mwalk: computed M.\n");
    5477     idString(M, "M");
    5478 #endif
     5533//#ifdef CHECK_IDEAL_MWALK
     5534    if(printout > 2)
     5535    {
     5536      idString(M, "//** Mrwalk: M");
     5537    }
     5538//#endif
    54795539    //change the ring to baseRing
    54805540    rChangeCurrRing(baseRing);
     
    54865546    to = clock();
    54875547#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
     5548    // compute a representation of the generators of submod (M) with respect to those of mod (Gomega),
     5549    // where Gomega is a reduced Groebner basis w.r.t. the current ring
    54895550    F = MLifttwoIdeal(Gomega2, M1, G);
    54905551#ifdef TIME_TEST
    54915552    tlift = tlift + clock() - to;
    54925553#endif
    5493 #ifdef CHECK_IDEAL_MWALK
    5494     idString(F, "F");
    5495 #endif
     5554//#ifdef CHECK_IDEAL_MWALK
     5555    if(printout > 2)
     5556    {
     5557      idString(F, "//** Mrwalk: F");
     5558    }
     5559//#endif
    54965560    idDelete(&Gomega2);
    54975561    idDelete(&M1);
     
    55025566#ifdef TIME_TEST
    55035567    to = clock();
    5504 #endif
    5505     //G = kStd(F1,NULL,testHomog,NULL,NULL,0,0,NULL);
    5506 #ifdef TIME_TEST
    55075568    tstd = tstd + clock() - to;
    55085569#endif
    55095570    idSkipZeroes(G);
    5510 #ifdef CHECK_IDEAL_MWALK
    5511     idString(G, "G");
    5512 #endif
     5571//#ifdef CHECK_IDEAL_MWALK
     5572    if(printout > 2)
     5573    {
     5574      idString(G, "//** Mrwalk: G");
     5575    }
     5576//#endif
     5577    NEXT_VECTOR:
    55135578#ifdef TIME_TEST
    55145579    to = clock();
    55155580#endif
    5516     intvec* next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, pert_deg);//next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
     5581    next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
     5582    if(polylength > 0)
     5583    {
     5584      //there is a polynomial in Gomega with at least 3 monomials,
     5585      //low-dimensional facet of the cone
     5586      delete next_weight;
     5587      next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, pert_deg);
     5588    }
    55175589#ifdef TIME_TEST
    55185590    tnw = tnw + clock() - to;
    55195591#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
     5592//#ifdef PRINT_VECTORS
     5593    if(printout > 0)
     5594    {
     5595      MivString(curr_weight, target_weight, next_weight);
     5596    }
     5597//#endif
     5598    if(MivComp(next_weight, ivNull) == 1 || MivComp(target_weight,curr_weight) == 1 || endwalks == TRUE)
     5599    {
     5600//#ifdef CHECK_IDEAL_MWALK
     5601      if(printout > 0)
     5602      {
     5603        PrintS("\n//** Mrwalk: entering last cone.\n");
     5604      }
     5605//#endif
     5606
    55285607      Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
    55295608      if(target_M->length() == nV)
     
    55385617      Gomega1 = idrMoveR(Gomega, baseRing,currRing);
    55395618      idDelete(&Gomega);
    5540 #ifdef CHECK_IDEAL_MWALK
    5541       idString(Gomega1, "Gomega");
    5542 #endif
     5619//#ifdef CHECK_IDEAL_MWALK
     5620      if(printout > 1)
     5621      {
     5622        idString(Gomega1, "//** Mrwalk: Gomega");
     5623      }
     5624      PrintS("\n //** Mrwalk: kStd(Gomega)");
     5625//#endif
    55435626      M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL);
    5544 #ifdef CHECK_IDEAL_MWALK
    5545       idString(M,"M");
    5546 #endif
     5627//#ifdef CHECK_IDEAL_MWALK
     5628      if(printout > 1)
     5629      {
     5630        idString(M,"//** Mrwalk: M");
     5631      }
     5632//#endif
    55475633      rChangeCurrRing(baseRing);
    55485634      M1 =  idrMoveR(M, newRing,currRing);
     
    55505636      Gomega2 = idrMoveR(Gomega1, newRing,currRing);
    55515637      idDelete(&Gomega1);
     5638      //PrintS("\n //** Mrwalk: MLifttwoIdeal");
    55525639      F = MLifttwoIdeal(Gomega2, M1, G);
    5553 #ifdef CHECK_IDEAL_MWALK
    5554       idString(F,"F");
    5555 #endif
     5640//#ifdef CHECK_IDEAL_MWALK
     5641      if(printout > 2)
     5642      {
     5643        idString(F,"//** Mrwalk: F");
     5644      }
     5645//#endif
    55565646      idDelete(&Gomega2);
    55575647      idDelete(&M1);
     
    55655655      to = clock();
    55665656#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   //    }
     5657      /*PrintS("\n //**Mrwalk: Interreduce");
     5658      //interreduce the Groebner basis <G> w.r.t. currRing
     5659      G = kInterRedCC(G,NULL);*/
    55715660#ifdef TIME_TEST
    55725661      tred = tred + clock() - to;
     
    55755664      delete next_weight;
    55765665      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
     5666    }
     5667
    55845668    for(i=nV-1; i>=0; i--)
    55855669    {
    5586       (*tmp_weight)[i] = (*curr_weight)[i];
     5670      //(*tmp_weight)[i] = (*curr_weight)[i];
    55875671      (*curr_weight)[i] = (*next_weight)[i];
    55885672    }
    55895673    delete next_weight;
    55905674  }
     5675  baseRing = currRing;
    55915676  rChangeCurrRing(XXRing);
    55925677  ideal result = idrMoveR(G,baseRing,currRing);
    55935678  idDelete(&G);
    5594 /*#ifdef CHECK_IDEAL_MWALK
    5595   pDelete(&p);
    5596 #endif*/
    5597   delete tmp_weight;
     5679  si_opt_1 = save1; //set original options, e. g. option(RedSB)
     5680  //delete tmp_weight;
    55985681  delete ivNull;
    55995682  delete exivlp;
     
    56015684  delete last_omega;
    56025685#endif
     5686  Print("\n//** Mrwalk: Groebner Walk took %d steps.\n", nstep);
    56035687#ifdef TIME_TEST
    5604   Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep);
    56055688  TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep);
    5606   Print("\n//** Mwalk: Ergebnis.\n");
    56075689  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
    56085690  //Print("\n// Overflow_Error? (%d)\n", Overflow_Error);
     
    57515833// use kStd, if nP = 0, else call LastGB
    57525834ideal Mpwalk(ideal Go, int op_deg, int tp_deg,intvec* curr_weight,
    5753              intvec* target_weight, int nP)
     5835             intvec* target_weight, int nP, int reduction, int printout)
    57545836{
     5837  BITSET save1 = si_opt_1; // save current options
     5838  if(reduction == 0)
     5839  {
     5840    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     5841    //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     5842    //si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB));
     5843  }
    57555844  Set_Error(FALSE  );
    57565845  Overflow_Error = FALSE;
     
    57665855  nstep = 0;
    57675856  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;
     5857  BOOLEAN endwalks = FALSE;
     5858
     5859  ideal Gomega, M, F, FF, G, Gomega1, Gomega2, M1,F1,Eresult,ssG;
    57715860  ring newRing, oldRing, TargetRing;
    57725861  intvec* iv_M_dp;
     
    57905879  ring XXRing = currRing;
    57915880
    5792 
    57935881  to = clock();
    5794   /* perturbs the original vector */
     5882  // perturbs the original vector
    57955883  if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) := "dp"
    57965884  {
     
    58095897      DefRingPar(curr_weight);
    58105898    else
    5811       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 1
     5899      rChangeCurrRing(VMrDefault(curr_weight));
    58125900
    58135901    G = idrMoveR(Go, XXRing,currRing);
     
    58245912  ring HelpRing = currRing;
    58255913
    5826   /* perturbs the target weight vector */
     5914  // perturbs the target weight vector
    58275915  if(tp_deg > 1 && tp_deg <= nV)
    58285916  {
     
    58305918      DefRingPar(target_weight);
    58315919    else
    5832       rChangeCurrRing(VMrDefault(target_weight)); // Aenderung 2
     5920      rChangeCurrRing(VMrDefault(target_weight));
    58335921
    58345922    TargetRing = currRing;
     
    58375925    {
    58385926      iv_M_lp = MivMatrixOrderlp(nV);
    5839       //ivString(iv_M_lp, "iv_M_lp");
    5840       //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg);
    58415927      target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
    58425928    }
     
    58445930    {
    58455931      iv_M_lp = MivMatrixOrder(target_weight);
    5846       //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg);
    58475932      target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
    58485933    }
     
    58525937    G = idrMoveR(ssG, TargetRing,currRing);
    58535938  }
    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   */
     5939    if(printout > 0)
     5940    {
     5941      Print("\n//** Mpwalk: Perturbation Walk of degree (%d,%d):",op_deg,tp_deg);
     5942      ivString(curr_weight, "//** Mpwalk: new current weight");
     5943      ivString(target_weight, "//** Mpwalk: new target weight");
     5944    }
    58595945  while(1)
    58605946  {
    58615947    nstep ++;
    58625948    to = clock();
    5863     /* compute an initial form ideal of <G> w.r.t. the weight vector
    5864        "curr_weight" */
     5949    // compute an initial form ideal of <G> w.r.t. the weight vector
     5950    // "curr_weight"
    58655951    Gomega = MwalkInitialForm(G, curr_weight);
    5866 
     5952//#ifdef CHECK_IDEAL_MWALK
     5953    if(printout > 1)
     5954    {
     5955      idString(Gomega,"//** Mpwalk: Gomega");
     5956    }
     5957//#endif
     5958    if(reduction == 0 && nstep > 1)
     5959    {
     5960      // check whether weight vector is in the interior of the cone
     5961      while(1)
     5962      {
     5963        FF = middleOfCone(G,Gomega);
     5964        if(FF != NULL)
     5965        {
     5966          idDelete(&G);
     5967          G = idCopy(FF);
     5968          idDelete(&FF);
     5969          next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
     5970          if(printout > 0)
     5971          {
     5972            MivString(curr_weight, target_weight, next_weight);
     5973          }
     5974        }
     5975        else
     5976        {
     5977          break;
     5978        }
     5979        for(i=nV-1; i>=0; i--)
     5980        {
     5981          (*curr_weight)[i] = (*next_weight)[i];
     5982        }
     5983        Gomega = MwalkInitialForm(G, curr_weight);
     5984        if(printout > 1)
     5985        {
     5986          idString(Gomega,"//** Mpwalk: Gomega");
     5987        }
     5988      }
     5989    }
    58675990
    58685991#ifdef ENDWALKS
    5869     if(endwalks == 1){
     5992    if(endwalks == TRUE)
     5993    {
    58705994      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
    58715995      idElements(G, "G");
    5872       // idElements(Gomega, "Gw");
    58735996      headidString(G, "G");
    5874       //headidString(Gomega, "Gw");
    58755997    }
    58765998#endif
     
    58916013      DefRingPar(curr_weight);
    58926014    else
    5893       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 3
     6015      rChangeCurrRing(VMrDefault(curr_weight));
    58946016
    58956017    newRing = currRing;
     
    58976019
    58986020#ifdef ENDWALKS
    5899     if(endwalks==1)
     6021    if(endwalks==TRUE)
    59006022    {
    59016023      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     
    59126034    tim = clock();
    59136035    to = clock();
    5914     /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
     6036    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
    59156037#ifdef  BUCHBERGER_ALG
    59166038    M = MstdhomCC(Gomega1);
     
    59186040    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
    59196041    delete hilb_func;
    5920 #endif // BUCHBERGER_ALG
    5921 
    5922     if(endwalks == 1){
     6042#endif
     6043//#ifdef CHECK_IDEAL_MWALK
     6044      if(printout > 2)
     6045      {
     6046        idString(M,"//** Mpwalk: M");
     6047      }
     6048//#endif
     6049
     6050    if(endwalks == TRUE){
    59236051      xtstd = xtstd+clock()-to;
    59246052#ifdef ENDWALKS
     
    59306058      tstd=tstd+clock()-to;
    59316059
    5932     /* change the ring to oldRing */
     6060    // change the ring to oldRing
    59336061    rChangeCurrRing(oldRing);
    59346062    M1 =  idrMoveR(M, newRing,currRing);
    59356063    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
    5936 
    5937     //if(endwalks==1)  PrintS("\n// Lifting is working:..");
    59386064
    59396065    to=clock();
     
    59426068       Gomega is a reduced Groebner basis w.r.t. the current ring */
    59436069    F = MLifttwoIdeal(Gomega2, M1, G);
    5944     if(endwalks != 1)
     6070    if(endwalks == FALSE)
    59456071      tlift = tlift+clock()-to;
    59466072    else
    59476073      xtlift=clock()-to;
    59486074
     6075//#ifdef CHECK_IDEAL_MWALK
     6076    if(printout > 2)
     6077    {
     6078      idString(F,"//** Mpwalk: F");
     6079    }
     6080//#endif
     6081
    59496082    idDelete(&M1);
    59506083    idDelete(&Gomega2);
    59516084    idDelete(&G);
    59526085
    5953     /* change the ring to newRing */
     6086    // change the ring to newRing
    59546087    rChangeCurrRing(newRing);
    5955     F1 = idrMoveR(F, oldRing,currRing);
    5956 
    5957     //if(endwalks==1)PrintS("\n// InterRed is working now:");
     6088    if(reduction == 0)
     6089    {
     6090      G = idrMoveR(F,oldRing,currRing);
     6091    }
     6092    else
     6093    {
     6094      F1 = idrMoveR(F, oldRing,currRing);
     6095      if(printout > 2)
     6096      {
     6097        PrintS("\n //** Mpwalk: reduce the Groebner basis.\n");
     6098      }
     6099      to=clock();
     6100      G = kInterRedCC(F1, NULL);
     6101      if(endwalks == FALSE)
     6102        tred = tred+clock()-to;
     6103      else
     6104        xtred=clock()-to;
     6105      idDelete(&F1);
     6106    }
     6107    if(endwalks == TRUE)
     6108      break;
    59586109
    59596110    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 */
     6111    // compute a next weight vector
    59746112    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
    59756113    tnw=tnw+clock()-to;
    5976 #ifdef PRINT_VECTORS
    5977     MivString(curr_weight, target_weight, next_weight);
    5978 #endif
     6114//#ifdef PRINT_VECTORS
     6115    if(printout > 0)
     6116    {
     6117      MivString(curr_weight, target_weight, next_weight);
     6118    }
     6119//#endif
    59796120
    59806121    if(Overflow_Error == TRUE)
     
    59946135    }
    59956136    if(MivComp(next_weight, target_weight) == 1)
    5996       endwalks = 1;
     6137      endwalks = TRUE;
    59976138
    59986139    for(i=nV-1; i>=0; i--)
     
    60006141
    60016142    delete next_weight;
    6002   }//while
     6143  }//end of while-loop
    60036144
    60046145  if(tp_deg != 1)
     
    60146155        DefRingPar(orig_target);
    60156156      else
    6016         rChangeCurrRing(VMrDefault(orig_target)); //Aenderung
     6157        rChangeCurrRing(VMrDefault(orig_target));
    60176158
    60186159    TargetRing=currRing;
     
    60306171    if( ntestw != 1 || ntwC == 0)
    60316172    {
    6032       /*
    6033       if(ntestw != 1){
     6173      if(ntestw != 1 && printout >2)
     6174      {
    60346175        ivString(pert_target_vector, "tau");
    60356176        PrintS("\n// ** perturbed target vector doesn't stay in cone!!");
    60366177        Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
    6037         idElements(F1, "G");
    6038       }
    6039       */
     6178        //idElements(F1, "G");
     6179      }
    60406180      // LastGB is "better" than the kStd subroutine
    60416181      to=clock();
     
    60686208    Eresult = idrMoveR(G, newRing,currRing);
    60696209  }
     6210  si_opt_1 = save1; //set original options, e. g. option(RedSB)
    60706211  delete ivNull;
    60716212  if(tp_deg != 1)
     
    60826223             tnw+xtnw);
    60836224
    6084   Print("\n// pSetm_Error = (%d)", ErrorCheck());
    6085   Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep,  Overflow_Error);
    6086 #endif
     6225  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
     6226  //Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep,  Overflow_Error);
     6227#endif
     6228  Print("\n//** Mpwalk: Perturbation Walk took %d steps.\n", nstep);
     6229  return(Eresult);
     6230}
     6231
     6232/*******************************************************
     6233 * THE PERTURBATION WALK ALGORITHM WITH RANDOM ELEMENT *
     6234 *******************************************************/
     6235ideal Mprwalk(ideal Go, intvec* orig_M, intvec* target_M, int weight_rad,
     6236              int op_deg, int tp_deg, int nP, int reduction, int printout)
     6237{
     6238  BITSET save1 = si_opt_1; // save current options
     6239  if(reduction == 0)
     6240  {
     6241    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     6242    //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     6243    //si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB));
     6244  }
     6245  Set_Error(FALSE);
     6246  Overflow_Error = FALSE;
     6247  //Print("// pSetm_Error = (%d)", ErrorCheck());
     6248
     6249  clock_t  tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
     6250  xtextra=0;
     6251  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0;
     6252  tinput = clock();
     6253
     6254  clock_t tim;
     6255
     6256  nstep = 0;
     6257  int i, ntwC=1, ntestw=1, polylength, nV = currRing->N;
     6258  BOOLEAN endwalks = FALSE;
     6259
     6260  ideal Gomega, M, F, FF, G, Gomega1, Gomega2, M1,F1,Eresult,ssG;
     6261  ring newRing, oldRing, TargetRing;
     6262  intvec* iv_M_dp;
     6263  intvec* iv_M_lp;
     6264  intvec* exivlp = Mivlp(nV);
     6265  intvec* curr_weight = new intvec(nV);
     6266  intvec* target_weight = new intvec(nV);
     6267  for(i=0; i<nV; i++)
     6268  {
     6269    (*curr_weight)[i] = (*orig_M)[i];
     6270    (*target_weight)[i] = (*target_M)[i];
     6271  }
     6272  intvec* orig_target = target_weight;
     6273  intvec* pert_target_vector = target_weight;
     6274  intvec* ivNull = new intvec(nV);
     6275  intvec* iv_dp = MivUnit(nV);// define (1,1,...,1)
     6276#ifndef BUCHBERGER_ALG
     6277  intvec* hilb_func;
     6278#endif
     6279  intvec* next_weight;
     6280
     6281  // to avoid (1,0,...,0) as the target vector
     6282  intvec* last_omega = new intvec(nV);
     6283  for(i=nV-1; i>0; i--)
     6284    (*last_omega)[i] = 1;
     6285  (*last_omega)[0] = 10000;
     6286
     6287  ring XXRing = currRing;
     6288
     6289  to = clock();
     6290  // perturbs the original vector
     6291  if(orig_M->length() == nV)
     6292  {
     6293    if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) := "dp"
     6294    {
     6295      G = MstdCC(Go);
     6296      tostd = clock()-to;
     6297      if(op_deg != 1)
     6298      {
     6299        iv_M_dp = MivMatrixOrderdp(nV);
     6300        curr_weight = MPertVectors(G, iv_M_dp, op_deg);
     6301      }
     6302    }
     6303    else
     6304    {
     6305      //define ring order := (a(curr_weight),lp);
     6306      if (rParameter(currRing) != NULL)
     6307        DefRingPar(curr_weight);
     6308      else
     6309        rChangeCurrRing(VMrDefault(curr_weight));
     6310
     6311      G = idrMoveR(Go, XXRing,currRing);
     6312      G = MstdCC(G);
     6313      tostd = clock()-to;
     6314      if(op_deg != 1)
     6315      {
     6316        iv_M_dp = MivMatrixOrder(curr_weight);
     6317        curr_weight = MPertVectors(G, iv_M_dp, op_deg);
     6318      }
     6319    }
     6320  }
     6321  else
     6322  {
     6323    rChangeCurrRing(VMatrDefault(orig_M));
     6324    G = idrMoveR(Go, XXRing,currRing);
     6325    G = MstdCC(G);
     6326    tostd = clock()-to;
     6327    if(op_deg != 1)
     6328    {
     6329      curr_weight = MPertVectors(G, orig_M, op_deg);
     6330    }
     6331  }
     6332
     6333  delete iv_dp;
     6334  if(op_deg != 1) delete iv_M_dp;
     6335
     6336  ring HelpRing = currRing;
     6337
     6338  // perturbs the target weight vector
     6339  if(target_M->length() == nV)
     6340  {
     6341    if(tp_deg > 1 && tp_deg <= nV)
     6342    {
     6343      if (rParameter(currRing) != NULL)
     6344        DefRingPar(target_weight);
     6345      else
     6346        rChangeCurrRing(VMrDefault(target_weight));
     6347
     6348      TargetRing = currRing;
     6349      ssG = idrMoveR(G,HelpRing,currRing);
     6350      if(MivSame(target_weight, exivlp) == 1)
     6351      {
     6352        iv_M_lp = MivMatrixOrderlp(nV);
     6353        target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
     6354      }
     6355      else
     6356      {
     6357        iv_M_lp = MivMatrixOrder(target_weight);
     6358        target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
     6359      }
     6360      delete iv_M_lp;
     6361      pert_target_vector = target_weight;
     6362      rChangeCurrRing(HelpRing);
     6363      G = idrMoveR(ssG, TargetRing,currRing);
     6364    }
     6365  }
     6366  else
     6367  {
     6368    if(tp_deg > 1 && tp_deg <= nV)
     6369    {
     6370      rChangeCurrRing(VMatrDefault(target_M));
     6371      TargetRing = currRing;
     6372      ssG = idrMoveR(G,HelpRing,currRing);
     6373      target_weight = MPertVectors(ssG, target_M, tp_deg);
     6374    }
     6375  }
     6376  if(printout > 0)
     6377  {
     6378    Print("\n//** Mprwalk: Random Perturbation Walk of degree (%d,%d):",op_deg,tp_deg);
     6379    ivString(curr_weight, "//** Mprwalk: new current weight");
     6380    ivString(target_weight, "//** Mprwalk: new target weight");
     6381  }
     6382  while(1)
     6383  {
     6384    nstep ++;
     6385    to = clock();
     6386    // compute an initial form ideal of <G> w.r.t. the weight vector
     6387    // "curr_weight"
     6388    Gomega = MwalkInitialForm(G, curr_weight);
     6389    polylength = lengthpoly(Gomega);
     6390//#ifdef CHECK_IDEAL_MWALK
     6391    if(printout > 1)
     6392    {
     6393      idString(Gomega,"//** Mprwalk: Gomega");
     6394    }
     6395//#endif
     6396
     6397    if(reduction == 0 && nstep > 1)
     6398    {
     6399      // check whether weight vector is in the interior of the cone
     6400      while(1)
     6401      {
     6402        FF = middleOfCone(G,Gomega);
     6403        if(FF != NULL)
     6404        {
     6405          idDelete(&G);
     6406          G = idCopy(FF);
     6407          idDelete(&FF);
     6408          next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
     6409          if(printout > 0)
     6410          {
     6411            MivString(curr_weight, target_weight, next_weight);
     6412          }
     6413        }
     6414        else
     6415        {
     6416          break;
     6417        }
     6418        for(i=nV-1; i>=0; i--)
     6419        {
     6420          (*curr_weight)[i] = (*next_weight)[i];
     6421        }
     6422        Gomega = MwalkInitialForm(G, curr_weight);
     6423        if(printout > 1)
     6424        {
     6425          idString(Gomega,"//** Mprwalk: Gomega");
     6426        }
     6427      }
     6428    }
     6429
     6430#ifdef ENDWALKS
     6431    if(endwalks == TRUE)
     6432    {
     6433      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6434      idElements(G, "G");
     6435      headidString(G, "G");
     6436    }
     6437#endif
     6438
     6439    tif = tif + clock()-to;
     6440
     6441#ifndef  BUCHBERGER_ALG
     6442    if(isNolVector(curr_weight) == 0)
     6443      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
     6444    else
     6445      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
     6446#endif // BUCHBERGER_ALG
     6447
     6448    oldRing = currRing;
     6449
     6450    if(target_M->length() == nV)
     6451    {
     6452      // define a new ring with ordering "(a(curr_weight),lp)
     6453      if (rParameter(currRing) != NULL)
     6454        DefRingPar(curr_weight);
     6455      else
     6456        rChangeCurrRing(VMrDefault(curr_weight));
     6457    }
     6458    else
     6459    {
     6460      rChangeCurrRing(VMatrRefine(target_M,curr_weight));
     6461    }
     6462    newRing = currRing;
     6463    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
     6464#ifdef ENDWALKS
     6465    if(endwalks == TRUE)
     6466    {
     6467      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6468      idElements(Gomega1, "Gw");
     6469      headidString(Gomega1, "headGw");
     6470      PrintS("\n// compute a rGB of Gw:\n");
     6471
     6472#ifndef  BUCHBERGER_ALG
     6473      ivString(hilb_func, "w");
     6474#endif
     6475    }
     6476#endif
     6477
     6478    tim = clock();
     6479    to = clock();
     6480    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
     6481#ifdef  BUCHBERGER_ALG
     6482    M = MstdhomCC(Gomega1);
     6483#else
     6484    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
     6485    delete hilb_func;
     6486#endif
     6487//#ifdef CHECK_IDEAL_MWALK
     6488    if(printout > 2)
     6489    {
     6490      idString(M,"//** Mprwalk: M");
     6491    }
     6492//#endif
     6493
     6494    if(endwalks == TRUE)
     6495    {
     6496      xtstd = xtstd+clock()-to;
     6497#ifdef ENDWALKS
     6498      Print("\n// time for the last std(Gw)  = %.2f sec\n",
     6499            ((double) clock())/1000000 -((double)tim) /1000000);
     6500#endif
     6501    }
     6502    else
     6503      tstd=tstd+clock()-to;
     6504
     6505    /* change the ring to oldRing */
     6506    rChangeCurrRing(oldRing);
     6507    M1 =  idrMoveR(M, newRing,currRing);
     6508    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
     6509
     6510    to=clock();
     6511    /* compute a representation of the generators of submod (M)
     6512       with respect to those of mod (Gomega).
     6513       Gomega is a reduced Groebner basis w.r.t. the current ring */
     6514    F = MLifttwoIdeal(Gomega2, M1, G);
     6515    if(endwalks == FALSE)
     6516      tlift = tlift+clock()-to;
     6517    else
     6518      xtlift=clock()-to;
     6519
     6520//#ifdef CHECK_IDEAL_MWALK
     6521    if(printout > 2)
     6522    {
     6523      idString(F,"//** Mprwalk: F");
     6524    }
     6525//#endif
     6526
     6527    idDelete(&M1);
     6528    idDelete(&Gomega2);
     6529    idDelete(&G);
     6530
     6531    // change the ring to newRing
     6532    rChangeCurrRing(newRing);
     6533    if(reduction == 0)
     6534    {
     6535      G = idrMoveR(F,oldRing,currRing);
     6536    }
     6537    else
     6538    {
     6539      F1 = idrMoveR(F, oldRing,currRing);
     6540      if(printout > 2)
     6541      {
     6542        PrintS("\n //** Mprwalk: reduce the Groebner basis.\n");
     6543      }
     6544      to=clock();
     6545      G = kInterRedCC(F1, NULL);
     6546      if(endwalks == FALSE)
     6547        tred = tred+clock()-to;
     6548      else
     6549        xtred=clock()-to;
     6550      idDelete(&F1);
     6551    }
     6552
     6553    if(endwalks == TRUE)
     6554      break;
     6555
     6556    to=clock();
     6557
     6558    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
     6559    if(polylength > 0)
     6560    {
     6561      if(printout > 2)
     6562      {
     6563        Print("\n //**Mprwalk: there is a polynomial in Gomega with at least 3 monomials, low-dimensional facet of the cone.\n");
     6564      }
     6565      delete next_weight;
     6566      next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, op_deg);
     6567    }
     6568    tnw=tnw+clock()-to;
     6569//#ifdef PRINT_VECTORS
     6570    if(printout > 0)
     6571    {
     6572      MivString(curr_weight, target_weight, next_weight);
     6573    }
     6574//#endif
     6575
     6576    if(Overflow_Error == TRUE)
     6577    {
     6578      ntwC = 0;
     6579      //ntestomega = 1;
     6580      //Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6581      //idElements(G, "G");
     6582      delete next_weight;
     6583      goto FINISH_160302;
     6584    }
     6585    if(MivComp(next_weight, ivNull) == 1){
     6586      newRing = currRing;
     6587      delete next_weight;
     6588      //Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6589      break;
     6590    }
     6591    if(MivComp(next_weight, target_weight) == 1)
     6592      endwalks = TRUE;
     6593
     6594    for(i=nV-1; i>=0; i--)
     6595      (*curr_weight)[i] = (*next_weight)[i];
     6596
     6597    delete next_weight;
     6598  }//while
     6599
     6600  if(tp_deg != 1)
     6601  {
     6602    FINISH_160302:
     6603    if(target_M->length() == nV)
     6604    {
     6605      if(MivSame(orig_target, exivlp) == 1)
     6606        if (rParameter(currRing) != NULL)
     6607          DefRingParlp();
     6608        else
     6609          VMrDefaultlp();
     6610      else
     6611        if (rParameter(currRing) != NULL)
     6612          DefRingPar(orig_target);
     6613        else
     6614          rChangeCurrRing(VMrDefault(orig_target));
     6615    }
     6616    else
     6617    {
     6618      rChangeCurrRing(VMatrDefault(target_M));
     6619    }
     6620    TargetRing=currRing;
     6621    F1 = idrMoveR(G, newRing,currRing);
     6622#ifdef CHECK_IDEAL
     6623      headidString(G, "G");
     6624#endif
     6625
     6626    // check whether the pertubed target vector stays in the correct cone
     6627    if(ntwC != 0)
     6628    {
     6629      ntestw = test_w_in_ConeCC(F1, pert_target_vector);
     6630    }
     6631    if(ntestw != 1 || ntwC == 0)
     6632    {
     6633      if(ntestw != 1 && printout > 2)
     6634      {
     6635        ivString(pert_target_vector, "tau");
     6636        PrintS("\n// **Mprwalk: perturbed target vector doesn't stay in cone.");
     6637        Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6638        //idElements(F1, "G");
     6639      }
     6640      // LastGB is "better" than the kStd subroutine
     6641      to=clock();
     6642      ideal eF1;
     6643      if(nP == 0 || tp_deg == 1 || MivSame(orig_target, exivlp) != 1 || target_M->length() != nV)
     6644      {
     6645        if(printout > 2)
     6646        {
     6647          PrintS("\n// ** Mprwalk: Call \"std\" to compute a Groebner basis.\n");
     6648        }
     6649        eF1 = MstdCC(F1);
     6650        idDelete(&F1);
     6651      }
     6652      else
     6653      {
     6654        if(printout > 2)
     6655        {
     6656          PrintS("\n// **Mprwalk: Call \"LastGB\" to compute a Groebner basis.\n");
     6657        }
     6658        rChangeCurrRing(newRing);
     6659        ideal F2 = idrMoveR(F1, TargetRing,currRing);
     6660        eF1 = LastGB(F2, curr_weight, tp_deg-1);
     6661        F2=NULL;
     6662      }
     6663      xtextra=clock()-to;
     6664      ring exTargetRing = currRing;
     6665
     6666      rChangeCurrRing(XXRing);
     6667      Eresult = idrMoveR(eF1, exTargetRing,currRing);
     6668    }
     6669    else{
     6670      rChangeCurrRing(XXRing);
     6671      Eresult = idrMoveR(F1, TargetRing,currRing);
     6672    }
     6673  }
     6674  else {
     6675    rChangeCurrRing(XXRing);
     6676    Eresult = idrMoveR(G, newRing,currRing);
     6677  }
     6678  si_opt_1 = save1; //set original options, e. g. option(RedSB)
     6679  delete ivNull;
     6680  if(tp_deg != 1)
     6681    delete target_weight;
     6682
     6683  if(op_deg != 1 )
     6684    delete curr_weight;
     6685
     6686  delete exivlp;
     6687  delete last_omega;
     6688
     6689#ifdef TIME_TEST
     6690  TimeStringFractal(tinput, tostd, tif+xtif, tstd+xtstd,0, tlift+xtlift, tred+xtred,
     6691             tnw+xtnw);
     6692
     6693  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
     6694  //Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep,  Overflow_Error);
     6695#endif
     6696  Print("\n//** Mprwalk: Perturbation Walk took %d steps.\n", nstep);
    60876697  return(Eresult);
    60886698}
     
    61106720 * Perturb the start weight vector at the top level, i.e. nlev = 1     *
    61116721 ***********************************************************************/
    6112 static ideal rec_fractal_call(ideal G, int nlev, intvec* omtmp)
     6722static ideal rec_fractal_call(ideal G, int nlev, intvec* ivtarget,
     6723             int reduction, int printout)
    61136724{
    61146725  Overflow_Error =  FALSE;
     
    61186729  ring new_ring, testring;
    61196730  //ring extoRing;
    6120   ideal Gomega, Gomega1, Gomega2, F, F1, Gresult, Gresult1, G1, Gt;
     6731  ideal Gomega, Gomega1, Gomega2, FF, F, F1, Gresult, Gresult1, G1, Gt;
    61216732  int nwalks = 0;
    61226733  intvec* Mwlp;
     
    61276738  intvec* next_vect;
    61286739  intvec* omega2 = new intvec(nV);
    6129   intvec* altomega = new intvec(nV);
    6130 
     6740  intvec* omtmp = new intvec(nV);
     6741  //intvec* altomega = new intvec(nV);
     6742
     6743  for(i = nV -1; i>0; i--)
     6744  {
     6745    (*omtmp)[i] = (*ivtarget)[i];
     6746  }
    61316747  //BOOLEAN isnewtarget = FALSE;
    61326748
     
    61696785    NEXT_VECTOR_FRACTAL:
    61706786    to=clock();
    6171     /* determine the next border */
     6787    // determine the next border
    61726788    next_vect = MkInterRedNextWeight(omega,omega2,G);
    61736789    xtnw=xtnw+clock()-to;
    6174 #ifdef PRINT_VECTORS
    6175     MivString(omega, omega2, next_vect);
    6176 #endif
     6790
    61776791    oRing = currRing;
    61786792
    6179     /* We only perturb the current target vector at the recursion  level 1 */
     6793    // We only perturb the current target vector at the recursion level 1
    61806794    if(Xngleich == 0 && nlev == 1) //(ngleich == 0) important, e.g. ex2, ex3
    61816795      if (MivComp(next_vect, omega2) != 1)
    61826796      {
    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");
     6797        // to dispense with taking initial (and lifting/interreducing
     6798        // after the call of recursion
     6799        if(printout > 0)
     6800        {
     6801          Print("\n//** rec_fractal_call: Perturb the both vectors with degree %d.",nlev);
     6802          //idElements(G, "G");
     6803        }
    61876804
    61886805        Xngleich = 1;
    61896806        nlev +=1;
    61906807
    6191         if (rParameter(currRing) != NULL)
    6192           DefRingPar(omtmp);
     6808        if(ivtarget->length() == nV)
     6809        {
     6810          if (rParameter(currRing) != NULL)
     6811            DefRingPar(omtmp);
     6812          else
     6813            rChangeCurrRing(VMrDefault(omtmp));
     6814        }
    61936815        else
    6194           rChangeCurrRing(VMrDefault1(omtmp)); //Aenderung3
    6195 
     6816        {
     6817          rChangeCurrRing(VMatrDefault(ivtarget));
     6818        }
    61966819        testring = currRing;
    61976820        Gt = idrMoveR(G, oRing,currRing);
    61986821
    6199         /* perturb the original target vector w.r.t. the current GB */
    6200         delete Xtau;
    6201         Xtau = NewVectorlp(Gt);
     6822        // perturb the original target vector w.r.t. the current GB
     6823        if(ivtarget->length() == nV)
     6824        {
     6825          delete Xtau;
     6826          Xtau = NewVectorlp(Gt);
     6827        }
     6828        else
     6829        {
     6830          delete Xtau;
     6831          Xtau = Mfpertvector(Gt,ivtarget);
     6832        }
    62026833
    62036834        rChangeCurrRing(oRing);
    62046835        G = idrMoveR(Gt, testring,currRing);
    62056836
    6206         /* perturb the current vector w.r.t. the current GB */
     6837        // perturb the current vector w.r.t. the current GB
    62076838        Mwlp = MivWeightOrderlp(omega);
    62086839        Xsigma = Mfpertvector(G, Mwlp);
     
    62226853        next_vect = MkInterRedNextWeight(omega,omega2,G);
    62236854        xtnw=xtnw+clock()-to;
    6224 
    6225 #ifdef PRINT_VECTORS
     6855      }// end of (if MivComp(next_vect, omega2) == 1)
     6856
     6857//#ifdef PRINT_VECTORS
     6858      if(printout > 0)
     6859      {
    62266860        MivString(omega, omega2, next_vect);
    6227 #endif
    6228       }
    6229 
     6861      }
     6862//#endif
    62306863
    62316864    /* check whether the the computed vector is in the correct cone */
     
    62366869    {
    62376870      delete next_vect;
    6238       if (rParameter(currRing) != NULL)
    6239       {
    6240         DefRingPar(omtmp);
     6871      if(ivtarget->length() == nV)
     6872      {
     6873        if (rParameter(currRing) != NULL)
     6874          DefRingPar(omtmp);
     6875        else
     6876          rChangeCurrRing(VMrDefault(omtmp));
    62416877      }
    62426878      else
    62436879      {
    6244         rChangeCurrRing(VMrDefault1(omtmp)); // Aenderung4
     6880        rChangeCurrRing(VMatrDefault(ivtarget));
    62456881      }
    62466882#ifdef TEST_OVERFLOW
     
    62486884      Gt = NULL; return(Gt);
    62496885#endif
    6250 
    6251       //Print("\n\n// apply BB's alg. in ring r = %s;", rString(currRing));
     6886      if(printout > 0)
     6887      {
     6888        Print("\n//** rec_fractal_call: Overflow. Applying Buchberger's algorithm in ring r = %s;",
     6889              rString(currRing));
     6890      }
    62526891      to=clock();
    62536892      Gt = idrMoveR(G, oRing,currRing);
     
    62576896
    62586897      delete omega2;
    6259       delete altomega;
    6260 
    6261       //Print("\n// Leaving the %d-th recursion with %d steps", nlev, nwalks);
    6262       //Print("  ** Overflow_Error? (%d)", Overflow_Error);
     6898      //delete altomega;
     6899      if(printout > 0)
     6900      {
     6901        Print("\n//** rec_fractal_call: Overflow. Leaving the %d-th recursion with %d steps.\n",
     6902              nlev, nwalks);
     6903        //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     6904      }
     6905
    62636906      nnflow ++;
    62646907
    62656908      Overflow_Error = FALSE;
    62666909      return (G1);
    6267     }
     6910    }//end overflow-check
    62686911
    62696912
     
    62776920    if (MivComp(next_vect, XivNull) == 1)
    62786921    {
    6279       if (rParameter(currRing) != NULL)
    6280         DefRingPar(omtmp);
     6922      if(ivtarget->length() == nV)
     6923      {
     6924        if (rParameter(currRing) != NULL)
     6925          DefRingPar(omtmp);
     6926        else
     6927          rChangeCurrRing(VMrDefault(omtmp));
     6928      }
    62816929      else
    6282         rChangeCurrRing(VMrDefault1(omtmp)); //Aenderung5
     6930      {
     6931        rChangeCurrRing(VMatrDefault(ivtarget));
     6932      }
    62836933
    62846934      testring = currRing;
    62856935      Gt = idrMoveR(G, oRing,currRing);
    62866936
    6287       if(test_w_in_ConeCC(Gt, omega2) == 1) {
     6937      if(test_w_in_ConeCC(Gt, omega2) == 1)
     6938      {
     6939        delete omega2;
     6940        delete next_vect;
     6941        //delete altomega;
     6942        if(printout > 0)
     6943        {
     6944          Print("\n//** rec_fractal_call: Correct cone. Leaving the %d-th recursion with %d steps.\n",
     6945              nlev, nwalks);
     6946        }
     6947        return (Gt);
     6948      }
     6949      else
     6950      {
     6951        if(printout > 0)
     6952        {
     6953          Print("\n//** rec_fractal_call: Wrong cone. Tau' doesn't stay in the correct cone.\n");
     6954        }
     6955/*
     6956#ifndef  MSTDCC_FRACTAL
     6957        intvec* Xtautmp;
     6958        if(ivtarget->length() == nV)
     6959        {
     6960          Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp));
     6961        }
     6962        else
     6963        {
     6964          Xtautmp = Mfpertvector(Gt, ivtarget);
     6965        }
     6966#ifdef TEST_OVERFLOW
     6967      if(Overflow_Error == TRUE)
     6968      Gt = NULL; return(Gt);
     6969#endif
     6970
     6971        if(MivSame(Xtau, Xtautmp) == 1)
     6972        {
     6973          if(printout > 0)
     6974          {
     6975            Print("\n//** rec_fractal_call: Updated vectors are equal to the old vectors.\n");
     6976          }
     6977          delete Xtautmp;
     6978          goto FRACTAL_MSTDCC;
     6979        }
     6980
     6981        Xtau = Xtautmp;
     6982        Xtautmp = NULL;
     6983
     6984        for(i=nV-1; i>=0; i--)
     6985          (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i];
     6986
     6987        rChangeCurrRing(oRing);
     6988        G = idrMoveR(Gt, testring,currRing);
     6989
     6990        goto NEXT_VECTOR_FRACTAL;
     6991#endif
     6992*/
     6993      FRACTAL_MSTDCC:
     6994        if(printout > 0)
     6995        {
     6996          Print("\n//** rec_fractal_call: Wrong cone. Applying Buchberger's algorithm in ring = %s.\n",
     6997                rString(currRing));
     6998        }
     6999        to=clock();
     7000        G = MstdCC(Gt);
     7001        xtextra=xtextra+clock()-to;
     7002
     7003        oRing = currRing;
     7004
     7005        // update the original target vector w.r.t. the current GB
     7006        if(ivtarget->length() == nV)
     7007        {
     7008          if(MivSame(Xivinput, Xivlp) == 1)
     7009            if (rParameter(currRing) != NULL)
     7010              DefRingParlp();
     7011            else
     7012              VMrDefaultlp();
     7013          else
     7014            if (rParameter(currRing) != NULL)
     7015              DefRingPar(Xivinput);
     7016            else
     7017              rChangeCurrRing(VMrDefault(Xivinput));
     7018        }
     7019        else
     7020        {
     7021          rChangeCurrRing(VMatrRefine(ivtarget,Xivinput));
     7022        }
     7023        testring = currRing;
     7024        Gt = idrMoveR(G, oRing,currRing);
     7025
     7026        // perturb the original target vector w.r.t. the current GB
     7027        if(ivtarget->length() == nV)
     7028        {
     7029          delete Xtau;
     7030          Xtau = NewVectorlp(Gt);
     7031        }
     7032        else
     7033        {
     7034          delete Xtau;
     7035          Xtau = Mfpertvector(Gt,ivtarget);
     7036        }
     7037
     7038        rChangeCurrRing(oRing);
     7039        G = idrMoveR(Gt, testring,currRing);
     7040
     7041        delete omega2;
     7042        delete next_vect;
     7043        //delete altomega;
     7044        if(printout > 0)
     7045        {
     7046          Print("\n//** rec_fractal_call: Vectors updated. Leaving the %d-th recursion with %d steps.\n",
     7047              nlev, nwalks);
     7048          //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     7049        }
     7050        if(Overflow_Error == TRUE)
     7051          nnflow ++;
     7052
     7053        Overflow_Error = FALSE;
     7054        return(G);
     7055      }
     7056    }// end of (if next_vect==nullvector)
     7057
     7058    for(i=nV-1; i>=0; i--) {
     7059      //(*altomega)[i] = (*omega)[i];
     7060      (*omega)[i] = (*next_vect)[i];
     7061    }
     7062    delete next_vect;
     7063
     7064    to=clock();
     7065    /* Take the initial form of <G> w.r.t. omega */
     7066    Gomega = MwalkInitialForm(G, omega);
     7067    xtif=xtif+clock()-to;
     7068    if(printout > 1)
     7069    {
     7070      idString(Gomega,"//** rec_fractal_call: Gomega");
     7071    }
     7072
     7073    if(reduction == 0)
     7074    {
     7075      /* Check whether the intermediate weight vector lies in the interior of the cone.
     7076       * If so, only perform reductions. Otherwise apply Buchberger's algorithm. */
     7077      FF = middleOfCone(G,Gomega);
     7078      if( FF != NULL)
     7079      {
     7080        idDelete(&G);
     7081        G = idCopy(FF);
     7082        idDelete(&FF);
     7083        /* Compue next vector. */
     7084        goto NEXT_VECTOR_FRACTAL;
     7085      }
     7086    }
     7087
     7088#ifndef  BUCHBERGER_ALG
     7089    if(isNolVector(omega) == 0)
     7090      hilb_func = hFirstSeries(Gomega,NULL,NULL,omega,currRing);
     7091    else
     7092      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
     7093#endif
     7094
     7095    if(ivtarget->length() == nV)
     7096    {
     7097      if (rParameter(currRing) != NULL)
     7098        DefRingPar(omega);
     7099      else
     7100        rChangeCurrRing(VMrDefault(omega));
     7101    }
     7102    else
     7103    {
     7104      rChangeCurrRing(VMatrRefine(ivtarget,omega));
     7105    }
     7106    Gomega1 = idrMoveR(Gomega, oRing,currRing);
     7107
     7108    // Maximal recursion depth, to compute a red. GB
     7109    // Fractal walk with the alternative recursion
     7110    // alternative recursion
     7111    if(nlev == Xnlev || lengthpoly(Gomega1) == 0)
     7112    {
     7113      if(printout > 1)
     7114      {
     7115        Print("\n//** rec_fractal_call: Maximal recursion depth.\n");
     7116      }
     7117
     7118      to=clock();
     7119#ifdef  BUCHBERGER_ALG
     7120      Gresult = MstdhomCC(Gomega1);
     7121#else
     7122      Gresult =kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,omega);
     7123      delete hilb_func;
     7124#endif
     7125      xtstd=xtstd+clock()-to;
     7126    }
     7127    else
     7128    {
     7129      rChangeCurrRing(oRing);
     7130      Gomega1 = idrMoveR(Gomega1, oRing,currRing);
     7131      Gresult = rec_fractal_call(idCopy(Gomega1),nlev+1,omega,reduction,printout);
     7132    }
     7133    if(printout > 2)
     7134    {
     7135      idString(Gresult,"//** rec_fractal_call: M");
     7136    }
     7137    //convert a Groebner basis from a ring to another ring
     7138    new_ring = currRing;
     7139
     7140    rChangeCurrRing(oRing);
     7141    Gresult1 = idrMoveR(Gresult, new_ring,currRing);
     7142    Gomega2 = idrMoveR(Gomega1, new_ring,currRing);
     7143
     7144    to=clock();
     7145    // Lifting process
     7146    F = MLifttwoIdeal(Gomega2, Gresult1, G);
     7147    xtlift=xtlift+clock()-to;
     7148    if(printout > 2)
     7149    {
     7150      idString(F,"//** rec_fractal_call: F");
     7151    }
     7152    idDelete(&Gresult1);
     7153    idDelete(&Gomega2);
     7154    idDelete(&G);
     7155
     7156    rChangeCurrRing(new_ring);
     7157    //F1 = idrMoveR(F, oRing,currRing);
     7158    G = idrMoveR(F,oRing,currRing);
     7159    to=clock();
     7160    /* Interreduce G */
     7161   // G = kInterRedCC(F1, NULL);
     7162    xtred=xtred+clock()-to;
     7163    //idDelete(&F1);
     7164  }
     7165}
     7166
     7167/************************************************************************
     7168 * Perturb the start weight vector at the top level with random element *
     7169 ************************************************************************/
     7170static ideal rec_r_fractal_call(ideal G, int nlev, intvec* ivtarget,
     7171                int weight_rad, int reduction, int printout)
     7172{
     7173  Overflow_Error =  FALSE;
     7174  //Print("\n\n// Entering the %d-th recursion:", nlev);
     7175
     7176  int i, polylength, nV = currRing->N;
     7177  ring new_ring, testring;
     7178  //ring extoRing;
     7179  ideal Gomega, Gomega1, Gomega2, F, FF, F1, Gresult, Gresult1, G1, Gt;
     7180  int nwalks = 0;
     7181  intvec* Mwlp;
     7182#ifndef BUCHBERGER_ALG
     7183  intvec* hilb_func;
     7184#endif
     7185//  intvec* extXtau;
     7186  intvec* next_vect;
     7187  intvec* omega2 = new intvec(nV);
     7188  intvec* omtmp = new intvec(nV);
     7189  intvec* altomega = new intvec(nV);
     7190
     7191  //BOOLEAN isnewtarget = FALSE;
     7192
     7193  for(i = nV -1; i>0; i--)
     7194  {
     7195    (*omtmp)[i] = (*ivtarget)[i];
     7196  }
     7197  // to avoid (1,0,...,0) as the target vector (Hans)
     7198  intvec* last_omega = new intvec(nV);
     7199  for(i=nV-1; i>0; i--)
     7200    (*last_omega)[i] = 1;
     7201  (*last_omega)[0] = 10000;
     7202
     7203  intvec* omega = new intvec(nV);
     7204  for(i=0; i<nV; i++) {
     7205    if(Xsigma->length() == nV)
     7206      (*omega)[i] =  (*Xsigma)[i];
     7207    else
     7208      (*omega)[i] = (*Xsigma)[(nV*(nlev-1))+i];
     7209
     7210    (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i];
     7211  }
     7212
     7213   if(nlev == 1)  Xcall = 1;
     7214   else Xcall = 0;
     7215
     7216  ring oRing = currRing;
     7217
     7218  while(1)
     7219  {
     7220#ifdef FIRST_STEP_FRACTAL
     7221    /*
     7222    perturb the current weight vector only on the top level or
     7223    after perturbation of the both vectors, nlev = 2 as the top level
     7224    */
     7225    if((nlev == 1 && Xcall == 0) || (nlev == 2 && Xngleich == 1))
     7226      if(islengthpoly2(G) == 1)
     7227      {
     7228        Mwlp = MivWeightOrderlp(omega);
     7229        Xsigma = Mfpertvector(G, Mwlp);
     7230        delete Mwlp;
     7231        Overflow_Error = FALSE;
     7232      }
     7233#endif
     7234    nwalks ++;
     7235    NEXT_VECTOR_FRACTAL:
     7236    to=clock();
     7237    /* determine the next border */
     7238    next_vect = MkInterRedNextWeight(omega,omega2,G);
     7239    if(polylength > 0 && G->m[0] != NULL)
     7240    {
     7241      /*
     7242      there is a polynomial in Gomega with at least 3 monomials,
     7243      low-dimensional facet of the cone
     7244      */
     7245      delete next_vect;
     7246      next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,nlev);
     7247      if(isNegNolVector(next_vect)==1)
     7248      {
     7249        delete next_vect;
     7250        next_vect = MkInterRedNextWeight(omega,omega2,G);
     7251      }
     7252    }
     7253    xtnw=xtnw+clock()-to;
     7254
     7255    oRing = currRing;
     7256
     7257    // We only perturb the current target vector at the recursion  level 1
     7258    if(Xngleich == 0 && nlev == 1) //(ngleich == 0) important, e.g. ex2, ex3
     7259      if (MivComp(next_vect, omega2) == 1)
     7260      {
     7261        // to dispense with taking initials and lifting/interreducing
     7262        // after the call of recursion.
     7263        if(printout > 0)
     7264        {
     7265          Print("\n//** rec_r_fractal_call: Perturb both vectors with degree %d.",nlev);
     7266          //idElements(G, "G");
     7267        }
     7268        Xngleich = 1;
     7269        nlev +=1;
     7270        if(ivtarget->length() == nV)
     7271        {
     7272          if (rParameter(currRing) != NULL)
     7273            DefRingPar(omtmp);
     7274          else
     7275            rChangeCurrRing(VMrDefault(omtmp));
     7276        }
     7277        else
     7278        {
     7279          rChangeCurrRing(VMatrDefault(ivtarget));
     7280        }
     7281        testring = currRing;
     7282        Gt = idrMoveR(G, oRing,currRing);
     7283
     7284        // perturb the original target vector w.r.t. the current GB
     7285        if(ivtarget->length() == nV)
     7286        {
     7287          delete Xtau;
     7288          Xtau = NewVectorlp(Gt);
     7289        }
     7290        else
     7291        {
     7292          delete Xtau;
     7293          Xtau = Mfpertvector(Gt,ivtarget);
     7294        }
     7295
     7296        rChangeCurrRing(oRing);
     7297        G = idrMoveR(Gt,testring,currRing);
     7298
     7299        // perturb the current vector w.r.t. the current GB
     7300        Mwlp = MivWeightOrderlp(omega);
     7301        if(ivtarget->length() > nV)
     7302        {
     7303          delete Mwlp;
     7304          Mwlp = MivMatrixOrderRefine(omega,ivtarget);
     7305        }
     7306        Xsigma = Mfpertvector(G, Mwlp);
     7307        delete Mwlp;
     7308
     7309        for(i=nV-1; i>=0; i--) {
     7310          (*omega2)[i] = (*Xtau)[nV+i];
     7311          (*omega)[i] = (*Xsigma)[nV+i];
     7312        }
     7313
     7314        delete next_vect;
     7315        to=clock();
     7316
     7317        /*
     7318        to avoid the value of Overflow_Error that occur in
     7319        Mfpertvector
     7320        */
     7321        Overflow_Error = FALSE;
     7322
     7323        next_vect = MkInterRedNextWeight(omega,omega2,G);
     7324        if(G->m[0] != NULL && polylength > 0)
     7325        {
     7326          /*
     7327          there is a polynomial in Gomega with at least 3 monomials,
     7328          low-dimensional facet of the cone
     7329          */
     7330          delete next_vect;
     7331          next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,nlev);
     7332          if(isNegNolVector(next_vect)==1)
     7333          {
     7334            delete next_vect;
     7335            next_vect = MkInterRedNextWeight(omega,omega2,G);
     7336          }
     7337        }
     7338        xtnw=xtnw+clock()-to;
     7339      }
     7340//#ifdef PRINT_VECTORS
     7341      if(printout > 0)
     7342      {
     7343        MivString(omega, omega2, next_vect);
     7344      }
     7345//#endif
     7346
     7347    /* check whether the the computed vector is in the correct cone
     7348       If no, the reduced GB of an omega-homogeneous ideal will be
     7349       computed by Buchberger algorithm and stop this recursion step*/
     7350    //if(test_w_in_ConeCC(G, next_vect) != 1) //e.g. Example s7, cyc6
     7351    if(Overflow_Error == TRUE)
     7352    {
     7353      delete next_vect;
     7354      if(ivtarget->length() == nV)
     7355      {
     7356        if (rParameter(currRing) != NULL)
     7357        {
     7358          DefRingPar(omtmp);
     7359        }
     7360        else
     7361        {
     7362          rChangeCurrRing(VMrDefault(omtmp));
     7363        }
     7364      }
     7365      else
     7366      {
     7367        rChangeCurrRing(VMatrDefault(ivtarget));
     7368      }
     7369#ifdef TEST_OVERFLOW
     7370      Gt = idrMoveR(G, oRing,currRing);
     7371      Gt = NULL;
     7372      return(Gt);
     7373#endif
     7374      if(printout > 0)
     7375      {
     7376        Print("\n//** rec_r_fractal_call: applying Buchberger's algorithm in ring r = %s;",
     7377              rString(currRing));
     7378      }
     7379      to=clock();
     7380      Gt = idrMoveR(G, oRing,currRing);
     7381      G1 = MstdCC(Gt);
     7382      xtextra=xtextra+clock()-to;
     7383      Gt = NULL;
     7384
     7385      delete omega2;
     7386      delete altomega;
     7387      if(printout > 0)
     7388      {
     7389        Print("\n//** rec_r_fractal_call: Leaving the %d-th recursion with %d steps.\n",
     7390              nlev, nwalks);
     7391        //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     7392      }
     7393      nnflow ++;
     7394      Overflow_Error = FALSE;
     7395      return (G1);
     7396    }
     7397    /*
     7398       If the perturbed target vector stays in the correct cone,
     7399       return the current Groebner basis.
     7400       Otherwise, return the Groebner basis computed with Buchberger's
     7401       algorithm.
     7402       Then we update the perturbed target vectors w.r.t. this GB.
     7403    */
     7404    if (MivComp(next_vect, XivNull) == 1)
     7405    {
     7406      // The computed vector is equal to the origin vector,
     7407      // because t is not defined
     7408      if(ivtarget->length() == nV)
     7409      {
     7410        if (rParameter(currRing) != NULL)
     7411          DefRingPar(omtmp);
     7412        else
     7413          rChangeCurrRing(VMrDefault(omtmp));
     7414      }
     7415      else
     7416      {
     7417        rChangeCurrRing(VMatrDefault(ivtarget));
     7418      }
     7419      testring = currRing;
     7420      Gt = idrMoveR(G, oRing,currRing);
     7421
     7422      if(test_w_in_ConeCC(Gt, omega2) == 1)
     7423      {
    62887424        delete omega2;
    62897425        delete next_vect;
    62907426        delete altomega;
    6291         //Print("\n// Leaving the %d-th recursion with %d steps ",nlev, nwalks);
    6292         //Print(" ** Overflow_Error? (%d)", Overflow_Error);
    6293 
     7427        if(printout > 0)
     7428        {
     7429          Print("\n//** rec_r_fractal_call: Leaving the %d-th recursion with %d steps.\n",
     7430                nlev, nwalks);
     7431          //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     7432        }
    62947433        return (Gt);
    62957434      }
    62967435      else
    6297       {
    6298         //ivString(omega2, "tau'");
    6299         //Print("\n//  tau' doesn't stay in the correct cone!!");
     7436      {
     7437        if(printout > 0)
     7438        {
     7439          Print("\n//** rec_r_fractal_call: target weight doesn't stay in the correct cone.\n");
     7440        }
    63007441
    63017442#ifndef  MSTDCC_FRACTAL
    6302         //07.08.03
    63037443        //ivString(Xtau, "old Xtau");
    6304         intvec* Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp));
     7444        intvec* Xtautmp;
     7445        if(ivtarget->length() == nV)
     7446        {
     7447          Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp));
     7448        }
     7449        else
     7450        {
     7451          Xtautmp = Mfpertvector(Gt, ivtarget);
     7452        }
    63057453#ifdef TEST_OVERFLOW
    63067454      if(Overflow_Error == TRUE)
     
    63307478
    63317479      FRACTAL_MSTDCC:
    6332         //Print("\n//  apply BB-Alg in ring = %s;", rString(currRing));
     7480        if(printout > 0)
     7481        {
     7482          Print("\n//** rec_r_fractal_call: apply Buchberger's algorithm in ring = %s.\n",
     7483                rString(currRing));
     7484        }
    63337485        to=clock();
    63347486        G = MstdCC(Gt);
     
    63387490
    63397491        // update the original target vector w.r.t. the current GB
    6340         if(MivSame(Xivinput, Xivlp) == 1)
    6341           if (rParameter(currRing) != NULL)
    6342             DefRingParlp();
     7492        if(ivtarget->length() == nV)
     7493        {
     7494          if(MivSame(Xivinput, Xivlp) == 1)
     7495            if (rParameter(currRing) != NULL)
     7496              DefRingParlp();
     7497            else
     7498              VMrDefaultlp();
    63437499          else
    6344             VMrDefaultlp();
     7500            if (rParameter(currRing) != NULL)
     7501              DefRingPar(Xivinput);
     7502            else
     7503              rChangeCurrRing(VMrDefault(Xivinput));
     7504        }
    63457505        else
    6346           if (rParameter(currRing) != NULL)
    6347             DefRingPar(Xivinput);
    6348           else
    6349             rChangeCurrRing(VMrDefault1(Xivinput)); //Aenderung6
    6350 
     7506        {
     7507          rChangeCurrRing(VMatrRefine(ivtarget,Xivinput));
     7508        }
    63517509        testring = currRing;
    63527510        Gt = idrMoveR(G, oRing,currRing);
    63537511
    6354         delete Xtau;
    6355         Xtau = NewVectorlp(Gt);
     7512        // perturb the original target vector w.r.t. the current GB
     7513        if(ivtarget->length() == nV)
     7514        {
     7515          delete Xtau;
     7516          Xtau = NewVectorlp(Gt);
     7517        }
     7518        else
     7519        {
     7520          delete Xtau;
     7521          Xtau = Mfpertvector(Gt,ivtarget);
     7522        }
    63567523
    63577524        rChangeCurrRing(oRing);
     
    63617528        delete next_vect;
    63627529        delete altomega;
    6363         /*
    6364           Print("\n// Leaving the %d-th recursion with %d steps,", nlev,nwalks);
    6365           Print(" ** Overflow_Error? (%d)", Overflow_Error);
    6366         */
     7530        if(printout > 0)
     7531        {
     7532          Print("\n//** rec_r_fractal_call: Leaving the %d-th recursion with %d steps.\n",
     7533                nlev,nwalks);
     7534          //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     7535        }
    63677536        if(Overflow_Error == TRUE)
    63687537          nnflow ++;
     
    63717540        return(G);
    63727541      }
    6373     }
    6374 
    6375     for(i=nV-1; i>=0; i--) {
     7542    } //end of if(MivComp(next_vect, XivNull) == 1)
     7543
     7544    for(i=nV-1; i>=0; i--)
     7545    {
    63767546      (*altomega)[i] = (*omega)[i];
    63777547      (*omega)[i] = (*next_vect)[i];
     
    63807550
    63817551    to=clock();
    6382     /* Take the initial form of <G> w.r.t. omega */
     7552    // Take the initial form of <G> w.r.t. omega
    63837553    Gomega = MwalkInitialForm(G, omega);
    63847554    xtif=xtif+clock()-to;
     7555    //polylength = 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise
     7556    polylength = lengthpoly(Gomega);
     7557    if(printout > 1)
     7558    {
     7559      idString(Gomega,"//** rec_r_fractal_call: Gomega");
     7560    }
     7561
     7562    if(reduction == 0)
     7563    {
     7564      /* Check whether the intermediate weight vector lies in the interior of the cone.
     7565       * If so, only perform reductions. Otherwise apply Buchberger's algorithm. */
     7566      FF = middleOfCone(G,Gomega);
     7567      if( FF != NULL)
     7568      {
     7569        idDelete(&G);
     7570        G = idCopy(FF);
     7571        idDelete(&FF);
     7572        /* Compue next vector. */
     7573        goto NEXT_VECTOR_FRACTAL;
     7574      }
     7575    }
    63857576
    63867577#ifndef  BUCHBERGER_ALG
     
    63897580    else
    63907581      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
    6391 #endif // BUCHBERGER_ALG
    6392 
    6393     if (rParameter(currRing) != NULL)
    6394       DefRingPar(omega);
     7582#endif
     7583    if(ivtarget->length() == nV)
     7584    {
     7585      if (rParameter(currRing) != NULL)
     7586        DefRingPar(omega);
     7587      else
     7588        rChangeCurrRing(VMrDefault(omega));
     7589    }
    63957590    else
    6396       rChangeCurrRing(VMrDefault1(omega)); //Aenderung7
    6397 
     7591    {
     7592      rChangeCurrRing(VMatrRefine(ivtarget,omega));
     7593    }
    63987594    Gomega1 = idrMoveR(Gomega, oRing,currRing);
    63997595
    6400     /* Maximal recursion depth, to compute a red. GB */
    6401     /* Fractal walk with the alternative recursion */
    6402     /* alternative recursion */
    6403     // if(nlev == nV || lengthpoly(Gomega1) == 0)
     7596    // Maximal recursion depth, to compute a red. GB
     7597    // Fractal walk with the alternative recursion
     7598    // alternative recursion
    64047599    if(nlev == Xnlev || lengthpoly(Gomega1) == 0)
    6405       //if(nlev == nV) // blind recursion
    6406     {
    6407       /*
    6408       if(Xnlev != nV)
    6409       {
    6410         Print("\n// ** Xnlev = %d", Xnlev);
    6411         ivString(Xtau, "Xtau");
    6412       }
    6413       */
     7600    {
    64147601      to=clock();
    64157602#ifdef  BUCHBERGER_ALG
     
    64187605      Gresult =kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,omega);
    64197606      delete hilb_func;
    6420 #endif // BUCHBERGER_ALG
     7607#endif
    64217608      xtstd=xtstd+clock()-to;
    64227609    }
    6423     else {
     7610    else
     7611    {
    64247612      rChangeCurrRing(oRing);
    64257613      Gomega1 = idrMoveR(Gomega1, oRing,currRing);
    6426       Gresult = rec_fractal_call(idCopy(Gomega1),nlev+1,omega);
    6427     }
    6428 
    6429     //convert a Groebner basis from a ring to another ring,
     7614      Gresult = rec_r_fractal_call(idCopy(Gomega1),nlev+1,omega,weight_rad,reduction,printout);
     7615    }
     7616    if(printout > 2)
     7617    {
     7618      idString(Gresult,"//** rec_r_fractal_call: M");
     7619    }
     7620    //convert a Groebner basis from a ring to another ring
    64307621    new_ring = currRing;
    64317622
     
    64357626
    64367627    to=clock();
    6437     /* Lifting process */
     7628    // Lifting process
    64387629    F = MLifttwoIdeal(Gomega2, Gresult1, G);
    64397630    xtlift=xtlift+clock()-to;
     7631
     7632    if(printout > 2)
     7633    {
     7634      idString(F,"//** rec_r_fractal_call: F");
     7635    }
     7636
    64407637    idDelete(&Gresult1);
    64417638    idDelete(&Gomega2);
     
    64437640
    64447641    rChangeCurrRing(new_ring);
    6445     F1 = idrMoveR(F, oRing,currRing);
    6446 
     7642    //F1 = idrMoveR(F, oRing,currRing);
     7643    G = idrMoveR(F,oRing,currRing);
    64477644    to=clock();
    6448     /* Interreduce G */
    6449     G = kInterRedCC(F1, NULL);
     7645    // Interreduce G
     7646    //G = kInterRedCC(F1, NULL);
    64507647    xtred=xtred+clock()-to;
    6451     idDelete(&F1);
     7648    //idDelete(&F1);
    64527649  }
    64537650}
    6454 
    6455 /************************************************************************
    6456  * Perturb the start weight vector at the top level with random element *
    6457  ************************************************************************/
    6458 static ideal rec_r_fractal_call(ideal G, int nlev, intvec* omtmp, int weight_rad)
    6459 {
    6460   Overflow_Error =  FALSE;
    6461   //Print("\n\n// Entering the %d-th recursion:", nlev);
    6462 
    6463   int i, nV = currRing->N;
    6464   ring new_ring, testring;
    6465   //ring extoRing;
    6466   ideal Gomega, Gomega1, Gomega2, F, F1, Gresult, Gresult1, G1, Gt;
    6467   int nwalks = 0;
    6468   intvec* Mwlp;
    6469 #ifndef BUCHBERGER_ALG
    6470   intvec* hilb_func;
    6471 #endif
    6472 //  intvec* extXtau;
    6473   intvec* next_vect;
    6474   intvec* omega2 = new intvec(nV);
    6475   intvec* altomega = new intvec(nV);
    6476 
    6477   //BOOLEAN isnewtarget = FALSE;
    6478 
    6479   // to avoid (1,0,...,0) as the target vector (Hans)
    6480   intvec* last_omega = new intvec(nV);
    6481   for(i=nV-1; i>0; i--)
    6482     (*last_omega)[i] = 1;
    6483   (*last_omega)[0] = 10000;
    6484 
    6485   intvec* omega = new intvec(nV);
    6486   for(i=0; i<nV; i++) {
    6487     if(Xsigma->length() == nV)
    6488       (*omega)[i] =  (*Xsigma)[i];
    6489     else
    6490       (*omega)[i] = (*Xsigma)[(nV*(nlev-1))+i];
    6491 
    6492     (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i];
    6493   }
    6494 
    6495    if(nlev == 1)  Xcall = 1;
    6496    else Xcall = 0;
    6497 
    6498   ring oRing = currRing;
    6499 
    6500   while(1)
    6501   {
    6502 #ifdef FIRST_STEP_FRACTAL
    6503     // perturb the current weight vector only on the top level or
    6504     // after perturbation of the both vectors, nlev = 2 as the top level
    6505     if((nlev == 1 && Xcall == 0) || (nlev == 2 && Xngleich == 1))
    6506       if(islengthpoly2(G) == 1)
    6507       {
    6508         Mwlp = MivWeightOrderlp(omega);
    6509         Xsigma = Mfpertvector(G, Mwlp);
    6510         delete Mwlp;
    6511         Overflow_Error = FALSE;
    6512       }
    6513 #endif
    6514     nwalks ++;
    6515     NEXT_VECTOR_FRACTAL:
    6516     to=clock();
    6517     /* determine the next border */
    6518     next_vect = MWalkRandomNextWeight(G, omega,omega2, weight_rad, 1+nlev);
    6519     //next_vect = MkInterRedNextWeight(omega,omega2,G);
    6520     xtnw=xtnw+clock()-to;
    6521 #ifdef PRINT_VECTORS
    6522     MivString(omega, omega2, next_vect);
    6523 #endif
    6524     oRing = currRing;
    6525 
    6526     /* We only perturb the current target vector at the recursion  level 1 */
    6527     if(Xngleich == 0 && nlev == 1) //(ngleich == 0) important, e.g. ex2, ex3
    6528       if (MivComp(next_vect, omega2) == 1)
    6529       {
    6530         /* to dispense with taking initial (and lifting/interreducing
    6531            after the call of recursion */
    6532         //Print("\n\n// ** Perturb the both vectors with degree %d with",nlev);
    6533         //idElements(G, "G");
    6534 
    6535         Xngleich = 1;
    6536         nlev +=1;
    6537 
    6538         if (rParameter(currRing) != NULL)
    6539           DefRingPar(omtmp);
    6540         else
    6541           rChangeCurrRing(VMrDefault1(omtmp)); //Aenderung3
    6542 
    6543         testring = currRing;
    6544         Gt = idrMoveR(G, oRing,currRing);
    6545 
    6546         /* perturb the original target vector w.r.t. the current GB */
    6547         delete Xtau;
    6548         Xtau = NewVectorlp(Gt);
    6549 
    6550         rChangeCurrRing(oRing);
    6551         G = idrMoveR(Gt, testring,currRing);
    6552 
    6553         /* perturb the current vector w.r.t. the current GB */
    6554         Mwlp = MivWeightOrderlp(omega);
    6555         Xsigma = Mfpertvector(G, Mwlp);
    6556         delete Mwlp;
    6557 
    6558         for(i=nV-1; i>=0; i--) {
    6559           (*omega2)[i] = (*Xtau)[nV+i];
    6560           (*omega)[i] = (*Xsigma)[nV+i];
    6561         }
    6562 
    6563         delete next_vect;
    6564         to=clock();
    6565 
    6566         /* to avoid the value of Overflow_Error that occur in Mfpertvector*/
    6567         Overflow_Error = FALSE;
    6568 
    6569         next_vect = MkInterRedNextWeight(omega,omega2,G);
    6570         xtnw=xtnw+clock()-to;
    6571 
    6572 #ifdef PRINT_VECTORS
    6573         MivString(omega, omega2, next_vect);
    6574 #endif
    6575       }
    6576 
    6577 
    6578     /* check whether the the computed vector is in the correct cone */
    6579     /* If no, the reduced GB of an omega-homogeneous ideal will be
    6580        computed by Buchberger algorithm and stop this recursion step*/
    6581     //if(test_w_in_ConeCC(G, next_vect) != 1) //e.g. Example s7, cyc6
    6582     if(Overflow_Error == TRUE)
    6583     {
    6584       delete next_vect;
    6585       if (rParameter(currRing) != NULL)
    6586       {
    6587         DefRingPar(omtmp);
    6588       }
    6589       else
    6590       {
    6591         rChangeCurrRing(VMrDefault1(omtmp)); // Aenderung4
    6592       }
    6593 #ifdef TEST_OVERFLOW
    6594       Gt = idrMoveR(G, oRing,currRing);
    6595       Gt = NULL; return(Gt);
    6596 #endif
    6597 
    6598       //Print("\n\n// apply BB's alg. in ring r = %s;", rString(currRing));
    6599       to=clock();
    6600       Gt = idrMoveR(G, oRing,currRing);
    6601       G1 = MstdCC(Gt);
    6602       xtextra=xtextra+clock()-to;
    6603       Gt = NULL;
    6604 
    6605       delete omega2;
    6606       delete altomega;
    6607 
    6608       //Print("\n// Leaving the %d-th recursion with %d steps", nlev, nwalks);
    6609       //Print("  ** Overflow_Error? (%d)", Overflow_Error);
    6610       nnflow ++;
    6611 
    6612       Overflow_Error = FALSE;
    6613       return (G1);
    6614     }
    6615 
    6616 
    6617     /* If the perturbed target vector stays in the correct cone,
    6618        return the current GB,
    6619        otherwise, return the computed  GB by the Buchberger-algorithm.
    6620        Then we update the perturbed target vectors w.r.t. this GB. */
    6621 
    6622     /* the computed vector is equal to the origin vector, since
    6623        t is not defined */
    6624     if (MivComp(next_vect, XivNull) == 1)
    6625     {
    6626       if (rParameter(currRing) != NULL)
    6627         DefRingPar(omtmp);
    6628       else
    6629         rChangeCurrRing(VMrDefault1(omtmp)); //Aenderung5
    6630 
    6631       testring = currRing;
    6632       Gt = idrMoveR(G, oRing,currRing);
    6633 
    6634       if(test_w_in_ConeCC(Gt, omega2) == 1) {
    6635         delete omega2;
    6636         delete next_vect;
    6637         delete altomega;
    6638         //Print("\n// Leaving the %d-th recursion with %d steps ",nlev, nwalks);
    6639         //Print(" ** Overflow_Error? (%d)", Overflow_Error);
    6640 
    6641         return (Gt);
    6642       }
    6643       else
    6644       {
    6645         //ivString(omega2, "tau'");
    6646         //Print("\n//  tau' doesn't stay in the correct cone!!");
    6647 
    6648 #ifndef  MSTDCC_FRACTAL
    6649         //07.08.03
    6650         //ivString(Xtau, "old Xtau");
    6651         intvec* Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp));
    6652 #ifdef TEST_OVERFLOW
    6653       if(Overflow_Error == TRUE)
    6654       Gt = NULL; return(Gt);
    6655 #endif
    6656 
    6657         if(MivSame(Xtau, Xtautmp) == 1)
    6658         {
    6659           //PrintS("\n// Update vectors are equal to the old vectors!!");
    6660           delete Xtautmp;
    6661           goto FRACTAL_MSTDCC;
    6662         }
    6663 
    6664         Xtau = Xtautmp;
    6665         Xtautmp = NULL;
    6666         //ivString(Xtau, "new  Xtau");
    6667 
    6668         for(i=nV-1; i>=0; i--)
    6669           (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i];
    6670 
    6671         //Print("\n//  ring tau = %s;", rString(currRing));
    6672         rChangeCurrRing(oRing);
    6673         G = idrMoveR(Gt, testring,currRing);
    6674 
    6675         goto NEXT_VECTOR_FRACTAL;
    6676 #endif
    6677 
    6678       FRACTAL_MSTDCC:
    6679         //Print("\n//  apply BB-Alg in ring = %s;", rString(currRing));
    6680         to=clock();
    6681         G = MstdCC(Gt);
    6682         xtextra=xtextra+clock()-to;
    6683 
    6684         oRing = currRing;
    6685 
    6686         // update the original target vector w.r.t. the current GB
    6687         if(MivSame(Xivinput, Xivlp) == 1)
    6688           if (rParameter(currRing) != NULL)
    6689             DefRingParlp();
    6690           else
    6691             VMrDefaultlp();
    6692         else
    6693           if (rParameter(currRing) != NULL)
    6694             DefRingPar(Xivinput);
    6695           else
    6696             rChangeCurrRing(VMrDefault1(Xivinput)); //Aenderung6
    6697 
    6698         testring = currRing;
    6699         Gt = idrMoveR(G, oRing,currRing);
    6700 
    6701         delete Xtau;
    6702         Xtau = NewVectorlp(Gt);
    6703 
    6704         rChangeCurrRing(oRing);
    6705         G = idrMoveR(Gt, testring,currRing);
    6706 
    6707         delete omega2;
    6708         delete next_vect;
    6709         delete altomega;
    6710         /*
    6711           Print("\n// Leaving the %d-th recursion with %d steps,", nlev,nwalks);
    6712           Print(" ** Overflow_Error? (%d)", Overflow_Error);
    6713         */
    6714         if(Overflow_Error == TRUE)
    6715           nnflow ++;
    6716 
    6717         Overflow_Error = FALSE;
    6718         return(G);