Changeset bbfc4a8 in git


Ignore:
Timestamp:
Jul 2, 2015, 7:32:53 PM (9 years ago)
Author:
Stephan Oberfranz <oberfran@…>
Branches:
(u'spielwiese', '4a9821a93ffdc22a6696668bd4f6b8c9de3e6c5f')
Children:
8d7ca03d4444afafa6c2111db95d65216b93ad46
Parents:
500e4b69e2df574c6ae8023f32cbdfbdbd0d5e8ef34cd445520f5220554ca71a6324d2a0ecd2c7f6
Message:
update

Merge branch 'spielwiese' of github.com:Singular/Sources into spielwiese
Files:
3 added
14 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/grwalk.lib

    rf34cd44 rbbfc4a8  
    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
    rf34cd44 rbbfc4a8  
    1616
    1717PROCEDURES:
    18 modWalk(ideal,int,int[,int,int,int,int]);        standard basis conversion of I using modular methods (chinese remainder)
     18
     19modWalk(I,#);                   standard basis conversion of I by Groebner Walk using modular methods
     20modrWalk(I,radius,pertdeg,#);   standard basis conversion of I by Random Walk using modular methods
     21modfWalk(I,#);                  standard basis conversion of I by Fractal Walk using modular methods
     22modfrWalk(I,radius,#);          standard basis conversion of I by Random Fractal Walk using modular methods
    1923";
    2024
    21 LIB "poly.lib";
    22 LIB "ring.lib";
    23 LIB "parallel.lib";
    2425LIB "rwalk.lib";
    2526LIB "grwalk.lib";
    26 LIB "modstd.lib";
    27 
    28 
    29 ////////////////////////////////////////////////////////////////////////////////
    30 
    31 static proc modpWalk(def II, int p, int variant, list #)
    32 "USAGE:  modpWalk(I,p,#); I ideal, p integer, variant integer
    33 ASSUME:  If size(#) > 0, then
    34            #[1] is an intvec describing the current weight vector
    35            #[2] is an intvec describing the target weight vector
    36 RETURN:  ideal - a standard basis of I mod p, integer - p
    37 NOTE:    The procedure computes a standard basis of the ideal I modulo p and
    38          fetches the result to the basering.
    39 EXAMPLE: example modpWalk; shows an example
    40 "
    41 {
    42   option(redSB);
    43   int k,nvar@r;
    44   def R0 = basering;
    45   string ordstr_R0 = ordstr(R0);
    46   list rl = ringlist(R0);
    47   int sizerl = size(rl);
    48   int neg = 1 - attrib(R0,"global");
    49   if(typeof(II) == "ideal")
    50   {
    51     ideal I = II;
     27LIB "modular.lib";
     28
     29proc modWalk(ideal I, list #)
     30"USAGE:   modWalk(I, [, v, w]); I ideal, v intvec, w intvec
     31RETURN:   a standard basis of I
     32NOTE:     The procedure computes a standard basis of I (over the rational
     33          numbers) by using modular methods.
     34SEE ALSO: modular
     35EXAMPLE:  example modWalk; shows an example"
     36{
     37    /* read optional parameter */
     38    if (size(#) > 0) {
     39        if (size(#) == 1) {
     40            intvec w = #[1];
     41        }
     42        if (size(#) == 2) {
     43            intvec v = #[1];
     44            intvec w = #[2];
     45        }
     46        if (size(#) > 2 || typeof(#[1]) != "intvec") {
     47            ERROR("wrong optional parameter");
     48        }
     49    }
     50
     51    /* save options */
     52    intvec opt = option(get);
     53    option(redSB);
     54
     55    /* set additional parameters */
     56    int reduction = 1;
     57    int printout = 0;
     58
     59    /* call modular() */
     60    if (size(#) > 0) {
     61        I = modular("gwalk", list(I,reduction,printout,#));
     62    }
     63    else {
     64        I = modular("gwalk", list(I,reduction,printout));
     65    }
     66
     67    /* return the result */
     68    attrib(I, "isSB", 1);
     69    option(set, opt);
     70    return(I);
     71}
     72example
     73{
     74    "EXAMPLE:";
     75    echo = 2;
     76    ring R1 = 0, (x,y,z,t), dp;
     77    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
     78    I = std(I);
     79    ring R2 = 0, (x,y,z,t), lp;
     80    ideal I = fetch(R1, I);
     81    ideal J = modWalk(I);
     82    J;
     83}
     84
     85proc modrWalk(ideal I, int radius, int pertdeg, list #)
     86"USAGE:   modrWalk(I, radius, pertdeg[, v, w]);
     87          I ideal, radius int, pertdeg int, v intvec, w intvec
     88RETURN:   a standard basis of I
     89NOTE:     The procedure computes a standard basis of I (over the rational
     90          numbers) by using modular methods.
     91SEE ALSO: modular
     92EXAMPLE:  example modrWalk; shows an example"
     93{
     94    /* read optional parameter */
     95    if (size(#) > 0) {
     96        if (size(#) == 1) {
     97            intvec w = #[1];
     98        }
     99        if (size(#) == 2) {
     100            intvec v = #[1];
     101            intvec w = #[2];
     102        }
     103        if (size(#) > 2 || typeof(#[1]) != "intvec") {
     104            ERROR("wrong optional parameter");
     105        }
     106    }
     107
     108    /* save options */
     109    intvec opt = option(get);
     110    option(redSB);
     111
     112    /* set additional parameters */
     113    int reduction = 1;
     114    int printout = 0;
     115
     116    /* call modular() */
     117    if (size(#) > 0) {
     118        I = modular("rwalk", list(I,radius,pertdeg,reduction,printout,#));
     119    }
     120    else {
     121        I = modular("rwalk", list(I,radius,pertdeg,reduction,printout));
     122    }
     123
     124    /* return the result */
     125    attrib(I, "isSB", 1);
     126    option(set, opt);
     127    return(I);
     128}
     129example
     130{
     131    "EXAMPLE:";
     132    echo = 2;
     133    ring R1 = 0, (x,y,z,t), dp;
     134    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
     135    I = std(I);
     136    ring R2 = 0, (x,y,z,t), lp;
     137    ideal I = fetch(R1, I);
    52138    int radius = 2;
    53     int pert_deg = 2;
    54   }
    55   if(typeof(II) == "list" && typeof(II[1]) == "ideal")
    56   {
    57     ideal I = II[1];
    58     if(size(II) == 2)
    59     {
    60       int radius = II[2];
    61       int pert_deg = 2;
    62     }
    63     if(size(II) == 3)
    64     {
    65       int radius = II[2];
    66       int pert_deg = II[3];
    67     }
    68   }
    69   rl[1] = p;
    70   int h = homog(I);
    71   def @r = ring(rl);
    72   setring @r;
    73   ideal i = fetch(R0,I);
    74   string order;
    75   if(system("nblocks") <= 2)
    76   {
    77     if(find(ordstr_R0, "M") + find(ordstr_R0, "lp") + find(ordstr_R0, "rp") <= 0)
    78     {
    79       order = "simple";
    80     }
    81   }
    82 
    83 //-------------------------  make i homogeneous  -----------------------------
    84   if(!mixedTest() && !h)
    85   {
    86     if(!((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg))
    87     {
    88       if(!((order == "simple") || (sizerl > 4)))
    89       {
    90         list rl@r = ringlist(@r);
    91         nvar@r = nvars(@r);
    92         intvec w;
    93         for(k = 1; k <= nvar@r; k++)
    94         {
    95           w[k] = deg(var(k));
    96         }
    97         w[nvar@r + 1] = 1;
    98         rl@r[2][nvar@r + 1] = "homvar";
    99         rl@r[3][2][2] = w;
    100         def HomR = ring(rl@r);
    101         setring HomR;
    102         ideal i = imap(@r, i);
    103         i = homog(i, homvar);
    104       }
    105     }
    106   }
    107 
    108 //-------------------------  compute a standard basis mod p  -----------------------------
    109 
    110   if(variant == 1)
    111   {
    112     if(size(#)>0)
    113     {
    114       i = rwalk(i,radius,pert_deg,#);
    115      // rwalk(i,radius,pert_deg,#); std(i);
    116     }
    117     else
    118     {
    119       i = rwalk(i,radius,pert_deg);
    120     }
    121   }
    122   if(variant == 2)
    123   {
    124     if(size(#) == 2)
    125     {
    126       i = gwalk(i,#);
    127     }
    128     else
    129     {
    130       i = gwalk(i);
    131     }
    132   }
    133   if(variant == 3)
    134   {
    135     if(size(#) == 2)
    136     {
    137       i = frandwalk(i,radius,#);
    138     }
    139     else
    140     {
    141       i = frandwalk(i,radius);
    142     }
    143   }
    144   if(variant == 4)
    145   {
    146     if(size(#) == 2)
    147     {
    148       i=fwalk(i,#);
    149     }
    150     else
    151     {
    152       i=fwalk(i);
    153     }
    154   }
    155   if(variant == 5)
    156   {
    157     if(size(#) == 2)
    158     {
    159      i=prwalk(i,radius,pert_deg,pert_deg,#);
    160     }
    161     else
    162     {
    163       i=prwalk(i,radius,pert_deg,pert_deg);
    164     }
    165   }
    166   if(variant == 6)
    167   {
    168     if(size(#) == 2)
    169     {
    170       i=pwalk(i,pert_deg,pert_deg,#);
    171     }
    172     else
    173     {
    174       i=pwalk(i,pert_deg,pert_deg);
    175     }
    176   }
    177 
    178   if(!mixedTest() && !h)
    179   {
    180     if(!((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg))
    181     {
    182       if(!((order == "simple") || (sizerl > 4)))
    183       {
    184         i = subst(i, homvar, 1);
    185         i = simplify(i, 34);
    186         setring @r;
    187         i = imap(HomR, i);
    188         i = interred(i);
    189         kill HomR;
    190       }
    191     }
    192   }
    193   setring R0;
    194   return(list(fetch(@r,i),p));
    195 }
    196 example
    197 {
    198   "EXAMPLE:"; echo = 2;
    199   option(redSB);
    200 
    201   int p = 181;
    202   intvec a = 2,1,3,4;
    203   intvec b = 1,9,1,1;
    204   ring ra = 0,x(1..4),(a(a),lp);
    205   ideal I = std(cyclic(4));
    206   ring rb = 0,x(1..4),(a(b),lp);
    207   ideal I = imap(ra,I);
    208   modpWalk(I,p,1,a,b);
    209   std(I);
    210 }
    211 
    212 ////////////////////////////////////////////////////////////////////////////////
    213 
    214 proc modWalk(def II, int variant, list #)
    215 "USAGE:  modWalk(II); II ideal or list(ideal,int)
    216 ASSUME:  If variant =
    217 @*       - 1 the Random Walk algorithm with radius II[2] is applied
    218            to II[1] if II = list(ideal, int). It is applied to II with radius 2
    219            if II is an ideal.
    220 @*       - 2, the Groebner Walk algorithm is applied to II[1] or to II, respectively.
    221 @*       - 3, the Fractal Walk algorithm with random element is applied to II[1] or II,
    222            respectively.
    223 @*       - 4, the Fractal Walk algorithm is applied to II[1] or II, respectively.
    224 @*       - 5, the Perturbation Walk algorithm with random element is applied to II[1]
    225            or II, respectively, with radius II[3] and perturbation degree II[2].
    226 @*       - 6, the Perturbation Walk algorithm is applied to II[1] or II, respectively,
    227            with perturbation degree II[3].
    228          If size(#) > 0, then # contains either 1, 2 or 4 integers such that
    229 @*       - #[1] is the number of available processors for the computation,
    230 @*       - #[2] is an optional parameter for the exactness of the computation,
    231                 if #[2] = 1, the procedure computes a standard basis for sure,
    232 @*       - #[3] is the number of primes until the first lifting,
    233 @*       - #[4] is the constant number of primes between two liftings until
    234            the computation stops.
    235 RETURN:  a standard basis of I if no warning appears.
    236 NOTE:    The procedure converts a standard basis of I (over the rational
    237          numbers) from the ordering \"a(v),lp\", "dp\" or \"Dp\" to the ordering
    238          \"(a(w),lp\" or \"a(1,0,...,0),lp\" by using modular methods.
    239          By default the procedure computes a standard basis of I for sure, but
    240          if the optional parameter #[2] = 0, it computes a standard basis of I
    241          with high probability.
    242 EXAMPLE: example modWalk; shows an example
    243 "
    244 {
    245   int TT = timer;
    246   int RT = rtimer;
    247   int i,j,pTest,sizeTest,weighted,n1;
    248   bigint N;
    249 
    250   def R0 = basering;
    251   list rl = ringlist(R0);
    252   if((npars(R0) > 0) || (rl[1] > 0))
    253   {
    254     ERROR("Characteristic of basering should be zero, basering should have no parameters.");
    255   }
    256 
    257   if(typeof(II) == "ideal")
    258   {
    259     ideal I = II;
    260     kill II;
    261     list II;
    262     II[1] = I;
    263     II[2] = 2;
    264     II[3] = 2;
    265   }
    266   else
    267   {
    268     if(typeof(II) == "list" && typeof(II[1]) == "ideal")
    269     {
    270       ideal I = II[1];
    271       if(size(II) == 1)
    272       {
    273         II[2] = 2;
    274         II[3] = 2;
    275       }
    276       if(size(II) == 2)
    277       {
    278         II[3] = 2;
    279       }
    280 
    281     }
    282     else
    283     {
    284       ERROR("Unexpected type of input.");
    285     }
    286   }
    287 
    288 //--------------------  Initialize optional parameters  ------------------------
    289   n1 = system("--cpus");
    290   if(size(#) == 0)
    291   {
    292     int exactness = 1;
    293     int n2 = 10;
    294     int n3 = 10;
    295   }
    296   else
    297   {
    298     if(size(#) == 1)
    299     {
    300       if(typeof(#[1]) == "int")
    301       {
    302         if(#[1] < n1)
    303         {
    304           n1 = #[1];
    305         }
    306         int exactness = 1;
    307         if(n1 >= 10)
    308         {
    309           int n2 = n1 + 1;
    310           int n3 = n1;
    311         }
    312         else
    313         {
    314           int n2 = 10;
    315           int n3 = 10;
    316         }
    317       }
    318       else
    319       {
    320         ERROR("Unexpected type of input.");
    321       }
    322     }
    323     if(size(#) == 2)
    324     {
    325       if(typeof(#[1]) == "int" && typeof(#[2]) == "int")
    326       {
    327         if(#[1] < n1)
    328         {
    329           n1 = #[1];
    330         }
    331         int exactness = #[2];
    332         if(n1 >= 10)
    333         {
    334           int n2 = n1 + 1;
    335           int n3 = n1;
    336         }
    337         else
    338         {
    339           int n2 = 10;
    340           int n3 = 10;
    341         }
    342       }
    343       else
    344       {
    345         if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec")
    346         {
    347           intvec curr_weight = #[1];
    348           intvec target_weight = #[2];
    349           weighted = 1;
    350           int n2 = 10;
    351           int n3 = 10;
    352           int exactness = 1;
    353         }
    354         else
    355         {
    356           ERROR("Unexpected type of input.");
    357         }
    358       }
    359     }
    360     if(size(#) == 3)
    361     {
    362       if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int")
    363       {
    364         intvec curr_weight = #[1];
    365         intvec target_weight = #[2];
    366         weighted = 1;
    367         n1 = #[3];
    368         int n2 = 10;
    369         int n3 = 10;
    370         int exactness = 1;
    371       }
    372       else
    373       {
    374         ERROR("Unexpected type of input.");
    375       }
    376     }
    377     if(size(#) == 4)
    378     {
    379       if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int"
    380           && typeof(#[4]) == "int")
    381       {
    382         intvec curr_weight = #[1];
    383         intvec target_weight = #[2];
    384         weighted = 1;
    385         if(#[1] < n1)
    386         {
    387           n1 = #[3];
    388         }
    389         int exactness = #[4];
    390         if(n1 >= 10)
    391         {
    392           int n2 = n1 + 1;
    393           int n3 = n1;
    394         }
    395         else
    396         {
    397           int n2 = 10;
    398           int n3 = 10;
    399         }
    400       }
    401       else
    402       {
    403         if(typeof(#[1]) == "int" && typeof(#[2]) == "int" && typeof(#[3]) == "int" && typeof(#[4]) == "int")
    404         {
    405           if(#[1] < n1)
    406           {
    407             n1 = #[1];
    408           }
    409           int exactness = #[2];
    410           if(n1 >= #[3])
    411           {
    412             int n2 = n1 + 1;
    413           }
    414           else
    415           {
    416             int n2 = #[3];
    417           }
    418           if(n1 >= #[4])
    419           {
    420             int n3 = n1;
    421           }
    422           else
    423           {
    424             int n3 = #[4];
    425           }
    426         }
    427         else
    428         {
    429           ERROR("Unexpected type of input.");
    430         }
    431       }
    432     }
    433     if(size(#) == 6)
    434     {
    435       if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int" && typeof(#[4]) == "int" && typeof(#[5]) == "int" && typeof(#[6]) == "int")
    436       {
    437         intvec curr_weight = #[1];
    438         intvec target_weight = #[2];
    439         weighted = 1;
    440         if(#[3] < n1)
    441         {
    442           n1 = #[3];
    443         }
    444         int exactness = #[4];
    445         if(n1 >= #[5])
    446         {
    447           int n2 = n1 + 1;
    448         }
    449         else
    450         {
    451           int n2 = #[5];
    452         }
    453         if(n1 >= #[6])
    454         {
    455           int n3 = n1;
    456         }
    457         else
    458         {
    459           int n3 = #[6];
    460         }
    461       }
    462       else
    463       {
    464         ERROR("Expected list(intvec,intvec,int,int,int,int) as optional parameter list.");
    465       }
    466     }
    467     if(size(#) == 1 || size(#) == 5 || size(#) > 6)
    468     {
    469       ERROR("Expected 0,2,3,4 or 5 optional arguments.");
    470     }
    471   }
    472   if(printlevel >= 10)
    473   {
    474   "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)+", exactness = "+string(exactness);
    475   }
    476 
    477 //-------------------------  Save current options  -----------------------------
    478   intvec opt = option(get);
    479   //option(redSB);
    480 
    481 //--------------------  Initialize the list of primes  -------------------------
    482   int tt = timer;
    483   int rt = rtimer;
    484   int en = 2134567879;
    485   int an = 1000000000;
    486   intvec L = primeList(I,n2);
    487   if(n2 > 4)
    488   {
    489   //  L[5] = prime(random(an,en));
    490   }
    491   if(printlevel >= 10)
    492   {
    493     "CPU-time for primeList: "+string(timer-tt)+" seconds.";
    494     "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
    495   }
    496   int h = homog(I);
    497   list P,T1,T2,LL,Arguments,PP;
    498   ideal J,K,H;
    499 
    500 //-------------------  parallelized Groebner Walk in positive characteristic  --------------------
    501 
    502   if(weighted)
    503   {
    504     for(i=1; i<=size(L); i++)
    505     {
    506       Arguments[i] = list(II,L[i],variant,list(curr_weight,target_weight));
    507     }
    508   }
    509   else
    510   {
    511     for(i=1; i<=size(L); i++)
    512     {
    513       Arguments[i] = list(II,L[i],variant);
    514     }
    515   }
    516   P = parallelWaitAll("modpWalk",Arguments);
    517   for(i=1; i<=size(P); i++)
    518   {
    519     T1[i] = P[i][1];
    520     T2[i] = bigint(P[i][2]);
    521   }
    522 
    523   while(1)
    524   {
    525     LL = deleteUnluckyPrimes(T1,T2,h);
    526     T1 = LL[1];
    527     T2 = LL[2];
    528 //-------------------  Now all leading ideals are the same  --------------------
    529 //-------------------  Lift results to basering via farey  ---------------------
    530 
    531     tt = timer; rt = rtimer;
    532     N = T2[1];
    533     for(i=2; i<=size(T2); i++)
    534     {
    535       N = N*T2[i];
    536     }
    537     H = chinrem(T1,T2);
    538     J = parallelFarey(H,N,n1);
    539     //J=farey(H,N);
    540     if(printlevel >= 10)
    541     {
    542       "CPU-time for lifting-process is "+string(timer - tt)+" seconds.";
    543       "Real-time for lifting-process is "+string(rtimer - rt)+" seconds.";
    544     }
    545 
    546 //----------------  Test if we already have a standard basis of I --------------
    547 
    548     tt = timer; rt = rtimer;
    549     pTest = pTestSB(I,J,L,variant);
    550     //pTest = primeTestSB(I,J,L,variant);
    551     if(printlevel >= 10)
    552     {
    553       "CPU-time for pTest is "+string(timer - tt)+" seconds.";
    554       "Real-time for pTest is "+string(rtimer - rt)+" seconds.";
    555     }
    556     if(pTest)
    557     {
    558       if(printlevel >= 10)
    559       {
    560         "CPU-time for computation without final tests is "+string(timer - TT)+" seconds.";
    561         "Real-time for computation without final tests is "+string(rtimer - RT)+" seconds.";
    562       }
    563       attrib(J,"isSB",1);
    564       if(exactness == 0)
    565       {
    566         option(set, opt);
    567         return(J);
    568       }
    569       else
    570       {
    571         tt = timer;
    572         rt = rtimer;
    573         sizeTest = 1 - isIdealIncluded(I,J,n1);
    574         if(printlevel >= 10)
    575         {
    576           "CPU-time for checking if I subset <G> is "+string(timer - tt)+" seconds.";
    577           "Real-time for checking if I subset <G> is "+string(rtimer - rt)+" seconds.";
    578         }
    579         if(sizeTest == 0)
    580         {
    581           tt = timer;
    582           rt = rtimer;
    583           K = std(J);
    584           if(printlevel >= 10)
    585           {
    586             "CPU-time for last std-computation is "+string(timer - tt)+" seconds.";
    587             "Real-time for last std-computation is "+string(rtimer - rt)+" seconds.";
    588           }
    589           if(size(reduce(K,J)) == 0)
    590           {
    591             option(set, opt);
    592             return(J);
    593           }
    594         }
    595       }
    596     }
    597 //--------------  We do not already have a standard basis of I, therefore do the main computation for more primes  --------------
    598 
    599     T1 = H;
    600     T2 = N;
    601     j = size(L)+1;
    602     tt = timer; rt = rtimer;
    603     L = primeList(I,n3,L,n1);
    604     if(printlevel >= 10)
    605     {
    606       "CPU-time for primeList: "+string(timer-tt)+" seconds.";
    607       "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
    608     }
    609     Arguments = list();
    610     PP = list();
    611     if(weighted)
    612     {
    613       for(i=j; i<=size(L); i++)
    614       {
    615         //Arguments[i-j+1] = list(II,L[i],variant,list(curr_weight,target_weight));
    616         Arguments[size(Arguments)+1] = list(II,L[i],variant,list(curr_weight,target_weight));
    617       }
    618     }
    619     else
    620     {
    621       for(i=j; i<=size(L); i++)
    622       {
    623         //Arguments[i-j+1] = list(II,L[i],variant);
    624         Arguments[size(Arguments)+1] = list(II,L[i],variant);
    625       }
    626     }
    627     PP = parallelWaitAll("modpWalk",Arguments);
    628     if(printlevel >= 10)
    629     {
    630       "parallel modpWalk";
    631     }
    632     for(i=1; i<=size(PP); i++)
    633     {
    634       //P[size(P) + 1] = PP[i];
    635       T1[size(T1) + 1] = PP[i][1];
    636       T2[size(T2) + 1] = bigint(PP[i][2]);
    637     }
    638   }
    639   if(printlevel >= 10)
    640   {
    641     "CPU-time for computation with final tests is "+string(timer - TT)+" seconds.";
    642     "Real-time for computation with final tests is "+string(rtimer - RT)+" seconds.";
    643   }
    644 }
    645 
    646 example
    647 {
    648   "EXAMPLE:";
    649   echo = 2;
    650   ring R=0,(x,y,z),lp;
    651   ideal I=-x+y2z-z,xz+1,x2+y2-1;
    652   // I is a standard basis in dp
    653   ideal J = modWalk(I,1);
    654   J;
    655 }
    656 
    657 ////////////////////////////////////////////////////////////////////////////////
    658 static proc isIdealIncluded(ideal I, ideal J, int n1)
    659 "USAGE:  isIdealIncluded(I,J,int n1); I ideal, J ideal, n1 integer
    660 "
    661 {
    662   if(n1 > 1)
    663   {
    664     int k;
    665     list args,results;
    666     for(k=1; k<=size(I); k++)
    667     {
    668       args[k] = list(ideal(I[k]),J,1);
    669     }
    670     results = parallelWaitAll("reduce",args);
    671     for(k=1; k<=size(results); k++)
    672     {
    673       if(results[k] == 0)
    674       {
    675         return(1);
    676       }
    677     }
    678     return(0);
    679   }
    680   else
    681   {
    682     if(reduce(I,J,1) == 0)
    683     {
    684       return(1);
    685     }
    686     else
    687     {
    688       return(0);
    689     }
    690   }
    691 }
    692 
    693 ////////////////////////////////////////////////////////////////////////////////
    694 static proc parallelChinrem(list T1, list T2, int n1)
    695 "USAGE:  parallelChinrem(T1,T2); T1 list of ideals, T2 list of primes, n1 integer"
    696 {
    697   int i,j,k;
    698 
    699   ideal H,J;
    700 
    701   list arguments_chinrem,results_chinrem;
    702   for(i=1; i<=size(T1); i++)
    703   {
    704     J = ideal(T1[i]);
    705     attrib(J,"isSB",1);
    706     arguments_chinrem[size(arguments_chinrem)+1] = list(list(J),T2);
    707   }
    708   results_chinrem = parallelWaitAll("chinrem",arguments_chinrem);
    709     for(j=1; j <= size(results_chinrem); j++)
    710     {
    711       J = results_chinrem[j];
    712       attrib(J,"isSB",1);
    713       if(isIdealIncluded(J,H,n1) == 0)
    714       {
    715         if(H == 0)
    716         {
    717           H = J;
    718         }
    719         else
    720         {
    721           H = H,J;
    722         }
    723       }
    724     }
    725   return(H);
    726 }
    727 
    728 ////////////////////////////////////////////////////////////////////////////////
    729 static proc parallelFarey(ideal H, bigint N, int n1)
    730 "USAGE:  parallelFarey(H,N,n1); H ideal, N bigint, n1 integer
    731 "
    732 {
    733   int i,j;
    734   int ii = 1;
    735   list arguments_farey,results_farey;
    736   for(i=1; i<=size(H); i++)
    737   {
    738     for(j=1; j<=size(H[i]); j++)
    739     {
    740       arguments_farey[size(arguments_farey)+1] = list(H[i][j],N);
    741     }
    742   }
    743   results_farey = parallelWaitAll("farey",arguments_farey);
    744   ideal J,K;
    745   poly f_farey;
    746   while(ii<=size(results_farey))
    747   {
    748     for(i=1; i<=size(H); i++)
    749     {
    750       f_farey = 0;
    751       for(j=1; j<=size(H[i]); j++)
    752       {
    753         f_farey = f_farey + results_farey[ii][1];
    754         ii++;
    755       }
    756       K = ideal(f_farey);
    757       attrib(K,"isSB",1);
    758       attrib(J,"isSB",1);
    759       if(isIdealIncluded(K,J,n1) == 0)
    760       {
    761         if(J == 0)
    762         {
    763           J = K;
    764         }
    765         else
    766         {
    767           J = J,K;
    768         }
    769       }
    770     }
    771   }
    772   return(J);
    773 }
    774 //////////////////////////////////////////////////////////////////////////////////
    775 static proc primeTestSB(def II, ideal J, list L, int variant, list #)
    776 "USAGE:  primeTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int
    777 RETURN:  1 (resp. 0) if for a randomly chosen prime p that is not in L
    778          J mod p is (resp. is not) a standard basis of I mod p
    779 EXAMPLE: example primeTestSB; shows an example
    780 "
    781 {
    782 if(typeof(II) == "ideal")
    783   {
    784   ideal I = II;
    785   int radius = 2;
    786   }
    787 if(typeof(II) == "list")
    788   {
    789   ideal I = II[1];
    790   int radius = II[2];
    791   }
    792 
    793 int i,j,k,p;
    794 def R = basering;
    795 list r = ringlist(R);
    796 
    797 while(!j)
    798   {
    799   j = 1;
    800   p = prime(random(1000000000,2134567879));
    801   for(i = 1; i <= size(L); i++)
    802     {
    803     if(p == L[i])
    804       {
    805       j = 0;
    806       break;
    807       }
    808     }
    809   if(j)
    810     {
    811     for(i = 1; i <= ncols(I); i++)
    812       {
    813       for(k = 2; k <= size(I[i]); k++)
    814         {
    815         if((denominator(leadcoef(I[i][k])) mod p) == 0)
    816           {
    817           j = 0;
    818           break;
    819           }
    820         }
    821       if(!j)
    822         {
    823         break;
    824         }
    825       }
    826     }
    827   if(j)
    828     {
    829     if(!primeTest(I,p))
    830       {
    831       j = 0;
    832       }
    833     }
    834   }
    835 r[1] = p;
    836 def @R = ring(r);
    837 setring @R;
    838 ideal I = imap(R,I);
    839 ideal J = imap(R,J);
    840 attrib(J,"isSB",1);
    841 
    842 int t = timer;
    843 j = 1;
    844 if(isIncluded(I,J) == 0)
    845   {
    846   j = 0;
    847   }
    848 if(printlevel >= 11)
    849   {
    850   "isIncluded(I,J) takes "+string(timer - t)+" seconds";
    851   "j = "+string(j);
    852   }
    853 t = timer;
    854 if(j)
    855   {
    856   if(size(#) > 0)
    857     {
    858     ideal K = modpWalk(I,p,variant,#)[1];
    859     }
    860   else
    861     {
    862     ideal K = modpWalk(I,p,variant)[1];
    863     }
    864   t = timer;
    865   if(isIncluded(J,K) == 0)
    866     {
    867     j = 0;
    868     }
    869   if(printlevel >= 11)
    870     {
    871     "isIncluded(K,J) takes "+string(timer - t)+" seconds";
    872     "j = "+string(j);
    873     }
    874   }
    875 setring R;
    876 
    877 return(j);
    878 }
    879 example
    880 { "EXAMPLE:"; echo = 2;
    881    intvec L = 2,3,5;
    882    ring r = 0,(x,y,z),lp;
    883    ideal I = x+1,x+y+1;
    884    ideal J = x+1,y;
    885    primeTestSB(I,I,L,1);
    886    primeTestSB(I,J,L,1);
    887 }
    888 
    889 ////////////////////////////////////////////////////////////////////////////////
    890 static proc pTestSB(ideal I, ideal J, list L, int variant, list #)
    891 "USAGE:  pTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int
    892 RETURN:  1 (resp. 0) if for a randomly chosen prime p that is not in L
    893          J mod p is (resp. is not) a standard basis of I mod p
    894 EXAMPLE: example pTestSB; shows an example
    895 "
    896 {
    897    int i,j,k,p;
    898    def R = basering;
    899    list r = ringlist(R);
    900 
    901    while(!j)
    902    {
    903       j = 1;
    904       p = prime(random(1000000000,2134567879));
    905       for(i = 1; i <= size(L); i++)
    906       {
    907          if(p == L[i]) { j = 0; break; }
    908       }
    909       if(j)
    910       {
    911          for(i = 1; i <= ncols(I); i++)
    912          {
    913             for(k = 2; k <= size(I[i]); k++)
    914             {
    915                if((denominator(leadcoef(I[i][k])) mod p) == 0) { j = 0; break; }
    916             }
    917             if(!j){ break; }
    918          }
    919       }
    920       if(j)
    921       {
    922          if(!primeTest(I,p)) { j = 0; }
    923       }
    924    }
    925    r[1] = p;
    926    def @R = ring(r);
    927    setring @R;
    928    ideal I = imap(R,I);
    929    ideal J = imap(R,J);
    930    attrib(J,"isSB",1);
    931 
    932    int t = timer;
    933    j = 1;
    934    if(isIncluded(I,J) == 0) { j = 0; }
    935 
    936    if(printlevel >= 11)
    937    {
    938       "isIncluded(I,J) takes "+string(timer - t)+" seconds";
    939       "j = "+string(j);
    940    }
    941 
    942    t = timer;
    943    if(j)
    944    {
    945       if(size(#) > 0)
    946       {
    947          ideal K = modpStd(I,p,variant,#[1])[1];
    948       }
    949       else
    950       {
    951          ideal K = groebner(I);
    952       }
    953       t = timer;
    954       if(isIncluded(J,K) == 0) { j = 0; }
    955 
    956       if(printlevel >= 11)
    957       {
    958          "isIncluded(J,K) takes "+string(timer - t)+" seconds";
    959          "j = "+string(j);
    960       }
    961    }
    962    setring R;
    963    return(j);
    964 }
    965 example
    966 { "EXAMPLE:"; echo = 2;
    967    intvec L = 2,3,5;
    968    ring r = 0,(x,y,z),dp;
    969    ideal I = x+1,x+y+1;
    970    ideal J = x+1,y;
    971    pTestSB(I,I,L,2);
    972    pTestSB(I,J,L,2);
    973 }
    974 ////////////////////////////////////////////////////////////////////////////////
    975 static proc mixedTest()
    976 "USAGE:  mixedTest();
    977 RETURN:  1 if ordering of basering is mixed, 0 else
    978 EXAMPLE: example mixedTest(); shows an example
    979 "
    980 {
    981    int i,p,m;
    982    for(i = 1; i <= nvars(basering); i++)
    983    {
    984       if(var(i) > 1)
    985       {
    986          p++;
    987       }
    988       else
    989       {
    990          m++;
    991       }
    992    }
    993    if((p > 0) && (m > 0)) { return(1); }
    994    return(0);
    995 }
    996 example
    997 { "EXAMPLE:"; echo = 2;
    998    ring R1 = 0,(x,y,z),dp;
    999    mixedTest();
    1000    ring R2 = 31,(x(1..4),y(1..3)),(ds(4),lp(3));
    1001    mixedTest();
    1002    ring R3 = 181,x(1..9),(dp(5),lp(4));
    1003    mixedTest();
    1004 }
     139    int pertdeg = 3;
     140    ideal J = modrWalk(I,radius,pertdeg);
     141    J;
     142}
     143
     144proc modfWalk(ideal I, list #)
     145"USAGE:   modfWalk(I, [, v, w]); I ideal, v intvec, w intvec
     146RETURN:   a standard basis of I
     147NOTE:     The procedure computes a standard basis of I (over the rational
     148          numbers) by using modular methods.
     149SEE ALSO: modular
     150EXAMPLE:  example modfWalk; shows an example"
     151{
     152    /* read optional parameter */
     153    if (size(#) > 0) {
     154        if (size(#) == 1) {
     155            intvec w = #[1];
     156        }
     157        if (size(#) == 2) {
     158            intvec v = #[1];
     159            intvec w = #[2];
     160        }
     161        if (size(#) > 2 || typeof(#[1]) != "intvec") {
     162            ERROR("wrong optional parameter");
     163        }
     164    }
     165
     166    /* save options */
     167    intvec opt = option(get);
     168    option(redSB);
     169
     170    /* set additional parameters */
     171    int reduction = 1;
     172    int printout = 0;
     173
     174    /* call modular() */
     175    if (size(#) > 0) {
     176        I = modular("fwalk", list(I,reduction,printout,#));
     177    }
     178    else {
     179        I = modular("fwalk", list(I,reduction,printout));
     180    }
     181
     182    /* return the result */
     183    attrib(I, "isSB", 1);
     184    option(set, opt);
     185    return(I);
     186}
     187example
     188{
     189    "EXAMPLE:";
     190    echo = 2;
     191    ring R1 = 0, (x,y,z,t), dp;
     192    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
     193    I = std(I);
     194    ring R2 = 0, (x,y,z,t), lp;
     195    ideal I = fetch(R1, I);
     196    ideal J = modfWalk(I);
     197    J;
     198}
     199
     200proc modfrWalk(ideal I, int radius, list #)
     201"USAGE:   modfrWalk(I, radius [, v, w]); I ideal, radius int, v intvec, w intvec
     202RETURN:   a standard basis of I
     203NOTE:     The procedure computes a standard basis of I (over the rational
     204          numbers) by using modular methods.
     205SEE ALSO: modular
     206EXAMPLE:  example modfrWalk; shows an example"
     207{
     208    /* read optional parameter */
     209    if (size(#) > 0) {
     210        if (size(#) == 1) {
     211            intvec w = #[1];
     212        }
     213        if (size(#) == 2) {
     214            intvec v = #[1];
     215            intvec w = #[2];
     216        }
     217        if (size(#) > 2 || typeof(#[1]) != "intvec") {
     218            ERROR("wrong optional parameter");
     219        }
     220    }
     221
     222    /* save options */
     223    intvec opt = option(get);
     224    option(redSB);
     225
     226    /* set additional parameters */
     227    int reduction = 1;
     228    int printout = 0;
     229
     230    /* call modular() */
     231    if (size(#) > 0) {
     232        I = modular("frandwalk", list(I,radius,reduction,printout,#));
     233    }
     234    else {
     235        I = modular("frandwalk", list(I,radius,reduction,printout));
     236    }
     237
     238    /* return the result */
     239    attrib(I, "isSB", 1);
     240    option(set, opt);
     241    return(I);
     242}
     243example
     244{
     245    "EXAMPLE:";
     246    echo = 2;
     247    ring R1 = 0, (x,y,z,t), dp;
     248    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
     249    I = std(I);
     250    ring R2 = 0, (x,y,z,t), lp;
     251    ideal I = fetch(R1, I);
     252    int radius = 2;
     253    ideal J = modfrWalk(I,radius);
     254    J;
     255}
  • Singular/LIB/rwalk.lib

    • Property mode changed from 100644 to 100755
    rf34cd44 rbbfc4a8  
    1010rwalk(ideal,int,int[,intvec,intvec]);   standard basis of ideal via Random Walk algorithm
    1111rwalk(ideal,int[,intvec,intvec]);       standard basis of ideal via Random Perturbation Walk algorithm
    12 frwalk(ideal,int[,intvec,intvec]);      standard basis of ideal via Random Fractal Walk algorithm
     12frandwalk(ideal,int[,intvec,intvec]);      standard basis of ideal via Random Fractal Walk algorithm
    1313";
    1414
     
    141141 * Random Walk  *
    142142 ****************/
    143 proc rwalk(ideal Go, int radius, int pert_deg, list #)
     143proc rwalk(ideal Go, int radius, int pert_deg, int reduction, int printout, list #)
    144144"SYNTAX: rwalk(ideal i, int radius);
    145145         if size(#)>0 then rwalk(ideal i, int radius, intvec v, intvec w);
     146         intermediate Groebner bases are not reduced if reduction = 0
    146147TYPE:    ideal
    147148PURPOSE: compute the standard basis of the ideal, calculated via
     
    178179
    179180ideal G = fetch(xR, Go);
    180 G = system("Mrwalk", G, curr_weight, target_weight, radius, pert_deg, basering);
     181G = system("Mrwalk", G, curr_weight, target_weight, radius, pert_deg, reduction, printout);
    181182
    182183setring xR;
     
    196197  int radius = 1;
    197198  int perturb_deg = 2;
    198   rwalk(I,radius,perturb_deg);
     199  int reduction = 0;
     200  int printout = 1;
     201  rwalk(I,radius,perturb_deg,reduction,printout);
    199202}
    200203
     
    202205 * Perturbation Walk with random element *
    203206 *****************************************/
    204 proc prwalk(ideal Go, int radius, int o_pert_deg, int t_pert_deg, list #)
     207proc prwalk(ideal Go, int radius, int o_pert_deg, int t_pert_deg, int reduction, int printout, list #)
    205208"SYNTAX: rwalk(ideal i, int radius);
    206209         if size(#)>0 then rwalk(ideal i, int radius, intvec v, intvec w);
     
    227230  OSCTW= OrderStringalp_NP("al", #);
    228231  }
     232int nP = OSCTW[1];
    229233string ord_str = OSCTW[2];
    230234intvec curr_weight = OSCTW[3]; // original weight vector
     
    238242
    239243ideal G = fetch(xR, Go);
    240 G = system("Mprwalk", G, curr_weight, target_weight, radius, o_pert_deg, t_pert_deg, basering);
     244G = system("Mprwalk", G, curr_weight, target_weight, radius, o_pert_deg, t_pert_deg,
     245           nP, reduction, printout);
    241246
    242247setring xR;
     
    257262  int o_perturb_deg = 2;
    258263  int t_perturb_deg = 2;
    259   prwalk(I,radius,o_perturb_deg,t_perturb_deg);
     264  int reduction = 0;
     265  int printout = 2;
     266  prwalk(I,radius,o_perturb_deg,t_perturb_deg,reduction,printout);
    260267}
    261268
     
    263270 * Fractal Walk with random element *
    264271 ************************************/
    265 proc frandwalk(ideal Go, int radius, list #)
    266 "SYNTAX: frwalk(ideal i, int radius);
    267          frwalk(ideal i, int radius, intvec v, intvec w);
     272proc frandwalk(ideal Go, int radius, int reduction, int printout, list #)
     273"SYNTAX: frwalk(ideal i, int radius, int reduction, int printout);
     274         frwalk(ideal i, int radius, int reduction, int printout, intvec v, intvec w);
    268275TYPE:    ideal
    269276PURPOSE: compute the standard basis of the ideal, calculated via
     
    299306   ideal G = fetch(xR, Go);
    300307   int pert_deg = 2;
    301    G = system("Mfrwalk", G, curr_weight, target_weight, radius);
     308
     309   G = system("Mfrwalk", G, curr_weight, target_weight, radius, reduction, printout);
    302310
    303311   setring xR;
     
    314322    ring r = 0,(z,y,x), lp;
    315323    ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
    316     frandwalk(I,2);
    317 }
     324    int reduction = 0;
     325    frandwalk(I,2,0,1);
     326}
  • Singular/LIB/swalk.lib

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

    rf34cd44 rbbfc4a8  
    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

    r500e4b rbbfc4a8  
    8383      {
    8484        assume( nMap != NULL );
    85 
    8685        number a = nMap((number)data, preimage_r->cf, currRing->cf);
    8786        if (nCoeff_is_Extension(currRing->cf))
     
    146145        }
    147146      }
    148       else
    149       if ( (what==IMAP_CMD) || /*(*/ (what==FETCH_CMD) /*)*/) /* && (nMap!=nCopy)*/
     147      else if ((what==IMAP_CMD) || (what==FETCH_CMD))
    150148      {
    151149        for (i=R*C-1;i>=0;i--)
     
    156154        }
    157155      }
    158       else /* if(what==MAP_CMD) */
    159       {
     156      else /* (what==MAP_CMD) */
     157      {
     158        assume(what==MAP_CMD);
    160159        matrix s=mpNew(N,maMaxDeg_Ma((ideal)data,preimage_r));
    161160        for (i=R*C-1;i>=0;i--)
  • Singular/walk.cc

    • Property mode changed from 100644 to 100755
    rf34cd44 rbbfc4a8  
    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;
    4297   intvec* next_weight2;
     4360  assume(currRing != NULL && curr_weight != NULL &&
     4361         target_weight != NULL && G->m[0] != NULL);
     4362
     4363  intvec* next_weight1 = MkInterRedNextWeight(curr_weight,target_weight,G);
     4364  if(weight_rad == 0)
     4365  {
     4366    return(next_weight1);
     4367  }
     4368
     4369  int i,weight_norm,nV = currRing->N;
     4370  intvec* next_weight2 = new intvec(nV);
    42984371  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)
     4372  intvec* result = new intvec(nV);
     4373
     4374  ideal G_test, G_test2;
     4375
     4376  //compute a random next weight vector "next_weight2"
     4377  while(1)
     4378  {
     4379    weight_norm = 0;
     4380    while(weight_norm == 0)
     4381    {
     4382
     4383      for(i=0; i<nV; i++)
     4384      {
     4385        (*next_weight2)[i] = rand() % 60000 - 30000;
     4386        weight_norm = weight_norm + (*next_weight2)[i]*(*next_weight2)[i];
     4387      }
     4388      weight_norm = 1 + floor(sqrt(weight_norm));
     4389    }
     4390
     4391    for(i=0; i<nV; i++)
     4392    {
     4393      if((*next_weight2)[i] < 0)
     4394      {
     4395        (*next_weight2)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight2)[i]/weight_norm);
     4396      }
     4397      else
     4398      {
     4399        (*next_weight2)[i] = (*curr_weight)[i] + floor(weight_rad*(*next_weight2)[i]/weight_norm);
     4400      }
     4401    }
     4402   
     4403    if(test_w_in_ConeCC(G, next_weight2) == 1)
     4404    {
     4405      PrintS("\n Hier 1\n");
     4406      G_test2 = MwalkInitialForm(G, next_weight2);
     4407      PrintS("\n Hier 2\n");
     4408      if(maxlengthpoly(G_test2) <= 1)
     4409      {
     4410        next_weight2 = MkInterRedNextWeight(next_weight2,target_weight,G);
     4411      }
     4412      idDelete(&G_test2);
     4413
     4414      if(MivAbsMax(next_weight2)>1147483647)
    43154415      {
    43164416        for(i=0; i<nV; i++)
    43174417        {
    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];
     4418          (*next_weight22)[i] = (*next_weight2)[i];
    43214419        }
    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)
     4420        i = 0;
     4421        // reduce the size of the maximal entry of the vector
     4422        while(test_w_in_ConeCC(G,next_weight22) == 1)
    43284423        {
    4329           (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm);
     4424          (*next_weight2)[i] = (*next_weight22)[i];
     4425          i = MivAbsMaxArg(next_weight22);
     4426          (*next_weight22)[i] = floor(0.1*(*next_weight22)[i] + 0.5);
    43304427        }
    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);
    43424428        delete next_weight22;
    4343         break;
    4344       }
    4345     }
    4346     intvec* result = new intvec(nV);
    4347     ideal G_test = MwalkInitialForm(G, next_weight);
     4429      }
     4430      break;
     4431    }
     4432    weight_rad++;
     4433  }
     4434 
     4435  // compute "usual" next weight vector
     4436  intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight, G);
     4437  G_test = MwalkInitialForm(G, next_weight);
     4438  G_test2 = MwalkInitialForm(G, next_weight2);
     4439
     4440  // compare next weights
     4441  if(Overflow_Error == FALSE)
     4442  {
    43484443    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|
     4444    if(G_test1->m[0] != NULL && maxlengthpoly(G_test1) < maxlengthpoly(G_test))
     4445    {
     4446      if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test1))
    43554447      {
    43564448        for(i=0; i<nV; i++)
     
    43594451        }
    43604452      }
    4361       else // |G_test1| < |G_test|, |G_test1| < |G_test2|
     4453      else
    43624454      {
    43634455        for(i=0; i<nV; i++)
     
    43654457          (*result)[i] = (*next_weight1)[i];
    43664458        }
    4367       }
     4459      }   
    43684460    }
    43694461    else
    43704462    {
    4371       if(IDELEMS(G_test2) <= IDELEMS(G_test)) // |G_test2| <= |G_test| <= |G_test1|
     4463      if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test))
    43724464      {
    43734465        for(i=0; i<nV; i++)
     
    43764468        }
    43774469      }
    4378       else // |G_test| <= |G_test1|, |G_test| < |G_test2|
     4470      else
    43794471      {
    43804472        for(i=0; i<nV; i++)
     
    43844476      }
    43854477    }
    4386     delete next_weight;
    4387     delete next_weight1;
    4388     idDelete(&G_test);
    43894478    idDelete(&G_test1);
    4390     idDelete(&G_test2);
    4391     if(test_w_in_ConeCC(G, result) == 1)
    4392     {
    4393       delete next_weight2;
    4394       return result;
     4479  }
     4480  else
     4481  {
     4482    Overflow_Error = FALSE;
     4483    if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test))
     4484    {
     4485      for(i=1; i<nV; i++)
     4486      {
     4487        (*result)[i] = (*next_weight2)[i];
     4488      }
    43954489    }
    43964490    else
    43974491    {
    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     {
    44264492      for(i=0; i<nV; i++)
    44274493      {
    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       {
    45134494        (*result)[i] = (*next_weight)[i];
    45144495      }
    45154496    }
    45164497  }
     4498
    45174499  idDelete(&G_test);
    45184500  idDelete(&G_test2);
     
    45754557    else
    45764558    {
    4577       rChangeCurrRing(VMrDefault(orig_target_weight)); // Aenderung
     4559      rChangeCurrRing(VMrDefault(orig_target_weight));
    45784560    }
    45794561    TargetRing = currRing;
     
    46464628    else
    46474629    {
    4648       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung
     4630      rChangeCurrRing(VMrDefault(curr_weight));
    46494631    }
    46504632    newRing = currRing;
     
    47554737    else
    47564738    {
    4757       rChangeCurrRing(VMrDefault(orig_target_weight)); // Aenderung
     4739      rChangeCurrRing(VMrDefault(orig_target_weight));
    47584740    }
    47594741    F1 = idrMoveR(G, newRing,currRing);
     
    47864768      else
    47874769      {
    4788         rChangeCurrRing(VMrDefault(orig_target_weight)); // Aenderung
     4770        rChangeCurrRing(VMrDefault(orig_target_weight));
    47894771      }
    47904772    KSTD_Finish:
     
    48844866      tim = clock();
    48854867      /*
    4886         Print("\n// **** Grï¿œbnerwalk took %d steps and ", nwalk);
     4868        Print("\n// **** Groebnerwalk took %d steps and ", nwalk);
    48874869        PrintS("\n// **** call the rec. Pert. Walk to compute a red GB of:");
    48884870        idElements(Gomega, "G_omega");
     
    49144896      oldRing = currRing;
    49154897
    4916       /* create a new ring newRing */
     4898      // create a new ring newRing
    49174899       if (rParameter(currRing) != NULL)
    49184900       {
     
    49214903       else
    49224904       {
    4923          rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung
     4905         rChangeCurrRing(VMrDefault(curr_weight));
    49244906       }
    49254907      newRing = currRing;
     
    49474929      else
    49484930      {
    4949         rChangeCurrRing(VMrDefault(curr_weight));  //Aenderung
     4931        rChangeCurrRing(VMrDefault(curr_weight));
    49504932      }
    49514933      newRing = currRing;
     
    49594941      M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
    49604942      delete hilb_func;
    4961 #endif // BUCHBERGER_ALG
     4943#endif
    49624944      tstd = tstd + clock() - to;
    49634945
     
    49684950
    49694951      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.
     4952      // compute a representation of the generators of submod (M) with respect
     4953      // to those of mod (Gomega).
     4954      // Gomega is a reduced Groebner basis w.r.t. the current ring.
    49714955      F = MLifttwoIdeal(Gomega2, M1, G);
    49724956      tlift = tlift + clock() - to;
     
    50185002      else
    50195003      {
    5020         rChangeCurrRing(VMrDefault(target_weight)); // Aenderung
     5004        rChangeCurrRing(VMrDefault(target_weight));
    50215005      }
    50225006      F1 = idrMoveR(G, newRing,currRing);
     
    50655049 * THE GROEBNER WALK ALGORITHM *
    50665050 *******************************/
    5067 ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M, ring baseRing)
     5051ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M,
     5052            ring baseRing, int reduction, int printout)
    50685053{
    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));
     5054  // save current options
     5055  BITSET save1 = si_opt_1;
     5056  if(reduction == 0)
     5057  {
     5058    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     5059    si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     5060  }
    50735061  Set_Error(FALSE);
    50745062  Overflow_Error = FALSE;
     5063  //BOOLEAN endwalks = FALSE;
    50755064#ifdef TIME_TEST
    50765065  clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
     
    50805069#endif
    50815070  nstep=0;
    5082   int i,nwalk,endwalks = 0;
     5071  int i,nwalk;
    50835072  int nV = baseRing->N;
    50845073
    5085   ideal Gomega, M, F, Gomega1, Gomega2, M1; //, F1;
     5074  ideal Gomega, M, F, FF, Gomega1, Gomega2, M1;
    50865075  ring newRing;
    50875076  ring XXRing = baseRing;
     5077  ring targetRing;
    50885078  intvec* ivNull = new intvec(nV);
    50895079  intvec* curr_weight = new intvec(nV);
    50905080  intvec* target_weight = new intvec(nV);
    50915081  intvec* exivlp = Mivlp(nV);
     5082/*
    50925083  intvec* tmp_weight = new intvec(nV);
    50935084  for(i=0; i<nV; i++)
    50945085  {
    5095     (*tmp_weight)[i] = (*target_M)[i];
    5096   }
     5086    (*tmp_weight)[i] = (*orig_M)[i];
     5087  }
     5088*/
    50975089  for(i=0; i<nV; i++)
    50985090  {
     
    51115103#endif
    51125104  rComplete(currRing);
    5113 #ifdef CHECK_IDEAL_MWALK
    5114     idString(Go,"Go");
    5115 #endif
     5105//#ifdef CHECK_IDEAL_MWALK
     5106  if(printout > 2)
     5107  {
     5108    idString(Go,"//** Mwalk: Go");
     5109  }
     5110//#endif
     5111
     5112  if(target_M->length() == nV)
     5113  {
     5114   // define the target ring
     5115    targetRing = VMrDefault(target_weight);
     5116  }
     5117  else
     5118  {
     5119    targetRing = VMatrDefault(target_M);
     5120  }
     5121  if(orig_M->length() == nV)
     5122  {
     5123    // define a new ring with ordering "(a(curr_weight),lp)
     5124    newRing = VMrDefault(curr_weight);
     5125  }
     5126  else
     5127  {
     5128    newRing = VMatrDefault(orig_M);
     5129  }
     5130  rChangeCurrRing(newRing);
     5131
    51165132#ifdef TIME_TEST
    51175133  to = clock();
    51185134#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);
     5135
    51285136  ideal G = MstdCC(idrMoveR(Go,baseRing,currRing));
    5129   baseRing = currRing;
     5137
    51305138#ifdef TIME_TEST
    51315139  tostd = clock()-to;
    51325140#endif
    51335141
     5142  baseRing = currRing;
    51345143  nwalk = 0;
     5144
    51355145  while(1)
    51365146  {
    51375147    nwalk ++;
    51385148    nstep ++;
     5149
    51395150#ifdef TIME_TEST
    51405151    to = clock();
    51415152#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"
     5153    // compute an initial form ideal of <G> w.r.t. "curr_vector"
     5154    Gomega = MwalkInitialForm(G, curr_weight);
    51465155#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
     5156    tif = tif + clock()-to;
     5157#endif
     5158
     5159//#ifdef CHECK_IDEAL_MWALK
     5160    if(printout > 1)
     5161    {
     5162      idString(Gomega,"//** Mwalk: Gomega");
     5163    }
     5164//#endif
     5165
     5166    if(reduction == 0)
     5167    {
     5168      FF = middleOfCone(G,Gomega);
     5169      if( FF != NULL)
     5170      {
     5171        idDelete(&G);
     5172        G = idCopy(FF);
     5173        idDelete(&FF);
     5174        goto NEXT_VECTOR;
     5175      }
     5176    }
     5177
    51525178#ifndef  BUCHBERGER_ALG
    51535179    if(isNolVector(curr_weight) == 0)
    51545180    {
    51555181      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
    5156     }
     5182    }   
    51575183    else
    51585184    {
     
    51605186    }
    51615187#endif
     5188
    51625189    if(nwalk == 1)
    51635190    {
    51645191      if(orig_M->length() == nV)
    51655192      {
    5166         newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
     5193        // define a new ring with ordering "(a(curr_weight),lp)
     5194        newRing = VMrDefault(curr_weight);
    51675195      }
    51685196      else
     
    51755203     if(target_M->length() == nV)
    51765204     {
    5177        newRing = VMrRefine(curr_weight,target_weight); //define a new ring with ordering "(a(curr_weight),Wp(target_weight))"
     5205       //define a new ring with ordering "(a(curr_weight),lp)"
     5206       newRing = VMrDefault(curr_weight);
    51785207     }
    51795208     else
    51805209     {
     5210       //define a new ring with matrix ordering
    51815211       newRing = VMatrRefine(target_M,curr_weight);
    51825212     }
     
    51995229#endif
    52005230    idSkipZeroes(M);
    5201 #ifdef CHECK_IDEAL_MWALK
    5202     PrintS("\n//** Mwalk: computed M.\n");
    5203     idString(M, "M");
    5204 #endif
     5231//#ifdef CHECK_IDEAL_MWALK
     5232    if(printout > 2)
     5233    {
     5234      idString(M, "//** Mwalk: M");
     5235    }
     5236//#endif
    52055237    //change the ring to baseRing
    52065238    rChangeCurrRing(baseRing);
     
    52125244    to = clock();
    52135245#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
     5246    // compute a representation of the generators of submod (M) with respect to those of mod (Gomega),
     5247    // where Gomega is a reduced Groebner basis w.r.t. the current ring
    52155248    F = MLifttwoIdeal(Gomega2, M1, G);
     5249
    52165250#ifdef TIME_TEST
    52175251    tlift = tlift + clock() - to;
    52185252#endif
    5219 #ifdef CHECK_IDEAL_MWALK
    5220     idString(F, "F");
    5221 #endif
     5253//#ifdef CHECK_IDEAL_MWALK
     5254    if(printout > 2)
     5255    {
     5256      idString(F, "//** Mwalk: F");
     5257    }
     5258//#endif
    52225259    idDelete(&Gomega2);
    52235260    idDelete(&M1);
     5261
    52245262    rChangeCurrRing(newRing); // change the ring to newRing
    52255263    G = idrMoveR(F,baseRing,currRing);
    52265264    idDelete(&F);
     5265    idSkipZeroes(G);
     5266
     5267//#ifdef CHECK_IDEAL_MWALK
     5268    if(printout > 2)
     5269    {
     5270      idString(G, "//** Mwalk: G");
     5271    }
     5272//#endif
     5273
     5274    rChangeCurrRing(targetRing);
     5275    G = idrMoveR(G,newRing,currRing);
     5276    // test whether target cone is reached
     5277    if(reduction !=0 && test_w_in_ConeCC(G,curr_weight) == 1)
     5278    {
     5279      baseRing = currRing;
     5280      break;
     5281      //endwalks = TRUE;
     5282    }
     5283
     5284    rChangeCurrRing(newRing);
     5285    G = idrMoveR(G,targetRing,currRing);
    52275286    baseRing = currRing;
     5287
     5288/*
    52285289#ifdef TIME_TEST
    52295290    to = clock();
    52305291#endif
    5231     //G = kStd(F1,NULL,testHomog,NULL,NULL,0,0,NULL);
     5292
    52325293#ifdef TIME_TEST
    52335294    tstd = tstd + clock() - to;
    52345295#endif
    5235     idSkipZeroes(G);
    5236 #ifdef CHECK_IDEAL_MWALK
    5237     idString(G, "G");
    5238 #endif
     5296*/
     5297
     5298
    52395299#ifdef TIME_TEST
    52405300    to = clock();
    52415301#endif
     5302    NEXT_VECTOR:
    52425303    intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
    52435304#ifdef TIME_TEST
    52445305    tnw = tnw + clock() - to;
    52455306#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
     5307//#ifdef PRINT_VECTORS
     5308    if(printout > 0)
     5309    {
     5310      MivString(curr_weight, target_weight, next_weight);
     5311    }
     5312//#endif
     5313    if(MivComp(target_weight,curr_weight) == 1)// || endwalks == TRUE)
     5314    {/*
     5315//#ifdef CHECK_IDEAL_MWALK
     5316      if(printout > 0)
     5317      {
     5318        PrintS("\n//** Mwalk: entering last cone.\n");
     5319      }
     5320//#endif
     5321
    52545322      Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
    52555323      if(target_M->length() == nV)
     
    52645332      Gomega1 = idrMoveR(Gomega, baseRing,currRing);
    52655333      idDelete(&Gomega);
    5266 #ifdef CHECK_IDEAL_MWALK
    5267       idString(Gomega1, "Gomega");
    5268 #endif
     5334//#ifdef CHECK_IDEAL_MWALK
     5335      if(printout > 1)
     5336      {
     5337        idString(Gomega1, "//** Mwalk: Gomega");
     5338      }
     5339      PrintS("\n //** Mwalk: kStd(Gomega)");
     5340//#endif
    52695341      M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL);
    5270 #ifdef CHECK_IDEAL_MWALK
    5271       idString(M,"M");
    5272 #endif
     5342//#ifdef CHECK_IDEAL_MWALK
     5343      if(printout > 1)
     5344      {
     5345        idString(M,"//** Mwalk: M");
     5346      }
     5347//#endif
    52735348      rChangeCurrRing(baseRing);
    52745349      M1 =  idrMoveR(M, newRing,currRing);
     
    52765351      Gomega2 = idrMoveR(Gomega1, newRing,currRing);
    52775352      idDelete(&Gomega1);
     5353      //PrintS("\n //** Mwalk: MLifttwoIdeal");
    52785354      F = MLifttwoIdeal(Gomega2, M1, G);
    5279 #ifdef CHECK_IDEAL_MWALK
    5280       idString(F,"F");
    5281 #endif
     5355//#ifdef CHECK_IDEAL_MWALK
     5356      if(printout > 2)
     5357      {
     5358        idString(F,"//** Mwalk: F");
     5359      }
     5360//#endif
    52825361      idDelete(&Gomega2);
    52835362      idDelete(&M1);
     
    52915370      to = clock();
    52925371#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   //    }
     5372      //PrintS("\n //**Mwalk: Interreduce");
     5373      //interreduce the Groebner basis <G> w.r.t. currRing
     5374      //G = kInterRedCC(G,NULL);
    52975375#ifdef TIME_TEST
    52985376      tred = tred + clock() - to;
    52995377#endif
    53005378      idSkipZeroes(G);
    5301       delete next_weight;
     5379      delete next_weight; */
    53025380      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
     5381    }
     5382
    53105383    for(i=nV-1; i>=0; i--)
    53115384    {
    5312       (*tmp_weight)[i] = (*curr_weight)[i];
     5385      //(*tmp_weight)[i] = (*curr_weight)[i];
    53135386      (*curr_weight)[i] = (*next_weight)[i];
    53145387    }
     
    53175390  rChangeCurrRing(XXRing);
    53185391  ideal result = idrMoveR(G,baseRing,currRing);
     5392  idDelete(&Go);
    53195393  idDelete(&G);
    5320 /*#ifdef CHECK_IDEAL_MWALK
    5321   pDelete(&p);
    5322 #endif*/
    5323   delete tmp_weight;
     5394  //delete tmp_weight;
    53245395  delete ivNull;
    53255396  delete exivlp;
     
    53285399#endif
    53295400#ifdef TIME_TEST
    5330   Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep);
    53315401  TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep);
    5332   Print("\n//** Mwalk: Ergebnis.\n");
    53335402  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
    53345403  //Print("\n// Overflow_Error? (%d)\n", Overflow_Error);
    53355404#endif
     5405  Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep);
    53365406  return(result);
    53375407}
    53385408
    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)
     5409// THE RANDOM WALK ALGORITHM
     5410ideal Mrwalk(ideal Go, intvec* orig_M, intvec* target_M, int weight_rad, int pert_deg,
     5411             int reduction, int printout)
    53425412{
    53435413  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));
     5414  if(reduction == 0)
     5415  {
     5416    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     5417    si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     5418  }
    53475419  Set_Error(FALSE);
    53485420  Overflow_Error = FALSE;
     5421  BOOLEAN endwalks = FALSE;
    53495422#ifdef TIME_TEST
    53505423  clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
     
    53545427#endif
    53555428  nstep=0;
    5356   int i,nwalk,endwalks = 0;
    5357   int nV = baseRing->N;
    5358 
    5359   ideal Gomega, M, F, Gomega1, Gomega2, M1; //, F1;
     5429  int i,polylength,nwalk;
     5430  int nV = currRing->N;
     5431
     5432  ideal Gomega, M, F,FF, Gomega1, Gomega2, M1;
    53605433  ring newRing;
    5361   ring XXRing = baseRing;
     5434  ring targetRing;
     5435  ring baseRing = currRing;
     5436  ring XXRing = currRing;
    53625437  intvec* ivNull = new intvec(nV);
    53635438  intvec* curr_weight = new intvec(nV);
    53645439  intvec* target_weight = new intvec(nV);
    5365   intvec* exivlp = Mivlp(nV);
    5366   intvec* tmp_weight = new intvec(nV);
    5367   for(i=0; i<nV; i++)
    5368   {
    5369     (*tmp_weight)[i] = (*target_M)[i];
    5370   }
     5440  intvec* next_weight= new intvec(nV);
     5441
    53715442  for(i=0; i<nV; i++)
    53725443  {
     
    53855456#endif
    53865457  rComplete(currRing);
    5387 #ifdef CHECK_IDEAL_MWALK
    5388     idString(Go,"Go");
    5389 #endif
    53905458#ifdef TIME_TEST
    53915459  to = clock();
    53925460#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       }
     5461
     5462  if(target_M->length() == nV)
     5463  {
     5464   // define the target ring
     5465    targetRing = VMrDefault(target_weight);
     5466  }
     5467  else
     5468  {
     5469    targetRing = VMatrDefault(target_M);
     5470  }
     5471  if(orig_M->length() == nV)
     5472  {
     5473    newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
     5474  }
     5475  else
     5476  {
     5477    newRing = VMatrDefault(orig_M);
     5478  }
    54015479  rChangeCurrRing(newRing);
    54025480  ideal G = MstdCC(idrMoveR(Go,baseRing,currRing));
     
    54075485
    54085486  nwalk = 0;
     5487
     5488#ifdef TIME_TEST
     5489  to = clock();
     5490#endif
     5491  Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
     5492#ifdef TIME_TEST
     5493  tif = tif + clock()-to; //time for computing initial form ideal
     5494#endif
     5495
    54095496  while(1)
    54105497  {
    54115498    nwalk ++;
    54125499    nstep ++;
    5413 #ifdef TIME_TEST
    5414     to = clock();
    5415 #endif
    5416 #ifdef CHECK_IDEAL_MWALK
    5417     idString(G,"G");
    5418 #endif
    5419     Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
    5420 #ifdef TIME_TEST
    5421     tif = tif + clock()-to; //time for computing initial form ideal
    5422 #endif
    5423 #ifdef CHECK_IDEAL_MWALK
    5424     idString(Gomega,"Gomega");
    5425 #endif
     5500    if(printout > 1)
     5501    {
     5502      idString(Gomega,"//** Mrwalk: Gomega");
     5503    }
     5504    if(reduction == 0)
     5505    {
     5506      FF = middleOfCone(G,Gomega);
     5507      if(FF != NULL)
     5508      {
     5509        idDelete(&G);
     5510        G = idCopy(FF);
     5511        idDelete(&FF);
     5512       
     5513        goto NEXT_VECTOR;
     5514      }
     5515    }
    54265516#ifndef  BUCHBERGER_ALG
    54275517    if(isNolVector(curr_weight) == 0)
    54285518    {
    54295519      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
    5430     }
     5520    }   
    54315521    else
    54325522    {
     
    54495539     if(target_M->length() == nV)
    54505540     {
    5451        newRing = VMrRefine(curr_weight,target_weight); //define a new ring with ordering "(a(curr_weight),Wp(target_weight))"
     5541       newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
    54525542     }
    54535543     else
     
    54735563#endif
    54745564    idSkipZeroes(M);
    5475 #ifdef CHECK_IDEAL_MWALK
    5476     PrintS("\n//** Mwalk: computed M.\n");
    5477     idString(M, "M");
    5478 #endif
     5565//#ifdef CHECK_IDEAL_MWALK
     5566    if(printout > 2)
     5567    {
     5568      idString(M, "//** Mrwalk: M");
     5569    }
     5570//#endif
    54795571    //change the ring to baseRing
    54805572    rChangeCurrRing(baseRing);
     
    54865578    to = clock();
    54875579#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
     5580    // compute a representation of the generators of submod (M) with respect to those of mod (Gomega),
     5581    // where Gomega is a reduced Groebner basis w.r.t. the current ring
    54895582    F = MLifttwoIdeal(Gomega2, M1, G);
    54905583#ifdef TIME_TEST
    54915584    tlift = tlift + clock() - to;
    54925585#endif
    5493 #ifdef CHECK_IDEAL_MWALK
    5494     idString(F, "F");
    5495 #endif
     5586//#ifdef CHECK_IDEAL_MWALK
     5587    if(printout > 2)
     5588    {
     5589      idString(F, "//** Mrwalk: F");
     5590    }
     5591//#endif
    54965592    idDelete(&Gomega2);
    54975593    idDelete(&M1);
     
    55025598#ifdef TIME_TEST
    55035599    to = clock();
    5504 #endif
    5505     //G = kStd(F1,NULL,testHomog,NULL,NULL,0,0,NULL);
    5506 #ifdef TIME_TEST
    55075600    tstd = tstd + clock() - to;
    55085601#endif
    55095602    idSkipZeroes(G);
    5510 #ifdef CHECK_IDEAL_MWALK
    5511     idString(G, "G");
    5512 #endif
     5603//#ifdef CHECK_IDEAL_MWALK
     5604    if(printout > 2)
     5605    {
     5606      idString(G, "//** Mrwalk: G");
     5607    }
     5608//#endif
     5609
     5610    rChangeCurrRing(targetRing);
     5611    G = idrMoveR(G,newRing,currRing);
     5612
     5613    // test whether target cone is reached
     5614    if(reduction !=0 && test_w_in_ConeCC(G,curr_weight) == 1)
     5615    {
     5616      baseRing = currRing;
     5617      break;
     5618    }
     5619
     5620    rChangeCurrRing(newRing);
     5621    G = idrMoveR(G,targetRing,currRing);
     5622    baseRing = currRing;
     5623
     5624    NEXT_VECTOR:
    55135625#ifdef TIME_TEST
    55145626    to = clock();
    55155627#endif
    5516     intvec* next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, pert_deg);//next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
     5628    next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
    55175629#ifdef TIME_TEST
    55185630    tnw = tnw + clock() - to;
    55195631#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
    5528       Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
    5529       if(target_M->length() == nV)
    5530       {
    5531         newRing = VMrDefault(target_weight); // define a new ring with ordering "(a(curr_weight),lp)
    5532       }
    5533       else
    5534       {
    5535         newRing = VMatrDefault(target_M);
    5536       }
    5537       rChangeCurrRing(newRing);
    5538       Gomega1 = idrMoveR(Gomega, baseRing,currRing);
     5632
     5633#ifdef TIME_TEST
     5634    to = clock();
     5635#endif
     5636    Gomega = MwalkInitialForm(G, next_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
     5637#ifdef TIME_TEST
     5638    tif = tif + clock()-to; //time for computing initial form ideal
     5639#endif
     5640
     5641    //polylength = 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise
     5642    polylength = lengthpoly(Gomega);
     5643    if(polylength > 0)
     5644    {
     5645      //there is a polynomial in Gomega with at least 3 monomials,
     5646      //low-dimensional facet of the cone
     5647      delete next_weight;
     5648#ifdef TIME_TEST
     5649    to = clock();
     5650#endif
     5651      next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, pert_deg);
     5652#ifdef TIME_TEST
     5653    tnw = tnw + clock() - to;
     5654#endif
    55395655      idDelete(&Gomega);
    5540 #ifdef CHECK_IDEAL_MWALK
    5541       idString(Gomega1, "Gomega");
    5542 #endif
    5543       M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL);
    5544 #ifdef CHECK_IDEAL_MWALK
    5545       idString(M,"M");
    5546 #endif
    5547       rChangeCurrRing(baseRing);
    5548       M1 =  idrMoveR(M, newRing,currRing);
    5549       idDelete(&M);
    5550       Gomega2 = idrMoveR(Gomega1, newRing,currRing);
    5551       idDelete(&Gomega1);
    5552       F = MLifttwoIdeal(Gomega2, M1, G);
    5553 #ifdef CHECK_IDEAL_MWALK
    5554       idString(F,"F");
    5555 #endif
    5556       idDelete(&Gomega2);
    5557       idDelete(&M1);
    5558       rChangeCurrRing(newRing); // change the ring to newRing
    5559       G = idrMoveR(F,baseRing,currRing);
    5560       idDelete(&F);
     5656#ifdef TIME_TEST
     5657    to = clock();
     5658#endif
     5659      Gomega = MwalkInitialForm(G, next_weight);
     5660#ifdef TIME_TEST
     5661    tif = tif + clock()-to; //time for computing initial form ideal
     5662#endif
     5663    }
     5664
     5665    // test whether target weight vector is reached
     5666    if(MivComp(next_weight, ivNull) == 1 || MivComp(target_weight,curr_weight) == 1)
     5667    {
    55615668      baseRing = currRing;
    5562       si_opt_1 = save1; //set original options, e. g. option(RedSB)
    5563       idSkipZeroes(G);
    5564 #ifdef TIME_TEST
    5565       to = clock();
    5566 #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   //    }
    5571 #ifdef TIME_TEST
    5572       tred = tred + clock() - to;
    5573 #endif
    5574       idSkipZeroes(G);
    55755669      delete next_weight;
    55765670      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
     5671    }
     5672
     5673//#ifdef PRINT_VECTORS
     5674    if(printout > 0)
     5675    {
     5676      MivString(curr_weight, target_weight, next_weight);
     5677    }
     5678//#endif
     5679
    55845680    for(i=nV-1; i>=0; i--)
    55855681    {
    5586       (*tmp_weight)[i] = (*curr_weight)[i];
    55875682      (*curr_weight)[i] = (*next_weight)[i];
    55885683    }
    55895684    delete next_weight;
    55905685  }
     5686  baseRing = currRing;
    55915687  rChangeCurrRing(XXRing);
    55925688  ideal result = idrMoveR(G,baseRing,currRing);
    55935689  idDelete(&G);
    5594 /*#ifdef CHECK_IDEAL_MWALK
    5595   pDelete(&p);
    5596 #endif*/
    5597   delete tmp_weight;
     5690  si_opt_1 = save1; //set original options, e. g. option(RedSB)
    55985691  delete ivNull;
    5599   delete exivlp;
    56005692#ifndef BUCHBERGER_ALG
    56015693  delete last_omega;
    56025694#endif
     5695  Print("\n//** Mrwalk: Groebner Walk took %d steps.\n", nstep);
    56035696#ifdef TIME_TEST
    5604   Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep);
    56055697  TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep);
    5606   Print("\n//** Mwalk: Ergebnis.\n");
    56075698  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
    56085699  //Print("\n// Overflow_Error? (%d)\n", Overflow_Error);
     
    57515842// use kStd, if nP = 0, else call LastGB
    57525843ideal Mpwalk(ideal Go, int op_deg, int tp_deg,intvec* curr_weight,
    5753              intvec* target_weight, int nP)
     5844             intvec* target_weight, int nP, int reduction, int printout)
    57545845{
     5846  BITSET save1 = si_opt_1; // save current options
     5847  if(reduction == 0)
     5848  {
     5849    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     5850    //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     5851    //si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB));
     5852  }
    57555853  Set_Error(FALSE  );
    57565854  Overflow_Error = FALSE;
     
    57665864  nstep = 0;
    57675865  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;
     5866  BOOLEAN endwalks = FALSE;
     5867
     5868  ideal Gomega, M, F, FF, G, Gomega1, Gomega2, M1,F1,Eresult,ssG;
    57715869  ring newRing, oldRing, TargetRing;
    57725870  intvec* iv_M_dp;
     
    57905888  ring XXRing = currRing;
    57915889
    5792 
    57935890  to = clock();
    5794   /* perturbs the original vector */
     5891  // perturbs the original vector
    57955892  if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) := "dp"
    57965893  {
     
    58095906      DefRingPar(curr_weight);
    58105907    else
    5811       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 1
     5908      rChangeCurrRing(VMrDefault(curr_weight));
    58125909
    58135910    G = idrMoveR(Go, XXRing,currRing);
     
    58245921  ring HelpRing = currRing;
    58255922
    5826   /* perturbs the target weight vector */
     5923  // perturbs the target weight vector
    58275924  if(tp_deg > 1 && tp_deg <= nV)
    58285925  {
     
    58305927      DefRingPar(target_weight);
    58315928    else
    5832       rChangeCurrRing(VMrDefault(target_weight)); // Aenderung 2
     5929      rChangeCurrRing(VMrDefault(target_weight));
    58335930
    58345931    TargetRing = currRing;
     
    58375934    {
    58385935      iv_M_lp = MivMatrixOrderlp(nV);
    5839       //ivString(iv_M_lp, "iv_M_lp");
    5840       //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg);
    58415936      target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
    58425937    }
     
    58445939    {
    58455940      iv_M_lp = MivMatrixOrder(target_weight);
    5846       //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg);
    58475941      target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
    58485942    }
     
    58525946    G = idrMoveR(ssG, TargetRing,currRing);
    58535947  }
    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   */
     5948    if(printout > 0)
     5949    {
     5950      Print("\n//** Mpwalk: Perturbation Walk of degree (%d,%d):",op_deg,tp_deg);
     5951      ivString(curr_weight, "//** Mpwalk: new current weight");
     5952      ivString(target_weight, "//** Mpwalk: new target weight");
     5953    }
    58595954  while(1)
    58605955  {
    58615956    nstep ++;
    58625957    to = clock();
    5863     /* compute an initial form ideal of <G> w.r.t. the weight vector
    5864        "curr_weight" */
     5958    // compute an initial form ideal of <G> w.r.t. the weight vector
     5959    // "curr_weight"
    58655960    Gomega = MwalkInitialForm(G, curr_weight);
    5866 
     5961//#ifdef CHECK_IDEAL_MWALK
     5962    if(printout > 1)
     5963    {
     5964      idString(Gomega,"//** Mpwalk: Gomega");
     5965    }
     5966//#endif
     5967    if(reduction == 0 && nstep > 1)
     5968    {
     5969/*
     5970      // check whether weight vector is in the interior of the cone
     5971      while(1)
     5972      {
     5973        FF = middleOfCone(G,Gomega);
     5974        if(FF != NULL)
     5975        {
     5976          idDelete(&G);
     5977          G = idCopy(FF);
     5978          idDelete(&FF);
     5979          next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
     5980          if(printout > 0)
     5981          {
     5982            MivString(curr_weight, target_weight, next_weight);
     5983          }
     5984        }
     5985        else
     5986        {
     5987          break;
     5988        }
     5989        for(i=nV-1; i>=0; i--)
     5990        {
     5991          (*curr_weight)[i] = (*next_weight)[i];
     5992        }
     5993        Gomega = MwalkInitialForm(G, curr_weight);
     5994        if(printout > 1)
     5995        {
     5996          idString(Gomega,"//** Mpwalk: Gomega");
     5997        }
     5998      }
     5999*/
     6000      FF = middleOfCone(G,Gomega);
     6001      if(FF != NULL)
     6002      {
     6003        idDelete(&G);
     6004        G = idCopy(FF);
     6005        idDelete(&FF);
     6006        goto NEXT_VECTOR;
     6007      }
     6008    }
    58676009
    58686010#ifdef ENDWALKS
    5869     if(endwalks == 1){
     6011    if(endwalks == TRUE)
     6012    {
    58706013      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
    58716014      idElements(G, "G");
    5872       // idElements(Gomega, "Gw");
    58736015      headidString(G, "G");
    5874       //headidString(Gomega, "Gw");
    58756016    }
    58766017#endif
     
    58916032      DefRingPar(curr_weight);
    58926033    else
    5893       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 3
     6034      rChangeCurrRing(VMrDefault(curr_weight));
    58946035
    58956036    newRing = currRing;
     
    58976038
    58986039#ifdef ENDWALKS
    5899     if(endwalks==1)
     6040    if(endwalks==TRUE)
    59006041    {
    59016042      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     
    59126053    tim = clock();
    59136054    to = clock();
    5914     /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
     6055    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
    59156056#ifdef  BUCHBERGER_ALG
    59166057    M = MstdhomCC(Gomega1);
     
    59186059    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
    59196060    delete hilb_func;
    5920 #endif // BUCHBERGER_ALG
    5921 
    5922     if(endwalks == 1){
     6061#endif
     6062//#ifdef CHECK_IDEAL_MWALK
     6063      if(printout > 2)
     6064      {
     6065        idString(M,"//** Mpwalk: M");
     6066      }
     6067//#endif
     6068
     6069    if(endwalks == TRUE){
    59236070      xtstd = xtstd+clock()-to;
    59246071#ifdef ENDWALKS
     
    59306077      tstd=tstd+clock()-to;
    59316078
    5932     /* change the ring to oldRing */
     6079    // change the ring to oldRing
    59336080    rChangeCurrRing(oldRing);
    59346081    M1 =  idrMoveR(M, newRing,currRing);
    59356082    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
    5936 
    5937     //if(endwalks==1)  PrintS("\n// Lifting is working:..");
    59386083
    59396084    to=clock();
     
    59426087       Gomega is a reduced Groebner basis w.r.t. the current ring */
    59436088    F = MLifttwoIdeal(Gomega2, M1, G);
    5944     if(endwalks != 1)
     6089    if(endwalks == FALSE)
    59456090      tlift = tlift+clock()-to;
    59466091    else
    59476092      xtlift=clock()-to;
    59486093
     6094//#ifdef CHECK_IDEAL_MWALK
     6095    if(printout > 2)
     6096    {
     6097      idString(F,"//** Mpwalk: F");
     6098    }
     6099//#endif
     6100
    59496101    idDelete(&M1);
    59506102    idDelete(&Gomega2);
    59516103    idDelete(&G);
    59526104
    5953     /* change the ring to newRing */
     6105    // change the ring to newRing
    59546106    rChangeCurrRing(newRing);
    5955     F1 = idrMoveR(F, oldRing,currRing);
    5956 
    5957     //if(endwalks==1)PrintS("\n// InterRed is working now:");
    5958 
     6107    if(reduction == 0)
     6108    {
     6109      G = idrMoveR(F,oldRing,currRing);
     6110    }
     6111    else
     6112    {
     6113      F1 = idrMoveR(F, oldRing,currRing);
     6114      if(printout > 2)
     6115      {
     6116        PrintS("\n //** Mpwalk: reduce the Groebner basis.\n");
     6117      }
     6118      to=clock();
     6119      G = kInterRedCC(F1, NULL);
     6120      if(endwalks == FALSE)
     6121        tred = tred+clock()-to;
     6122      else
     6123        xtred=clock()-to;
     6124      idDelete(&F1);
     6125    }
     6126    if(endwalks == TRUE)
     6127      break;
     6128
     6129    NEXT_VECTOR:
    59596130    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 */
     6131    // compute a next weight vector
    59746132    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
    59756133    tnw=tnw+clock()-to;
    5976 #ifdef PRINT_VECTORS
    5977     MivString(curr_weight, target_weight, next_weight);
    5978 #endif
     6134//#ifdef PRINT_VECTORS
     6135    if(printout > 0)
     6136    {
     6137      MivString(curr_weight, target_weight, next_weight);
     6138    }
     6139//#endif
    59796140
    59806141    if(Overflow_Error == TRUE)
     
    59946155    }
    59956156    if(MivComp(next_weight, target_weight) == 1)
    5996       endwalks = 1;
     6157      endwalks = TRUE;
    59976158
    59986159    for(i=nV-1; i>=0; i--)
     
    60006161
    60016162    delete next_weight;
    6002   }//while
     6163  }//end of while-loop
    60036164
    60046165  if(tp_deg != 1)
     
    60146175        DefRingPar(orig_target);
    60156176      else
    6016         rChangeCurrRing(VMrDefault(orig_target)); //Aenderung
     6177        rChangeCurrRing(VMrDefault(orig_target));
    60176178
    60186179    TargetRing=currRing;
     
    60306191    if( ntestw != 1 || ntwC == 0)
    60316192    {
    6032       /*
    6033       if(ntestw != 1){
     6193      if(ntestw != 1 && printout >2)
     6194      {
    60346195        ivString(pert_target_vector, "tau");
    60356196        PrintS("\n// ** perturbed target vector doesn't stay in cone!!");
    60366197        Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
    6037         idElements(F1, "G");
    6038       }
    6039       */
     6198        //idElements(F1, "G");
     6199      }
    60406200      // LastGB is "better" than the kStd subroutine
    60416201      to=clock();
     
    60686228    Eresult = idrMoveR(G, newRing,currRing);
    60696229  }
     6230  si_opt_1 = save1; //set original options, e. g. option(RedSB)
    60706231  delete ivNull;
    60716232  if(tp_deg != 1)
     
    60826243             tnw+xtnw);
    60836244
    6084   Print("\n// pSetm_Error = (%d)", ErrorCheck());
    6085   Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep,  Overflow_Error);
    6086 #endif
     6245  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
     6246  //Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep,  Overflow_Error);
     6247#endif
     6248  Print("\n//** Mpwalk: Perturbation Walk took %d steps.\n", nstep);
     6249  return(Eresult);
     6250}
     6251
     6252/*******************************************************
     6253 * THE PERTURBATION WALK ALGORITHM WITH RANDOM ELEMENT *
     6254 *******************************************************/
     6255ideal Mprwalk(ideal Go, intvec* orig_M, intvec* target_M, int weight_rad,
     6256              int op_deg, int tp_deg, int nP, int reduction, int printout)
     6257{
     6258  BITSET save1 = si_opt_1; // save current options
     6259  if(reduction == 0)
     6260  {
     6261    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     6262    si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     6263  }
     6264  Set_Error(FALSE);
     6265  Overflow_Error = FALSE;
     6266  //Print("// pSetm_Error = (%d)", ErrorCheck());
     6267
     6268  clock_t  tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
     6269  xtextra=0;
     6270  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0;
     6271  tinput = clock();
     6272
     6273  clock_t tim;
     6274
     6275  nstep = 0;
     6276  int i, ntwC=1, ntestw=1, polylength, nV = currRing->N;
     6277  BOOLEAN endwalks = FALSE;
     6278
     6279  ideal Gomega, M, F, FF, G, Gomega1, Gomega2, M1,F1,Eresult,ssG;
     6280  ring newRing, oldRing, TargetRing;
     6281  intvec* iv_M_dp;
     6282  intvec* iv_M_lp;
     6283  intvec* exivlp = Mivlp(nV);
     6284  intvec* curr_weight = new intvec(nV);
     6285  intvec* target_weight = new intvec(nV);
     6286  for(i=0; i<nV; i++)
     6287  {
     6288    (*curr_weight)[i] = (*orig_M)[i];
     6289    (*target_weight)[i] = (*target_M)[i];
     6290  }
     6291  intvec* orig_target = target_weight;
     6292  intvec* pert_target_vector = target_weight;
     6293  intvec* ivNull = new intvec(nV);
     6294  intvec* iv_dp = MivUnit(nV);// define (1,1,...,1)
     6295#ifndef BUCHBERGER_ALG
     6296  intvec* hilb_func;
     6297#endif
     6298  intvec* next_weight;
     6299
     6300  // to avoid (1,0,...,0) as the target vector
     6301  intvec* last_omega = new intvec(nV);
     6302  for(i=nV-1; i>0; i--)
     6303    (*last_omega)[i] = 1;
     6304  (*last_omega)[0] = 10000;
     6305
     6306  ring XXRing = currRing;
     6307
     6308  to = clock();
     6309  // perturbs the original vector
     6310  if(orig_M->length() == nV)
     6311  {
     6312    if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) := "dp"
     6313    {
     6314      G = MstdCC(Go);
     6315      tostd = clock()-to;
     6316      if(op_deg != 1)
     6317      {
     6318        iv_M_dp = MivMatrixOrderdp(nV);
     6319        curr_weight = MPertVectors(G, iv_M_dp, op_deg);
     6320      }
     6321    }
     6322    else
     6323    {
     6324      //define ring order := (a(curr_weight),lp);
     6325      if (rParameter(currRing) != NULL)
     6326        DefRingPar(curr_weight);
     6327      else
     6328        rChangeCurrRing(VMrDefault(curr_weight));
     6329
     6330      G = idrMoveR(Go, XXRing,currRing);
     6331      G = MstdCC(G);
     6332      tostd = clock()-to;
     6333      if(op_deg != 1)
     6334      {
     6335        iv_M_dp = MivMatrixOrder(curr_weight);
     6336        curr_weight = MPertVectors(G, iv_M_dp, op_deg);
     6337      }
     6338    }
     6339  }
     6340  else
     6341  {
     6342    rChangeCurrRing(VMatrDefault(orig_M));
     6343    G = idrMoveR(Go, XXRing,currRing);
     6344    G = MstdCC(G);
     6345    tostd = clock()-to;
     6346    if(op_deg != 1)
     6347    {
     6348      curr_weight = MPertVectors(G, orig_M, op_deg);
     6349    }
     6350  }
     6351
     6352  delete iv_dp;
     6353  if(op_deg != 1) delete iv_M_dp;
     6354
     6355  ring HelpRing = currRing;
     6356
     6357  // perturbs the target weight vector
     6358  if(target_M->length() == nV)
     6359  {
     6360    if(tp_deg > 1 && tp_deg <= nV)
     6361    {
     6362      if (rParameter(currRing) != NULL)
     6363        DefRingPar(target_weight);
     6364      else
     6365        rChangeCurrRing(VMrDefault(target_weight));
     6366
     6367      TargetRing = currRing;
     6368      ssG = idrMoveR(G,HelpRing,currRing);
     6369      if(MivSame(target_weight, exivlp) == 1)
     6370      {
     6371        iv_M_lp = MivMatrixOrderlp(nV);
     6372        target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
     6373      }
     6374      else
     6375      {
     6376        iv_M_lp = MivMatrixOrder(target_weight);
     6377        target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
     6378      }
     6379      delete iv_M_lp;
     6380      pert_target_vector = target_weight;
     6381      rChangeCurrRing(HelpRing);
     6382      G = idrMoveR(ssG, TargetRing,currRing);
     6383    }
     6384  }
     6385  else
     6386  {
     6387    if(tp_deg > 1 && tp_deg <= nV)
     6388    {
     6389      rChangeCurrRing(VMatrDefault(target_M));
     6390      TargetRing = currRing;
     6391      ssG = idrMoveR(G,HelpRing,currRing);
     6392      target_weight = MPertVectors(ssG, target_M, tp_deg);
     6393    }
     6394  }
     6395  if(printout > 0)
     6396  {
     6397    Print("\n//** Mprwalk: Random Perturbation Walk of degree (%d,%d):",op_deg,tp_deg);
     6398    ivString(curr_weight, "//** Mprwalk: new current weight");
     6399    ivString(target_weight, "//** Mprwalk: new target weight");
     6400  }
     6401
     6402#ifdef TIME_TEST
     6403  to = clock();
     6404#endif
     6405  Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
     6406#ifdef TIME_TEST
     6407  tif = tif + clock()-to; //time for computing initial form ideal
     6408#endif
     6409
     6410  while(1)
     6411  {
     6412    nstep ++;
     6413    to = clock();
     6414    // compute an initial form ideal of G w.r.t. the weight vector "curr_weight"
     6415/*    Gomega = MwalkInitialForm(G, curr_weight);
     6416    polylength = lengthpoly(Gomega);
     6417*/
     6418//#ifdef CHECK_IDEAL_MWALK
     6419    if(printout > 1)
     6420    {
     6421      idString(Gomega,"//** Mprwalk: Gomega");
     6422    }
     6423//#endif
     6424
     6425    if(reduction == 0 && nstep > 1)
     6426    {
     6427/*
     6428      // check whether weight vector is in the interior of the cone
     6429      while(1)
     6430      {
     6431        FF = middleOfCone(G,Gomega);
     6432        if(FF != NULL)
     6433        {
     6434          idDelete(&G);
     6435          G = idCopy(FF);
     6436          idDelete(&FF);
     6437          next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
     6438          if(printout > 0)
     6439          {
     6440            MivString(curr_weight, target_weight, next_weight);
     6441          }
     6442        }
     6443        else
     6444        {
     6445          break;
     6446        }
     6447        for(i=nV-1; i>=0; i--)
     6448        {
     6449          (*curr_weight)[i] = (*next_weight)[i];
     6450        }
     6451        Gomega = MwalkInitialForm(G, curr_weight);
     6452        if(printout > 1)
     6453        {
     6454          idString(Gomega,"//** Mprwalk: Gomega");
     6455        }
     6456      }
     6457*/
     6458      FF = middleOfCone(G,Gomega);
     6459      if(FF != NULL)
     6460      {
     6461        idDelete(&G);
     6462        G = idCopy(FF);
     6463        idDelete(&FF);
     6464        goto NEXT_VECTOR;
     6465      }
     6466    }
     6467
     6468#ifdef ENDWALKS
     6469    if(endwalks == TRUE)
     6470    {
     6471      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6472      idElements(G, "G");
     6473      headidString(G, "G");
     6474    }
     6475#endif
     6476
     6477#ifndef  BUCHBERGER_ALG
     6478    if(isNolVector(curr_weight) == 0)
     6479      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
     6480    else
     6481      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
     6482#endif // BUCHBERGER_ALG
     6483
     6484    oldRing = currRing;
     6485
     6486    if(target_M->length() == nV)
     6487    {
     6488      // define a new ring with ordering "(a(curr_weight),lp)
     6489      if (rParameter(currRing) != NULL)
     6490        DefRingPar(curr_weight);
     6491      else
     6492        rChangeCurrRing(VMrDefault(curr_weight));
     6493    }
     6494    else
     6495    {
     6496      rChangeCurrRing(VMatrRefine(target_M,curr_weight));
     6497    }
     6498    newRing = currRing;
     6499    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
     6500#ifdef ENDWALKS
     6501    if(endwalks == TRUE)
     6502    {
     6503      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6504      idElements(Gomega1, "Gw");
     6505      headidString(Gomega1, "headGw");
     6506      PrintS("\n// compute a rGB of Gw:\n");
     6507
     6508#ifndef  BUCHBERGER_ALG
     6509      ivString(hilb_func, "w");
     6510#endif
     6511    }
     6512#endif
     6513
     6514    tim = clock();
     6515    to = clock();
     6516    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
     6517#ifdef  BUCHBERGER_ALG
     6518    M = MstdhomCC(Gomega1);
     6519#else
     6520    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
     6521    delete hilb_func;
     6522#endif
     6523//#ifdef CHECK_IDEAL_MWALK
     6524    if(printout > 2)
     6525    {
     6526      idString(M,"//** Mprwalk: M");
     6527    }
     6528//#endif
     6529
     6530    if(endwalks == TRUE)
     6531    {
     6532      xtstd = xtstd+clock()-to;
     6533#ifdef ENDWALKS
     6534      Print("\n// time for the last std(Gw)  = %.2f sec\n",
     6535            ((double) clock())/1000000 -((double)tim) /1000000);
     6536#endif
     6537    }
     6538    else
     6539      tstd=tstd+clock()-to;
     6540
     6541    /* change the ring to oldRing */
     6542    rChangeCurrRing(oldRing);
     6543    M1 =  idrMoveR(M, newRing,currRing);
     6544    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
     6545
     6546    to=clock();
     6547    /* compute a representation of the generators of submod (M)
     6548       with respect to those of mod (Gomega).
     6549       Gomega is a reduced Groebner basis w.r.t. the current ring */
     6550    F = MLifttwoIdeal(Gomega2, M1, G);
     6551    if(endwalks == FALSE)
     6552      tlift = tlift+clock()-to;
     6553    else
     6554      xtlift=clock()-to;
     6555
     6556//#ifdef CHECK_IDEAL_MWALK
     6557    if(printout > 2)
     6558    {
     6559      idString(F,"//** Mprwalk: F");
     6560    }
     6561//#endif
     6562
     6563    idDelete(&M1);
     6564    idDelete(&Gomega2);
     6565    idDelete(&G);
     6566
     6567    // change the ring to newRing
     6568    rChangeCurrRing(newRing);
     6569    if(reduction == 0)
     6570    {
     6571      G = idrMoveR(F,oldRing,currRing);
     6572    }
     6573    else
     6574    {
     6575      F1 = idrMoveR(F, oldRing,currRing);
     6576      if(printout > 2)
     6577      {
     6578        PrintS("\n //** Mprwalk: reduce the Groebner basis.\n");
     6579      }
     6580      to=clock();
     6581      G = kInterRedCC(F1, NULL);
     6582      if(endwalks == FALSE)
     6583        tred = tred+clock()-to;
     6584      else
     6585        xtred=clock()-to;
     6586      idDelete(&F1);
     6587    }
     6588
     6589    if(endwalks == TRUE)
     6590      break;
     6591
     6592Print("\n Next weight");
     6593    NEXT_VECTOR:
     6594#ifdef TIME_TEST
     6595    to = clock();
     6596#endif
     6597    next_weight = next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
     6598#ifdef TIME_TEST
     6599    tnw = tnw + clock() - to;
     6600#endif
     6601
     6602#ifdef TIME_TEST
     6603    to = clock();
     6604#endif
     6605    Gomega = MwalkInitialForm(G, next_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
     6606#ifdef TIME_TEST
     6607    tif = tif + clock()-to; //time for computing initial form ideal
     6608#endif
     6609
     6610    //polylength = 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise
     6611    polylength = lengthpoly(Gomega);
     6612    if(polylength > 0)
     6613    {
     6614      Print("\n there is a polynomial in Gomega with at least 3 monomials");
     6615      //low-dimensional facet of the cone
     6616      delete next_weight;
     6617#ifdef TIME_TEST
     6618      to = clock();
     6619#endif
     6620      next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, op_deg);
     6621#ifdef TIME_TEST
     6622      tnw = tnw + clock() - to;
     6623#endif
     6624      idDelete(&Gomega);
     6625#ifdef TIME_TEST
     6626      to = clock();
     6627#endif
     6628      Gomega = MwalkInitialForm(G, next_weight);
     6629#ifdef TIME_TEST
     6630      tif = tif + clock()-to; //time for computing initial form ideal
     6631#endif
     6632    }
     6633
     6634/*
     6635    to=clock();
     6636    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
     6637    if(polylength > 0)
     6638    {
     6639      if(printout > 2)
     6640      {
     6641        Print("\n //**Mprwalk: there is a polynomial in Gomega with at least 3 monomials, low-dimensional facet of the cone.\n");
     6642      }
     6643      delete next_weight;
     6644      next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, op_deg);
     6645    }
     6646    tnw=tnw+clock()-to;
     6647*/
     6648//#ifdef PRINT_VECTORS
     6649    if(printout > 0)
     6650    {
     6651      MivString(curr_weight, target_weight, next_weight);
     6652    }
     6653//#endif
     6654
     6655    if(Overflow_Error == TRUE)
     6656    {
     6657      ntwC = 0;
     6658      //Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6659      //idElements(G, "G");
     6660      delete next_weight;
     6661      goto FINISH_160302;
     6662    }
     6663    if(MivComp(next_weight, ivNull) == 1){
     6664      newRing = currRing;
     6665      delete next_weight;
     6666      //Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6667      break;
     6668    }
     6669    if(MivComp(next_weight, target_weight) == 1)
     6670      endwalks = TRUE;
     6671
     6672    for(i=nV-1; i>=0; i--)
     6673      (*curr_weight)[i] = (*next_weight)[i];
     6674
     6675    delete next_weight;
     6676  }//while
     6677
     6678  if(tp_deg != 1)
     6679  {
     6680    FINISH_160302:
     6681    if(target_M->length() == nV)
     6682    {
     6683      if(MivSame(orig_target, exivlp) == 1)
     6684        if (rParameter(currRing) != NULL)
     6685          DefRingParlp();
     6686        else
     6687          VMrDefaultlp();
     6688      else
     6689        if (rParameter(currRing) != NULL)
     6690          DefRingPar(orig_target);
     6691        else
     6692          rChangeCurrRing(VMrDefault(orig_target));
     6693    }
     6694    else
     6695    {
     6696      rChangeCurrRing(VMatrDefault(target_M));
     6697    }
     6698    TargetRing=currRing;
     6699    F1 = idrMoveR(G, newRing,currRing);
     6700#ifdef CHECK_IDEAL
     6701      headidString(G, "G");
     6702#endif
     6703
     6704    // check whether the pertubed target vector stays in the correct cone
     6705    if(ntwC != 0)
     6706    {
     6707      ntestw = test_w_in_ConeCC(F1, pert_target_vector);
     6708    }
     6709    if(ntestw != 1 || ntwC == 0)
     6710    {
     6711      if(ntestw != 1 && printout > 2)
     6712      {
     6713        ivString(pert_target_vector, "tau");
     6714        PrintS("\n// **Mprwalk: perturbed target vector doesn't stay in cone.");
     6715        Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
     6716        //idElements(F1, "G");
     6717      }
     6718      // LastGB is "better" than the kStd subroutine
     6719      to=clock();
     6720      ideal eF1;
     6721      if(nP == 0 || tp_deg == 1 || MivSame(orig_target, exivlp) != 1 || target_M->length() != nV)
     6722      {
     6723        if(printout > 2)
     6724        {
     6725          PrintS("\n// ** Mprwalk: Call \"std\" to compute a Groebner basis.\n");
     6726        }
     6727        eF1 = MstdCC(F1);
     6728        idDelete(&F1);
     6729      }
     6730      else
     6731      {
     6732        if(printout > 2)
     6733        {
     6734          PrintS("\n// **Mprwalk: Call \"LastGB\" to compute a Groebner basis.\n");
     6735        }
     6736        rChangeCurrRing(newRing);
     6737        ideal F2 = idrMoveR(F1, TargetRing,currRing);
     6738        eF1 = LastGB(F2, curr_weight, tp_deg-1);
     6739        F2=NULL;
     6740      }
     6741      xtextra=clock()-to;
     6742      ring exTargetRing = currRing;
     6743
     6744      rChangeCurrRing(XXRing);
     6745      Eresult = idrMoveR(eF1, exTargetRing,currRing);
     6746    }
     6747    else
     6748    {
     6749      rChangeCurrRing(XXRing);
     6750      Eresult = idrMoveR(F1, TargetRing,currRing);
     6751    }
     6752  }
     6753  else
     6754  {
     6755    rChangeCurrRing(XXRing);
     6756    Eresult = idrMoveR(G, newRing,currRing);
     6757  }
     6758  si_opt_1 = save1; //set original options, e. g. option(RedSB)
     6759  delete ivNull;
     6760  if(tp_deg != 1)
     6761    delete target_weight;
     6762
     6763  if(op_deg != 1 )
     6764    delete curr_weight;
     6765
     6766  delete exivlp;
     6767  delete last_omega;
     6768
     6769#ifdef TIME_TEST
     6770  TimeStringFractal(tinput, tostd, tif+xtif, tstd+xtstd,0, tlift+xtlift, tred+xtred,
     6771             tnw+xtnw);
     6772
     6773  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
     6774  //Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep,  Overflow_Error);
     6775#endif
     6776  Print("\n//** Mprwalk: Perturbation Walk took %d steps.\n", nstep);
    60876777  return(Eresult);
    60886778}
     
    61106800 * Perturb the start weight vector at the top level, i.e. nlev = 1     *
    61116801 ***********************************************************************/
    6112 static ideal rec_fractal_call(ideal G, int nlev, intvec* omtmp)
     6802static ideal rec_fractal_call(ideal G, int nlev, intvec* ivtarget,
     6803             int reduction, int printout)
    61136804{
    61146805  Overflow_Error =  FALSE;
    6115   //Print("\n\n// Entering the %d-th recursion:", nlev);
    6116 
     6806  if(printout >0)
     6807  {
     6808    Print("\n\n// Entering the %d-th recursion:", nlev);
     6809  }
    61176810  int i, nV = currRing->N;
    61186811  ring new_ring, testring;
    61196812  //ring extoRing;
    6120   ideal Gomega, Gomega1, Gomega2, F, F1, Gresult, Gresult1, G1, Gt;
     6813  ideal Gomega, Gomega1, Gomega2, FF, F, F1, Gresult, Gresult1, G1, Gt;
    61216814  int nwalks = 0;
    61226815  intvec* Mwlp;
     
    61276820  intvec* next_vect;
    61286821  intvec* omega2 = new intvec(nV);
    6129   intvec* altomega = new intvec(nV);
    6130 
     6822  intvec* omtmp = new intvec(nV);
     6823  //intvec* altomega = new intvec(nV);
     6824
     6825  for(i = nV -1; i>0; i--)
     6826  {
     6827    (*omtmp)[i] = (*ivtarget)[i];
     6828  }
    61316829  //BOOLEAN isnewtarget = FALSE;
    61326830
     
    61696867    NEXT_VECTOR_FRACTAL:
    61706868    to=clock();
    6171     /* determine the next border */
     6869    // determine the next border
    61726870    next_vect = MkInterRedNextWeight(omega,omega2,G);
    61736871    xtnw=xtnw+clock()-to;
    6174 #ifdef PRINT_VECTORS
    6175     MivString(omega, omega2, next_vect);
    6176 #endif
     6872
    61776873    oRing = currRing;
    61786874
    6179     /* We only perturb the current target vector at the recursion  level 1 */
     6875    // We only perturb the current target vector at the recursion level 1
    61806876    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");
     6877      if (MivComp(next_vect, omega2) == 1)
     6878      {
     6879        // to dispense with taking initial (and lifting/interreducing
     6880        // after the call of recursion
     6881        if(printout > 0)
     6882        {
     6883          Print("\n//** rec_fractal_call: Perturb the both vectors with degree %d.",nlev);
     6884          //idElements(G, "G");
     6885        }
    61876886
    61886887        Xngleich = 1;
    61896888        nlev +=1;
    61906889
    6191         if (rParameter(currRing) != NULL)
    6192           DefRingPar(omtmp);
     6890        if(ivtarget->length() == nV)
     6891        {
     6892          if (rParameter(currRing) != NULL)
     6893            DefRingPar(omtmp);
     6894          else
     6895            rChangeCurrRing(VMrDefault(omtmp));
     6896        }
    61936897        else
    6194           rChangeCurrRing(VMrDefault1(omtmp)); //Aenderung3
    6195 
     6898        {
     6899          rChangeCurrRing(VMatrDefault(ivtarget));
     6900        }
    61966901        testring = currRing;
    61976902        Gt = idrMoveR(G, oRing,currRing);
    61986903
    6199         /* perturb the original target vector w.r.t. the current GB */
    6200         delete Xtau;
    6201         Xtau = NewVectorlp(Gt);
     6904        // perturb the original target vector w.r.t. the current GB
     6905        if(ivtarget->length() == nV)
     6906        {
     6907          delete Xtau;
     6908          Xtau = NewVectorlp(Gt);
     6909        }
     6910        else
     6911        {
     6912          delete Xtau;
     6913          Xtau = Mfpertvector(Gt,ivtarget);
     6914        }
    62026915
    62036916        rChangeCurrRing(oRing);
    62046917        G = idrMoveR(Gt, testring,currRing);
    62056918
    6206         /* perturb the current vector w.r.t. the current GB */
     6919        // perturb the current vector w.r.t. the current GB
    62076920        Mwlp = MivWeightOrderlp(omega);
    62086921        Xsigma = Mfpertvector(G, Mwlp);
     
    62176930        to=clock();
    62186931
    6219         /* to avoid the value of Overflow_Error that occur in Mfpertvector*/
     6932        // to avoid the value of Overflow_Error that occur in Mfpertvector
    62206933        Overflow_Error = FALSE;
    6221 
    62226934        next_vect = MkInterRedNextWeight(omega,omega2,G);
    62236935        xtnw=xtnw+clock()-to;
    6224 
    6225 #ifdef PRINT_VECTORS
     6936      }// end of (if MivComp(next_vect, omega2) == 1)
     6937
     6938//#ifdef PRINT_VECTORS
     6939      if(printout > 0)
     6940      {
    62266941        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)
     6942      }
     6943//#endif
     6944
     6945    // check whether the the computed vector is in the correct cone.
     6946    // If no, compute the reduced Groebner basis of an omega-homogeneous
     6947    // ideal with Buchberger's algorithm and stop this recursion step
     6948    if(Overflow_Error == TRUE || test_w_in_ConeCC(G, next_vect) != 1)  //e.g. Example s7, cyc6
    62366949    {
    62376950      delete next_vect;
    6238       if (rParameter(currRing) != NULL)
    6239       {
    6240         DefRingPar(omtmp);
     6951      if(ivtarget->length() == nV)
     6952      {
     6953        if (rParameter(currRing) != NULL)
     6954          DefRingPar(omtmp);
     6955        else
     6956          rChangeCurrRing(VMrDefault(omtmp));
    62416957      }
    62426958      else
    62436959      {
    6244         rChangeCurrRing(VMrDefault1(omtmp)); // Aenderung4
     6960        rChangeCurrRing(VMatrDefault(ivtarget));
    62456961      }
    62466962#ifdef TEST_OVERFLOW
     
    62486964      Gt = NULL; return(Gt);
    62496965#endif
    6250 
    6251       //Print("\n\n// apply BB's alg. in ring r = %s;", rString(currRing));
     6966      if(printout > 0)
     6967      {
     6968        Print("\n//** rec_fractal_call: Applying Buchberger's algorithm in ring r = %s;",
     6969              rString(currRing));
     6970      }
    62526971      to=clock();
    62536972      Gt = idrMoveR(G, oRing,currRing);
     
    62576976
    62586977      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);
     6978      //delete altomega;
     6979      if(printout > 0)
     6980      {
     6981        Print("\n//** rec_fractal_call: Overflow. Leaving the %d-th recursion with %d steps.\n",
     6982              nlev, nwalks);
     6983        //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     6984      }
     6985
    62636986      nnflow ++;
    6264 
    62656987      Overflow_Error = FALSE;
    62666988      return (G1);
     
    62776999    if (MivComp(next_vect, XivNull) == 1)
    62787000    {
    6279       if (rParameter(currRing) != NULL)
    6280         DefRingPar(omtmp);
     7001      if(ivtarget->length() == nV)
     7002      {
     7003        if (rParameter(currRing) != NULL)
     7004          DefRingPar(omtmp);
     7005        else
     7006          rChangeCurrRing(VMrDefault(omtmp));
     7007      }
    62817008      else
    6282         rChangeCurrRing(VMrDefault1(omtmp)); //Aenderung5
     7009      {
     7010        rChangeCurrRing(VMatrDefault(ivtarget));
     7011      }
    62837012
    62847013      testring = currRing;
    62857014      Gt = idrMoveR(G, oRing,currRing);
    62867015
    6287       if(test_w_in_ConeCC(Gt, omega2) == 1) {
     7016      if(test_w_in_ConeCC(Gt, omega2) == 1)
     7017      {
     7018        delete omega2;
     7019        delete next_vect;
     7020        //delete altomega;
     7021        if(printout > 0)
     7022        {
     7023          Print("\n//** rec_fractal_call: Correct cone. Leaving the %d-th recursion with %d steps.\n",
     7024              nlev, nwalks);
     7025        }
     7026        return (Gt);
     7027      }
     7028      else
     7029      {
     7030        if(printout > 0)
     7031        {
     7032          Print("\n//** rec_fractal_call: Wrong cone. Tau doesn't stay in the correct cone.\n");
     7033        }
     7034
     7035#ifndef  MSTDCC_FRACTAL
     7036        intvec* Xtautmp;
     7037        if(ivtarget->length() == nV)
     7038        {
     7039          Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp));
     7040        }
     7041        else
     7042        {
     7043          Xtautmp = Mfpertvector(Gt, ivtarget);
     7044        }
     7045#ifdef TEST_OVERFLOW
     7046      if(Overflow_Error == TRUE)
     7047      Gt = NULL; return(Gt);
     7048#endif
     7049
     7050        if(MivSame(Xtau, Xtautmp) == 1)
     7051        {
     7052          if(printout > 0)
     7053          {
     7054            Print("\n//** rec_fractal_call: Updated vectors are equal to the old vectors.\n");
     7055          }
     7056          delete Xtautmp;
     7057          goto FRACTAL_MSTDCC;
     7058        }
     7059
     7060        Xtau = Xtautmp;
     7061        Xtautmp = NULL;
     7062
     7063        for(i=nV-1; i>=0; i--)
     7064          (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i];
     7065
     7066        rChangeCurrRing(oRing);
     7067        G = idrMoveR(Gt, testring,currRing);
     7068
     7069        goto NEXT_VECTOR_FRACTAL;
     7070#endif
     7071
     7072      FRACTAL_MSTDCC:
     7073        if(printout > 0)
     7074        {
     7075          Print("\n//** rec_fractal_call: Wrong cone. Applying Buchberger's algorithm in ring = %s.\n",
     7076                rString(currRing));
     7077        }
     7078        to=clock();
     7079        G = MstdCC(Gt);
     7080        xtextra=xtextra+clock()-to;
     7081
     7082        oRing = currRing;
     7083
     7084        // update the original target vector w.r.t. the current GB
     7085        if(ivtarget->length() == nV)
     7086        {
     7087          if(MivSame(Xivinput, Xivlp) == 1)
     7088            if (rParameter(currRing) != NULL)
     7089              DefRingParlp();
     7090            else
     7091              VMrDefaultlp();
     7092          else
     7093            if (rParameter(currRing) != NULL)
     7094              DefRingPar(Xivinput);
     7095            else
     7096              rChangeCurrRing(VMrDefault(Xivinput));
     7097        }
     7098        else
     7099        {
     7100          rChangeCurrRing(VMatrRefine(ivtarget,Xivinput));
     7101        }
     7102        testring = currRing;
     7103        Gt = idrMoveR(G, oRing,currRing);
     7104
     7105        // perturb the original target vector w.r.t. the current GB
     7106        if(ivtarget->length() == nV)
     7107        {
     7108          delete Xtau;
     7109          Xtau = NewVectorlp(Gt);
     7110        }
     7111        else
     7112        {
     7113          delete Xtau;
     7114          Xtau = Mfpertvector(Gt,ivtarget);
     7115        }
     7116
     7117        rChangeCurrRing(oRing);
     7118        G = idrMoveR(Gt, testring,currRing);
     7119
     7120        delete omega2;
     7121        delete next_vect;
     7122        //delete altomega;
     7123        if(printout > 0)
     7124        {
     7125          Print("\n//** rec_fractal_call: Vectors updated. Leaving the %d-th recursion with %d steps.\n",
     7126              nlev, nwalks);
     7127          //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     7128        }
     7129        if(Overflow_Error == TRUE)
     7130          nnflow ++;
     7131
     7132        Overflow_Error = FALSE;
     7133        return(G);
     7134      }
     7135    }// end of (if next_vect==nullvector)
     7136
     7137    for(i=nV-1; i>=0; i--) {
     7138      //(*altomega)[i] = (*omega)[i];
     7139      (*omega)[i] = (*next_vect)[i];
     7140    }
     7141    delete next_vect;
     7142
     7143    to=clock();
     7144    // Take the initial form of <G> w.r.t. omega
     7145    Gomega = MwalkInitialForm(G, omega);
     7146    xtif=xtif+clock()-to;
     7147    if(printout > 1)
     7148    {
     7149      idString(Gomega,"//** rec_fractal_call: Gomega");
     7150    }
     7151
     7152    if(reduction == 0)
     7153    {
     7154      // Check whether the intermediate weight vector lies in the interior of the cone.
     7155      // If so, only perform reductions. Otherwise apply Buchberger's algorithm.
     7156      FF = middleOfCone(G,Gomega);
     7157      if( FF != NULL)
     7158      {
     7159        idDelete(&G);
     7160        G = idCopy(FF);
     7161        idDelete(&FF);
     7162        // Compue next vector.
     7163        goto NEXT_VECTOR_FRACTAL;
     7164      }
     7165    }
     7166
     7167#ifndef  BUCHBERGER_ALG
     7168    if(isNolVector(omega) == 0)
     7169      hilb_func = hFirstSeries(Gomega,NULL,NULL,omega,currRing);
     7170    else
     7171      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
     7172#endif
     7173
     7174    if(ivtarget->length() == nV)
     7175    {
     7176      if (rParameter(currRing) != NULL)
     7177        DefRingPar(omega);
     7178      else
     7179        rChangeCurrRing(VMrDefault(omega));
     7180    }
     7181    else
     7182    {
     7183      rChangeCurrRing(VMatrRefine(ivtarget,omega));
     7184    }
     7185    Gomega1 = idrMoveR(Gomega, oRing,currRing);
     7186
     7187    // Maximal recursion depth, to compute a red. GB
     7188    // Fractal walk with the alternative recursion
     7189    // alternative recursion
     7190    if(nlev == Xnlev || lengthpoly(Gomega1) == 0)
     7191    {
     7192      if(printout > 1)
     7193      {
     7194        Print("\n//** rec_fractal_call: Maximal recursion depth.\n");
     7195      }
     7196
     7197      to=clock();
     7198#ifdef  BUCHBERGER_ALG
     7199      Gresult = MstdhomCC(Gomega1);
     7200#else
     7201      Gresult =kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,omega);
     7202      delete hilb_func;
     7203#endif
     7204      xtstd=xtstd+clock()-to;
     7205    }
     7206    else
     7207    {
     7208      rChangeCurrRing(oRing);
     7209      Gomega1 = idrMoveR(Gomega1, oRing,currRing);
     7210      Gresult = rec_fractal_call(idCopy(Gomega1),nlev+1,omega,reduction,printout);
     7211    }
     7212    if(printout > 2)
     7213    {
     7214      idString(Gresult,"//** rec_fractal_call: M");
     7215    }
     7216    //convert a Groebner basis from a ring to another ring
     7217    new_ring = currRing;
     7218
     7219    rChangeCurrRing(oRing);
     7220    Gresult1 = idrMoveR(Gresult, new_ring,currRing);
     7221    Gomega2 = idrMoveR(Gomega1, new_ring,currRing);
     7222
     7223    to=clock();
     7224    // Lifting process
     7225    F = MLifttwoIdeal(Gomega2, Gresult1, G);
     7226    xtlift=xtlift+clock()-to;
     7227    if(printout > 2)
     7228    {
     7229      idString(F,"//** rec_fractal_call: F");
     7230    }
     7231    idDelete(&Gresult1);
     7232    idDelete(&Gomega2);
     7233    idDelete(&G);
     7234
     7235    rChangeCurrRing(new_ring);
     7236    //F1 = idrMoveR(F, oRing,currRing);
     7237    G = idrMoveR(F,oRing,currRing);
     7238    to=clock();
     7239    // Interreduce G
     7240    // G = kInterRedCC(F1, NULL);
     7241    xtred=xtred+clock()-to;
     7242    //idDelete(&F1);
     7243  }
     7244}
     7245
     7246/************************************************************************
     7247 * Perturb the start weight vector at the top level with random element *
     7248 ************************************************************************/
     7249static ideal rec_r_fractal_call(ideal G, int nlev, intvec* ivtarget,
     7250                int weight_rad, int reduction, int printout)
     7251{
     7252  Overflow_Error =  FALSE;
     7253  //Print("\n\n// Entering the %d-th recursion:", nlev);
     7254
     7255  int i, polylength, nV = currRing->N;
     7256  ring new_ring, testring;
     7257  //ring extoRing;
     7258  ideal Gomega, Gomega1, Gomega2, F, FF, F1, Gresult, Gresult1, G1, Gt;
     7259  int nwalks = 0;
     7260  intvec* Mwlp;
     7261#ifndef BUCHBERGER_ALG
     7262  intvec* hilb_func;
     7263#endif
     7264//  intvec* extXtau;
     7265  intvec* next_vect;
     7266  intvec* omega2 = new intvec(nV);
     7267  intvec* omtmp = new intvec(nV);
     7268  intvec* altomega = new intvec(nV);
     7269
     7270  //BOOLEAN isnewtarget = FALSE;
     7271
     7272  for(i = nV -1; i>0; i--)
     7273  {
     7274    (*omtmp)[i] = (*ivtarget)[i];
     7275  }
     7276  // to avoid (1,0,...,0) as the target vector (Hans)
     7277  intvec* last_omega = new intvec(nV);
     7278  for(i=nV-1; i>0; i--)
     7279    (*last_omega)[i] = 1;
     7280  (*last_omega)[0] = 10000;
     7281
     7282  intvec* omega = new intvec(nV);
     7283  for(i=0; i<nV; i++) {
     7284    if(Xsigma->length() == nV)
     7285      (*omega)[i] =  (*Xsigma)[i];
     7286    else
     7287      (*omega)[i] = (*Xsigma)[(nV*(nlev-1))+i];
     7288
     7289    (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i];
     7290  }
     7291
     7292   if(nlev == 1)  Xcall = 1;
     7293   else Xcall = 0;
     7294
     7295  ring oRing = currRing;
     7296
     7297  while(1)
     7298  {
     7299#ifdef FIRST_STEP_FRACTAL
     7300    /*
     7301    perturb the current weight vector only on the top level or
     7302    after perturbation of the both vectors, nlev = 2 as the top level
     7303    */
     7304    if((nlev == 1 && Xcall == 0) || (nlev == 2 && Xngleich == 1))
     7305      if(islengthpoly2(G) == 1)
     7306      {
     7307        Mwlp = MivWeightOrderlp(omega);
     7308        Xsigma = Mfpertvector(G, Mwlp);
     7309        delete Mwlp;
     7310        Overflow_Error = FALSE;
     7311      }
     7312#endif
     7313    nwalks ++;
     7314    NEXT_VECTOR_FRACTAL:
     7315    to=clock();
     7316    /* determine the next border */
     7317    next_vect = MkInterRedNextWeight(omega,omega2,G);
     7318    if(polylength > 0 && G->m[0] != NULL)
     7319    {
     7320      if(printout > 0)
     7321      {
     7322        PrintS("\n**// rec_r_fractal_call: there is a polynomial in Gomega with at least 3 monomials, low-dimensional facet of the cone.\n");
     7323      }
     7324      delete next_vect;
     7325      next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,nlev);
     7326      if(isNegNolVector(next_vect) == 1)
     7327      {
     7328        delete next_vect;
     7329        next_vect = MkInterRedNextWeight(omega,omega2,G);
     7330      }
     7331    }
     7332    xtnw=xtnw+clock()-to;
     7333
     7334    oRing = currRing;
     7335
     7336    // We only perturb the current target vector at the recursion  level 1
     7337    if(Xngleich == 0 && nlev == 1) //(ngleich == 0) important, e.g. ex2, ex3
     7338      if (MivComp(next_vect, omega2) == 1)
     7339      {
     7340        // to dispense with taking initials and lifting/interreducing
     7341        // after the call of recursion.
     7342        if(printout > 0)
     7343        {
     7344          Print("\n//** rec_r_fractal_call: Perturb both vectors with degree %d.",nlev);
     7345          //idElements(G, "G");
     7346        }
     7347        Xngleich = 1;
     7348        nlev +=1;
     7349        if(ivtarget->length() == nV)
     7350        {
     7351          if (rParameter(currRing) != NULL)
     7352            DefRingPar(omtmp);
     7353          else
     7354            rChangeCurrRing(VMrDefault(omtmp));
     7355        }
     7356        else
     7357        {
     7358          rChangeCurrRing(VMatrDefault(ivtarget));
     7359        }
     7360        testring = currRing;
     7361        Gt = idrMoveR(G, oRing,currRing);
     7362
     7363        // perturb the original target vector w.r.t. the current GB
     7364        if(ivtarget->length() == nV)
     7365        {
     7366          delete Xtau;
     7367          Xtau = NewVectorlp(Gt);
     7368        }
     7369        else
     7370        {
     7371          delete Xtau;
     7372          Xtau = Mfpertvector(Gt,ivtarget);
     7373        }
     7374
     7375        rChangeCurrRing(oRing);
     7376        G = idrMoveR(Gt,testring,currRing);
     7377
     7378        // perturb the current vector w.r.t. the current GB
     7379        Mwlp = MivWeightOrderlp(omega);
     7380        if(ivtarget->length() > nV)
     7381        {
     7382          delete Mwlp;
     7383          Mwlp = MivMatrixOrderRefine(omega,ivtarget);
     7384        }
     7385        Xsigma = Mfpertvector(G, Mwlp);
     7386        delete Mwlp;
     7387
     7388        for(i=nV-1; i>=0; i--) {
     7389          (*omega2)[i] = (*Xtau)[nV+i];
     7390          (*omega)[i] = (*Xsigma)[nV+i];
     7391        }
     7392
     7393        delete next_vect;
     7394        to=clock();
     7395
     7396        /*
     7397        to avoid the value of Overflow_Error that occur in
     7398        Mfpertvector
     7399        */
     7400        Overflow_Error = FALSE;
     7401
     7402        next_vect = MkInterRedNextWeight(omega,omega2,G);
     7403        if(G->m[0] != NULL && polylength > 0)
     7404        {
     7405         
     7406          PrintS("//** rec_r_fractal_call: there is a polynomial in Gomega with at least 3 monomials, low-dimensional facet of the cone");
     7407         
     7408          delete next_vect;
     7409          next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,nlev);
     7410          if(isNegNolVector(next_vect) == 1)
     7411          {
     7412            delete next_vect;
     7413            next_vect = MkInterRedNextWeight(omega,omega2,G);
     7414          }
     7415        }
     7416        xtnw=xtnw+clock()-to;
     7417      }
     7418//#ifdef PRINT_VECTORS
     7419      if(printout > 0)
     7420      {
     7421        MivString(omega, omega2, next_vect);
     7422      }
     7423//#endif
     7424
     7425    /* check whether the the computed vector is in the correct cone
     7426       If no, the reduced GB of an omega-homogeneous ideal will be
     7427       computed by Buchberger algorithm and stop this recursion step*/
     7428    //if(test_w_in_ConeCC(G, next_vect) != 1) //e.g. Example s7, cyc6
     7429    if(Overflow_Error == TRUE || test_w_in_ConeCC(G,next_vect) != 1)
     7430    {
     7431      delete next_vect;
     7432      if(ivtarget->length() == nV)
     7433      {
     7434        if (rParameter(currRing) != NULL)
     7435        {
     7436          DefRingPar(omtmp);
     7437        }
     7438        else
     7439        {
     7440          rChangeCurrRing(VMrDefault(omtmp));
     7441        }
     7442      }
     7443      else
     7444      {
     7445        rChangeCurrRing(VMatrDefault(ivtarget));
     7446      }
     7447#ifdef TEST_OVERFLOW
     7448      Gt = idrMoveR(G, oRing,currRing);
     7449      Gt = NULL;
     7450      return(Gt);
     7451#endif
     7452      if(printout > 0)
     7453      {
     7454        Print("\n//** rec_r_fractal_call: applying Buchberger's algorithm in ring r = %s;",
     7455              rString(currRing));
     7456      }
     7457      to=clock();
     7458      Gt = idrMoveR(G, oRing,currRing);
     7459      G1 = MstdCC(Gt);
     7460      xtextra=xtextra+clock()-to;
     7461      Gt = NULL;
     7462
     7463      delete omega2;
     7464      delete altomega;
     7465      if(printout > 0)
     7466      {
     7467        Print("\n//** rec_r_fractal_call: Leaving the %d-th recursion with %d steps.\n",
     7468              nlev, nwalks);
     7469        //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     7470      }
     7471      nnflow ++;
     7472      Overflow_Error = FALSE;
     7473      return (G1);
     7474    }
     7475    /*
     7476       If the perturbed target vector stays in the correct cone,
     7477       return the current Groebner basis.
     7478       Otherwise, return the Groebner basis computed with Buchberger's
     7479       algorithm.
     7480       Then we update the perturbed target vectors w.r.t. this GB.
     7481    */
     7482    if (MivComp(next_vect, XivNull) == 1)
     7483    {
     7484      // The computed vector is equal to the origin vector,
     7485      // because t is not defined
     7486      if(ivtarget->length() == nV)
     7487      {
     7488        if (rParameter(currRing) != NULL)
     7489          DefRingPar(omtmp);
     7490        else
     7491          rChangeCurrRing(VMrDefault(omtmp));
     7492      }
     7493      else
     7494      {
     7495        rChangeCurrRing(VMatrDefault(ivtarget));
     7496      }
     7497      testring = currRing;
     7498      Gt = idrMoveR(G, oRing,currRing);
     7499
     7500      if(test_w_in_ConeCC(Gt, omega2) == 1)
     7501      {
    62887502        delete omega2;
    62897503        delete next_vect;
    62907504        delete altomega;
    6291         //Print("\n// Leaving the %d-th recursion with %d steps ",nlev, nwalks);
    6292         //Print(" ** Overflow_Error? (%d)", Overflow_Error);
    6293 
     7505        if(printout > 0)
     7506        {
     7507          Print("\n//** rec_r_fractal_call: Leaving the %d-th recursion with %d steps.\n",
     7508                nlev, nwalks);
     7509          //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     7510        }
    62947511        return (Gt);
    62957512      }
    62967513      else
    6297       {
    6298         //ivString(omega2, "tau'");
    6299         //Print("\n//  tau' doesn't stay in the correct cone!!");
     7514      {
     7515        if(printout > 0)
     7516        {
     7517          Print("\n//** rec_r_fractal_call: target weight doesn't stay in the correct cone.\n");
     7518        }
    63007519
    63017520#ifndef  MSTDCC_FRACTAL
    6302         //07.08.03
    63037521        //ivString(Xtau, "old Xtau");
    6304         intvec* Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp));
     7522        intvec* Xtautmp;
     7523        if(ivtarget->length() == nV)
     7524        {
     7525          Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp));
     7526        }
     7527        else
     7528        {
     7529          Xtautmp = Mfpertvector(Gt, ivtarget);
     7530        }
    63057531#ifdef TEST_OVERFLOW
    63067532      if(Overflow_Error == TRUE)
     
    63307556
    63317557      FRACTAL_MSTDCC:
    6332         //Print("\n//  apply BB-Alg in ring = %s;", rString(currRing));
     7558        if(printout > 0)
     7559        {
     7560          Print("\n//** rec_r_fractal_call: apply Buchberger's algorithm in ring = %s.\n",
     7561                rString(currRing));
     7562        }
    63337563        to=clock();
    63347564        G = MstdCC(Gt);
     
    63387568
    63397569        // update the original target vector w.r.t. the current GB
    6340         if(MivSame(Xivinput, Xivlp) == 1)
    6341           if (rParameter(currRing) != NULL)
    6342             DefRingParlp();
     7570        if(ivtarget->length() == nV)
     7571        {
     7572          if(MivSame(Xivinput, Xivlp) == 1)
     7573            if (rParameter(currRing) != NULL)
     7574              DefRingParlp();
     7575            else
     7576              VMrDefaultlp();
    63437577          else
    6344             VMrDefaultlp();
     7578            if (rParameter(currRing) != NULL)
     7579              DefRingPar(Xivinput);
     7580            else
     7581              rChangeCurrRing(VMrDefault(Xivinput));
     7582        }
    63457583        else
    6346           if (rParameter(currRing) != NULL)
    6347             DefRingPar(Xivinput);
    6348           else
    6349             rChangeCurrRing(VMrDefault1(Xivinput)); //Aenderung6
    6350 
     7584        {
     7585          rChangeCurrRing(VMatrRefine(ivtarget,Xivinput));
     7586        }
    63517587        testring = currRing;
    63527588        Gt = idrMoveR(G, oRing,currRing);
    63537589
    6354         delete Xtau;
    6355         Xtau = NewVectorlp(Gt);
     7590        // perturb the original target vector w.r.t. the current GB
     7591        if(ivtarget->length() == nV)
     7592        {
     7593          delete Xtau;
     7594          Xtau = NewVectorlp(Gt);
     7595        }
     7596        else
     7597        {
     7598          delete Xtau;
     7599          Xtau = Mfpertvector(Gt,ivtarget);
     7600        }
    63567601
    63577602        rChangeCurrRing(oRing);
     
    63617606        delete next_vect;
    63627607        delete altomega;
    6363         /*
    6364           Print("\n// Leaving the %d-th recursion with %d steps,", nlev,nwalks);
    6365           Print(" ** Overflow_Error? (%d)", Overflow_Error);
    6366         */
     7608        if(printout > 0)
     7609        {
     7610          Print("\n//** rec_r_fractal_call: Leaving the %d-th recursion with %d steps.\n",
     7611                nlev,nwalks);
     7612          //Print(" ** Overflow_Error? (%d)", Overflow_Error);
     7613        }
    63677614        if(Overflow_Error == TRUE)
    63687615          nnflow ++;
     
    63717618        return(G);
    63727619      }
    6373     }
    6374 
    6375     for(i=nV-1; i>=0; i--) {
     7620    } //end of if(MivComp(next_vect, XivNull) == 1)
     7621
     7622    for(i=nV-1; i>=0; i--)
     7623    {
    63767624      (*altomega)[i] = (*omega)[i];
    63777625      (*omega)[i] = (*next_vect)[i];
     
    63807628
    63817629    to=clock();
    6382     /* Take the initial form of <G> w.r.t. omega */
     7630    // Take the initial form of <G> w.r.t. omega
    63837631    Gomega = MwalkInitialForm(G, omega);
    63847632    xtif=xtif+clock()-to;
     7633    //polylength = 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise
     7634    polylength = lengthpoly(Gomega);
     7635    if(printout > 1)
     7636    {
     7637      idString(Gomega,"//** rec_r_fractal_call: Gomega");
     7638    }
     7639
     7640    if(reduction == 0)
     7641    {
     7642      /* Check whether the intermediate weight vector lies in the interior of the cone.
     7643       * If so, only perform reductions. Otherwise apply Buchberger's algorithm. */
     7644      FF = middleOfCone(G,Gomega);
     7645      if( FF != NULL)
     7646      {
     7647        idDelete(&G);
     7648        G = idCopy(FF);
     7649        idDelete(&FF);
     7650        /* Compue next vector. */
     7651        goto NEXT_VECTOR_FRACTAL;
     7652      }
     7653    }
    63857654
    63867655#ifndef  BUCHBERGER_ALG
     
    63897658    else
    63907659      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
    6391 #endif // BUCHBERGER_ALG
    6392 
    6393     if (rParameter(currRing) != NULL)
    6394       DefRingPar(omega);
     7660#endif
     7661    if(ivtarget->length() == nV)
     7662    {
     7663      if (rParameter(currRing) != NULL)
     7664        DefRingPar(omega);
     7665      else
     7666        rChangeCurrRing(VMrDefault(omega));
     7667    }
    63957668    else
    6396       rChangeCurrRing(VMrDefault1(omega)); //Aenderung7
    6397 
     7669    {
     7670      rChangeCurrRing(VMatrRefine(ivtarget,omega));
     7671    }
    63987672    Gomega1 = idrMoveR(Gomega, oRing,currRing);
    63997673
    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)
     7674    // Maximal recursion depth, to compute a red. GB
     7675    // Fractal walk with the alternative recursion
     7676    // alternative recursion
    64047677    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       */
     7678    {
    64147679      to=clock();
    64157680#ifdef  BUCHBERGER_ALG
     
    64187683      Gresult =kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,omega);
    64197684      delete hilb_func;
    6420 #endif // BUCHBERGER_ALG
     7685#endif
    64217686      xtstd=xtstd+clock()-to;
    64227687    }
    6423     else {
     7688    else
     7689    {
    64247690      rChangeCurrRing(oRing);
    64257691      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,
     7692      Gresult = rec_r_fractal_call(idCopy(Gomega1),nlev+1,omega,weight_rad,reduction,printout);
     7693    }
     7694    if(printout > 2)
     7695    {
     7696      idString(Gresult,"//** rec_r_fractal_call: M");
     7697    }
     7698    //convert a Groebner basis from a ring to another ring
    64307699    new_ring = currRing;
    64317700
     
    64357704
    64367705    to=clock();
    6437     /* Lifting process */
     7706    // Lifting process
    64387707    F = MLifttwoIdeal(Gomega2, Gresult1, G);
    64397708    xtlift=xtlift+clock()-to;
     7709
     7710    if(printout > 2)
     7711    {
     7712      idString(F,"//** rec_r_fractal_call: F");
     7713    }
     7714
    64407715    idDelete(&Gresult1);
    64417716    idDelete(&Gomega2);
     
    64437718
    64447719    rChangeCurrRing(new_ring);
    6445     F1 = idrMoveR(F, oRing,currRing);
    6446 
     7720    //F1 = idrMoveR(F, oRing,currRing);
     7721    G = idrMoveR(F,oRing,currRing);
    64477722    to=clock();
    6448     /* Interreduce G */
    6449     G = kInterRedCC(F1, NULL);
     7723    // Interreduce G
     7724    //G = kInterRedCC(F1, NULL);
    64507725    xtred=xtred+clock()-to;
    6451     idDelete(&F1);
     7726    //idDelete(&F1);
    64527727  }
    64537728}
    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 
    68037729
    68047730
     
    68077733 *                                                                             *
    68087734 * 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"   *
     7735 * rec_fractal_call to compute the wanted Groebner basis.                      *
     7736 * At the main procedur we compute the reduced Groebner basis w.r.t. a "fast"  *
    68117737 * order, e.g. "dp" and a sequence of weight vectors which are row vectors     *
    68127738 * of a matrix. This matrix defines the given monomial order, e.g. "lp"        *
    68137739 *******************************************************************************/
    6814 ideal Mfwalk(ideal G, intvec* ivstart, intvec* ivtarget)
     7740ideal Mfwalk(ideal G, intvec* ivstart, intvec* ivtarget,
     7741             int reduction, int printout)
    68157742{
     7743  BITSET save1 = si_opt_1; // save current options
     7744  if(reduction == 0)
     7745  {
     7746    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     7747    //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     7748  }
    68167749  Set_Error(FALSE);
    68177750  Overflow_Error = FALSE;
     
    68487781      intvec* iv_dp = MivUnit(nV); // define (1,1,...,1)
    68497782      intvec* Mdp;
    6850 
    6851       if(MivSame(ivstart, iv_dp) != 1)
    6852         Mdp = MivWeightOrderdp(ivstart);
     7783      if(ivstart->length() == nV)
     7784      {
     7785        if(MivSame(ivstart, iv_dp) != 1)
     7786          Mdp = MivWeightOrderdp(ivstart);
     7787        else
     7788          Mdp = MivMatrixOrderdp(nV);
     7789      }
    68537790      else
    6854         Mdp = MivMatrixOrderdp(nV);
     7791      {
     7792        Mdp = ivstart;
     7793      }
    68557794
    68567795      Xsigma = Mfpertvector(I, Mdp);
     
    68697808  Xivlp = Mivlp(nV);
    68707809
    6871   if(MivComp(ivtarget, Xivlp)  != 1)
    6872   {
    6873     if (rParameter(currRing) != NULL)
    6874       DefRingPar(ivtarget);
     7810  if(ivtarget->length() == nV)
     7811  {
     7812    if(MivComp(ivtarget, Xivlp)  != 1)
     7813    {
     7814      if (rParameter(currRing) != NULL)
     7815        DefRingPar(ivtarget);
     7816      else
     7817        rChangeCurrRing(VMrDefault(ivtarget));
     7818
     7819      I1 = idrMoveR(I, oldRing,currRing);
     7820      Mlp = MivWeightOrderlp(ivtarget);
     7821      Xtau = Mfpertvector(I1, Mlp);
     7822    }
    68757823    else
    6876       rChangeCurrRing(VMrDefault1(ivtarget)); //Aenderung1
    6877 
    6878     I1 = idrMoveR(I, oldRing,currRing);
    6879     Mlp = MivWeightOrderlp(ivtarget);
    6880     Xtau = Mfpertvector(I1, Mlp);
     7824    {
     7825      if (rParameter(currRing) != NULL)
     7826        DefRingParlp();
     7827      else
     7828        VMrDefaultlp();
     7829
     7830      I1 = idrMoveR(I, oldRing,currRing);
     7831      Mlp =  MivMatrixOrderlp(nV);
     7832      Xtau = Mfpertvector(I1, Mlp);
     7833    }
    68817834  }
    68827835  else
    68837836  {
    6884     if (rParameter(currRing) != NULL)
    6885       DefRingParlp();
    6886     else
    6887       VMrDefaultlp();
    6888 
    6889     I1 = idrMoveR(I, oldRing,currRing);
    6890     Mlp =  MivMatrixOrderlp(nV);
     7837    rChangeCurrRing(VMatrDefault(ivtarget));
     7838    I1 = idrMoveR(I,oldRing,currRing);
     7839    Mlp =  ivtarget;
    68917840    Xtau = Mfpertvector(I1, Mlp);
    68927841  }
     
    68997848  id_Delete(&I, oldRing);
    69007849  ring tRing = currRing;
    6901 
    6902   if (rParameter(currRing) != NULL)
    6903     DefRingPar(ivstart);
     7850  if(ivtarget->length() == nV)
     7851  {
     7852    if (rParameter(currRing) != NULL)
     7853      DefRingPar(ivstart);
     7854    else
     7855      rChangeCurrRing(VMrDefault(ivstart));
     7856  }
    69047857  else
    6905     rChangeCurrRing(VMrDefault1(ivstart)); //Aenderung2
     7858  {
     7859    rChangeCurrRing(VMatrDefault(ivstart));
     7860  }
    69067861
    69077862  I = idrMoveR(I1,tRing,currRing);
     
    69147869  ring helpRing = currRing;
    69157870
    6916   J = rec_fractal_call(J, 1, ivtarget);
     7871  J = rec_fractal_call(J,1,ivtarget,reduction,printout);
    69177872
    69187873  rChangeCurrRing(oldRing);
     
    69207875  idSkipZeroes(resF);
    69217876
     7877  si_opt_1 = save1; //set original options, e. g. option(RedSB)
    69227878  delete Xivlp;
    69237879  delete Xsigma;
     
    69387894}
    69397895
    6940 ideal Mfrwalk(ideal G, intvec* ivstart, intvec* ivtarget,int weight_rad)
     7896/*******************************************************************************
     7897 * The implementation of the fractal walk algorithm with random element        *
     7898 *                                                                             *
     7899 * The main procedur Mfwalk calls the recursive Subroutine                     *
     7900 * rec_r_fractal_call to compute the wanted Groebner basis.                    *
     7901 * At the main procedure we compute the reduced Groebner basis w.r.t. a "fast" *
     7902 * order, e.g. "dp" and a sequence of weight vectors which are row vectors     *
     7903 * of a matrix. This matrix defines the given monomial order, e.g. "lp"        *
     7904 *******************************************************************************/
     7905ideal Mfrwalk(ideal G, intvec* ivstart, intvec* ivtarget,
     7906              int weight_rad, int reduction, int printout)
    69417907{
     7908  BITSET save1 = si_opt_1; // save current options
     7909  if(reduction == 0)
     7910  {
     7911    si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis
     7912    //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions
     7913  }
    69427914  Set_Error(FALSE);
    69437915  Overflow_Error = FALSE;
     
    69747946      intvec* iv_dp = MivUnit(nV); // define (1,1,...,1)
    69757947      intvec* Mdp;
    6976 
    6977       if(MivSame(ivstart, iv_dp) != 1)
    6978         Mdp = MivWeightOrderdp(ivstart);
     7948      if(ivstart->length() == nV)
     7949      {
     7950        if(MivSame(ivstart, iv_dp) != 1)
     7951          Mdp = MivWeightOrderdp(ivstart);
     7952        else
     7953          Mdp = MivMatrixOrderdp(nV);
     7954      }
    69797955      else
    6980         Mdp = MivMatrixOrderdp(nV);
     7956      {
     7957        Mdp = ivstart;
     7958      }
    69817959
    69827960      Xsigma = Mfpertvector(I, Mdp);
     
    69957973  Xivlp = Mivlp(nV);
    69967974
    6997   if(MivComp(ivtarget, Xivlp)  != 1)
    6998   {
    6999     if (rParameter(currRing) != NULL)
    7000       DefRingPar(ivtarget);
     7975  if(ivtarget->length() == nV)
     7976  {
     7977    if(MivComp(ivtarget, Xivlp)  != 1)
     7978    {
     7979      if (rParameter(currRing) != NULL)
     7980        DefRingPar(ivtarget);
     7981      else
     7982        rChangeCurrRing(VMrDefault(ivtarget));
     7983
     7984      I1 = idrMoveR(I, oldRing,currRing);
     7985      Mlp = MivWeightOrderlp(ivtarget);
     7986      Xtau = Mfpertvector(I1, Mlp);
     7987    }
    70017988    else
    7002       rChangeCurrRing(VMrDefault1(ivtarget)); //Aenderung1
    7003 
    7004     I1 = idrMoveR(I, oldRing,currRing);
    7005     Mlp = MivWeightOrderlp(ivtarget);
    7006     Xtau = Mfpertvector(I1, Mlp);
     7989    {
     7990      if (rParameter(currRing) != NULL)
     7991        DefRingParlp();
     7992      else
     7993        VMrDefaultlp();
     7994
     7995      I1 = idrMoveR(I, oldRing,currRing);
     7996      Mlp =  MivMatrixOrderlp(nV);
     7997      Xtau = Mfpertvector(I1, Mlp);
     7998    }
    70077999  }
    70088000  else
    70098001  {
    7010     if (rParameter(currRing) != NULL)
    7011       DefRingParlp();
    7012     else
    7013       VMrDefaultlp();
    7014 
    7015     I1 = idrMoveR(I, oldRing,currRing);
    7016     Mlp =  MivMatrixOrderlp(nV);
     8002    rChangeCurrRing(VMatrDefault(ivtarget));
     8003    I1 = idrMoveR(I,oldRing,currRing);
     8004    Mlp =  ivtarget;
    70178005    Xtau = Mfpertvector(I1, Mlp);
    70188006  }
     
    70258013  id_Delete(&I, oldRing);
    70268014  ring tRing = currRing;
    7027 
    7028   if (rParameter(currRing) != NULL)
    7029     DefRingPar(ivstart);
     8015  if(ivtarget->length() == nV)
     8016  {
     8017    if (rParameter(currRing) != NULL)
     8018      DefRingPar(ivstart);
     8019    else
     8020      rChangeCurrRing(VMrDefault(ivstart));
     8021  }
    70308022  else
    7031     rChangeCurrRing(VMrDefault1(ivstart)); //Aenderung2
     8023  {
     8024    rChangeCurrRing(VMatrDefault(ivstart));
     8025  }
    70328026
    70338027  I = idrMoveR(I1,tRing,currRing);
     
    70398033  ideal resF;
    70408034  ring helpRing = currRing;
    7041 //ideal G, int nlev, intvec* omtmp, int weight_rad)
    7042   J = rec_r_fractal_call(J, 1, ivtarget,weight_rad);
     8035
     8036  J = rec_r_fractal_call(J,1,ivtarget,weight_rad,reduction,printout);
    70438037
    70448038  rChangeCurrRing(oldRing);
     
    70468040  idSkipZeroes(resF);
    70478041
     8042  si_opt_1 = save1; //set original options, e. g. option(RedSB)
    70488043  delete Xivlp;
    70498044  delete Xsigma;
     
    71018096  intvec* hilb_func;
    71028097#endif
    7103   /* to avoid (1,0,...,0) as the target vector */
     8098  // to avoid (1,0,...,0) as the target vector
    71048099  intvec* last_omega = new intvec(nV);
    71058100  for(i=nV-1; i>0; i--)
     
    71168111
    71178112  to=clock();
    7118   /* compute a red. GB w.r.t. the help ring */
     8113  // compute a red. GB w.r.t. the help ring
    71198114  if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) = "dp"
    71208115    G = MstdCC(G);
     
    71258120      DefRingPar(curr_weight);
    71268121    else
    7127       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung 4
     8122      rChangeCurrRing(VMrDefault(curr_weight));
    71288123    G = idrMoveR(G, XXRing,currRing);
    71298124    G = MstdCC(G);
     
    71518146      DefRingPar(curr_weight);
    71528147    else
    7153       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 5
     8148      rChangeCurrRing(VMrDefault(curr_weight));
    71548149    to=clock();
    71558150    Gw = idrMoveR(G, exring,currRing);
     
    71868181      DefRingPar(curr_weight);
    71878182    else
    7188       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 6
     8183      rChangeCurrRing(VMrDefault(curr_weight));
    71898184
    71908185    newRing = currRing;
     
    72718266          DefRingPar(target_tmp);
    72728267        else
    7273           rChangeCurrRing(VMrDefault(target_tmp)); // Aenderung 8
     8268          rChangeCurrRing(VMrDefault(target_tmp));
    72748269
    72758270      lpRing = currRing;
     
    73318326          DefRingPar(target_tmp);
    73328327        else
    7333           rChangeCurrRing(VMrDefault(target_tmp)); //Aenderung 9
     8328          rChangeCurrRing(VMrDefault(target_tmp));
    73348329
    73358330      lpRing = currRing;
     
    75348529    else
    75358530    {
    7536       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung 10
     8531      rChangeCurrRing(VMrDefault(curr_weight));
    75378532    }
    75388533    G = idrMoveR(G, XXRing,currRing);
     
    75678562    else
    75688563    {
    7569       rChangeCurrRing(VMrDefault(curr_weight)); //Aenderung 11
     8564      rChangeCurrRing(VMrDefault(curr_weight));
    75708565    }
    75718566    to=clock();
     
    76108605    else
    76118606    {
    7612       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung 12
     8607      rChangeCurrRing(VMrDefault(curr_weight));
    76138608    }
    76148609    newRing = currRing;
     
    77798774        else
    77808775        {
    7781           rChangeCurrRing(VMrDefault(target_tmp)); // Aenderung 13
     8776          rChangeCurrRing(VMrDefault(target_tmp));
    77828777        }
    77838778      }
     
    78528847        else
    78538848        {
    7854           rChangeCurrRing(VMrDefault(target_tmp)); // Aenderung 14
     8849          rChangeCurrRing(VMrDefault(target_tmp));
    78558850        }
    78568851      }
     
    80959090    else
    80969091    {
    8097       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung 15
     9092      rChangeCurrRing(VMrDefault(curr_weight));
    80989093    }
    80999094    newRing = currRing;
     
    82379232
    82389233  //Print("\n// \"Mpwalk\" (1,%d) took %d steps and %.2f sec. Overflow_Error (%d)", tp_deg, nwalk, ((double) clock()-tinput)/1000000, nOverflow_Error);
    8239 
     9234  Print("\n// Mprwalk took %d steps. Ring= %s;\n", nwalk, rString(currRing));
    82409235  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;
    8401 #endif
    8402 #ifdef PRINT_VECTORS
    8403     MivString(curr_weight, target_weight, next_weight);
    8404 #endif
    8405     if(Overflow_Error == TRUE)
    8406     {
    8407       PrintS("\n//**Mprwalk: OVERFLOW: The computed vector does not stay in cone, the result may be wrong.\n");
    8408       delete next_weight;
    8409       break;
    8410     }
    8411 
    8412     if(test_w_in_ConeCC(G,target_weight) == 1 || MivComp(next_weight, ivNull) == 1)
    8413     {
    8414       delete next_weight;
    8415       break;
    8416     }
    8417     //update tmp_weight and curr_weight
    8418     for(i=nV-1; i>=0; i--)
    8419     {
    8420       (*tmp_weight)[i] = (*curr_weight)[i];
    8421       (*curr_weight)[i] = (*next_weight)[i];
    8422     }
    8423     delete next_weight;
    8424   } //end of while-loop
    8425   Print("\n// Mprwalk took %d steps. Ring= %s;\n", nwalk, rString(currRing));
    8426   idSkipZeroes(G);
    8427   si_opt_1 = save1; //set original options, e. g. option(RedSB)
    8428   baseRing = currRing;
    8429   rChangeCurrRing(XXRing);
    8430   ideal Res = idrMoveR(G,baseRing,currRing);
    8431   delete tmp_weight;
    8432   delete ivNull;
    8433   delete exivlp;
    8434 #ifndef BUCHBERGER_ALG
    8435   delete last_omega;
    8436 #endif
    8437 #ifdef TIME_TEST
    8438   TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep);
    8439 #endif
    8440   return(Res);
    84419236}
    84429237
     
    85179312        else
    85189313        {
    8519           rChangeCurrRing(VMrDefault(cw_tmp)); // Aenderung 16
     9314          rChangeCurrRing(VMrDefault(cw_tmp));
    85209315        }
    85219316        G = idrMoveR(Go, XXRing,currRing);
     
    85919386    else
    85929387    {
    8593       rChangeCurrRing(VMrDefault(curr_weight)); // Aenderung 17
     9388      rChangeCurrRing(VMrDefault(curr_weight));
    85949389    }
    85959390    newRing = currRing;
     
    86539448      else
    86549449      {
    8655         rChangeCurrRing(VMrDefault(target_weight)); // Aenderung 18
     9450        rChangeCurrRing(VMrDefault(target_weight));
    86569451      }
    86579452      F1 = idrMoveR(G, newRing,currRing);
  • Singular/walk.h

    • Property mode changed from 100644 to 100755
    rf34cd44 rbbfc4a8  
    5252
    5353/* Okt -- Nov'01 */
    54 // compute a Groebner basis of an ideal G w.r.t. lexicographic order
     54// compute a Groebner basis of an ideal G w.r.t. lexicographic order 
    5555//ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M);
    56 ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M, ring baseRing);
     56ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M, ring baseRing, int reduction, int printout);
    5757
    5858// random walk algorithm to compute a Groebner basis
    59 ideal Mrwalk(ideal Go, intvec* curr_weight, intvec* target_weight, int weight_rad, int pert_deg, ring baseRing);
     59
     60ideal Mrwalk(ideal Go, intvec* orig_M, intvec* target_M, int weight_rad, int pert_deg, int reduction, int printout);
    6061
    6162/* the perturbation walk algorithm */
    6263
    63 ideal Mpwalk(ideal Go, int op_deg, int tp_deg,intvec* curr_weight,intvec* target_weight, int nP);
     64ideal Mpwalk(ideal Go, int op_deg, int tp_deg,intvec* curr_weight,intvec* target_weight, int nP, int reduction, int printout);
    6465
    6566/* the perturbation walk algorithm with random element */
    66 
    67 ideal Mprwalk(ideal Go, intvec* curr_weight, intvec* target_weight, int weight_rad, int op_deg, int tp_deg, ring baseRing);
     67ideal Mprwalk(ideal Go, intvec* orig_M, intvec* target_M, int weight_rad, int op_deg, int tp_deg, int nP, int reduction, int printout);
    6868
    6969/* The fractal walk algorithm */
    70 ideal Mfwalk(ideal G, intvec* ivstart, intvec* ivtarget);
     70ideal Mfwalk(ideal G, intvec* ivstart, intvec* ivtarget, int reduction, int printout);
    7171
    7272/* The fractal walk algorithm with random element */
    73 ideal Mfrwalk(ideal G, intvec* ivstart, intvec* ivtarget,int weight_rad);
     73ideal Mfrwalk(ideal G, intvec* ivstart, intvec* ivtarget, int weight_rad, int reduction, int printout);
    7474
    7575/* Implement Tran's idea */
  • Tst/Short/ok_s.lst

    r500e4b rbbfc4a8  
    4444bug_tr707
    4545bug_tr709
     46bug_tr723
    4647bug_tr724
    4748bug_genus_etc
  • libpolys/coeffs/numbers.cc

    r500e4b rbbfc4a8  
    228228number ndCopyMap(number a, const coeffs aRing, const coeffs r)
    229229{
    230   assume( getCoeffType(r) == getCoeffType(aRing) );
     230  // aRing and r need not be the same, but must be the same representation
     231  assume(aRing->rep==r->rep);
    231232  if ( nCoeff_has_simple_Alloc(r) && nCoeff_has_simple_Alloc(aRing) )
    232233    return a;
  • libpolys/coeffs/rintegers.cc

    r500e4b rbbfc4a8  
    339339{
    340340  /* dst = currRing */
    341   if (nCoeff_is_Ring_Z(src) || nCoeff_is_Ring_ModN(src) || nCoeff_is_Ring_PtoM(src))
     341  /* dst = nrn */
     342  if ((src->rep==n_rep_gmp)
     343  && (nCoeff_is_Ring_Z(src) || nCoeff_is_Ring_ModN(src) || nCoeff_is_Ring_PtoM(src)))
     344  {
     345    return ndCopyMap; //nrzCopyMap;
     346  }
     347  if ((src->rep==n_rep_gap_gmp) /*&& nCoeff_is_Ring_Z(src)*/)
    342348  {
    343349    return ndCopyMap; //nrzCopyMap;
  • libpolys/polys/ext_fields/transext.cc

    r500e4b rbbfc4a8  
    193193  {
    194194    p_Test(den, ntRing);
    195 
    196195    if(p_IsConstant(den, ntRing) && (n_IsOne(pGetCoeff(den), ntCoeffs)))
    197196    {
     
    199198      return FALSE;
    200199    }
    201 
    202200    if( !n_GreaterZero(pGetCoeff(den), ntCoeffs) )
    203201    {
     
    205203      return FALSE;
    206204    }
    207 
    208205    // test that den is over integers!?
    209 
    210   } else
     206  }
     207  else
    211208  {  // num != NULL // den == NULL
    212 
    213209//    if( COM(t) != 0 )
    214210//    {
     
    343339  if (IS0(a)) return NULL;
    344340  fraction f = (fraction)a;
    345   poly g = p_Copy(NUM(f), ntRing);
    346   poly h = NULL; if (!DENIS1(f)) h = p_Copy(DEN(f), ntRing);
    347   fraction result = (fraction)omAllocBin(fractionObjectBin);
    348   NUM(result) = g;
    349   DEN(result) = h;
     341  poly g = NUM(f);
     342  poly h = NULL;
     343  h =DEN(f);
     344  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
     345  NUM(result) = p_Copy(g,cf->extRing);
     346  DEN(result) = p_Copy(h,cf->extRing);
    350347  COM(result) = COM(f);
    351348  ntTest((number)result);
     
    384381    if( !n_GreaterZero(g, ntCoeffs) )
    385382    {
    386       NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
     383      NUM (f) = p_Neg(NUM (f), ntRing);
    387384      g = n_InpNeg(g, ntCoeffs);
    388385    }
     
    393390    if( !n_IsOne(g, ntCoeffs) )
    394391    {
    395       DEN (f) = p_NSet(g, ntRing); // update COM(f)???
     392      DEN (f) = p_NSet(g, ntRing);
    396393      COM (f) ++;
    397394      assume( DEN (f) != NULL );
     
    12501247
    12511248  if (IS0(a)) return;
     1249  if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
    12521250  if (!simpleTestsHaveAlreadyBeenPerformed)
    12531251  {
    1254     if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
    12551252
    12561253    /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
     
    12601257      p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
    12611258      COM(f) = 0;
    1262       ntTest(a); // !!!!
     1259      ntTest(a);
    12631260      return;
    12641261    }
     
    14551452    //PrintS(" den=");p_wrp(DEN(a),ntRing);PrintLn();
    14561453    definiteGcdCancellation(a, cf, FALSE);
    1457     fraction f=(fraction)a;
    1458     if ((DEN(f)!=NULL)
    1459     &&(!n_GreaterZero(pGetCoeff(DEN(f)),ntCoeffs)))
    1460     {
    1461       NUM(f)=p_Neg(NUM(f),ntRing);
    1462       DEN(f)=p_Neg(DEN(f),ntRing);
    1463       a=(number)f;
     1454    if ((DEN((fraction)a)!=NULL)
     1455    &&(!n_GreaterZero(pGetCoeff(DEN((fraction)a)),ntCoeffs)))
     1456    {
     1457      NUM((fraction)a)=p_Neg(NUM((fraction)a),ntRing);
     1458      DEN((fraction)a)=p_Neg(DEN((fraction)a),ntRing);
    14641459    }
    14651460  }
     
    18431838  fraction f = (fraction)a;
    18441839  poly g = prMapR(NUM(f), nMap, rSrc, rDst);
     1840  /* g may contain summands with coeff 0 */
     1841  poly hh=g;
     1842  poly prev=NULL;
     1843  while(hh!=NULL)
     1844  {
     1845    if (n_IsZero(pGetCoeff(hh),rDst->cf))
     1846    {
     1847      if (prev==NULL)
     1848      {
     1849        g=p_LmFreeAndNext(g,rDst);
     1850        hh=g;
     1851      }
     1852      else
     1853      {
     1854        prev->next=p_LmFreeAndNext(prev->next,rDst);
     1855        hh=prev->next;
     1856      }
     1857    }
     1858    else
     1859    {
     1860      prev=hh;
     1861      pIter(hh);
     1862    }
     1863  }
     1864  if (g==NULL) return NULL;
    18451865
    18461866  poly h = NULL;
    18471867
    18481868  if (!DENIS1(f))
     1869  {
    18491870     h = prMapR(DEN(f), nMap, rSrc, rDst);
     1871     /* h may contain summands with coeff 0 */
     1872    hh=h;
     1873    prev=NULL;
     1874    while(hh!=NULL)
     1875    {
     1876      if (n_IsZero(pGetCoeff(hh),rDst->cf))
     1877      {
     1878        if (prev==NULL)
     1879        {
     1880          h=p_LmFreeAndNext(h,rDst);
     1881          hh=h;
     1882        }
     1883        else
     1884        {
     1885          prev->next=p_LmFreeAndNext(prev->next,rDst);
     1886          hh=prev->next;
     1887        }
     1888      }
     1889      else
     1890      {
     1891        prev=hh;
     1892        pIter(hh);
     1893      }
     1894    }
     1895    if (h==NULL) WerrorS("mapping to */0");
     1896  }
    18501897
    18511898  fraction result = (fraction)omAllocBin(fractionObjectBin);
  • libpolys/polys/monomials/p_polys.cc

    r500e4b rbbfc4a8  
    38473847    if( !DENIS1((fraction)z) )
    38483848    {
    3849       if (p_IsConstant(DEN((fraction)z),srcExtRing))
    3850       {
    3851         number n=pGetCoeff(DEN((fraction)z));
    3852         zz=p_Div_nn(zz,n,srcExtRing);
    3853         p_Normalize(zz,srcExtRing);
    3854       }
    3855       else
     3849      if (!p_IsConstant(DEN((fraction)z),srcExtRing))
    38563850        WarnS("Not defined: Cannot map a rational fraction and make a polynomial out of it! Ignoring the denumerator.");
    38573851    }
     
    38853879    qq = p_PermPoly(zz, par_perm-1, srcExtRing, dst, nMap, NULL, rVar (srcExtRing)-1);
    38863880
     3881  if(nCoeff_is_transExt(srcCf)
     3882  && (!DENIS1((fraction)z))
     3883  && p_IsConstant(DEN((fraction)z),srcExtRing))
     3884  {
     3885    number n=nMap(pGetCoeff(DEN((fraction)z)),srcExtRing->cf, dstCf);
     3886    qq=p_Div_nn(qq,n,dst);
     3887    n_Delete(&n,dstCf);
     3888    p_Normalize(qq,dst);
     3889  }
    38873890  p_Test (qq, dst);
    3888 
    3889 //       poly p_PermPoly (poly p, int * perm, const ring oldRing, const ring dst, nMapFunc nMap, int *par_perm, int OldPar)
    3890 
    3891 //  assume( FALSE );  WarnS("longalg missing 2");
    38923891
    38933892  return qq;
     
    39053904    PrintS("\np_PermPoly::p: "); p_Write(p, oldRing, oldRing); PrintLn();
    39063905#endif
    3907 
    39083906  const int OldpVariables = rVar(oldRing);
    39093907  poly result = NULL;
     
    39113909  poly aq = NULL; /* the map coefficient */
    39123910  poly qq; /* the mapped monomial */
    3913 
    39143911  assume(dst != NULL);
    39153912  assume(dst->cf != NULL);
    3916 
    39173913  while (p != NULL)
    39183914  {
     
    39223918      qq = p_Init(dst);
    39233919      assume( nMap != NULL );
    3924 
    39253920      number n = nMap(p_GetCoeff(p, oldRing), oldRing->cf, dst->cf);
    3926 
    39273921      n_Test (n,dst->cf);
    3928 
    39293922      if ( nCoeff_is_algExt(dst->cf) )
    39303923        n_Normalize(n, dst->cf);
    3931 
    39323924      p_GetCoeff(qq, dst) = n;// Note: n can be a ZERO!!!
    3933       // coef may be zero:
    3934 //      p_Test(qq, dst);
    39353925    }
    39363926    else
    39373927    {
    39383928      qq = p_One(dst);
    3939 
    39403929//      aq = naPermNumber(p_GetCoeff(p, oldRing), par_perm, OldPar, oldRing); // no dst???
    39413930//      poly    n_PermNumber(const number z, const int *par_perm, const int P, const ring src, const ring dst)
    39423931      aq = n_PermNumber(p_GetCoeff(p, oldRing), par_perm, OldPar, oldRing, dst);
    3943 
    39443932      p_Test(aq, dst);
    3945 
    39463933      if ( nCoeff_is_algExt(dst->cf) )
    39473934        p_Normalize(aq,dst);
    3948 
    39493935      if (aq == NULL)
    39503936        p_SetCoeff(qq, n_Init(0, dst->cf),dst); // Very dirty trick!!!
    3951 
    39523937      p_Test(aq, dst);
    39533938    }
    3954 
    39553939    if (rRing_has_Comp(dst))
    39563940       p_SetComp(qq, p_GetComp(p, oldRing), dst);
    3957 
    39583941    if ( n_IsZero(pGetCoeff(qq), dst->cf) )
    39593942    {
     
    39813964              assume( dst->cf->extRing == NULL );
    39823965              number ee = n_Param(1, dst);
    3983 
    39843966              number eee;
    39853967              n_Power(ee, e, &eee, dst->cf); //nfDelete(ee,dst);
    3986 
    39873968              ee = n_Mult(c, eee, dst->cf);
    39883969              //nfDelete(c,dst);nfDelete(eee,dst);
     
    39933974              const int par = -perm[i];
    39943975              assume( par > 0 );
    3995 
    39963976//              WarnS("longalg missing 3");
    39973977#if 1
    39983978              const coeffs C = dst->cf;
    39993979              assume( C != NULL );
    4000 
    40013980              const ring R = C->extRing;
    40023981              assume( R != NULL );
    4003 
    40043982              assume( par <= rVar(R) );
    4005 
    40063983              poly pcn; // = (number)c
    4007 
    40083984              assume( !n_IsZero(c, C) );
    4009 
    40103985              if( nCoeff_is_algExt(C) )
    40113986                 pcn = (poly) c;
    40123987               else //            nCoeff_is_transExt(C)
    40133988                 pcn = NUM((fraction)c);
    4014 
    40153989              if (pNext(pcn) == NULL) // c->z
    40163990                p_AddExp(pcn, -perm[i], e, R);
     
    40203994                p_SetExp(mmc, -perm[i], e, R);
    40213995                p_Setm(mmc, R);
    4022 
    40233996                number nnc;
    40243997                // convert back to a number: number nnc = mmc;
     
    40274000                else //            nCoeff_is_transExt(C)
    40284001                  nnc = ntInit(mmc, C);
    4029 
    40304002                p_GetCoeff(qq, dst) = n_Mult((number)c, nnc, C);
    40314003                n_Delete((number *)&c, C);
    40324004                n_Delete((number *)&nnc, C);
    40334005              }
    4034 
    40354006              mapped_to_par=1;
    40364007#endif
     
    40754046      if (aq!=NULL)
    40764047         qq=p_Mult_q(aq,qq,dst);
    4077 
    40784048      aq = qq;
    4079 
    40804049      while (pNext(aq) != NULL) pIter(aq);
    4081 
    40824050      if (result_last==NULL)
    40834051      {
     
    40964064    }
    40974065  }
    4098 
    40994066  result=p_SortAdd(result,dst);
    41004067#else
     
    41274094#endif
    41284095  p_Test(result,dst);
    4129 
    41304096#if 0
    41314097  p_Test(result,dst);
  • libpolys/polys/polys0.cc

    r500e4b rbbfc4a8  
    142142  }
    143143  p_Normalize(p,lmRing);
    144   p_Normalize(p,lmRing); /* Manual/absfact.tst */
     144  if ((n_GetChar(lmRing->cf) == 0)
     145  && (nCoeff_is_transExt(lmRing->cf)))
     146    p_Normalize(p,lmRing); /* Manual/absfact.tst */
    145147  if ((p_GetComp(p, lmRing) == 0) || (!lmRing->VectorOut))
    146148  {
Note: See TracChangeset for help on using the changeset viewer.