Changeset 38301d in git


Ignore:
Timestamp:
Jul 11, 2015, 10:21:42 AM (8 years ago)
Author:
Stephan Oberfranz <oberfran@…>
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
47362dcbb6b4fbb3b9d6cded63720e333b4c97b3
Parents:
d52203afadb18a3c05d4c08a2bc18bc327690cd9e9478bb30a6f3a2d77328b3805644a035ac3e934
Message:
update

Merge branch 'spielwiese' of github.com:Singular/Sources into spielwiese
Files:
6 added
1 deleted
42 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/crypto.lib

    rd52203 r38301d  
    6969  merkle_hellman_transformation(list knapsack, int key, int mod1                generates a hard knapsack for the  Merkle-Hellman Kryptosystem for a given easy knapsack , a multiplicator key and a modulus mod1
    7070  merkle_hellman_encryption(list knapsack, list message)                    encrypts a message with the Merkle-Hellman Kryptosystem, using a hard knapsack and a message encoded as binary list
    71   merkle_hellman_decryption(list knapsack, int key, int mod1, int message)          decrypts a message with the multiplicative Merkle-Hellman Kryptosystem, using the hard knapsack, the key, the modulus mod1 and the message encoded as integer   
     71  merkle_hellman_decryption(list knapsack, int key, int mod1, int message)          decrypts a message with the multiplicative Merkle-Hellman Kryptosystem, using the hard knapsack, the key, the modulus mod1 and the message encoded as integer
    7272  super_increasing_knapsack(int ksize)                            Creates the smallest super-increasing knapsack of given size ksize
    7373  h_increasing_knapsack(int ksize, int h)                            Creates the smallest h-increasing knapsack of given size ksize and h
     
    12721272   if((k mod 2)==0)
    12731273   {
    1274       resu=ellipticMult(N,a,b,P,k/2);
     1274      resu=ellipticMult(N,a,b,P,k div 2);
    12751275      return(ellipticAdd(N,a,b,resu,resu));
    12761276   }
  • Singular/LIB/gradedModules.lib

    rd52203 r38301d  
    3131    grsum(M,N)      direct sum of two graded modules coker(M) + coker(N)
    3232    grpower(M,p)    direct p-th power of graded module coker(M)
    33     grtranspose(M)  un-ordered graded transpose of map M       
     33    grtranspose(M)  un-ordered graded transpose of map M
    3434    grgens(M)       try to compute submodule generators of coker(M)
    3535    grpres(F)       presentation of submodule generated by columns of F
     
    5151    mappingcone(M,N)   mapping cone?
    5252    grlifting3(A,B)    RND! chain lifting? probably wrong one
    53     mappingcone3(A,B)  mapping cone3? 
     53    mappingcone3(A,B)  mapping cone3?
    5454    grrange(M)         get the row-weightings
    5555    grneg(A)           graded object given by -A
    5656    matrixpres(a)      matrix presentation of direct sum of Omega^{a[i]}(i)
    5757";
    58    
     58
    5959//    grisequal(A,B)  check whether A is exactly eqal to B? TODO: isomorphic!
    6060
     
    8686
    8787  v = -v; // NOTE: due to Mathematical meanings of Singular data
    88  
     88
    8989
    9090  int lst = v[1];
    91   int cnt = 1; 
     91  int cnt = 1;
    9292
    9393  string p = R;
     
    9595
    9696  int k, d;
    97   for (k = 2; k <= n; k++ ) 
     97  for (k = 2; k <= n; k++ )
    9898  {
    9999    d = v[k];
    100100    if( d == lst ) { cnt = cnt + 1; }
    101     else 
     101    else
    102102    {
    103103      if (cnt > 1){ p = p + "^" + string(cnt); }
     
    120120  def E = grtwist(2, 0);
    121121  def v = grrange(E); // grdeg(E);
    122   grsumstr("", v ); 
     122  grsumstr("", v );
    123123}
    124124
     
    180180  {
    181181    int n = size(N); ASSUME(0, n > 0);
    182    
     182
    183183    string msg1 = "";
    184184    if( size(R) >= 2 )
     
    186186      msg1 = msg1 + "(let R:="+R+")";
    187187      R = "R"; // !!!
    188     }   
     188    }
    189189    msg1 = msg1 + ": " ;
    190190
    191    
     191
    192192
    193193    list arr; arr[n] = list();
    194194    int exact = (1==1);
    195    
     195
    196196    int i = 1;
    197    
     197
    198198    ASSUME(1, grtest(N[i]));
    199    
     199
    200200    string dst = grsumstr(R, grrange(N[i]));
    201201    string src = grsumstr(R, grdeg(N[i]));
    202    
     202
    203203    arr[i] = list(dst,  src);
    204204
    205205    i = i + 1;
    206    
     206
    207207    while( i <= n )
    208208    {
    209209      ASSUME(1, grtest(N[i]));
    210      
     210
    211211      dst = grsumstr(R, grrange(N[i]));
    212      
     212
    213213      if( exact && (src != dst) )
    214214      {
     
    218218
    219219      src = grsumstr(R, grdeg(N[i]));
    220      
     220
    221221      arr[i] = list(dst,  src);
    222222
     
    236236        msg = msg + newline + arr[i][1] + " <-- "+o+"_" + string(i) + " --";
    237237      };
    238      
     238
    239239      msg =  msg + newline + arr[n][2];
    240       msg = msg + ", given by maps: "; 
     240      msg = msg + ", given by maps: ";
    241241    } else
    242242    {
    243243//      print(arr);
    244244
    245       msg = msg + "-object collection";     
     245      msg = msg + "-object collection";
    246246      o = "o";
    247      
     247
    248248//      for( i = 1; i <= n; i++ )
    249249//      {
    250250//        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)+"]): "; 
     251//      };
     252      msg = msg + ", given by the following maps (named here as "+o+"_[1 .. "+string(n)+"]): ";
    253253    }
    254254
     
    257257
    258258    for( i = 1; i <= size(N); i++ )
    259     { 
     259    {
    260260      print( o+"_" + string(i) + " :" );
    261       grview( N[i] ); 
     261      grview( N[i] );
    262262    };
    263263
     
    268268
    269269  ASSUME(1, grtest(N) );
    270  
     270
    271271  msg = msg + " homomorphism";
    272272  if( size(R) >= 2 )
     
    280280  intvec gr = grrange(N); // grading weights?
    281281  string dst = grsumstr(R, gr);
    282  
     282
    283283  intvec G = grdeg(N);
    284284  string src = grsumstr(R, G);
    285  
     285
    286286  if( ncols(N) == 0 )
    287287  {
    288288    src = "0";
    289   } 
     289  }
    290290
    291291  lst = msg;
    292292
    293   if( (size(lst) + size(dst) + size(src) + 4) > 80 ) 
    294   { 
     293  if( (size(lst) + size(dst) + size(src) + 4) > 80 )
     294  {
    295295    if( (size(lst) + size(dst)) > 80 ) { msg = msg + newline; lst = ""; }
    296296
     
    312312//  lst = lst + ", given by ";
    313313
    314  
     314
    315315  int nc = ncols(N); int nr = nrows(N);
    316316
     
    323323  {
    324324    ASSUME(0, nc > 0);
    325    
     325
    326326    matrix M = module(N);
    327    
     327
    328328    int r,c;
    329329    int d = 1; // number of extra cols/rows for extra info around the central degree(N) block in D
     
    361361      }
    362362
    363     } else 
    364     { 
     363    } else
     364    {
    365365      msg = msg + "a matrix";
    366366    }
    367367
    368     print(msg + ", with degrees: " );   
    369     draw(D, d); // print it nicely!   
     368    print(msg + ", with degrees: " );
     369    draw(D, d); // print it nicely!
    370370  }
    371371}
     
    462462"
    463463{
    464   ASSUME(1, grtest(M) ); 
    465  
     464  ASSUME(1, grtest(M) );
     465
    466466  if ( typeof(attrib(M, "degHomog")) == "intvec" )
    467467  {
     
    469469
    470470    if( size(t) == 0 ){ return (t); } // ZERO!
    471    
     471
    472472    ASSUME(2, ncols(M) == size(t) );
    473473    return (t);
     
    629629  "Input degrees: "; grview(I);
    630630
    631   def RR = grres(I, 0, 1); list L = RR; 
     631  def RR = grres(I, 0, 1); list L = RR;
    632632
    633633  " = Non-minimal betti numbers: "; print(betti(L, 0), "betti");
     
    830830
    831831  ring r=32003,(x,y,z),dp;
    832  
     832
    833833  grview( grtwists ( intvec(-4, 1, 6 )) );
    834834
     
    843843"
    844844{
    845   ASSUME(0, a > 0); 
     845  ASSUME(0, a > 0);
    846846  def Z = grtwists( intvec(d:a) ); // will set the rank as well
    847847//  ASSUME(2, grisequal(Z, grpower( grshift(grzero(), d), a ) )); // optional check
     
    931931
    932932  def SS = grobj(S, c, dc);
    933  
     933
    934934  ASSUME(0, size( grrange(SS) ) == (size(a) + size(b)) );
    935935  ASSUME(0, size( grdeg(SS) ) == (size(da) + size(db)) );
    936936  ASSUME(0, ncols( SS ) == size( grdeg(SS) ) );
    937937  ASSUME(0, nrows( SS ) == size( grrange(SS) ) );
    938  
     938
    939939  return(SS);
    940940}
     
    984984{
    985985  ASSUME(1, grtest(M) );
    986            
    987            
     986
     987
    988988
    989989  intvec a = grrange(M);
     
    993993  {
    994994    "!! Warning: shifting '0 <- 0' leaves it as it unchanged!";
    995     return (M);   
    996   }
    997  
     995    return (M);
     996  }
     997
    998998  a = a - intvec(d:size(a));
    999999  t = t - intvec(d:size(t));
     
    10561056{
    10571057  ASSUME(0, size(w) >= nrows(A) );
    1058  
     1058
    10591059  module M = module(A);
    10601060
     
    10631063
    10641064  intvec @ww = 0:0;
    1065  
     1065
    10661066  if( size(#) > 0 )
    10671067  {
    10681068    ASSUME(0, typeof(#[1]) == "intvec" );
    1069    
     1069
    10701070    @ww = intvec( #[1] );
    10711071
     
    10771077      }
    10781078    }
    1079    
     1079
    10801080    ASSUME(0, size(@ww) == ncols(M) );
    10811081  }
     
    10911091        M = freemodule(0);
    10921092      }
    1093      
     1093
    10941094      attrib( M, "rank", size(w) );
    10951095      attrib( M, "isHomog", w );
    1096      
     1096
    10971097//      ASSUME(0, /* PROBLEM WITH ZERO COLUMNS / THEIR DEGREES! */ (ncols(M) == 0) );
    10981098    }
    10991099  }
    1100  
     1100
    11011101//  type(@ww);  type(M);
    1102  
     1102
    11031103  ASSUME(0, size(@ww) == ncols(M) ); // ?!
    1104  
     1104
    11051105  attrib(M, "degHomog", @ww);
    11061106
     
    11321132  grview( grobj( freemodule(0), intvec(1,2,3) ) );
    11331133
    1134   matrix z1[0][3]; grview( grobj( z1, 0:0, intvec(1,2,3) ) ); 
     1134  matrix z1[0][3]; grview( grobj( z1, 0:0, intvec(1,2,3) ) );
    11351135  grview( grobj( freemodule(0), 0:0, intvec(1,2,3) ) );
    11361136
    11371137  matrix z0[0][0]; grview( grobj( z0, 0:0 ) );
    11381138  grview( grobj( freemodule(0), 0:0 ) );
    1139  
    1140  
     1139
     1140
    11411141
    11421142}
     
    11451145"USAGE:  grtest(M[,b]), anyting M, optionally int b
    11461146RETURN:  1 if M is a valid graded object, 0 otherwise
    1147 PURPOSE: validate a graded object. Print an invalid object message if b is not given 
     1147PURPOSE: validate a graded object. Print an invalid object message if b is not given
    11481148NOTE: M should be an ideal or module or matrix, with weighting attribute
    11491149   'isHomog' and optionally total graded degrees attribute 'degHomog'.
     
    11691169  if ( nrows(N) != size(gr) )
    11701170  {
    1171     if(b) { "   ? grtest: Input has wrong number of rows!"; };   
     1171    if(b) { "   ? grtest: Input has wrong number of rows!"; };
    11721172    return (0);
    11731173  };
     
    11771177    return(1);
    11781178  }
    1179  
     1179
    11801180//  if( attrib(N, "rank") != size(gr) ){ return (0); } // wrong rank :(
    11811181
     
    11891189      return (1);
    11901190    }
    1191    
     1191
    11921192    if ( ncols(N) != size(T) )
    11931193    {
    1194       if(b) { "   ? grtest: Input has wrong number of cols!"; };   
     1194      if(b) { "   ? grtest: Input has wrong number of cols!"; };
    11951195      return (0);
    11961196    };
    1197    
     1197
    11981198    int k = nvars(basering) + 1; // index of mod. column in the leadexp
    11991199
     
    12591259{
    12601260  ASSUME(0, d >= 0 );
    1261  
     1261
    12621262  if( d == 0 ) { return (A); }
    1263  
     1263
    12641264  if( ncols(A) == 0 )
    12651265  {
    12661266    matrix B[nrows(A) + d][0];
    1267     return (B);   
    1268   }
    1269  
     1267    return (B);
     1268  }
     1269
    12701270  module T; T[d] = 0;
    12711271  T = T, module(transpose(A));
     
    12811281  matrix m[0][2];
    12821282  type( align(m, 3) );
    1283 } 
     1283}
    12841284
    12851285
     
    13021302  module A = grobj( module([x+y, x, 0, 0], [0, x+y, y, 0]), intvec(0,0,0,1) );
    13031303  grview(A);
    1304  
     1304
    13051305  module B = grgroebner(A);
    13061306  grview(B);
     
    13241324  ring r=32003,(x,y,z),dp;
    13251325
    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) ) ); 
     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) ) );
    13271327  grview(A);
    1328  
     1328
    13291329  grview(grsyz(A));
    1330  
    1331   module X = grgroebner( grobj( module([x]), intvec(2) ) ); 
     1330
     1331  module X = grgroebner( grobj( module([x]), intvec(2) ) );
    13321332  grview(X);
    13331333
    13341334  // syzygy module should be zero!
    13351335  grview(grsyz(X));
    1336  
    1337  
     1336
     1337
    13381338}
    13391339
     
    13511351  intvec a = grdeg(A);
    13521352  intvec b = grrange(B);
    1353  
     1353
    13541354  ASSUME(0, (size(a) == size(b)) && (a == b));  // == for intvec :(
    13551355
     
    13611361  ring r=32003,(x,y,z),dp;
    13621362
    1363   module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); 
     1363  module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) );
    13641364  grview(A);
    1365  
     1365
    13661366  A = grgroebner(A);
    13671367  grview(A);
    1368  
     1368
    13691369  module B = grsyz(A);
    13701370  grview(B);
    13711371  print(B);
    1372  
     1372
    13731373  module D = grprod( A, B );
    13741374  grview(D);
     
    13841384"USAGE:  grres(M, l[, b]), graded object M, int l, int b
    13851385RETURN:  graded resolution = list of graded objects
    1386 PURPOSE: compute graded resolution of M (of length l) and minimise it if b was given                             
     1386PURPOSE: compute graded resolution of M (of length l) and minimise it if b was given
    13871387EXAMPLE: example grres; shows an example
    13881388"
     
    13921392
    13931393  intvec v = grrange(A);
    1394  
     1394
    13951395  int b = (size(#) > 0);
    13961396  if(b) { list r = res(A, l, #[1]); } else { list r = res(A, l); }
    13971397
    13981398  l = size(r);
    1399  
     1399
    14001400  int i;
    1401  
     1401
    14021402  for ( i = 1; i <= l; i++ )
    14031403  {
     
    14111411    r[i] = grobj(r[i], v); v = grdeg(r[i]);
    14121412  }
    1413  
     1413
    14141414  i = i-1;
    14151415
     
    14211421  ring r=32003,(x,y,z),dp;
    14221422
    1423   module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); 
     1423  module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) );
    14241424  grview(A);
    1425  
     1425
    14261426  module B = grgroebner(A);
    14271427  grview(B);
     
    14381438RETURN:  graded object
    14391439PURPOSE: graded transpose of M
    1440 NOTE:    no reordering is performend by this procedure   
     1440NOTE:    no reordering is performend by this procedure
    14411441EXAMPLE: example grtranspose; shows an example
    14421442"
     
    14571457
    14581458  module K = grsyz( L ); grview(K);
    1459  
     1459
    14601460
    14611461  // Corner cases: 0 <- 0!
     
    14631463  grview( grtranspose( Z ) );
    14641464
    1465  
     1465
    14661466  // Corner cases: * <- 0
    14671467  matrix M1[3][0];
    1468  
     1468
    14691469  module Z1 = grobj( M1, intvec(-1, 0, 1) ); grview(Z1);
    14701470  grview( grtranspose( Z1 ) );
    1471  
    1472  
     1471
     1472
    14731473  // Corner cases: 0 <- *
    14741474  matrix M2[0][3];
     
    14761476  module Z2 = grobj( M2, 0:0, intvec(-1, 0, 1) ); grview(Z2);
    14771477  grview( grtranspose( Z2 ) );
    1478  
     1478
    14791479}
    14801480
     
    14821482proc grgens(def M)
    14831483"USAGE:   grgens(M), graded object M (map)
    1484 RETURN:  graded object 
     1484RETURN:  graded object
    14851485PURPOSE: try compute graded generators of coker(M) and return them as columns
    14861486         of a graded map.
     
    14921492
    14931493  module N = grtranspose( grsyz( grtranspose(M) ) );
    1494  
     1494
    14951495//  ASSUME(3, grisequal( grgroebner(M), grgroebner( grpres( N ) ) ) ); // FIXME: not always true!?
    1496  
     1496
    14971497  return ( N );
    14981498}
     
    15051505
    15061506  module N = grgens(M);
    1507  
     1507
    15081508  grview( N ); print(N); // fine == M
    15091509
    15101510
    1511   module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) ); 
     1511  module A = grobj( module([x+y, x, 0, 3], [0, x+y, y, 2], [y, y, z, 1]), intvec(0,0,0,1) );
    15121512
    15131513  A = grgroebner(A); grview(A);
     
    15181518
    15191519  grview( grgens( grzero() ) );
    1520  
     1520
    15211521}
    15221522
     
    15341534
    15351535//  ASSUME(3, grisequal( M, grgens( N ) ) );
    1536  
     1536
    15371537  return ( N );
    15381538}
     
    15451545  grview(A);
    15461546
    1547   "graded transpose: "; def B = grtranspose(A); grview( B ); print(B); 
     1547  "graded transpose: "; def B = grtranspose(A); grview( B ); print(B);
    15481548
    15491549  "... syzygy: "; def C = grsyz(B); grview(C);
     
    15531553  "... and back to presentation: "; def E = grsyz( D ); grview(E); print(E);
    15541554
    1555   def F = grgens( E ); grview(F); print(F); 
     1555  def F = grgens( E ); grview(F); print(F);
    15561556  def G = grpres( F ); grview(G); print(G);
    15571557
    15581558
    15591559  def M = grtwists( intvec(-2, 0, 4, 4) ); grview(M);
    1560  
     1560
    15611561  def N = grgens(M); grview( N ); print(N);
    1562  
    1563   def L = grpres( N ); grview( L ); print(L); 
     1562
     1563  def L = grpres( N ); grview( L ); print(L);
    15641564}
    15651565
     
    16521652  module N;
    16531653
    1654   if(i>n) 
     1654  if(i>n)
    16551655    { // no middle part
    16561656      if(a[1]>0)
     
    16591659
    16601660          if(a[n+1]>0)
    1661             { N=grsum(N,grtwist(a[n+1],-1));}       
    1662         } 
     1661            { N=grsum(N,grtwist(a[n+1],-1));}
     1662        }
    16631663      else
    16641664        { N=grtwist(a[n+1],-1);}
    1665      
    1666       return (N); // grorder(N)); 
     1665
     1666      return (N); // grorder(N));
    16671667    }
    1668   else // i <= n: middle part is present, a_i != 0 
     1668  else // i <= n: middle part is present, a_i != 0
    16691669    { // a = a1  ... |  i:2, a_2 ..... i: n, a_n | .... i: n+1a_(n+1)
    16701670      j = i - 1;
     
    16971697
    16981698
    1699       return ((N)); //      return (grorder(N)); 
     1699      return ((N)); //      return (grorder(N));
    17001700    }
    17011701}
     
    17041704  ring r = 32003,(x(0..4)),dp;
    17051705
    1706   def N1 = KeneshlouMatrixPresentation(intvec(2,0,0,0,0)); 
     1706  def N1 = KeneshlouMatrixPresentation(intvec(2,0,0,0,0));
    17071707  grview(N1);
    17081708
     
    17381738  ASSUME(1, grtest(B));
    17391739  ASSUME(0, grrange(A)==grrange(B));
    1740  
     1740
    17411741  intvec v = grrange(A);
    17421742  intvec w=grdeg(A),grdeg(B);
     
    17461746{ "EXAMPLE:"; echo = 2;
    17471747  ring r;
    1748  
     1748
    17491749  module R=grobj(module([x,y,z]),intvec(0:3));
    17501750  grview(R);
     
    17711771//  matrix U;
    17721772  matrix L =lift(A,B/*,U*/);  //  module(B*U) = module(matrix(A)*L)
    1773  
     1773
    17741774  return(grobj(L, grdeg(A), grdeg(B)));
    17751775}
     
    17821782
    17831783
    1784   module D=grobj(module([y,0,z],[x2+y2,z,0]),intvec(0,1,0)); 
     1784  module D=grobj(module([y,0,z],[x2+y2,z,0]),intvec(0,1,0));
    17851785  grview(D);
    17861786
     
    17881788  grview(G);
    17891789
    1790   ASSUME(0, grisequal( grprod(D, G), P) ); 
     1790  ASSUME(0, grisequal( grprod(D, G), P) );
    17911791}
    17921792
     
    18071807
    18081808  module Z = grobj(freemodule(0),intvec(0:0),intvec(0:0));
    1809  
     1809
    18101810  grrange(Z);
    18111811  grdeg(Z);
    1812  
     1812
    18131813  grview(Z);
    18141814
     
    18381838       grtranspose( grprod( N,  alpha1 ) ) )
    18391839       ) ); // alpha0!
    1840    
     1840
    18411841}
    18421842example
     
    18991899
    19001900  ASSUME(0, t >= 2);
    1901  
     1901
    19021902  list P;
    19031903
    19041904  "t: ", t;
    1905  
     1905
    19061906  P[1]= grrndmap( rM[1], rN[1] ); // alpha1
    19071907
    1908   if(t==2){return(P[1]);} 
    1909    
     1908  if(t==2){return(P[1]);}
     1909
    19101910  for(k=2; k<=t; k++)
    19111911  {
    1912     P[k] = grlift( grprod(P[k-1],rM[k]), rN[k] ); 
    1913      grview(P[k]); 
    1914    
    1915   }
    1916      
     1912    P[k] = grlift( grprod(P[k-1],rM[k]), rN[k] );
     1913     grview(P[k]);
     1914
     1915  }
     1916
    19171917  return(P);
    1918    
     1918
    19191919}
    19201920example
     
    19231923  ring r=32003,(x,y,z),dp;
    19241924
    1925   module P=grobj(module([xy,0,xz]),intvec(0,1,0)); 
     1925  module P=grobj(module([xy,0,xz]),intvec(0,1,0));
    19261926  grview(P);
    19271927
    1928   module D=grobj(module([y,0,z],[x2+y2,z,0]),intvec(0,1,0)); 
     1928  module D=grobj(module([y,0,z],[x2+y2,z,0]),intvec(0,1,0));
    19291929  grview(D);
    19301930
     
    19371937  module D=grobj(module([y,0,z],[x2+y2,z,0], [z3, xy, xy2]),intvec(0,1,0));
    19381938  D = grgroebner(D);
    1939   grview( grres(D, 0)); 
     1939  grview( grres(D, 0));
    19401940
    19411941  def G=grlifting(D, D);
     
    19511951  def M = KeneshlouMatrixPresentation(intvec(0,0,1,0));
    19521952  grview( grres(M, 0) );
    1953  
     1953
    19541954//   module N = grshift(kos[3], 1); // psi, Syz_2(K(1))
    19551955  def N = KeneshlouMatrixPresentation(intvec(0,1,0,0));
     
    19661966proc mappingcone(M,N)
    19671967"USAGE: mappingcone(M,N), M,N graded objects
    1968 RETURN: chain complex (as a list) 
     1968RETURN: chain complex (as a list)
    19691969PURPOSE: construct a free resolution of the cokernel of a random map between Img(M), and Img(N).
    19701970EXAMPLE: example mappingcone; shows an example
     
    20262026
    20272027// correct
    2028 proc grrndmap(def S, def D, list #) 
     2028proc grrndmap(def S, def D, list #)
    20292029"USAGE: grrndmap(S,D), graded objects S and D
    20302030RETURN: graded object
     
    20922092// 0<---M<----F0<----F1<----F2<----F3<----
    20932093//              |p1   |p2
    2094 // 
     2094//
    20952095// 0<---N<----G0<----G1<----G2<----G3<----
    20962096//                B(g1)      g2     g3
    2097 // 
     2097//
    20982098proc grlifting2(A,B)
    20992099"USAGE: grlifting2(A,B), graded objects A and B (matrices defining maps)
    2100 RETURN: map of chain complexes (as a list) 
     2100RETURN: map of chain complexes (as a list)
    21012101PURPOSE: construct a map of chain complexes between free resolution of
    21022102M=coker(A) and N=coker(B).
     
    21272127  P[1]=grrndmap2(A,B);
    21282128
    2129   // A(or B)=0 
     2129  // A(or B)=0
    21302130  if(t==1){return(P[1])};
    2131  
     2131
    21322132  for(k=2;k<=t;k++)
    21332133  {
     
    21612161def T=grlifting2(DD,PP); T;
    21622162
    2163 // def Z=grlifting2(P,D); Z; // WRONG!!! 
    2164              
     2163// def Z=grlifting2(P,D); Z; // WRONG!!!
     2164
    21652165}
    21662166
     
    21712171proc mappingcone2(A,B)
    21722172"USAGE: mappingcone2(A,B), graded objects A and B (matrices defining maps)
    2173 RETURN: chain complex (as a list) 
     2173RETURN: chain complex (as a list)
    21742174PURPOSE: construct the free resolution of a cokernel of a random map between M=coker(A), and N=coker(B)
    21752175EXAMPLE: example mappingcone2;
     
    22342234*/
    22352235
    2236 proc grlifting3(A,B) 
     2236proc grlifting3(A,B)
    22372237"TODO: grlifting4 was newer and had more documentation than this proc, but was removed... Please verify and update!
    22382238"
     
    22472247  list rN = grres(B,0,1);
    22482248  print( betti(rN), "betti");
    2249  
     2249
    22502250  int i,j,k;
    22512251
     
    22602260  }
    22612261  int t=min(i,j);
    2262  
     2262
    22632263  list P;
    22642264
    22652265  "t: ", t;
    22662266//  grview(rM[t]);  grview(rN[t]);
    2267  
     2267
    22682268  P[t]= grrndmap2(rM[t],rN[t]);
    22692269  grview(P[t]);
     
    23122312//def I=KeneshlouMatrixPresentation(intvec(2,3,0,6,2));
    23132313//def J=KeneshlouMatrixPresentation(intvec(4,0,1,2,1));
    2314 //def N=grlifting3(I,J); grview(N); 
     2314//def N=grlifting3(I,J); grview(N);
    23152315}
    23162316
     
    23182318"USAGE: grneg(A), graded object A
    23192319RETURN: graded object
    2320 PURPOSE: graded map defined by -A 
     2320PURPOSE: graded map defined by -A
    23212321EXAMPLE: example grneg; shows an example
    23222322"
     
    23412341proc mappingcone3(A,B)
    23422342"USAGE: mappingcone3(A,B), graded objects A and B (matrices defining maps)
    2343 RETURN: chain complex (as a list) 
     2343RETURN: chain complex (as a list)
    23442344PURPOSE: construct a free resolution of the cokernel of a random map between M=coker(A), and N=coker(B)
    23452345EXAMPLE: example mappingcone3; shows an example
     
    23762376
    23772377    T[i]=grtranspose(D);
    2378    
     2378
    23792379    kill A, B, C, D, v, w, r, s, zero;
    23802380  }
     
    23942394def F=grlifting3(A,T); grview(F);
    23952395
    2396 // BUG in the proc 
     2396// BUG in the proc
    23972397def G=mappingcone3(A,T); grview(G);
    23982398
     
    24032403ideal P=groebner(flatten(U[2]));
    24042404resolution L=mres(P,0);
    2405 print(betti(L),"betti");   
     2405print(betti(L),"betti");
    24062406*/
    24072407
     
    24212421def I=KeneshlouMatrixPresentation(intvec(2,3,0,6,2));
    24222422def J=KeneshlouMatrixPresentation(intvec(4,0,1,2,1));
    2423 // def N=grlifting3(I,J); 
     2423// def N=grlifting3(I,J);
    24242424// 2nd module does not lie in the first:
    24252425// def NN=mappingcone3(I,J); // ????????
     
    24582458  module N;
    24592459
    2460   if(i>n) 
     2460  if(i>n)
    24612461    { // no middle part
    24622462      if(a[1]>0)
     
    24652465
    24662466          if(a[n+1]>0)
    2467             { N=grsum(N,grtwist(a[n+1],0));}         
    2468         } 
     2467            { N=grsum(N,grtwist(a[n+1],0));}
     2468        }
    24692469      else
    24702470        { N=grtwist(a[n+1],0);}
    2471      
    2472       return (N); // grorder(N)); 
     2471
     2472      return (N); // grorder(N));
    24732473    }
    24742474
    2475 else // i <= n: middle part is present, a_i != 0 
     2475else // i <= n: middle part is present, a_i != 0
    24762476    { // a = a1  ... |  i:2, a_2 ..... i: n, a_n | .... i: n+1a_(n+1)
    24772477      module I = maxideal(1);
     
    25052505
    25062506
    2507       return ((N)); //      return (grorder(N)); 
     2507      return ((N)); //      return (grorder(N));
    25082508    }
    25092509}
     
    25182518grview(S);
    25192519
    2520 def N1 = matrixpres(intvec(2,0,0,0,0)); 
     2520def N1 = matrixpres(intvec(2,0,0,0,0));
    25212521grview(N1);
    25222522
  • Singular/LIB/grwalk.lib

    re9478b r38301d  
    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
    re9478b r38301d  
    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/nfmodstd.lib

    rd52203 r38301d  
    524524    I;
    525525}
     526
     527//////////////////////////////////////////////////////////////////////////////
     528
     529/*
     530Benchmark Problems from
     531
     532Boku, Decker, Fieker, Steenpass: Groebner Bases over Algebraic Number
     533Fields.
     534
     535// 1
     536ring R = (0,a), (x,y,z), dp;
     537minpoly = (a^2+1);
     538poly f1 = (a+8)*x^2*y^2+5*x*y^3+(-a+3)*x^3*z
     539          +x^2*y*z;
     540poly f2 = x^5+2*y^3*z^2+13*y^2*z^3+5*y*z^4;
     541poly f3 = 8*x^3+(a+12)*y^3+x*z^2+3;
     542poly f4 = (-a+7)*x^2*y^4+y^3*z^3+18*y^3*z^2;
     543ideal I1 = f1,f2,f3,f4;
     544
     545// 2
     546ring R = (0,a), (x,y,z), dp;
     547minpoly = (a^5+a^2+2);
     548poly f1 = 2*x*y^4*z^2+(a-1)*x^2*y^3*z
     549          +(2*a)*x*y*z^2+7*y^3+(7*a+1);
     550poly f2 = 2*x^2*y^4*z+(a)*x^2*y*z^2-x*y^2*z^2
     551          +(2*a+3)*x^2*y*z-12*x+(12*a)*y;
     552poly f3 = (2*a)*y^5*z+x^2*y^2*z-x*y^3*z
     553          +(-a)*x*y^3+y^4+2*y^2*z;
     554poly f4 = (3*a)*x*y^4*z^3+(a+1)*x^2*y^2*z
     555          -x*y^3*z+4*y^3*z^2+(3*a)*x*y*z^3
     556          +4*z^2-x+(a)*y;
     557ideal I2 = f1,f2,f3,f4;
     558
     559// 3a
     560ring R = (0,a), (v,w,x,y,z), dp;
     561minpoly = (a^7-7*a+3);
     562poly f1 = (a)*v+(a-1)*w+x+(a+2)*y+z;
     563poly f2 = v*w+(a-1)*w*x+(a+2)*v*y+x*y+(a)*y*z;
     564poly f3 = (a)*v*w*x+(a+5)*w*x*y+(a)*v*w*z
     565          +(a+2)*v*y*z+(a)*x*y*z;
     566poly f4 = (a-11)*v*w*x*y+(a+5)*v*w*x*z
     567          +(a)*v*w*y*z+(a)*v*x*y*z
     568          +(a)*w*x*y*z;
     569poly f5 = (a+3)*v*w*x*y*z+(a+23);
     570ideal I3a = f1,f2,f3,f4,f5;
     571
     572// 3b
     573ring R = (0,a), (u,v,w,x,y,z), dp;
     574minpoly = (a^7-7*a+3);
     575poly f1 = (a)*u+(a+2)*v+w+x+y+z;
     576poly f2 = u*v+v*w+w*x+x*y+(a+3)*u*z+y*z;
     577poly f3 = u*v*w+v*w*x+(a+1)*w*x*y+u*v*z+u*y*z
     578          +x*y*z;
     579poly f4 = (a-1)*u*v*w*x+v*w*x*y+u*v*w*z
     580          +u*v*y*z+u*x*y*z+w*x*y*z;
     581poly f5 = u*v*w*x*y+(a+1)*u*v*w*x*z+u*v*w*y*z
     582          +u*v*x*y*z+u*w*x*y*z+v*w*x*y*z;
     583poly f6 = u*v*w*x*y*z+(-a+2);
     584ideal I3b = f1,f2,f3,f4,f5,f6;
     585
     586// 4
     587ring R = (0,a), (w,x,y,z), dp;
     588minpoly = (a^6+a^5+a^4+a^3+a^2+a+1);
     589poly f1 = (a+5)*w^3*x^2*y+(a-3)*w^2*x^3*y
     590          +(a+7)*w*x^2*y^2;
     591poly f2 = (a)*w^5+(a+3)*w*x^2*y^2
     592          +(a^2+11)*x^2*y^2*z;
     593poly f3 = (a+7)*w^3+12*x^3+4*w*x*y+(a)*z^3;
     594poly f4 = 3*w^3+(a-4)*x^3+x*y^2;
     595ideal I4 = f1,f2,f3,f4;
     596
     597// 5
     598ring R = (0,a), (w,x,y,z), dp;
     599minpoly = (a^12-5*a^11+24*a^10-115*a^9+551*a^8
     600          -2640*a^7+12649*a^6-2640*a^5+551*a^4
     601          -115*a^3+24*a^2-5*a+1);
     602poly f1 = (2*a+3)*w*x^4*y^2+(a+1)*w^2*x^3*y*z
     603          +2*w*x*y^2*z^3+(7*a-1)*x^3*z^4;
     604poly f2 = 2*w^2*x^4*y+w^2*x*y^2*z^2
     605          +(-a)*w*x^2*y^2*z^2
     606          +(a+11)*w^2*x*y*z^3-12*w*z^6
     607          +12*x*z^6;
     608poly f3 = 2*x^5*y+w^2*x^2*y*z-w*x^3*y*z
     609          -w*x^3*z^2+(a)*x^4*z^2+2*x^2*y*z^3;
     610poly f4 = 3*w*x^4*y^3+w^2*x^2*y*z^3
     611          -w*x^3*y*z^3+(a+4)*x^3*y^2*z^3
     612          +3*w*x*y^3*z^3+(4*a)*y^2*z^6-w*z^7
     613          +x*z^7;
     614ideal I5 = f1,f2,f3,f4;
     615
     616// 6
     617ring R = (0,a), (u,v,w,x,y,z), dp;
     618minpoly = (a^2+5*a+1);
     619poly f1 = u+v+w+x+y+z+(a);
     620poly f2 = u*v+v*w+w*x+x*y+y*z+(a)*u+(a)*z;
     621poly f3 = u*v*w+v*w*x+w*x*y+x*y*z+(a)*u*v
     622          +(a)*u*z+(a)*y*z;
     623poly f4 = u*v*w*x+v*w*x*y+w*x*y*z+(a)*u*v*w
     624          +(a)*u*v*z+(a)*u*y*z+(a)*x*y*z;
     625poly f5 = u*v*w*x*y+v*w*x*y*z+(a)*u*v*w*x
     626          +(a)*u*v*w*z+(a)*u*v*y*z+(a)*u*x*y*z
     627          +(a)*w*x*y*z;
     628poly f6 = u*v*w*x*y*z+(a)*u*v*w*x*y
     629          +(a)*u*v*w*x*z+(a)*u*v*w*y*z
     630          +(a)*u*v*x*y*z+(a)*u*w*x*y*z
     631          +(a)*v*w*x*y*z;
     632poly f7 = (a)*u*v*w*x*y*z-1;
     633ideal I6 = f1,f2,f3,f4,f5,f6,f7;
     634
     635// 7
     636ring R = (0,a), (w,x,y,z), dp;
     637minpoly = (a^8-16*a^7+19*a^6-a^5-5*a^4+13*a^3
     638          -9*a^2+13*a+17);
     639poly f1 = (-a^2-1)*x^2*y+2*w*x*z-2*w
     640          +(a^2+1)*y;
     641poly f2 = (a^3-a-3)*w^3*y+4*w*x^2*y+4*w^2*x*z
     642          +2*x^3*z+(a)*w^2-10*x^2+4*w*y-10*x*z
     643          +(2*a^2+a);
     644poly f3 = (a^2+a+11)*x*y*z+w*z^2-w-2*y;
     645poly f4 = -w*y^3+4*x*y^2*z+4*w*y*z^2+2*x*z^3
     646          +(2*a^3+a^2)*w*y+4*y^2-10*x*z-10*z^2
     647          +(3*a^2+5);
     648ideal I7 = f1,f2,f3,f4;
     649
     650// 8
     651ring R = (0,a), (t,u,v,w,x,y,z), dp;
     652minpoly = (a^7+10*a^5+5*a^3+10*a+1);
     653poly f1 = v*x+w*y-x*z-w-y;
     654poly f2 = v*w-u*x+x*y-w*z+v+x+z;
     655poly f3 = t*w-w^2+x^2-t;
     656poly f4 = (-a)*v^2-u*y+y^2-v*z-z^2+u;
     657poly f5 = t*v+v*w+(-a^2-a-5)*x*y-t*z+w*z+v+x+z
     658          +(a+1);
     659poly f6 = t*u+u*w+(-a-11)*v*x-t*y+w*y-x*z-t-u
     660          +w+y;
     661poly f7 = w^2*y^3-w*x*y^3+x^2*y^3+w^2*y^2*z
     662          -w*x*y^2*z+x^2*y^2*z+w^2*y*z^2
     663          -w*x*y*z^2+x^2*y*z^2+w^2*z^3-w*x*z^3
     664          +x^2*z^3;
     665poly f8 = t^2*u^3+t^2*u^2*v+t^2*u*v^2+t^2*v^3
     666          -t*u^3*x-t*u^2*v*x-t*u*v^2*x-t*v^3*x
     667          +u^3*x^2+u^2*v*x^2+u*v^2*x^2
     668          +v^3*x^2;
     669ideal I8 = f1,f2,f3,f4,f5,f6,f7,f8;
     670*/
  • Singular/LIB/primdec.lib

    rd52203 r38301d  
    807807        if((size(sact)==1)&&(m==2))
    808808        {
    809           l[2*i]=l[2*i-1];
    810           attrib(l[2*i],"isSB",1);
     809          l[2*i]=std(l[2*i-1],sact[1]);
    811810        }
    812811        if((size(sact)==1)&&(m>2))
     
    875874  int uytrewq;
    876875  int nva = nvars(basering);
    877   int @k,@s,@n,@k1,zz;
     876  int @k,@s,@n,@k1,@zz;
    878877  list primary,lres0,lres1,act,@lh,@wh;
    879878  map phi,psi,phi1,psi1;
     
    10561055  {
    10571056    poly randp;
    1058     for(zz=1;zz<nvars(basering);zz++)
     1057    for(@zz=1;@zz<nvars(basering);@zz++)
    10591058    {
    10601059      randp=randp
    1061               +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(zz);
     1060              +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(@zz);
    10621061    }
    10631062    randp=randp+var(nvars(basering));
     
    10691068    if (size(primary[2*@k])==0)
    10701069    {
    1071       for(zz=1;zz<size(primary[2*@k-1])-1;zz++)
     1070      for(@zz=1;@zz<size(primary[2*@k-1])-1;@zz++)
    10721071      {
    10731072        attrib(primary[2*@k-1],"isSB",1);
    1074         if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][zz]))
     1073        if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][@zz]))
    10751074        {
    10761075          primary[2*@k]=primary[2*@k-1];
     
    11001099        if(deg(lead(primary[2*@k-1][@n]))==1)
    11011100        {
    1102           for(zz=1;zz<=nva;zz++)
     1101          for(@zz=1;@zz<=nva;@zz++)
    11031102          {
    1104             if(lead(primary[2*@k-1][@n])/var(zz)!=0)
     1103            if(lead(primary[2*@k-1][@n])/var(@zz)!=0)
    11051104            {
    1106               jmap1[zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]
     1105              jmap1[@zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]
    11071106                   +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]);
    1108               jmap2[zz]=primary[2*@k-1][@n];
    1109               @qht[@n]=var(zz);
     1107              jmap2[@zz]=primary[2*@k-1][@n];
     1108              @qht[@n]=var(@zz);
    11101109            }
    11111110          }
     
    90249023  int uytrewq;
    90259024  int nva = nvars(basering);
    9026   int @k,@s,@n,@k1,zz;
     9025  int @k,@s,@n,@k1,@zz;
    90279026  list primary,lres0,lres1,act,@lh,@wh;
    90289027  map phi,psi,phi1,psi1;
     
    92079206  {
    92089207    poly randp;
    9209     for(zz=1;zz<nvars(basering);zz++)
     9208    for(@zz=1;@zz<nvars(basering);@zz++)
    92109209    {
    92119210      randp=randp
    9212               +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(zz);
     9211              +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(@zz);
    92139212    }
    92149213    randp=randp+var(nvars(basering));
     
    92209219    if (size(primary[2*@k])==0)
    92219220    {
    9222       for(zz=1;zz<size(primary[2*@k-1])-1;zz++)
     9221      for(@zz=1;@zz<size(primary[2*@k-1])-1;@zz++)
    92239222      {
    92249223        attrib(primary[2*@k-1],"isSB",1);
    9225         if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][zz]))
     9224        if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][@zz]))
    92269225        {
    92279226          primary[2*@k]=primary[2*@k-1];
     
    92519250        if(deg(lead(primary[2*@k-1][@n]))==1)
    92529251        {
    9253           for(zz=1;zz<=nva;zz++)
     9252          for(@zz=1;@zz<=nva;@zz++)
    92549253          {
    9255             if(lead(primary[2*@k-1][@n])/var(zz)!=0)
     9254            if(lead(primary[2*@k-1][@n])/var(@zz)!=0)
    92569255            {
    9257               jmap1[zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]
     9256              jmap1[@zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]
    92589257                   +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]);
    9259               jmap2[zz]=primary[2*@k-1][@n];
    9260               @qht[@n]=var(zz);
     9258              jmap2[@zz]=primary[2*@k-1][@n];
     9259              @qht[@n]=var(@zz);
    92619260            }
    92629261          }
  • Singular/LIB/rwalk.lib

    • Property mode changed from 100644 to 100755
    re9478b r38301d  
    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/surf.lib

    rd52203 r38301d  
    274274
    275275  string surf_call; i = 0;
    276  
     276
    277277  if (isWindows())
    278278  {
     
    303303  else
    304304  {
    305     surf_call = "surfer";   
     305    surf_call = "surfer";
    306306    surf_call = surf_call + " " + l + " >/dev/null 2>&1";
    307307    "Close window to exit from `surfer`.";
    308308    i = system("sh", surf_call);
    309    
    310     if ( (i != 0) && isMacOSX() ) 
     309
     310    if ( (i != 0) && isMacOSX() )
    311311    {
    312312      "*!* Sorry: calling `surfer` failed ['"+surf_call+"']" + newline
    313313      + " (The shell returned the error code " + string(i) + "." + newline
    314       + "But since we are on Mac OS X, let us try to open SURFER.app instead..." + newline 
     314      + "But since we are on Mac OS X, let us try to open SURFER.app instead..." + newline
    315315      + "Appropriate SURFER.app is available for instance at http://www.mathematik.uni-kl.de/~motsak/files/SURFER.dmg";
    316      
     316
    317317      // fallback, will only work if SURFER.app is available (e.g. in /Applications)
    318318      // get SURFER.app e.g. from http://www.mathematik.uni-kl.de/~motsak/files/SURFER.dmg
    319319      // note that the newer (Java-based) variant of Surfer may not support command line usage yet :(
    320      
     320
    321321      surf_call = "open -a SURFER -W --args -t -s";
    322322      surf_call = surf_call + " " + l + " >/dev/null 2>&1";
     
    324324      i = system("sh", surf_call);
    325325    }
    326    
    327    
    328    
     326
     327
     328
    329329  }
    330330  system("sh", "/bin/rm " + l);
     
    381381{
    382382  string s = system("uname");
    383  
     383
    384384  for (int i = 1; i <= size(s)-2; i = i + 1)
    385385  {
  • Singular/LIB/swalk.lib

    • Property mode changed from 100644 to 100755
  • Singular/Makefile.am

    rd52203 r38301d  
    9898   mmalloc.h \
    9999   mod_lib.h \
    100    omSingularConfig.h \
    101100   links/ndbm.h \
    102101   newstruct.h \
  • Singular/dyn_modules/singmathic/singmathic.cc

    rd52203 r38301d  
    468468  result->rtyp=NONE;
    469469
    470   if (arg == NULL || arg->next != NULL || 
     470  if (arg == NULL || arg->next != NULL ||
    471471  ((arg->Typ() != IDEAL_CMD) &&(arg->Typ() != MODUL_CMD))){
    472472    WerrorS("Syntax: mathicgb(<ideal>/<module>)");
  • Singular/dyn_modules/syzextra/mod_main.cc

    rd52203 r38301d  
    8282    for (int j=0; j<l; j++)
    8383      if (id->m[j] != NULL && p_GetComp(id->m[j], r) > 0)
    84         return TRUE;   
     84        return TRUE;
    8585
    8686    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!
     
    9191
    9292
    93    
     93
    9494
    9595static inline void NoReturn(leftv& res)
     
    15501550  const int iLimit = r->typ[pos].data.is.limit;
    15511551  const ideal F = r->typ[pos].data.is.F;
    1552  
     1552
    15531553  ideal FF = id_Copy(F, r);
    15541554
  • Singular/extra.cc

    rd52203 r38301d  
    37943794    }
    37953795    else
    3796 #endif   
     3796#endif
    37973797/*==================== Error =================*/
    37983798      Werror( "(extended) system(\"%s\",...) %s", sys_cmd, feNotImplemented );
  • Singular/maps_ip.cc

    rd52203 r38301d  
    7373      if (P!=0)
    7474      {
    75 //        WerrorS("Sorry 'napPermNumber' was lost in the refactoring process (due to Frank): needs to be fixed");
    76 //        return TRUE;
    77 #if 1
    7875// poly n_PermNumber(const number z, const int *par_perm, const int OldPar, const ring src, const ring dst);
    7976        res->data= (void *) n_PermNumber((number)data, par_perm, P, preimage_r, currRing);
    80 #endif
    8177        res->rtyp=POLY_CMD;
    8278        if (nCoeff_is_algExt(currRing->cf))
     
    8783      {
    8884        assume( nMap != NULL );
    89 
    9085        number a = nMap((number)data, preimage_r->cf, currRing->cf);
    9186        if (nCoeff_is_Extension(currRing->cf))
     
    150145        }
    151146      }
    152       else
    153       if ( (what==IMAP_CMD) || /*(*/ (what==FETCH_CMD) /*)*/) /* && (nMap!=nCopy)*/
     147      else if ((what==IMAP_CMD) || (what==FETCH_CMD))
    154148      {
    155149        for (i=R*C-1;i>=0;i--)
     
    160154        }
    161155      }
    162       else /* if(what==MAP_CMD) */
    163       {
     156      else /* (what==MAP_CMD) */
     157      {
     158        assume(what==MAP_CMD);
    164159        matrix s=mpNew(N,maMaxDeg_Ma((ideal)data,preimage_r));
    165160        for (i=R*C-1;i>=0;i--)
     
    170165        idDelete((ideal *)&s);
    171166      }
    172       if (nCoeff_is_Extension(currRing->cf))
     167      if (nCoeff_is_algExt(currRing->cf))
    173168      {
    174169        for (i=R*C-1;i>=0;i--)
  • Singular/misc_ip.cc

    rd52203 r38301d  
    716716  } while (v!=NULL);
    717717
    718 #ifdef OM_SINGULAR_CONFIG_H
    719718   // set global variable to show memory usage
    720719  extern int om_sing_opt_show_mem;
    721720  if (BVERBOSE(V_SHOW_MEM)) om_sing_opt_show_mem = 1;
    722721  else om_sing_opt_show_mem = 0;
    723 #endif
    724722
    725723  return FALSE;
     
    10951093#endif   // HAVE_SIMPLEIPC
    10961094    fe_reset_input_mode();
     1095    monitor(NULL,0);
    10971096#ifdef PAGE_TEST
    10981097    mmEndStat();
  • Singular/test.cc

    rd52203 r38301d  
    194194#include <Singular/links/ndbm.h>
    195195#include <Singular/newstruct.h>
    196 #include <Singular/omSingularConfig.h>
    197196#include <Singular/pcv.h>
    198197#include <Singular/links/pipeLink.h>
  • Singular/walk.cc

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