Changeset 3e1c75 in git


Ignore:
Timestamp:
Jun 17, 2015, 2:09:30 PM (8 years ago)
Author:
Stephan Oberfranz <oberfran@…>
Branches:
(u'spielwiese', 'a7324b6e0b44a1a8ed3fa4d9ca3e2ff210ddd52c')
Children:
35aef31280be0abfb1b00fad6f744c7141e32d12
Parents:
269fc4d27717409c0438dbf31005bd82968fc09d66bb58083828927467d2480a460e91aec1dc789d
Message:
update

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

Legend:

Unmodified
Added
Removed
  • Singular/LIB/grwalk.lib

    r66bb58 r3e1c75  
    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
    r66bb58 r3e1c75  
    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/poly.lib

    r269fc4 r3e1c75  
    895895{
    896896  if (f==0) { return(number(1)); }
    897   return(leadcoef(f)/leadcoef(cleardenom(f)));
     897  if (attrib(basering,"ring_cf")
     898  && ((typeof(f)=="poly")||(typeof(f)=="vector")))
     899  {
     900    number c=leadcoef(f);
     901    while((c!=1)&&(f!=0))
     902    {
     903      c=gcd(c,leadcoef(f));
     904      f=f-lead(f);
     905    }
     906    return(c);
     907  }
     908  else
     909  {
     910    return(leadcoef(f)/leadcoef(cleardenom(f)));
     911  }
    898912}
    899913example
  • Singular/LIB/rwalk.lib

    • Property mode changed from 100644 to 100755
    r66bb58 r3e1c75  
    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/standard.lib

    r269fc4 r3e1c75  
    954954
    955955 //------------------ classify the possible settings ---------------------
    956   string algorithm;       //possibilities: std, slimgb, stdorslimgb
     956  string algorithm;       //possibilities: std, slimgb, stdorslimgb, mathicgb
    957957  string conversion;      //possibilities: hilb, fglm, hilborfglm, no
    958958  string partovar;        //possibilities: yes, no
     
    961961
    962962  //define algorithm:
     963  if( (was_minpoly == 0) && (npars_P == 0) && (was_qring == 0) && (attrib (P,"global") == 1) && (char(P) > 0) && (size(BRlist)<=4) )
     964  {
     965    if( defined(Singmathic) )
     966    {
     967      algorithm = "mathicgb"; // make it default for any appropriate setting... if mathicgb is available...
     968    } else
     969    {
     970      if( p_opt && find(method,"mathicgb")  ) { "Sorry Singmathic::mathicgb is not available!"; }
     971    }
     972  }
    963973  if( find(method,"std") && !find(method,"slimgb") )
    964974  {
     
    10181028        (order=="simple" && (method==",par2var" && npars_P==0 )) ||
    10191029         (conversion=="no" && partovar=="no" &&
    1020            (algorithm=="std" || algorithm=="slimgb" ||
    1021             (find(method,"std") && find(method,"slimgb")) ) ) )
     1030           (algorithm=="std" || algorithm=="slimgb" || algorithm=="mathicgb" ||
     1031            (find(method,"std") && find(method,"slimgb"))
     1032           )
     1033         )
     1034     )
    10221035  {
    10231036    direct = "yes";
     
    10401053  //BRlist (=ringlist of basering) > 4 if the basering is non-commutative
    10411054//---------------------------- direct methods -----------------------------
     1055  if ( algorithm=="mathicgb" )
     1056  {
     1057    if (p_opt) { algorithm + " in " + string(P); }
     1058    return( mathicgb(i) );
     1059  }
    10421060  if ( direct == "yes" )
    10431061  {
  • Singular/LIB/swalk.lib

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

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

    r269fc4 r3e1c75  
    8080#endif
    8181        res->rtyp=POLY_CMD;
    82         if (nCoeff_is_Extension(currRing->cf))
     82        if (nCoeff_is_algExt(currRing->cf))
    8383          res->data=(void *)p_MinPolyNormalize((poly)res->data, currRing);
    8484        pTest((poly) res->data);
     
    8989
    9090        number a = nMap((number)data, preimage_r->cf, currRing->cf);
    91 
    92 
    9391        if (nCoeff_is_Extension(currRing->cf))
    9492        {
    95           n_Normalize(a, currRing->cf); // ???
     93          n_Normalize(a, currRing->cf);
    9694/*
    9795          number a = (number)res->data;
     
    270268        n_Delete(&d, currRing); d = NULL;
    271269      }
    272        else if (!p_IsConstant((poly)NUM((fraction)d), R))
     270      else if (!p_IsConstant((poly)NUM((fraction)d), R))
    273271      {
    274272        WarnS("ignoring denominators of coefficients...");
  • Singular/walk.cc

    • Property mode changed from 100644 to 100755
    <
    r66bb58 r3e1c75  
    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);
     7203
     7204    rChangeCurrRing(new_ring);
     7205    //F1 = idrMoveR(F, oRing,currRing);
     7206    G = idrMoveR(F,oRing,currRing);
     7207    to=clock();
     7208    // Interreduce G
     7209    // G = kInterRedCC(F1, NULL);
     7210    xtred=xtred+clock()-to;
     7211    //idDelete(&F1);
     7212  }
     7213}
     7214
     7215/************************************************************************
     7216 * Perturb the start weight vector at the top level with random element *
     7217 ************************************************************************/
     7218static ideal rec_r_fractal_call(ideal G, int nlev, intvec* ivtarget,
     7219                int weight_rad, int reduction, int printout)
     7220{
     7221  Overflow_Error =  FALSE;
     7222  //Print("\n\n// Entering the %d-th recursion:", nlev);
     7223
     7224  int i, polylength, nV = currRing->N;
     7225  ring new_ring, testring;
     7226  //ring extoRing;
     7227  ideal Gomega, Gomega1, Gomega2, F, FF, F1, Gresult, Gresult1, G1, Gt;
     7228  int nwalks = 0;
     7229  intvec* Mwlp;
     7230#ifndef BUCHBERGER_ALG
     7231  intvec* hilb_func;
     7232#endif
     7233//  intvec* extXtau;
     7234  intvec* next_vect;
     7235  intvec* omega2 = new intvec(nV);
     7236  intvec* omtmp = new intvec(nV);
     7237  intvec* altomega = new intvec(nV);
     7238
     7239  //BOOLEAN isnewtarget = FALSE;
     7240
     7241  for(i = nV -1; i>0; i--)
     7242  {
     7243    (*omtmp)[i] = (*ivtarget)[i];
     7244  }
     7245  // to avoid (1,0,...,0) as the target vector (Hans)
     7246  intvec* last_omega = new intvec(nV);
     7247  for(i=nV-1; i>0; i--)
     7248    (*last_omega)[i] = 1;
     7249  (*last_omega)[0] = 10000;
     7250
     7251  intvec* omega = new intvec(nV);
     7252  for(i=0; i<nV; i++) {
     7253    if(Xsigma->length() == nV)
     7254      (*omega)[i] =  (*Xsigma)[i];
     7255    else
     7256      (*omega)[i] = (*Xsigma)[(nV*(nlev-1))+i];
     7257
     7258    (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i];
     7259  }
     7260
     7261   if(nlev == 1)  Xcall = 1;
     7262   else Xcall = 0;
     7263
     7264  ring oRing = currRing;
     7265
     7266  while(1)
     7267  {
     7268#ifdef FIRST_STEP_FRACTAL
     7269    /*
     7270    perturb the current weight vector only on the top level or
     7271    after perturbation of the both vectors, nlev = 2 as the top level
     7272    */
     7273    if((nlev == 1 && Xcall == 0) || (nlev == 2 && Xngleich == 1))
     7274      if(islengthpoly2(G) == 1)
     7275      {
     7276        Mwlp = MivWeightOrderlp(omega);
     7277        Xsigma = Mfpertvector(G, Mwlp);
     7278        delete Mwlp;
     7279        Overflow_Error = FALSE;
     7280      }
     7281#endif
     7282    nwalks ++;
     7283    NEXT_VECTOR_FRACTAL:
     7284    to=clock();
     7285    /* determine the next border */
     7286    next_vect = MkInterRedNextWeight(omega,omega2,G);
     7287    if(polylength > 0 && G->m[0] != NULL)
     7288    {
     7289     
     7290      PrintS("\n**// rec_r_fractal_call: there is a polynomial in Gomega with at least 3 monomials, low-dimensional facet of the cone.\n");
     7291     
     7292      delete next_vect;
     7293      next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,nlev);
     7294      if(isNegNolVector(next_vect)==1)
     7295      {
     7296        delete next_vect;
     7297        next_vect = MkInterRedNextWeight(omega,omega2,G);
     7298      }
     7299    }
     7300    xtnw=xtnw+clock()-to;
     7301
     7302    oRing = currRing;
     7303
     7304    // We only perturb the current target vector at the recursion  level 1
     7305    if(Xngleich == 0 && nlev == 1) //(ngleich == 0) important, e.g. ex2, ex3
     7306      if (MivComp(next_vect, omega2) == 1)
     7307      {
     7308        // to dispense with taking initials and lifting/interreducing
     7309        // after the call of recursion.
     7310        if(printout > 0)
     7311        {
     7312          Print("\n//** rec_r_fractal_call: Perturb both vectors with degree %d.",nlev);
     7313          //idElements(G, "G");
     7314        }
     7315        Xngleich = 1;
     7316        nlev +=1;
     7317        if(ivtarget->length() == nV)
     7318        {
     7319          if (rParameter(currRing) != NULL)
     7320            DefRingPar(omtmp);
     7321          else
     7322            rChangeCurrRing(VMrDefault(omtmp));
     7323        }
     7324        else
     7325        {
     7326          rChangeCurrRing(VMatrDefault(ivtarget));
     7327        }
     7328        testring = currRing;
     7329        Gt = idrMoveR(G, oRing,currRing);
     7330
     7331        // perturb the original target vector w.r.t. the current GB
     7332        if(ivtarget->length() == nV)
     7333        {
     7334          delete Xtau;
     7335          Xtau = NewVectorlp(Gt);
     7336        }
     7337        else
     7338        {
     7339          delete Xtau;
     7340          Xtau = Mfpertvector(Gt,ivtarget);
     7341        }
     7342
     7343        rChangeCurrRing(oRing);
     7344        G = idrMoveR(Gt,testring,currRing);
     7345
     7346        // perturb the current vector w.r.t. the current GB
     7347        Mwlp = MivWeightOrderlp(omega);
     7348        if(ivtarget->length() > nV)
     7349        {
     7350          delete Mwlp;
     7351          Mwlp = MivMatrixOrderRefine(omega,ivtarget);
     7352        }
     7353        Xsigma = Mfpertvector(G, Mwlp);
     7354        delete Mwlp;
     7355
     7356        for(i=nV-1; i>=0; i--) {
     7357          (*omega2)[i] = (*Xtau)[nV+i];
     7358          (*omega)[i] = (*Xsigma)[nV+i];
     7359        }
     7360
     7361        delete next_vect;
     7362        to=clock();
     7363
     7364        /*
     7365        to avoid the value of Overflow_Error that occur in
     7366        Mfpertvector
     7367        */
     7368        Overflow_Error = FALSE;
     7369
     7370        next_vect = MkInterRedNextWeight(omega,omega2,G);
     7371        if(G->m[0] != NULL && polylength > 0)
     7372        {
     7373         
     7374          PrintS("//** rec_r_fractal_call: there is a polynomial in Gomega with at least 3 monomials, low-dimensional facet of the cone");
     7375         
     7376          delete next_vect;
     7377          next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,nlev);
     7378          if(isNegNolVector(next_vect)==1)
     7379          {
     7380            delete next_vect;
     7381            next_vect = MkInterRedNextWeight(omega,omega2,G);
     7382          }
     7383        }
     7384        xtnw=xtnw+clock()-to;
     7385      }
     7386//#ifdef PRINT_VECTORS
     7387      if(printout > 0)
     7388      {
     7389        MivString(omega, omega2, next_vect);
     7390      }
     7391//#endif
     7392
     7393    /* check whether the the computed vector is in the correct cone
     7394       If no, the reduced GB of an omega-homogeneous ideal will be
     7395       computed by Buchberger algorithm and stop this recursion step*/
     7396    //if(test_w_in_ConeCC(G, next_vect) != 1) //e.g. Example s7, cyc6
     7397    if(Overflow_Error == TRUE || test_w_in_ConeCC(G,next_vect) != 1)
     7398    {
     7399      delete next_vect;
     7400      if(ivtarget->length() == nV)
     7401      {
     7402        if (rParameter(currRing) != NULL)
     7403        {
     7404          DefRingPar(omtmp);
     7405        }
     7406        else
     7407        {
     7408          rChangeCurrRing(VMrDefault(omtmp));
     7409        }
     7410      }
     7411      else
     7412      {
     7413        rChangeCurrRing(VMatrDefault(ivtarget));
     7414      }
     7415#ifdef TEST_OVERFLOW
     7416      Gt = idrMoveR(G, oRing,currRing);
     7417      Gt = NULL;
     7418      return(Gt);
     7419#endif
     7420      if(printout > 0)
     7421      {
     7422        Print("\n//** rec_r_fractal_call: applying Buchberger's algorithm in ring r = %s;",
     7423              rString(currRing));
     7424      }
     7425      to=clock();
     7426      Gt = idrMoveR(G, oRing,currRing);
     7427      G1 = MstdCC(Gt);
     7428      xtextra=xtextra+clock()-to;
     7429      Gt = NULL;
     7430
     7431      delete omega2;
     7432      delete altomega;
     7433      if(printout > 0)
     7434      {
     7435        Print("\n//** rec_r_fractal_call: Leaving the %d-th recursion with %d steps.\n",
     7436              nlev, nwalks);
     7437        //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     7438      }
     7439      nnflow ++;
     7440      Overflow_Error = FALSE;
     7441      return (G1);
     7442    }
     7443    /*
     7444       If the perturbed target vector stays in the correct cone,
     7445       return the current Groebner basis.
     7446       Otherwise, return the Groebner basis computed with Buchberger's
     7447       algorithm.
     7448       Then we update the perturbed target vectors w.r.t. this GB.
     7449    */
     7450    if (MivComp(next_vect, XivNull) == 1)
     7451    {
     7452      // The computed vector is equal to the origin vector,
     7453      // because t is not defined
     7454      if(ivtarget->length() == nV)
     7455      {
     7456        if (rParameter(currRing) != NULL)
     7457          DefRingPar(omtmp);
     7458        else
     7459          rChangeCurrRing(VMrDefault(omtmp));
     7460      }
     7461      else
     7462      {
     7463        rChangeCurrRing(VMatrDefault(ivtarget));
     7464      }
     7465      testring = currRing;
     7466      Gt = idrMoveR(G, oRing,currRing);
     7467
     7468      if(test_w_in_ConeCC(Gt, omega2) == 1)
     7469      {
    62887470        delete omega2;
    62897471        delete next_vect;
    62907472        delete altomega;
    6291         //Print("\n// Leaving the %d-th recursion with %d steps ",nlev, nwalks);
    6292         //Print(" ** Overflow_Error? (%d)", Overflow_Error);
    6293 
     7473        if(printout > 0)
     7474        {
     7475          Print("\n//** rec_r_fractal_call: Leaving the %d-th recursion with %d steps.\n",
     7476                nlev, nwalks);
     7477          //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     7478        }
    62947479        return (Gt);
    62957480      }
    62967481      else
    6297       {
    6298         //ivString(omega2, "tau'");
    6299         //Print("\n//  tau' doesn't stay in the correct cone!!");
     7482      {
     7483        if(printout > 0)
     7484        {
     7485          Print("\n//** rec_r_fractal_call: target weight doesn't stay in the correct cone.\n");
     7486        }
    63007487
    63017488#ifndef  MSTDCC_FRACTAL
    6302         //07.08.03
    63037489        //ivString(Xtau, "old Xtau");
    6304         intvec* Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp));
     7490        intvec* Xtautmp;
     7491        if(ivtarget->length() == nV)
     7492        {
     7493          Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp));
     7494        }
     7495        else
     7496        {
     7497          Xtautmp = Mfpertvector(Gt, ivtarget);
     7498        }
    63057499#ifdef TEST_OVERFLOW
    63067500      if(Overflow_Error == TRUE)
     
    63307524
    63317525      FRACTAL_MSTDCC:
    6332         //Print("\n//  apply BB-Alg in ring = %s;", rString(currRing));
     7526        if(printout > 0)
     7527        {
     7528          Print("\n//** rec_r_fractal_call: apply Buchberger's algorithm in ring = %s.\n",
     7529                rString(currRing));
     7530        }
    63337531        to=clock();
    63347532        G = MstdCC(Gt);
     
    63387536
    63397537        // update the original target vector w.r.t. the current GB
    6340         if(MivSame(Xivinput, Xivlp) == 1)
    6341           if (rParameter(currRing) != NULL)
    6342             DefRingParlp();
     7538        if(ivtarget->length() == nV)
     7539        {
     7540          if(MivSame(Xivinput, Xivlp) == 1)
     7541            if (rParameter(currRing) != NULL)
     7542              DefRingParlp();
     7543            else
     7544              VMrDefaultlp();
    63437545          else
    6344             VMrDefaultlp();
     7546            if (rParameter(currRing) != NULL)
     7547              DefRingPar(Xivinput);
     7548            else
     7549              rChangeCurrRing(VMrDefault(Xivinput));
     7550        }
    63457551        else
    6346           if (rParameter(currRing) != NULL)
    6347             DefRingPar(Xivinput);
    6348           else
    6349             rChangeCurrRing(VMrDefault1(Xivinput)); //Aenderung6
    6350 
     7552        {
     7553          rChangeCurrRing(VMatrRefine(ivtarget,Xivinput));
     7554        }
    63517555        testring = currRing;
    63527556        Gt = idrMoveR(G, oRing,currRing);
    63537557
    6354         delete Xtau;
    6355         Xtau = NewVectorlp(Gt);
     7558        // perturb the original target vector w.r.t. the current GB
     7559        if(ivtarget->length() == nV)
     7560        {
     7561          delete Xtau;
     7562          Xtau = NewVectorlp(Gt);
     7563        }
     7564        else
     7565        {
     7566          delete Xtau;
     7567          Xtau = Mfpertvector(Gt,ivtarget);
     7568        }
    63567569
    63577570        rChangeCurrRing(oRing);
     
    63617574        delete next_vect;
    63627575        delete altomega;
    6363         /*
    6364           Print("\n// Leaving the %d-th recursion with %d steps,", nlev,nwalks);
    6365           Print(" ** Overflow_Error? (%d)", Overflow_Error);
    6366         */
     7576        if(printout > 0)
     7577        {
     7578          Print("\n//** rec_r_fractal_call: Leaving the %d-th recursion with %d steps.\n",
     7579                nlev,nwalks);
     7580          //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     7581        }
    63677582        if(Overflow_Error == TRUE)
    63687583          nnflow ++;
     
    63717586        return(G);
    63727587      }
    6373     }
    6374 
    6375     for(i=nV-1; i>=0; i--) {
     7588    } //end of if(MivComp(next_vect, XivNull) == 1)
     7589
     7590    for(i=nV-1; i>=0; i--)
     7591    {
    63767592      (*altomega)[i] = (*omega)[i];
    63777593      (*omega)[i] = (*next_vect)[i];
     
    63807596
    63817597    to=clock();
    6382     /* Take the initial form of <G> w.r.t. omega */
     7598    // Take the initial form of <G> w.r.t. omega
    63837599    Gomega = MwalkInitialForm(G, omega);
    63847600    xtif=xtif+clock()-to;
     7601    //polylength = 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise
     7602    polylength = lengthpoly(Gomega);
     7603    if(printout > 1)
     7604    {
     7605      idString(Gomega,"//** rec_r_fractal_call: Gomega");
     7606    }
     7607
     7608    if(reduction == 0)
     7609    {
     7610      /* Check whether the intermediate weight vector lies in the interior of the cone.
     7611       * If so, only perform reductions. Otherwise apply Buchberger's algorithm. */
     7612      FF = middleOfCone(G,Gomega);
     7613      if( FF != NULL)
     7614      {
     7615        idDelete(&G);
     7616        G = idCopy(FF);
     7617        idDelete(&FF);
     7618        /* Compue next vector. */
     7619        goto NEXT_VECTOR_FRACTAL;
     7620      }
     7621    }
    63857622
    63867623#ifndef  BUCHBERGER_ALG
     
    63897626    else
    63907627      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
    6391 #endif // BUCHBERGER_ALG
    6392 
    6393     if (rParameter(currRing) != NULL)
    6394       DefRingPar(omega);
     7628#endif
     7629    if(ivtarget->length() == nV)
     7630    {
     7631      if (rParameter(currRing) != NULL)
     7632        DefRingPar(omega);
     7633      else
     7634        rChangeCurrRing(VMrDefault(omega));
     7635    }
    63957636    else
    6396       rChangeCurrRing(VMrDefault1(omega)); //Aenderung7
    6397 
     7637    {
     7638      rChangeCurrRing(VMatrRefine(ivtarget,omega));
     7639    }
    63987640    Gomega1 = idrMoveR(Gomega, oRing,currRing);
    63997641
    6400     /* Maximal recursion depth, to compute a red. GB */
    6401     /* Fractal walk with the alternative recursion */
    6402     /* alternative recursion */
    6403     // if(nlev == nV || lengthpoly(Gomega1) == 0)
     7642    // Maximal recursion depth, to compute a red. GB
     7643    // Fractal walk with the alternative recursion
     7644    // alternative recursion
    64047645    if(nlev == Xnlev || lengthpoly(Gomega1) == 0)
    6405       //if(nlev == nV) // blind recursion
    6406     {
    6407       /*
    6408       if(Xnlev != nV)
    6409       {
    6410         Print("\n// ** Xnlev = %d", Xnlev);
    6411         ivString(Xtau, "Xtau");
    6412       }
    6413       */
     7646    {
    64147647      to=clock();
    64157648#ifdef  BUCHBERGER_ALG
     
    64187651      Gresult =kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,omega);
    64197652      delete hilb_func;
    6420 #endif // BUCHBERGER_ALG
     7653#endif
    64217654      xtstd=xtstd+clock()-to;
    64227655    }
    6423     else {
     7656    else
     7657    {
    64247658      rChangeCurrRing(oRing);
    64257659      Gomega1 = idrMoveR(Gomega1, oRing,currRing);
    6426       Gresult = rec_fractal_call(idCopy(Gomega1),nlev+1,omega);
    6427     }
    6428 
    6429     //convert a Groebner basis from a ring to another ring,
     7660      Gresult = rec_r_fractal_call(idCopy(Gomega1),nlev+1,omega,weight_rad,reduction,printout);
     7661    }
     7662    if(printout > 2)
     7663    {
     7664      idString(Gresult,"//** rec_r_fractal_call: M");
     7665    }
     7666    //convert a Groebner basis from a ring to another ring
    64307667    new_ring = currRing;
    64317668
     
    64357672
    64367673    to=clock();
    6437     /* Lifting process */
     7674    // Lifting process
    64387675    F = MLifttwoIdeal(Gomega2, Gresult1, G);
    64397676    xtlift=xtlift+clock()-to;
     7677
     7678    if(printout > 2)
     7679    {
     7680      idString(F,"//** rec_r_fractal_call: F");
     7681    }
     7682
    64407683    idDelete(&Gresult1);
    64417684    idDelete(&Gomega2);
     
    64437686
    64447687    rChangeCurrRing(new_ring);
    6445     F1 = idrMoveR(F, oRing,currRing);
    6446 
     7688    //F1 = idrMoveR(F, oRing,currRing);
     7689    G = idrMoveR(F,oRing,currRing);
    64477690    to=clock();
    6448     /* Interreduce G */
    6449     G = kInterRedCC(F1, NULL);
     7691    // Interreduce G
     7692    //G = kInterRedCC(F1, NULL);
    64507693    xtred=xtred+clock()-to;
    6451     idDelete(&F1);
     7694    //idDelete(&F1);
    64527695  }
    64537696}
    6454 
    6455 /************************************************************************
    6456  * Perturb the start weight vector at the top level with random element *
    6457  ************************************************************************/
    6458 static ideal rec_r_fractal_call(ideal G, int nlev, intvec* omtmp, int weight_rad)
    6459 {
    6460   Overflow_Error =  FALSE;
    6461   //Print("\n\n// Entering the %d-th recursion:", nlev);
    6462 
    6463   int i, nV = currRing->N;
    6464   ring new_ring, testring;
    6465   //ring extoRing;
    6466   ideal Gomega, Gomega1, Gomega2, F, F1, Gresult, Gresult1, G1, Gt;
    6467   int nwalks = 0;
    6468   intvec* Mwlp;
    6469 #ifndef BUCHBERGER_ALG
    6470   intvec* hilb_func;
    6471 #endif
    6472 //  intvec* extXtau;
    6473   intvec* next_vect;
    6474   intvec* omega2 = new intvec(nV);
    6475   intvec* altomega = new intvec(nV);
    6476 
    6477   //BOOLEAN isnewtarget = FALSE;
    6478 
    6479   // to avoid (1,0,...,0) as the target vector (Hans)
    6480   intvec* last_omega = new intvec(nV);
    6481   for(i=nV-1; i>0; i--)
    6482     (*last_omega)[i] = 1;
    6483   (*last_omega)[0] = 10000;
    6484 
    6485   intvec* omega = new intvec(nV);
    6486   for(i=0; i<nV; i++) {
    6487     if(Xsigma->length() == nV)
    6488       (*omega)[i] =  (*Xsigma)[i];
    6489     else
    6490       (*omega)[i] = (*Xsigma)[(nV*(nlev-1))+i];
    6491 
    6492     (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i];
    6493   }
    6494 
    6495    if(nlev == 1)  Xcall = 1;
    6496    else Xcall = 0;
    6497 
    6498   ring oRing = currRing;
    6499 
    6500   while(1)
    6501   {
    6502 #ifdef FIRST_STEP_FRACTAL
    6503     // perturb the current weight vector only on the top level or
    6504     // after perturbation of the both vectors, nlev = 2 as the top level
    6505     if((nlev == 1 && Xcall == 0) || (nlev == 2 && Xngleich == 1))
    6506       if(islengthpoly2(G) == 1)
    6507       {
    6508         Mwlp = MivWeightOrderlp(omega);
    6509         Xsigma = Mfpertvector(G, Mwlp);
    6510         delete Mwlp;
    6511         Overflow_Error = FALSE;
    6512       }
    6513 #endif
    6514     nwalks ++;
    6515     NEXT_VECTOR_FRACTAL:
    6516     to=clock();
    6517     /* determine the next border */
    6518     next_vect = MWalkRandomNextWeight(G, omega,omega2, weight_rad, 1+nlev);
    6519     //next_vect = MkInterRedNextWeight(omega,omega2,G);
    6520     xtnw=xtnw+clock()-to;
    6521 #ifdef PRINT_VECTORS
    6522     MivString(omega, omega2, next_vect);
    6523 #endif
    6524     oRing = currRing;
    6525 
    6526     /* We only perturb the current target vector at the recursion  level 1 */
    6527     if(Xngleich == 0 && nlev == 1) //(ngleich == 0) important, e.g. ex2, ex3
    6528       if (MivComp(next_vect, omega2) == 1)
    6529       {
    6530         /* to dispense with taking initial (and lifting/interreducing
    6531            after the call of recursion */
    6532         //Print("\n\n// ** Perturb the both vectors with degree %d with",nlev);
    6533         //idElements(G, "G");
    6534 
    6535         Xngleich = 1;
    6536         nlev +=1;
    6537 
    6538         if (rParameter(currRing) != NULL)
    6539           DefRingPar(omtmp);
    6540         else
    6541           rChangeCurrRing(VMrDefault1(omtmp)); //Aenderung3
    6542 
    6543         testring = currRing;
    6544         Gt = idrMoveR(G, oRing,currRing);
    6545 
    6546         /* perturb the original target vector w.r.t. the current GB */
    6547         delete Xtau;
    6548         Xtau = NewVectorlp(Gt);
    6549 
    6550         rChangeCurrRing(oRing);
    6551         G = idrMoveR(Gt, testring,currRing);
    6552 
    6553         /* perturb the current vector w.r.t. the current GB */
    6554         Mwlp = MivWeightOrderlp(omega);
    6555         Xsigma = Mfpertvector(G, Mwlp);
    6556         delete Mwlp;
    6557 
    6558         for(i=nV-1; i>=0; i--) {
    6559           (*omega2)[i] = (*Xtau)[nV+i];
    6560           (*omega)[i] = (*Xsigma)[nV+i];
    6561         }
    6562 
    6563         delete next_vect;
    6564         to=clock();
    6565 
    6566         /* to avoid the value of Overflow_Error that occur in Mfpertvector*/
    6567         Overflow_Error = FALSE;
    6568 
    6569         next_vect = MkInterRedNextWeight(omega,omega2,G);
    6570         xtnw=xtnw+clock()-to;
    6571 
    6572 #ifdef PRINT_VECTORS
    6573         MivString(omega, omega2, next_vect);
    6574 #endif
    6575       }
    6576 
    6577 
    6578     /* check whether the the computed vector is in the correct cone */
    6579     /* If no, the reduced GB of an omega-homogeneous ideal will be
    6580        computed by Buchberger algorithm and stop this recursion step*/
    6581     //if(test_w_in_ConeCC(G, next_vect) != 1) //e.g. Example s7, cyc6
    6582     if(Overflow_Error == TRUE)
    6583     {
    6584       delete next_vect;
    6585       if (rParameter(currRing) != NULL)
    6586       {
    6587         DefRingPar(omtmp);
    6588       }
    6589       else
    6590       {
    6591         rChangeCurrRing(VMrDefault1(omtmp)); // Aenderung4
    6592       }
    6593 #ifdef TEST_OVERFLOW
    6594       Gt = idrMoveR(G, oRing,currRing);
    6595       Gt = NULL; return(Gt);
    6596 #endif
    6597 
    6598       //Print("\n\n// apply BB's alg. in ring r = %s;", rString(currRing));
    6599       to=clock();
    6600       Gt = idrMoveR(G, oRing,currRing);
    6601       G1 = MstdCC(Gt);
    6602       xtextra=xtextra+clock()-to;
    6603       Gt = NULL;
    6604 
    6605       delete omega2;
    6606       delete altomega;
    6607 
    6608       //Print("\n// Leaving the %d-th recursion with %d steps", nlev, nwalks);
    6609       //Print("  ** Overflow_Error? (%d)", Overflow_Error);
    6610       nnflow ++;
    6611 
    6612       Overflow_Error = FALSE;
    6613       return (G1);
    6614     }
    6615 
    6616 
    6617     /* If the perturbed target vector stays in the correct cone,
    6618        return the current GB,
    6619        otherwise, return the computed  GB by the Buchberger-algorithm.
    6620        Then we update the perturbed target vectors w.r.t. this GB. */
    6621 
    6622     /* the computed vector is equal to the origin vector, since
    6623        t is not defined */
    6624     if (MivComp(next_vect, XivNull) == 1)
    6625     {
    6626       if (rParameter(currRing) != NULL)
    6627         DefRingPar(omtmp);
    6628       else
    6629         rChangeCurrRing(VMrDefault1(omtmp)); //Aenderung5
    6630 
    6631       testring = currRing;
    6632       Gt = idrMoveR(G, oRing,currRing);
    6633 
    6634       if(test_w_in_ConeCC(Gt, omega2) == 1) {
    6635         delete omega2;
    6636         delete next_vect;
    6637         delete altomega;
    6638         //Print("\n// Leaving the %d-th recursion with %d steps ",nlev, nwalks);
    6639         //Print(" ** Overflow_Error? (%d)", Overflow_Error);
    6640 
    6641         return (Gt);
    6642       }
    6643       else
    6644       {
    6645         //ivString(omega2, "tau'");
    6646         //Print("\n//  tau' doesn't stay in the correct cone!!");
    6647 
    6648 #ifndef  MSTDCC_FRACTAL
    6649         //07.08.03
    6650         //ivString(Xtau, "old Xtau");
    6651         intvec* Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp));
    6652 #ifdef TEST_OVERFLOW
    6653       if(Overflow_Error == TRUE)
    6654       Gt = NULL; return(Gt);
    6655 #endif
    6656 
    6657         if(MivSame(Xtau, Xtautmp) == 1)
    6658         {
    6659           //PrintS("\n// Update vectors are equal to the old vectors!!");
    6660           delete Xtautmp;
    6661           goto FRACTAL_MSTDCC;
    6662         }
    6663 
    6664         Xtau = Xtautmp;
    6665         Xtautmp = NULL;
    6666         //ivString(Xtau, "new  Xtau");
    6667 
    6668         for(i=nV-1; i>=0; i--)
    6669           (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i];
    6670 
    6671         //Print("\n//  ring tau = %s;", rString(currRing));
    6672         rChangeCurrRing(oRing);
    6673         G = idrMoveR(Gt, testring,currRing);
    6674 
    6675         goto NEXT_VECTOR_FRACTAL;
    6676 #endif
    6677 
    6678       FRACTAL_MSTDCC:
    6679         //Print("\n//  apply BB-Alg in ring = %s;", rString(currRing));
    6680         to=clock();
    6681         G = MstdCC(Gt);
    6682         xtextra=xtextra+clock()-to;
    6683 
    6684         oRing = currRing;
    6685 
    6686         // update the original target vector w.r.t. the current GB
    6687         if(MivSame(Xivinput, Xivlp) == 1)
    6688           if (rParameter(currRing) != NULL)
    6689             DefRingParlp();
    6690           else
    6691             VMrDefaultlp();
    6692         else
    6693           if (rParameter(currRing) != NULL)
    6694             DefRingPar(Xivinput);
    6695           else
    6696             rChangeCurrRing(VMrDefault1(Xivinput)); //Aenderung6
    6697 
    6698         testring = currRing;
    6699         Gt = idrMoveR(G, oRing,currRing);
    6700 
    6701         delete Xtau;
    6702         Xtau = NewVectorlp(Gt);
    6703 
    6704         rChangeCurrRing(oRing);
    6705         G = idrMoveR(Gt, testring,currRing);
    6706 
    6707         delete omega2;
    6708         delete next_vect;
    6709         delete altomega;
    6710         /*
    6711           Print("\n// Leaving the %d-th recursion with %d steps,", nlev,nwalks);
    6712           Print(" ** Overflow_Error? (%d)", Overflow_Error);
    6713         */
    6714         if(Overflow_Error == TRUE)
    6715           nnflow ++;
    6716 
    6717         Overflow_Error = FALSE;
    6718         return(G);
    6719       }
    6720     }
    6721 
    6722     for(i=nV-1; i>=0; i--) {
    6723       (*altomega)[i] = (*omega)[i];
    6724       (*omega)[i] = (*next_vect)[i];
    6725     }
    6726     delete next_vect;
    6727 
    6728     to=clock();
    6729     /* Take the initial form of <G> w.r.t. omega */
    6730     Gomega = MwalkInitialForm(G, omega);
    6731     xtif=xtif+clock()-to;
    6732 
    6733 #ifndef  BUCHBERGER_ALG
    6734     if(isNolVector(omega) == 0)
    6735       hilb_func = hFirstSeries(Gomega,NULL,NULL,omega,currRing);
    6736     else
    6737       hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
    6738 #endif // BUCHBERGER_ALG
    6739 
    6740     if (rParameter(currRing) != NULL)
    6741       DefRingPar(omega);
    6742     else
    6743       rChangeCurrRing(VMrDefault1(omega)); //Aenderung7
    6744 
    6745     Gomega1 = idrMoveR(Gomega, oRing,currRing);
    6746 
    6747     /* Maximal recursion depth, to compute a red. GB */
    6748     /* Fractal walk with the alternative recursion */
    6749     /* alternative recursion */
    6750     // if(nlev == nV || lengthpoly(Gomega1) == 0)
    6751     if(nlev == Xnlev || lengthpoly(Gomega1) == 0)
    6752       //if(nlev == nV) // blind recursion
    6753     {
    6754       /*
    6755       if(Xnlev != nV)
    6756       {
    6757         Print("\n// ** Xnlev = %d", Xnlev);
    6758         ivString(Xtau, "Xtau");
    6759       }
    6760       */
    6761       to=clock();
    6762 #ifdef  BUCHBERGER_ALG
    6763       Gresult = MstdhomCC(Gomega1);
    6764 #else
    6765       Gresult =kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,omega);
    6766       delete hilb_func;
    6767 #endif // BUCHBERGER_ALG
    6768       xtstd=xtstd+clock()-to;
    6769     }
    6770     else {
    6771       rChangeCurrRing(oRing);
    6772       Gomega1 = idrMoveR(Gomega1, oRing,currRing);
    6773       Gresult = rec_fractal_call(idCopy(Gomega1),nlev+1,omega);
    6774     }
    6775 
    6776     //convert a Groebner basis from a ring to another ring,
    6777     new_ring = currRing;
    6778 
    6779     rChangeCurrRing(oRing);
    6780     Gresult1 = idrMoveR(Gresult, new_ring,currRing);
    6781     Gomega2 = idrMoveR(Gomega1, new_ring,currRing);
    6782 
    6783     to=clock();
    6784     /* Lifting process */
    6785     F = MLifttwoIdeal(Gomega2, Gresult1, G);
    6786     xtlift=xtlift+clock()-to;
    6787     idDelete(&Gresult1);
    6788     idDelete(&Gomega2);
    6789     idDelete(&G);
    6790 
    6791     rChangeCurrRing(new_ring);
    6792     F1 = idrMoveR(F, oRing,currRing);
    6793 
    6794     to=clock();
    6795     /* Interreduce G */
    6796     G = kInterRedCC(F1, NULL);
    6797     xtred=xtred+clock()-to;
    6798     idDelete(&F1);
    6799   }
    6800 }
    6801 
    6802 
    68037697
    68047698
     
    68077701 *                                                                             *
    68087702 * The main procedur Mfwalk calls the recursive Subroutine                     *
    6809  * rec_fractal_call to compute the wanted Grï¿œbner basis.                       *
    6810  * At the main procedur we compute the reduced Grï¿œbner basis w.r.t. a "fast"   *
     7703 * rec_fractal_call to compute the wanted Groebner basis.                      *
     7704 * At the main procedur we compute the reduced Groebner basis w.r.t. a "fast"  *
    68117705 * order, e.g. "dp" and a sequence of weight vectors which are row vectors     *
    68127706 * of a matrix. This matrix defines the given monomial order, e.g. "lp"        *
    68137707 *******************************************************************************/
    6814 ideal Mfwalk(ideal G, intvec* ivstart, intvec* ivtarget)
     7708ideal Mfwalk(ideal G, intvec* ivstart, intvec* ivtarget,
     7709             int reduction, int printout)
    68157710{
     7711  BITSET save1 = si_opt_1; // save current options
     7712  if(reduction == 0)
     7713  {
     7714    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     7715    //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     7716  }
    68167717  Set_Error(FALSE);
    68177718  Overflow_Error = FALSE;
     
    68487749      intvec* iv_dp = MivUnit(nV); // define (1,1,...,1)
    68497750      intvec* Mdp;
    6850 
    6851       if(MivSame(ivstart, iv_dp) != 1)
    6852         Mdp = MivWeightOrderdp(ivstart);
     7751      if(ivstart->length() == nV)
     7752      {
     7753        if(MivSame(ivstart, iv_dp) != 1)
     7754          Mdp = MivWeightOrderdp(ivstart);
     7755        else
     7756          Mdp = MivMatrixOrderdp(nV);
     7757      }
    68537758      else
    6854         Mdp = MivMatrixOrderdp(nV);
     7759      {
     7760        Mdp = ivstart;
     7761      }
    68557762
    68567763      Xsigma = Mfpertvector(I, Mdp);
     
    68697776  Xivlp = Mivlp(nV);
    68707777
    6871   if(MivComp(ivtarget, Xivlp)  != 1)
    6872   {
    6873     if (rParameter(currRing) != NULL)
    6874       DefRingPar(ivtarget);
     7778  if(ivtarget->length() == nV)
     7779  {
     7780    if(MivComp(ivtarget, Xivlp)  != 1)
     7781    {
     7782      if (rParameter(currRing) != NULL)
     7783        DefRingPar(ivtarget);
     7784      else
     7785        rChangeCurrRing(VMrDefault(ivtarget));
     7786
     7787      I1 = idrMoveR(I, oldRing,currRing);
     7788      Mlp = MivWeightOrderlp(ivtarget);
     7789      Xtau = Mfpertvector(I1, Mlp);
     7790    }
    68757791    else
    6876       rChangeCurrRing(VMrDefault1(ivtarget)); //Aenderung1
    6877 
    6878     I1 = idrMoveR(I, oldRing,currRing);
    6879     Mlp = MivWeightOrderlp(ivtarget);
    6880     Xtau = Mfpertvector(I1, Mlp);
     7792    {
     7793      if (rParameter(currRing) != NULL)
     7794        DefRingParlp();
     7795      else
     7796        VMrDefaultlp();
     7797
     7798      I1 = idrMoveR(I, oldRing,currRing);
     7799      Mlp =  MivMatrixOrderlp(nV);
     7800      Xtau = Mfpertvector(I1, Mlp);
     7801    }
    68817802  }
    68827803  else
    68837804  {
    6884     if (rParameter(currRing) != NULL)
    6885       DefRingParlp();
    6886     else
    6887       VMrDefaultlp();
    6888 
    6889     I1 = idrMoveR(I, oldRing,currRing);
    6890     Mlp =  MivMatrixOrderlp(nV);
     7805    rChangeCurrRing(VMatrDefault(ivtarget));
     7806    I1 = idrMoveR(I,oldRing,currRing);
     7807    Mlp =  ivtarget;
    68917808    Xtau = Mfpertvector(I1, Mlp);
    68927809  }
     
    68997816  id_Delete(&I, oldRing);
    69007817  ring tRing = currRing;
    6901 
    6902   if (rParameter(currRing) != NULL)
    6903     DefRingPar(ivstart);
     7818  if(ivtarget->length() == nV)
     7819  {
     7820    if (rParameter(currRing) != NULL)
     7821      DefRingPar(ivstart);
     7822    else
     7823      rChangeCurrRing(VMrDefault(ivstart));
     7824  }
    69047825  else
    6905     rChangeCurrRing(VMrDefault1(ivstart)); //Aenderung2
     7826  {
     7827    rChangeCurrRing(VMatrDefault(ivstart));
     7828  }
    69067829
    69077830  I = idrMoveR(I1,tRing,currRing);
     
    69147837  ring helpRing = currRing;
    69157838
    6916   J = rec_fractal_call(J, 1, ivtarget);
     7839  J = rec_fractal_call(J,1,ivtarget,reduction,printout);
    69177840
    69187841  rChangeCurrRing(oldRing);
     
    69207843  idSkipZeroes(resF);
    69217844
     7845  si_opt_1 = save1; //set original options, e. g. option(RedSB)
    69227846  delete Xivlp;
    69237847  delete Xsigma;
     
    69387862}
    69397863
    6940 ideal Mfrwalk(ideal G, intvec* ivstart, intvec* ivtarget,int weight_rad)
     7864/*******************************************************************************
     7865 * The implementation of the fractal walk algorithm with random element        *
     7866 *                                                                             *
     7867 * The main procedur Mfwalk calls the recursive Subroutine                     *
     7868 * rec_r_fractal_call to compute the wanted Groebner basis.                    *
     7869 * At the main procedure we compute the reduced Groebner basis w.r.t. a "fast" *
     7870 * order, e.g. "dp" and a sequence of weight vectors which are row vectors     *
     7871 * of a matrix. This matrix defines the given monomial order, e.g. "lp"        *
     7872 *******************************************************************************/
     7873ideal Mfrwalk(ideal G, intvec* ivstart, intvec* ivtarget,
     7874              int weight_rad, int reduction, int printout)
    69417875{
     7876  BITSET save1 = si_opt_1; // save current options
     7877  if(reduction == 0)
     7878  {
     7879    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     7880    //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     7881  }
    69427882  Set_Error(FALSE);
    69437883  Overflow_Error = FALSE;
     
    69747914      intvec* iv_dp = MivUnit(nV); // define (1,1,...,1)
    69757915      intvec* Mdp;
    6976 
    6977       if(MivSame(ivstart, iv_dp) != 1)
    6978         Mdp = MivWeightOrderdp(ivstart);
     7916      if(ivstart->length() == nV)
     7917      {
     7918        if(MivSame(ivstart, iv_dp) != 1)
     7919          Mdp = MivWeightOrderdp(ivstart);
     7920        else
     7921          Mdp = MivMatrixOrderdp(nV);
     7922      }
    69797923      else
    6980         Mdp = MivMatrixOrderdp(nV);
     7924      {
     7925        Mdp = ivstart;
     7926      }
    69817927
    69827928      Xsigma = Mfpertvector(I, Mdp);
     
    69957941  Xivlp = Mivlp(nV);
    69967942
    6997   if(MivComp(ivtarget, Xivlp)  != 1)
    6998   {
    6999     if (rParameter(currRing) != NULL)
    7000       DefRingPar(ivtarget);
     7943  if(ivtarget->length() == nV)
     7944  {
     7945    if(MivComp(ivtarget, Xivlp)  != 1)
     7946    {
     7947      if (rParameter(currRing) != NULL)
     7948        DefRingPar(ivtarget);
     7949      else
     7950        rChangeCurrRing(VMrDefault(ivtarget));
     7951
     7952      I1 = idrMoveR(I, oldRing,currRing);
     7953      Mlp = MivWeightOrderlp(ivtarget);
     7954      Xtau = Mfpertvector(I1, Mlp);
     7955    }
    70017956    else
    7002       rChangeCurrRing(VMrDefault1(ivtarget)); //Aenderung1
    7003 
    7004     I1 = idrMoveR(I, oldRing,currRing);
    7005     Mlp = MivWeightOrderlp(ivtarget);
    7006     Xtau = Mfpertvector(I1, Mlp);
     7957    {
     7958      if (rParameter(currRing) != NULL)
     7959        DefRingParlp();
     7960      else
     7961        VMrDefaultlp();
     7962
     7963      I1 = idrMoveR(I, oldRing,currRing);
     7964      Mlp =  MivMatrixOrderlp(nV);
     7965      Xtau = Mfpertvector(I1, Mlp);
     7966    }
    70077967  }
    70087968  else
    70097969  {
    7010     if (rParameter(currRing) != NULL)
    7011       DefRingParlp();
    7012     else
    7013       VMrDefaultlp();
    7014 
    7015     I1 = idrMoveR(I, oldRing,currRing);
    7016     Mlp =  MivMatrixOrderlp(nV);
     7970    rChangeCurrRing(VMatrDefault(ivtarget));
     7971    I1 = idrMoveR(I,oldRing,currRing);
     7972    Mlp =  ivtarget;
    70177973    Xtau = Mfpertvector(I1, Mlp);
    70187974  }
     
    70257981  id_Delete(&I, oldRing);
    70267982  ring tRing = currRing;
    7027 
    7028   if (rParameter(currRing) != NULL)
    7029     DefRingPar(ivstart);
     7983  if(ivtarget->length() == nV)
     7984  {
     7985    if (rParameter(currRing) != NULL)
     7986      DefRingPar(ivstart);
     7987    else
     7988      rChangeCurrRing(VMrDefault(ivstart));
     7989  }
    70307990  else
    7031     rChangeCurrRing(VMrDefault1(ivstart)); //Aenderung2
     7991  {
     7992    rChangeCurrRing(VMatrDefault(ivstart));
     7993  }
    70327994
    70337995  I = idrMoveR(I1,tRing,currRing);
     
    70398001  ideal resF;
    70408002  ring helpRing = currRing;
    7041 //ideal G, int nlev, intvec* omtmp, int weight_rad)
    7042   J = rec_r_fractal_call(J, 1, ivtarget,weight_rad);
     8003
     8004  J = rec_r_fractal_call(J,1,ivtarget,weight_rad,reduction,printout);
    70438005
    70448006  rChangeCurrRing(oldRing);
     
    70468008  idSkipZeroes(resF);
    70478009
     8010  si_opt_1 = save1; //set original options, e. g. option(RedSB)
    70488011  delete Xivlp;
    70498012  delete Xsigma;
     
    71018064  intvec* hilb_func;
    71028065#endif
    7103   /* to avoid (1,0,...,0) as the target vector */
     8066  // to avoid (1,0,...,0) as the target vector
    71048067  intvec* last_omega = new intvec(nV);
    71058068  for(i=nV-1; i>0; i--)
     
    71168079
    71178080  to=clock();
    7118   /* compute a red. GB w.r.t. the help ring */
     8081  // compute a red. GB w.r.t. the help ring
    71198082  if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) = "dp"
    71208083    G = MstdCC(G);
     
    71258088      DefRingPar(curr_weight);
    71268089    else
    7127       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung 4
     8090      rChangeCurrRing(VMrDefault(curr_weight));
    71288091    G = idrMoveR(G, XXRing,currRing);
    71298092    G = MstdCC(G);
     
    71518114      DefRingPar(curr_weight);
    71528115    else
    7153       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 5
     8116      rChangeCurrRing(VMrDefault(curr_weight));
    71548117    to=clock();
    71558118    Gw = idrMoveR(G, exring,currRing);
     
    71868149      DefRingPar(curr_weight);
    71878150    else
    7188       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 6
     8151      rChangeCurrRing(VMrDefault(curr_weight));
    71898152
    71908153    newRing = currRing;
     
    72718234          DefRingPar(target_tmp);
    72728235        else
    7273           rChangeCurrRing(VMrDefault(target_tmp)); // Aenderung 8
     8236          rChangeCurrRing(VMrDefault(target_tmp));
    72748237
    72758238      lpRing = currRing;
     
    73318294          DefRingPar(target_tmp);
    73328295        else
    7333           rChangeCurrRing(VMrDefault(target_tmp)); //Aenderung 9
     8296          rChangeCurrRing(VMrDefault(target_tmp));
    73348297
    73358298      lpRing = currRing;
     
    75348497    else
    75358498    {
    7536       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung 10
     8499      rChangeCurrRing(VMrDefault(curr_weight));
    75378500    }
    75388501    G = idrMoveR(G, XXRing,currRing);
     
    75678530    else
    75688531    {
    7569       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 11
     8532      rChangeCurrRing(VMrDefault(curr_weight));
    75708533    }
    75718534    to=clock();
     
    76108573    else
    76118574    {
    7612       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung 12
     8575      rChangeCurrRing(VMrDefault(curr_weight));
    76138576    }
    76148577    newRing = currRing;
     
    77798742        else
    77808743        {
    7781           rChangeCurrRing(VMrDefault(target_tmp)); // Aenderung 13
     8744          rChangeCurrRing(VMrDefault(target_tmp));
    77828745        }
    77838746      }
     
    78528815        else
    78538816        {
    7854           rChangeCurrRing(VMrDefault(target_tmp)); // Aenderung 14
     8817          rChangeCurrRing(VMrDefault(target_tmp));
    78558818        }
    78568819      }
     
    80959058    else
    80969059    {
    8097       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung 15
     9060      rChangeCurrRing(VMrDefault(curr_weight));
    80989061    }
    80999062    newRing = currRing;
     
    82379200
    82389201  //Print("\n// \"Mpwalk\" (1,%d) took %d steps and %.2f sec. Overflow_Error (%d)", tp_deg, nwalk, ((double) clock()-tinput)/1000000, nOverflow_Error);
    8239 
     9202  Print("\n// Mprwalk took %d steps. Ring= %s;\n", nwalk, rString(currRing));
    82409203  return(result);
    8241 }
    8242 
    8243 /*******************************************************
    8244  * THE PERTURBATION WALK ALGORITHM WITH RANDOM ELEMENT *
    8245  *******************************************************/
    8246 ideal Mprwalk(ideal Go, intvec* curr_weight, intvec* target_weight, int weight_rad, int op_deg, int tp_deg, ring baseRing)
    8247 {
    8248   BITSET save1 = si_opt_1; // save current options
    8249   si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
    8250   Set_Error(FALSE);
    8251   Overflow_Error = FALSE;
    8252 #ifdef TIME_TEST
    8253   clock_t tinput=0, tostd=0, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
    8254   xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0;
    8255   tinput = clock();
    8256   clock_t tim;
    8257 #endif
    8258   int i,nwalk,nV = baseRing->N;
    8259 
    8260   ideal G, Gomega, M, F, Gomega1, Gomega2, M1;
    8261   ring newRing;
    8262   ring XXRing = baseRing;
    8263   intvec* exivlp = Mivlp(nV);
    8264   intvec* orig_target = target_weight;
    8265   intvec* pert_target_vector = target_weight;
    8266   intvec* ivNull = new intvec(nV);
    8267   intvec* tmp_weight = new intvec(nV);
    8268 #ifdef CHECK_IDEAL_MWALK
    8269   poly p;
    8270 #endif
    8271   for(i=0; i<nV; i++)
    8272   {
    8273     (*tmp_weight)[i] = (*curr_weight)[i];
    8274   }
    8275 #ifndef BUCHBERGER_ALG
    8276   intvec* hilb_func;
    8277    // to avoid (1,0,...,0) as the target vector
    8278   intvec* last_omega = new intvec(nV);
    8279   for(i=0 i<nV; i++)
    8280   {
    8281     (*last_omega)[i] = 1;
    8282   }
    8283   (*last_omega)[0] = 10000;
    8284 #endif
    8285   baseRing = currRing;
    8286   newRing = VMrDefault(curr_weight);
    8287   rChangeCurrRing(newRing);
    8288   G = idrMoveR(Go,baseRing,currRing);
    8289 #ifdef TIME_TEST
    8290   to = clock();
    8291 #endif
    8292   G = kStd(G,NULL,testHomog,NULL,NULL,0,0,NULL);
    8293   idSkipZeroes(G);
    8294 #ifdef TIME_TEST
    8295   tostd = tostd + to - clock();
    8296 #endif
    8297 #ifdef CHECK_IDEAL_MWALK
    8298   idString(G,"G");
    8299 #endif
    8300   if(op_deg >1)
    8301   {
    8302     if(MivComp(curr_weight,MivUnit(nV)) == 1) //ring order is "dp"
    8303     {
    8304       curr_weight = MPertVectors(G, MivMatrixOrderdp(nV), op_deg);
    8305     }
    8306     else //ring order is not "dp"
    8307     {
    8308       curr_weight = MPertVectors(G, MivMatrixOrder(curr_weight), op_deg);
    8309     }
    8310   }
    8311   baseRing = currRing;
    8312   if(tp_deg > 1 && tp_deg <= nV)
    8313   {
    8314     pert_target_vector = target_weight;
    8315   }
    8316 #ifdef CHECK_IDEAL_MWALK
    8317   ivString(curr_weight, "new curr_weight");
    8318   ivString(target_weight, "new target_weight");
    8319 #endif
    8320   nwalk = 0;
    8321   while(1)
    8322   {
    8323     nwalk ++;
    8324 #ifdef TIME_TEST
    8325     to = clock();
    8326 #endif
    8327     Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
    8328 #ifdef TIME_TEST
    8329     tif = tif + clock()-to; //time for computing initial form ideal
    8330 #endif
    8331 #ifdef CHECK_IDEAL_MWALK
    8332     idString(Gomega,"Gomega");
    8333 #endif
    8334 #ifndef  BUCHBERGER_ALG
    8335     if(isNolVector(curr_weight) == 0)
    8336     {
    8337       hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
    8338     }
    8339     else
    8340     {
    8341       hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
    8342     }
    8343 #endif
    8344     if(nwalk == 1)
    8345     {
    8346       newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
    8347     }
    8348     else
    8349     {
    8350       newRing = VMrRefine(curr_weight,target_weight); //define a new ring with ordering "(a(curr_weight),Wp(target_weight))"
    8351     }
    8352     rChangeCurrRing(newRing);
    8353     Gomega1 = idrMoveR(Gomega, baseRing,currRing);
    8354     idDelete(&Gomega);
    8355     // compute a Groebner basis of <Gomega> w.r.t. "newRing"
    8356 #ifdef TIME_TEST
    8357     to = clock();
    8358 #endif
    8359 #ifndef  BUCHBERGER_ALG
    8360     M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
    8361     delete hilb_func;
    8362 #else
    8363     M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL);
    8364 #endif
    8365     idSkipZeroes(M);
    8366 #ifdef TIME_TEST
    8367     tstd = tstd + clock() - to;
    8368 #endif
    8369 #ifdef CHECK_IDEAL_MWALK
    8370     idString(M, "M");
    8371 #endif
    8372     //change the ring to baseRing
    8373     rChangeCurrRing(baseRing);
    8374     M1 =  idrMoveR(M, newRing,currRing);
    8375     idDelete(&M);
    8376     Gomega2 = idrMoveR(Gomega1, newRing,currRing);
    8377     idDelete(&Gomega1);
    8378     to = clock();
    8379     // 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
    8380     F = MLifttwoIdeal(Gomega2, M1, G);
    8381     idSkipZeroes(F);
    8382 #ifdef TIME_TEST
    8383     tlift = tlift + clock() - to;
    8384 #endif
    8385 #ifdef CHECK_IDEAL_MWALK
    8386     idString(F,"F");
    8387 #endif
    8388     rChangeCurrRing(newRing); // change the ring to newRing
    8389     G = idrMoveR(F,baseRing,currRing);
    8390     idDelete(&F);
    8391     baseRing = currRing; // set baseRing equal to newRing
    8392 #ifdef CHECK_IDEAL_MWALK
    8393     idString(G,"G");
    8394 #endif
    8395 #ifdef TIME_TEST
    8396     to = clock();
    8397 #endif
    8398     intvec* next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, op_deg);
    8399 #ifdef TIME_TEST
    8400     tnw = tnw + clock() - to;